Inference of compressed Potts graphical models

TL;DR

This paper introduces a regularization approach for inferring Potts graphical models that reduces complexity by grouping less frequent states, improving computational efficiency without sacrificing accuracy for high-frequency symbols.

Contribution

The study proposes a double regularization scheme combining color compression and sparsity, validated with two inference algorithms on synthetic and biological data.

Findings

Color compression preserves high-frequency symbol accuracy.

Significant reduction in parameters and computational time.

Effective application to protein family sequence data.

Abstract

We consider the problem of inferring a graphical Potts model on a population of variables, with a non-uniform number of Potts colors (symbols) across variables. This inverse Potts problem generally involves the inference of a large number of parameters, often larger than the number of available data, and, hence, requires the introduction of regularization. We study here a double regularization scheme, in which the number of colors available to each variable is reduced, and interaction networks are made sparse. To achieve this color compression scheme, only Potts states with large empirical frequency (exceeding some threshold) are explicitly modeled on each site, while the others are grouped into a single state. We benchmark the performances of this mixed regularization approach, with two inference algorithms, the Adaptive Cluster Expansion (ACE) and the PseudoLikelihood Maximization (PLM) on synthetic data obtained by sampling disordered Potts models on an Erdos-Renyi random graphs. We show in particular that color compression does not affect the quality of reconstruction of the parameters corresponding to high-frequency symbols, while drastically reducing the number of the other parameters and thus the computational time. Our procedure is also applied to multi-sequence alignments of protein families, with similar results.

Tables

Table 1. Table 1: Inferred models obtained varying the expansion thresholds for the reference data (ER05,B=1000 configurations), of Fig.  2 , with compression frequency f 0 = 0.01 subscript 𝑓 0 0.01 f_{0}=0.01 . The Table gives the cluster inclusion threshold t 𝑡 t , the number N 2 subscript 𝑁 2 N_{2} of 2-site clusters, the maximal relative error ϵ m ​ a ​ x subscript italic-ϵ 𝑚 𝑎 𝑥 \epsilon_{max} , the regularized cross entropy S c λ subscript superscript 𝑆 𝜆 𝑐 S^{\lambda}_{c} and the cross entropy S c subscript 𝑆 𝑐 S_{c} obtained with the cluster expansion at different thresholds t 𝑡 t . The entropy of the model having generated the data is S = 50.3 𝑆 50.3 S=50.3 . Possible solutions ( ϵ max < 5 ) subscript italic-ϵ max 5 (\epsilon_{\rm max}<5) are indicated by asterisks. The optimal threshold determined by the spACE procedure is 3.6 × 10 − 2 3.6 superscript 10 2 3.6\times 10^{-2} .

t	${ϵ}_{max}$	$N_2$	$S_c^{\lambda }$	$S_c$
1	17	0	56	56
$9.6\times {10}^{-2}$	${4.9}^{\ast }$	29	51.7	51.2
$3.6\times {10}^{-2}$	${3.9}^{\ast }$	55	50.6	49.9
$2.2\times {10}^{-2}$	34	90	49.7	48.8
$3.4\times {10}^{-5}$	$1^{\ast }$	1225	42.3	39

Table 2. Table 2: Comparison between importance sampling (AIS) and ACE methods to obtain the logarithm of the partition function Z 𝑍 Z , the entropy S 𝑆 S , and the Kullback-Leibler divergence (KL), at the different sparsity threshold t 𝑡 t for the reference model, ER05, data sampling: B=1000 and color compression f o = 0.01 subscript 𝑓 𝑜 0.01 f_{o}=0.01 , the optimal threshold determined by the spACE procedure is 3.6 × 10 − 2 3.6 superscript 10 2 3.6\times 10^{-2} . Fluctuations of the above values in AIS estimation over repetition of MC sampling are of the order 5 10 − 3 superscript 510 3 5\,10^{-3} for the sparse models and 5 10 − 2 superscript 510 2 5\,10^{-2} for the fully connected model.

cluster inclusion threshold t	logZ(AIS)	logZ(ACE)	S(AIS)	S(ACE)	KL(AIS)	KL(ACE)
1	28.9	28.6	56.9	56	6.4	6
$9.6\times {10}^{-2}$	32	31.9	52.6	51.4	2.2	2.3
$3.6\times {10}^{-2}$	32.5	31.8	51.8	50.7	1.5	1.5
$2.2\times {10}^{-2}$	34.2	32.8	52.9	50.9	2.5	3
$3.4\times {10}^{-5}$	36	26	57.8	49.1	7.4	5

Table 3. Table 3: KL divergences between inferred and empirical distributions on the reference data sets at different B 𝐵 B and with no color compression. The value obtained by weakly regularized PLM γ J = 1 B subscript 𝛾 𝐽 1 𝐵 \gamma_{J}=\frac{1}{B} is compared to the ones obtained by the standards strongly regularized PLM inference and weakly regularized spACE inference, used in Fig.  3 .

	PLM	PLM	ACE
B	(${\gamma }_J=\frac{1}{B}$,${\gamma }_h=\frac{0.002}{B}$)	(${\gamma }_J=\frac{50}{B}$,${\gamma }_h=\frac{0.1}{B}$)	(${\gamma }_J=\frac{1}{B}$,${\gamma }_h=\frac{0.01}{B}$)
${10}^3$	12.97	3.80	1.35
${10}^4$	2.85	1.28	0.37
${10}^5$	2.10	0.30	0.18

Table 4. Table 6: Difference between the fields fixed by regularization and the one computed with the pseudo-count for the unseen Potts variables in the approximation of independent sites respectively for the fields regularization γ h = 0.01 / B subscript 𝛾 ℎ 0.01 𝐵 \gamma_{h}=0.01/B used in ACE and γ h = 0.1 / B subscript 𝛾 ℎ 0.1 𝐵 \gamma_{h}=0.1/B used in PLM

$B$	$h_u^{\Gamma \left(ACE\right)}$	$h_u^p$	$\Delta h_u$
${10}^2$	-6.6	-6.9	0.3
${10}^3$	-8.6	- 9.2	0.5
${10}^4$	-10.7	-11.5	0.8
${10}^5$	-12.9	-13.8	0.9

Equations

$$E(\mathbf{a}|\mathbf{J})=-\sum _{i=1}^N{h}_i(a_i)-\sum _{i=1}^{N-1}\sum _{j=i+1}^N{J}_{ij}(a_i,a_j)$$

$$P(\mathbf{a}|\mathbf{J})=\frac{\exp \left(-E(\mathbf{a}|\mathbf{J})\right)}{Z(\mathbf{J})},$$

$$\begin{array}{rl}\sum _{a_1,\dots a_N}\delta (a_i,a)P(a_1\dots ,a_N|\mathbf{J})& =f_i(a),\\ \sum _{a_1,\dots a_N}\delta (a_i,a)\delta (a_j,b)P(a_1\dots ,a_N|\mathbf{J})& =f_{ij}(a,b),\end{array}$$

$$\begin{array}{rl}S(\mathbf{J}|\mathbf{f})& =\log Z(\mathbf{J})-\sum _{i=1}^N\sum _{a=1}^q{h}_i(a)f_i^s(a)-\sum _{i=1}^{N-1}\sum _{j=i+1}^N\sum _{a=1}^q\sum _{b=1}^q{J}_{ij}(a,b)f_{ij}^s(a,b),\end{array}$$

$$\Delta S={\gamma }_h\sum _{i=1}^N\sum _{a=1}^q{h}_i(a{)}^2+{\gamma }_J\sum _{i=1}^{N-1}\sum _{j=i+1}^N\sum _{a=1}^q\sum _{b=1}^q{J}_{ij}(a,b{)}^2.$$

$$h_i(a)$$

$$J_{ij}(a,b)$$

$$J_{ij}(a,c_j)=J_{ij}(c_i,b)=h_i(c_i)=0,$$

$$h_i(a)$$

$$J_{ij}(a,b)$$

$$\sum _{c=1}^q{J}_{ij}(a,c)=\sum _{c=1}^q{J}_{ij}(c,a)=\sum _{c=1}^q{h}_i(c)=0,$$

$$h_i(a)$$

$$J_{ij}(a,b)$$

$${\gamma }_J\sum _{b=1}^{q_j}{J}_{ij}(a,b)={\gamma }_h{h}_i(a),{\gamma }_J\sum _{a=1}^{q_i}{J}_{ij}(a,b)={\gamma }_h{h}_j(b),\sum _{a=1}^{q_i}{h}_i(a)=0.$$

$$f_i(a)>f_0,$$

$$h_i(a^{\prime })=h_i(q_i+1)+\log \left(\frac{f_i(a^{\prime })}{f_i(q_i+1)}\right),J_{ij}(a^{\prime },b)=J_{ij}(q_i+1,b),$$

$${ϵ}_{\mathrm{p}}=\frac{1}{Nq}\sqrt{\sum _{i,a}{\left(\frac{\Delta f_i(a)}{{\sigma }_i(a)}\right)}^2}.$$

$${ϵ}_{\mathrm{c}}=\frac{2}{N(N-1)q^2}\sqrt{\sum _{i<j,a,b}{\left(\frac{\Delta c_{i,j}(a,b)}{{\sigma }_{i,j}^c(a,b)}\right)}^2},$$

$${ϵ}_{\max }=\underset{\{i,j,a,b\}}{\max }\frac{1}{\sqrt{2\log (M)}}\left(\frac{|\Delta f_i(a)|}{{\sigma }_i(a)},\frac{|\Delta f_{ij}(a,b)|}{{\sigma }_{ij}(ab)}\right),$$

$$D(P_{{\mathbf{J}}^{real}}||P_{{\mathbf{J}}^{inf}})$$

$$\Delta f=\sqrt{\frac{{\sum }_i{\sum }_a{\left(f_i^{gen}(a)-f_i^{test}(a)\right)}^2}{{\sum }_i{q}_i}},\Delta \text{corr}=\sqrt{\frac{{\sum }_{ij}{\sum }_{ab}(co{v}_{ij}^{test}(a,b)-co{v}_{ij}^{gen}(a,b){)}^2}{{\sum }_{ij}{q}_i\cdot q_j}}.$$

$$\Delta h=\sqrt{\frac{{\sum }_i{\sum }_a{\left(h_i^{inf}(a)-h_i^{real}(a)\right)}^2}{{\sum }_i{q}_i}},\Delta J=\sqrt{\frac{{\sum }_{ij}{\sum }_{ab}(J_{ij}^{real}(a,b)-J_{ij}^{inf}(a,b){)}^2}{{\sum }_{ij}{q}_i\cdot q_j}};$$

$$S(\mathbf{J}|\mathbf{f})=\sum _{\Gamma }\Delta S_{\Gamma },$$

$$\Delta S_{\Gamma }=S_{\Gamma }-\sum _{{\Gamma }^{\prime }\subset \Gamma }\Delta S_{{\Gamma }^{\prime }}.$$

$$\mathbf{J}(t)=\sum _k\sum _{\Gamma \in L_k^{\prime }(t)}\Delta {\mathbf{J}}_{\Gamma },\Delta {\mathbf{J}}_{\Gamma }={\mathbf{J}}_{\Gamma }-\sum _{{\Gamma }^{\prime }\subset \Gamma }\Delta {\mathbf{J}}_{{\Gamma }^{\prime }}.$$

$$S_B=-\sum _{\mathbf{\sigma }}{P}_{{\mathbf{J}}^B}(\mathbf{\sigma })\log P_{{\mathbf{J}}^B}(\mathbf{\sigma }),$$

$$P_{{\mathbf{J}}^B}(\mathbf{\sigma })=\frac{\exp \left({\sum }_{i=1}^N{h}_i^B{\sigma }_i+{\sum }_{k<l}^N{J}_{kl}^B{\sigma }_k{\sigma }_l\right)}{{\mathcal{Z}}_B}.$$

$$D(P_{{\mathbf{J}}^{true}}||P_{{\mathbf{J}}^B})$$

$$=-S_{true}-\sum _{\mathbf{\sigma }}{P}_{{\mathbf{J}}^{true}}(\mathbf{\sigma })\left\{\sum _i{h}_i^B{\sigma }_i+\sum _{k<l}{J}_{kl}^B{\sigma }_k{\sigma }_l-\log {\mathcal{Z}}^B\right\}$$

$$=-S_{true}+\log {\mathcal{Z}}^B-\sum _{\mathbf{\sigma }}{P}_{{\mathbf{J}}^{true}}(\mathbf{\sigma })\left\{\sum _i{h}_i^B{\sigma }_i+\sum _{k<l}{J}_{kl}^B{\sigma }_k{\sigma }_l\right\}.$$

$$\log {\mathcal{Z}}^B=S_B+\sum _{\mathbf{\sigma }}{P}_{{\mathbf{J}}^B}(\mathbf{\sigma })\left\{\sum _i{h}_i^B{\sigma }_i+\sum _{k<l}{J}_{kl}^B{\sigma }_k{\sigma }_l\right\}.$$

$$D(P_{{\mathbf{J}}^{true}}||P_{{\mathbf{J}}^B})=$$

$$+\sum _{\mathbf{\sigma }}{P}_{{\mathbf{J}}^B}(\mathbf{\sigma })\left\{\sum _i{h}_i^B{\sigma }_i+\sum _{k<l}{J}_{kl}^B{\sigma }_k{\sigma }_l\right\}.$$

Full text

I Introduction

Graphical models are an important tool to model dependencies and infer effective interactions among a population of variables from data. We hereafter refer to this approach as the inverse Ising problem ackley1985learning ; cocco2012adaptive ; nguyen2017inverse in the case of binary variables, and as the inverse Potts problem in the more general case of of multi-categorial variableswu1982potts . Applications include inferring functional couplings among a set of neurons from the recording of their activity cocco2017functional ; schneidman2006weak ; cocco2009neuronal , between birds in flocks bialek2012statistical and among amino acids in sequences that belong to the same protein family ekeberg2013improved . Over the last decade, it was shown that describing these protein families by Potts models, whose parameters were learned from the corresponding sequence alignment could provide information on the protein structure burger2010disentangling ; weigt2009identification ; marks2011protein ; sulkowska2012genomics ; hopf2012three ; ekeberg2013improved ; sulkowska2012genomics ; nugent2012accurate ; hopf2012three ; de2013emerging ; jones2012psicov ; morcos2014coevolutionary , help predict fitness variations following mutations nguyen2017inverse ; ferguson2013translating ; mann2014fitness ; figliuzzi2016coevolutionary ; hopf2017mutation and design new working proteins of the same family cocco2018inverse ; socolich2005evolutionary ; russ2019 . Given the computational untractability of achieving exact solutions, different effective methods have been proposed to infer the Potts parameters from sequence data, including Gaussian approximation with different priors jones2012psicov ; baldassi2014fast ; burger2010disentangling ; pensar2019high , message passing weigt2009identification , PseudoLikelihood Maximization (PLM)ekeberg2013improved ; ekeberg2014fast ; vuffray2019efficient , minimum probability flow sohl2011new , and the Adaptive Cluster Expansion (ACE) method cocco2011adaptive ; barton2016adaptive .

Even if Potts models only include pairwise interactions, the number of parameters to be inferred is huge. For an $N$-site protein, where each site can be one out of the 20 natural amino acids or an extra symbol standing for a site insertion or deletion, the number of independent parameters is $20\cdot N+20^{2}\cdot N(N-1)/2$. This gives about $10^{6}$ parameters for $N=100$ and almost $10^{8}$ parameters for $N=500$, while protein sequence alignments typically include few thousands to few tens thousand sequences. Moreover, amino-acid frequencies, and, hence, sampling quality may vary substantially from site to site, making it impossible to accurately reconstruct the complete set of Potts parameters. To overcome the problem of undersampling regularization terms are generally included. Standard $L_{2}$-regularization helps constraining parameter values, but does not change their number. $L_{1}$ and $L_{0}$-based regularization, on the contrary, may effectively remove many interaction parameters associated to low (in absolute value) connected correlations.

In this paper a simple procedure to reduce the number $q$ of Potts parameters is described and analyzed. In physics Potts states are often referred to as colors, so we call this state reduction procedure color compression. Our goal is to infer a compressed Potts model, where the number of states $q_{i}\leq q$ depends on the site $i$. The basic idea is to group together rarely observed states on each site, defined as those below a given frequency threshold $f_{0}$. This way, the number of Potts states $q_{i}$ on each site $i$ is variable, leading to the reduced number of parameters $M_{cc}=\sum_{i}^{N}q_{i}+\sum_{i<j}^{N}q_{i}\,q_{j}$. Color compression is therefore equivalent to an $L_{0}$ regularization on the number of inferred parameters associated to the Potts states. Slightly different schemes are based on grouping colors according to their entropy contributions to the site variability barton2016adaptive or to their mutual information haldane2018coevolutionary or compression to a fixed number of colors anton2018radi , and comes to a similar outcome. As expected color compression can help in limiting the computational time, in avoiding overfitting and, in a more theoretical framework, in understanding the intrinsic dimensionality of the problem, by distinguishing the parameters that can be reliably inferred from those fixed through regularization only. Color compression was already used by some of us in the Adaptive Cluster Expansion (ACE) algorithm barton2016adaptive in order to reduce the computational time and have a simpler inferred model for the analysis of protein sequence data, but its performance has not been systematically tested up to now.

The aim of the present study is to benchmark the color compression procedure in a systematic way as a function of the compression frequency threshold $f_{0}$ and of the sampling depth $B$: Correlations between variables that are grouped in the color compression procedure are lost. We expect therefore that compression is useful in the inference at low $B$ and low $f_{0}$ to discard correlations between poorly sampled colors, affected by large statistical fluctuations. Conversely, information on the correlated structure of the data will be lost for large $B$ and $f_{0}$. Consequences of compression on the performances of inferred models are here investigated for different regularization types and strengths, for two particular inference procedures: PLM with large or small $L_{2}$ regularization parameters, and the ACE procedure with a fully connected or sparse interaction graph. We then introduce a procedure to recover, after inference of the compressed Potts model, a full Potts model over all possible $q$ states, which we refer to as color decompression. Decompression is necessary to compare the inferred parameters to the ground truth (when available), to compare the quality of the inferences for different color compression strengths (and therefore different $q_{i}$ values), and, more generally, when the model is used to predict the behavior of poorly sampled variables in the original data set.

The first part of the paper briefly sketches the methodological background and the inference algorithms (Section II). In Section III the procedures of color compression and decompression are introduced. We then assess the performances of the procedures on synthetic data generated from Potts model on random graphs in Section IV. In Section V we show an illustrative example on fitness prediction for real proteins, to verify that the results obtained on synthetic data model translate to real cases. Section VI shows the gain our color compression procedure provides in terms of computational time. Some conclusion and perspectives are presented in Sec. VII.

II Reminder on Inference and Algorithms

II.1 Inverse Potts Problem

The Potts model describes a system of $N$ interacting sites, each assuming one of $q$ possible Potts states (or colors). The probability distribution of each color on each site is controlled by a set of parameters that can be divided into local fields $h_{i}(a_{i})$, depending only on one site $i$ and its color $a_{i}$, and pairwise couplings $J_{ij}(a_{i},a_{j})$, depending on the pair of sites $i,j$ and the two Potts states $a_{i}$, $a_{j}$. An energy value is associated to each system configuration ${\mathbf{a}}=a_{1},\ldots a_{N}$,

[TABLE]

and, consequently, a probability

[TABLE]

where $Z({\bf J})=\sum_{{\mathbf{a}}}\exp\left(-E({\mathbf{a}}|{\bf J})\right)$ is the partition function and ensures that all probabilities sum to one. For simplicity, here we label the set of fields and couplings as ${\bf J}$.

Given a sample of configurations, one may be interested in inferring back the model from which these samples were generated, or at least a model reproducing the statistical properties of such configurations, such as the one- and two-site frequencies, $f_{i}(a)$ and $f_{ij}(a,b)$. In general the Potts model defined above is the simplest, or maximum entropy (jaynes1982rationale, ), probabilistic model capable of reproducing the observed frequencies. In the present case we know by construction that the Potts model is not only the simplest model to fit the data, but also the real model from which the sample was generated. To reproduce the statistics of the data, the parameters $h_{i}(a)$ and $J_{ij}(a,b)$ must be chosen such that site averages and correlations in the model match those in the data, i.e.,

[TABLE]

where $\delta(a_{i},a)$ is the Kronecker delta function, which is one if the symbol $a_{i}$ at site $i$ is equal to $a$ and zero otherwise. Finding the parameters $h_{i}(a)$, $J_{ij}(a,b)$ that satisfy Eq. 3 constitutes the inverse Potts problem.

II.2 Cross-entropy and regularization

Formally, the solution to the inverse Potts problem is the set of fields and couplings that maximize the average log-likelihood or, equivalently, that minimize the cross-entropy between the data and the model. This cross-entropy can be written as

[TABLE]

where, for simplicity, we indicate the set of single and pairwise frequencies as ${\bf f}$ and the set of fields and couplings as ${\bf J}$.

To guarantee that the minimization of the cross-entropy is a well defined problem, a regularization term $\Delta S$ is added, which, n the Bayesian formulation, corresponds to prior knowledge over the parameters $\bf J$. A Gaussian prior distribution, also referred to as $L_{2}$-regularization, is a usual choice:

[TABLE]

The regularization parameters $\gamma_{J}$ and $\gamma_{h}$ are related to the prior variances of fields ($\sigma^{2}_{h}$) and couplings ($\sigma^{2}_{J}$) through $\gamma_{h}=1/(B\sigma^{2}_{h})$, and $\gamma_{J}=1/(B\sigma^{2}_{J})$, where $B$ is the number of configurations in the sample. When the penalties are relatively weak ($\gamma\sim\mathcal{O}(1/B)$), this regularization can be thought of as a weakly informative prior gelman2008weakly whose main purpose is to prevent pathologies in the inference.

II.3 Gauge invariance

The $N\cdot q$ frequencies $f_{i}(a)$ and $\frac{1}{2}N(N-1)q^{2}$ correlations $f_{ij}(a,b)$, with $i<j$, are related to each other: the former sum up to 1, while the latter have the frequencies as marginals. Therefore, not all constraints in Eq. 3 are independent and multiple sets of parameters give the same probability distribution. In the language of physics this over-parameterization of the model is referred to as gauge invariance and the choice of one particular parameter set among the equivalent ones as gauge choice. This gauge invariance reduces the number of free parameters in the Potts model to $q-1$ fields for each site and $(q-1)^{2}$ couplings for each pair of sites.

In particular, we can reparameterize the model without changing the probabilities by an arbitrary transformation of the form:

[TABLE]

for any $K_{ij}(a)$, $H_{i}$ and $\kappa_{ij}$. In the so called lattice-gas gauge, this freedom is used to define a gauge state $c_{i}$ at each site such that

[TABLE]

for all states $a,b$ and sites $i,j$. The couplings and fields are transformed as follows:

[TABLE]

Two common gauge states are the most and the least frequent states of each site, defining respectively the consensus gauge and the least-frequent gauge. In protein analysis, the gauge state is often fixed to the amino acid present at site $i$ in a reference sequence, called wild-type sequence. An alternative choice is the so-called zero-sum gauge, in which

[TABLE]

for all states $a$ and all variables $i,j$. Fields and couplings can also be put in the so-called zero-sum gauge through

[TABLE]

where $g(\cdot)$ denotes the uniform average of $g(a)$ over all states $a$ at fixed position.

Note that, while all observables such as the moments of the distribution are invariant with respect to the gauge choice, the fields and the couplings are not. Arbitrary functions of the couplings and fields, such as the commonly-used Frobenius norm of the couplings, are also not generally gauge invariant. If not explicitly stated, the comparisons shown in this paper are performed in the consensus gauge, but the choice of the gauge for the inference and for the analysis of the inferred network can be different. The gauge chosen during the inference will be further discussed in section II.4 in the descriptions of ACE and PLM.

II.4 Algorithms

The presence of the partition function $Z$ in Eq. 4 precludes direct numerical minimization of the cross-entropy when the system size is large, since this requires summing over all $\prod_{i=1}^{N}q_{i}$ possible configurations of the system. However many approximate solutions have been proposed to tackle this issue. We briefly recall two of these methods to respectively approximate the cross-entropy or the log-likelihood: the Adaptive Cluster Expansion (ACE) and PseudoLikelihood Maximization (PLM).

II.4.1 Adaptive cluster expansion (ACE)

The cross-entropy (Eq. 4) can be exactly decomposed as a sum of cross-entropy contributions, calculated recursively (see Appendix VIII.1 and VIII.2). The adaptive cluster expansion cocco2011adaptive ; cocco2012adaptive ; barton2016adaptive is based on the idea of summing up cluster contributions based on their importance as quantified by their absolute contribution to the cross entropy. To this end an inclusion threshold parameter $t$ is introduced and only clusters with cross-entropy contributions larger than the threshold $t$ are included. The inclusion threshold $t$ is then progressively decreased to include more and more clusters in the summation. The expansion is usually stopped when the frequencies and correlations of the inferred model reproduce the empirical ones to within the statistical error bars due to finite sampling. The inference routine which has been used in this paper is publicly available at https://github.com/johnbarton/ACE . For an input sample of size $B$, the regularization parameters are set to $\gamma_{J}=1/B$ and $\gamma_{h}=0.01/B$, corresponding to a variance of the prior distribution of couplings of order 1 and a variance of fields of order 100.

II.4.2 Pseudo-likelihood maximization (PLM)

The idea behind Pseudo-Likelihood Maximization is to approximate the full likelihood of the data given the model (or equivalently the full cross-entropy (Eq. 4) by the site-by-site maximization of the conditional probability of observing one state at a site, given the observed states on the other sites. This approximation makes the problem tractable, and it also makes possible to parallelize the computation for the different sites. Pseudolikelihood is a consistent estimator of the likelihood in the limit of infinite input data. For this study, a version of the asymmetric pseudolikelihood maximization ekeberg2013improved ; ekeberg2014fast capable of working with a site-dependent number of Potts states has been implemented adapting the public code by M. Ekenberg and E. Aurell at https://github.com/magnusekeberg/plmDCA.

Unlike ACE, the networks inferred by PLM with $L_{2}$-regularization are always fully connected. As empirically shown in protein sequence analysis ekeberg2013improved ; hopf2017mutation ; cocco2018inverse and in theoretical analyses barton2014large , large regularization is needed in the presence of fully connected networks to avoid overfitting and thus to improve contact and fitness predictions. We have tested different regularization strengths, see Sec. 3, and fixed $\gamma_{J}=N/B$, $\gamma_{h}=0.1/B$ for system with $N$ variables and input sampling of size $B$.

With PLM gauge invariance is automatically broken. The inference is performed in the gauge that minimizes the $L_{2}$-regularization:

[TABLE]

As for ACE, the PLM fields and couplings are subsequently transformed to the consensus gauge for comparison.

III Regularizations

III.1 Removing variable states

III.1.1 Color compression

So far we have described (Eq. 4 and 5) how to infer the parameters of a Potts model where the number of states $q$ is the same at all sites, but it is easy to generalize this procedure to Potts models in which the number of states $q_{i}$ depends on the site $i$. This situation naturally arises due to sampling: states with very small probabilities are rarely observed. For instance, in multiple sequence alignments of protein families, for the large majority of sites. only a subset of the full $q=21$ possible amino acids are observed. In the color compression procedure, for each site $i$, we model explicitly only the $q_{i}$ states observed with a frequency $f_{i}(a)$ larger than a frequency cutoff threshold $f_{0}$

[TABLE]

and we group together the remaining $q-q_{i}$ low frequency states into a single one. The frequency of the grouped/compressed Potts state $q_{i}+1$ is then the total frequency of the states that have been grouped together: $f_{i}(q_{i}+1)\equiv\sum_{a^{\prime}=q_{i}+1}^{q}f_{i}(a^{\prime})$.

III.1.2 Artificial data sets on Erdős-Rényi random graphs

To benchmark color compressed inference we have generated synthetic data from a Potts model with $N=50$ variables carrying $q=10$ Potts states and interacting on an Erdős-Rényi random graph. To build the interaction graph, each edge in the network is included with probability $0.05$, giving a mean connectivity 2.5, with a maximum connectivity equal to 7. An exemple of contact map is shown in Supplementary Fig. 16. The field and couplings parametes on interacting sites are selected from Gaussian distributions of mean $\mu=0$ and standard deviations $\sigma_{J}^{2}=1$ and $\sigma_{h}^{2}=5.$ Therefore if $i$ and $j$ interact $J_{ij}$ is a $10\times 10$ matrix whose elements are chosen independently according to the above distributions, and each element of the matrix is zero when the sites do not interact.

We have generated 10 independent realizations of such ER models (networks of interactions and sets of fields and couplings). For each realization, $B=5\;10^{2},10^{3},10^{4}$, or $10^{5}$ configurations are generated by Markov Chain Monte-Carlo sampling. The number of available data, $B\times N$, can be compared to the number of parameters to be inferred, $M_{tot}=qN+q^{2}N(N-1)/2\simeq 1.2\;10^{5}$. We have, for the previously listed values of $B$, $B\times N=2.5\;10^{4},5\;10^{4},5\;10^{5}$, or $5\;10^{6}$, the first two being in heavy undersampling conditions, the third in scarce sampling, and only the last one being relatively well sampled.

III.1.3 Reduction in the number of parameters in compressed data: comparison between ER data with $q=10$ and sequence data with $q=21$

Tanks to color compression the number of colors per site, $q_{i}$, and thus the number of parameters $M_{cc}=\sum_{i}^{N}q_{i}+\sum_{i<j}^{N}q_{i}\,q_{j}$ are reduced with respect to the full number of colors per site, $q$, and the total number of parameters to be inferred $M_{tot}$. Fig. 1(left) illustrates how such reduction is achieved as a function of the compression threshold $f_{0}$ and for different sampling depths $B$ for the artificial data. Increasing $f_{0}$ the inferred model goes from the full-color $\langle q_{i}\rangle=q$ Potts model to the Ising model with only two colors $\langle q_{i}\rangle=2$ per sites indicating only the presence or absence of the most probable color on each site. The number of parameters to be inferred is divided by 50 at high compression.

Fig. 1 (right) shows a similar behavior for parameter reduction as a function of $f_{0}$ and sampling size $B$ in protein sequence data, for three families that will be further described and analyzed in Sec. V. In the latter case the number of Potts states is the number of possible amino-acid types (plus the gap symbol), $q=21$, and the number of parameters is decreased by 100 times at high compression. While protein sequence data are the direct application of color compression that we will consider in Sec. V, we have generated synthetic data from a Potts model with $q=10$ instead of $q=21$ states to speed up the numerical analysis performed varying sampling depth, compression threshold and inference methods; especially for the ACE algorithm, indeed, as it will be shown in Sec. VI, the computational cost of inferring a complete $q=21$ Potts model may become very large. Moreover the aim of the present study is to investigate the interplay between sampling and effect of compression, and we do not expect qualitative changes between $q=10$ and $q=21$: we expect that the compression threshold $f_{0}$, at which the performances of the inferred model worsen with respect to the full Potts model does not depend on the maximal numbers of Potts state in the original model but rather on their occurrence in the sampling. At large compression threshold, we expect that model performance deteriorates because highly frequent, well-sampled colors are grouped and their correlations are lost.

III.1.4 Color decompression

Once the restricted Potts model is inferred, we need to recover the complete model with $q$ states at each site in order to compare it to the ground truth. To this aim, we have to determine the parameters for the states that were grouped on site $i$. We use the following procedure: For each grouped state $a^{\prime}$, the fields and the couplings are estimated through

[TABLE]

where $q_{i}+1$ refers to the unique Potts state for the grouped symbols. This procedure allows us to correctly recover the local fields parameters reproducing the frequencies $f_{i}(a^{\prime})$ for the grouped symbols, while a common coupling parameters is assigned to all the grouped symbols. Then, we associate fields and couplings to states that are never observed in the sampling, hereafter referred to as unseen states, through a natural extension of the procedure described above. We assign a pseudocount frequency $f_{i}(a")=\alpha/B$ to the unseen symbols $(a")$. Herafter, we have fixed $\alpha=0.1$, in the expected range $0<\alpha<1$. When the grouped state is not present while unseens are present, we use the least probable state on the given site ($l_{i}$) instead of the grouped state in the denominator of Eq. 12, by replacing $q_{i}+1$ by $l_{i}.$

The choice of associating the unseen states to the grouped state or the least probable state is both simple and effective. Indeed, it follows the gauge choice, and yields fields with lower values than for the observed states. A detailed comparison of the effects of the pseudo-count above and of the standard $L_{2}$ regularization in the simple case of an independent model is presented in Appendix VIII.5.

III.1.5 Gauge used in the ACE inference

ACE inference has been implemented in the lattice-gas gauge (Eq. 6) barton2016adaptive . It is important to notice that the choice of the gauge symbol may have some effect on the performance of the inference procedure because the regularization term is not gauge invariant. For abundant data or in the limit of large compression, Potts states are well sampled and the choice of the gauge symbol is largely irrelevant. However, for few data or in the absence of color compression, or at small compression, the best performance is obtained by choosing as gauge symbol the least-probable Potts state on each site. In this way, all the fields and couplings corresponding to at least one poorly sampled state are fixed to zero, and have therefore null statistical variances by construction. Fields and couplings are then translated to the gauge in which the most probable symbol is chosen as gauge-symbol (consensus gauge) to perform the comparisons described in the next sections using Eq. 7. This consensus gauge is best for comparison because the statistics of the consensus symbols (on all sites) are the easiest ones to measure accurately.

III.2 Removing interactions

We also regularize the inference procedure by limiting the number of non-zero interactions. Sparsification of the interaction network is sometimes achieved through $L_{1}$ regularization of the couplings. Hereafter, we show that the inclusion threshold of the ACE inference procedure defined in Section II.4.1 plays a similar role, while not affecting the amplitude of non-zero couplings.

III.2.1 Role of ACE inclusion threshold: sparse versus dense inferred graphs

Fig. 2 shows the behavior of the ACE algorithm as a function of the inclusion threshold $t$ for one particular graph, hereafter called ER05, with $B=1000$ sampled configurations, analyzed with a color compression of $f_{0}=0.01$. This representative data set will be our reference case. For each threshold $t$ used to select clusters in the ACE expansion, the model frequencies $\langle\delta(a_{i},a)\rangle$ and $\langle\delta(a_{i},a)\delta(b_{i},b)\rangle$ calculated by Monte-Carlo simulation are compared to the data frequencies $f_{i}(a)$ and $f_{ij}(a,b)$ (see Eq. 3).

As detailed in cocco2012adaptive ; barton2016adaptive , to monitor the ability of the inferred model’s ability to reproduce the measured frequencies and correlations while avoiding overfitting, we define a relative error that is the ratio between the deviations of the predicted observables from the data, $\Delta f_{i}(a)=\langle\delta(a_{i},a)\rangle-f_{i}(a)$ and $\Delta f_{ij}(a,b)=\langle\delta(a_{i},a)\delta(b_{i},b)\rangle-f_{ij}(a,b)$, and the expected statistical fluctuations due to finite sampling, $\sigma_{i}(a)=\sqrt{f_{i}(a)(1-f_{i}(a))/B}$ and $\sigma_{ij}(a,b)=\sqrt{f_{ij}(a,b)(1-f_{ij}(a,b))/B}$. The relative error on frequencies is

[TABLE]

The relative error on connected correlations, $c_{ij}(a,b)=\langle\delta(a_{i},a)\delta(b_{i},b)\rangle-\langle\delta(a_{i},a)\rangle\langle\delta(b_{i},b)\rangle$, is

[TABLE]

where we estimate the standard deviation in the connected correlations as $\sigma^{c}_{i,j}(a,b)=\sigma_{ij}(a,b)+f_{j}(b)\sigma_{i}(a)+f_{i}(a)\sigma_{j}(b)$. Finally, the maximum relative error is

[TABLE]

where ${M}=Nq+(N(N-1)/2)q^{2}$ is the total number of one- and two-point correlations. As shown in Fig. 2 (top panel) the relative errors defined above have a nonmonotonic behavior as a function of the threshold, reaching relative minima that successfully reconstruct the data ($\epsilon_{\rm max}<5$) at multiple values (marked by asterisks) of the expansion threshold $t$, see Table 1. The regularized cross entropy, the total number of clusters included in the expansion, and their maximal size as a function of the cluster inclusion threshold $t$ are also shown Fig. 2.

The cluster inclusion threshold acts as an additional regularization. There are 3 plateaus in the regularized cross entropy as a function of the cluster inclusion threshold $t$ of Fig. 2: the first plateau corresponds to an independent model, the second one to a sparse interaction network, and the third one to a fully connected network. The number of edges present in the inferred graph of Fig. 2 is given by the number $N_{2}$ of 2-site clusters in the ACE expansion and is shown in Table 1 for the threshold corresponding to the minimal relative errors $\epsilon_{max}$. In particular the two relative minima better reproducing the data correspond to 2 different inferred networks. The minimum with $\epsilon_{\rm max}=3.9$ is at high threshold ($t=0.036$) and is characterized by a numbers of edges $N_{2}=55$ smaller than the total number $N\,(N-1)/2=1225$ of possible pairs. The inferred model is therefore a sparse graph, with a number of edges $N_{2}$ comparable with the number of edges of the model used to generate the data ($N_{0}=59$ for the model used to generate the data in Fig. 2). The second relative minimum with $\epsilon_{\rm max}=1$ is at low threshold $t=6.4\times 10^{-5}$ where the expansion includes the maximal number of 2 site clusters $N_{2}=N\times(N-1)/2$, corresponding to a fully connected graph. As can be guessed by the difference in the connectivity between the original and inferred model, and as we will better quantify in Section IV.1, the fully connected solution is overfitting the data.

III.2.2 spACE, a variant of ACE for sparse interaction networks inference

To force the ACE algorithm towards a sparse solution we introduced a new procedure in the cluster expansion (available at https://github.com/johnbarton/ACE),, which stops the algorithm at a maximal number $N_{2}^{max}$ of 2-site clusters, see Appendix VIII.2. This procedures imposes a prior knowledge on the sparsity of the interaction graph by giving an upper bound for the number of edges. The Erdős-Rényi random graph models used here to generate the data have an average connectivity of 2.5 neighbors per site, so we can use this prior knowledge to fix a maximal number of edges to $N_{2}^{max}=N\times 2.5\approx 125$. In practice to find the best sparse graph with a number of edges smaller than the prescribed value $N_{2}^{max}$ we spanned the threshold range in regular intervals 555The first interval is $I/\tau<t<1$, the second is for $1/\tau^{2}<t<1/\tau$ and so on. We have chosen here $\tau$ =3.4 (see command -trec on the GitHub site), at which the local minima of the absolute error and the corresponding t parameters are recorded (see Table 1) The parameters corresponding to the global minimum (under the sparsity conditions) are then chosen ($t=3.6\times 10^{-2}$ for Table 1). The spACE procedure is stable when changing $N_{2}^{max}$ around its prior fixed average value. In the following we have stopped the algorithm for $N_{2}^{max}=200$ (shown on Fig. 2), after checking that the value $N_{2}^{max}=100$ gave similar results. This procedure greatly reduces the computational time, which increases linearly with the number of computed clusters and grows exponentially, as $q^{K_{max}}$, with their maximal size $K_{max}$ (see Sec. VI and Appendix VIII.3): as shown in Fig. 2 $K_{max}=3$ at $N_{2}^{max}=200$ while $K_{max}=7$ for the fully connected graph. The number of inferred parameters $M_{cc}$ in Fig. 1 is further reduced by a factor $5\;10^{-2}$ for the best sparse model on ER data in Table 1. having couplings only on connected sites.

IV Benchmarking of color compression and decompression on synthetic data

To carry out an extensive analysis of the effects of the color compression introduced in Section III on the quality of the inference, we will apply it to artificial data generated by Potts models on Erdős-Rényi (ER) random graphs. The model and the generation of the data are described in Section III.1.2.

Once the data are obtained, we apply the compression schemes introduced in Section III with no color compression and with frequency cut-off $f_{0}=[0,10^{-4},10^{-3},3\cdot 10^{-3},10^{-2},3\cdot 10^{-2},10^{-1}]$. Note that all frequency thresholds $0<f_{0}\leqslant 1/B$ give the same color compression, so we infer the model only for the upper value in this range and thus the number of the tested frequency thresholds depends on $B$. Moreover, $f_{0}=0$ corresponds to removing from the inference only the unseen states.

Given the 10 realizations of the Erdős-Rényi model, the 4 sample sizes and the 5 to 8 (depending on the sampling) values of the frequency threshold define 280 data sets. For each of them, we have inferred the corresponding Potts parameters, both with the ACE and the PLM algorithms.

IV.1 Probability distributions

The Kullback-Leibler (KL) divergence measures how the inferred probability distribution of the possible configurations diverges from the empirical one (defined from the data samples), and can be computed as:

[TABLE]

where ${\mathbf{a}}=\{a_{1},\ldots a_{N}\}$ is a configuration and $\langle\cdot\rangle_{{\mathbf{a}}\>in\>real}$ indicates the average over the configurations generated by Markov Chain Monte Carlo (MCMC) from the real model. The first and the second lines are identical only when an infinite configuration sample is employed. Here, we estimate the average over $P_{\boldsymbol{J}^{real}}$ using an ensemble of 50,000 MCMC configurations sampled from the model.

As described before, the computation of the partition function $Z$ is far from being trivial, and was done in two ways. First, we used Annealed Importance Sampling (AIS) Salakhutdinov08 , starting from the independent-site model: All initial couplings were set to zero, while initial fields were computed as $h^{0}_{i}(a)=\log(f_{i}(a)+\alpha)-\log(f_{i}(cons_{i}))$ where $cons_{i}$ is the most common state at site $i$ and $\alpha=1/B$ is the smallest observed frequency used as regularization. Then a chain of models with increasing couplings (up to the inferred values) are thermalized and the ratios of their partition functions may be estimated. Secondly, the Kullback-Leibler divergence and the logarithm of the partition function can also be directly estimated by the ACE procedure (Table 1), see Appendix VIII.1. The KL divergences obtained directly from the ACE expansion and the ones obtained with importance sampling are very similar, as shown in Table 2, for the reference case in Fig. 2 at the optimal cluster inclusion threshold corresponding to a sparse inferred network. The values of the logarithm of the partition function, and of the entropy are also consistent between the two methods. In the following we will use the annealed importance sampling to calculate the KL divergence to compare results from PLM and ACE.

IV.1.1 KL Divergence for ACE models at different inclusion thresholds $t$

Table 2 displays the KL divergences for the reference data set and different cluster inclusion thresholds of Fig. 2 and in table 1 obtained both with importance sampling and the ACE expansion.

The fully-connected-graph model has a larger KL divergence and is therefore overfitting the data, while the sparse-graph model reproduces better the original model. All results for the ACE expansion presented in the following are obtained with the spACE procedure to infer a sparse graph. For the fully connected solution, due to overfitting, the cross entropy of Table 1 is not a good approximation to the entropy. Therefore, the estimate of the logarithm of the partition function and of the entropy given in table 2 are significantly different from the ones obtained by the AIS method.

IV.1.2 KL Divergence as a function of the sampling depth B and the color compression threshold $f_{0}$

Fig. 3 shows the mean over the ten ER realizations of the KL divergence between the real and the inferred distributions for various sampling depths and compression parameters for both ACE (left) and PLM (right). As expected, the KL divergence decreases for bigger samples, becoming very close to zero for $B=10^{5}$. The same happens for the standard deviations over the 10 realizations. ACE gives smaller KL divergences with respect to PLM, showing that the sparsity imposed in the spACE procedure gives a model reproducing better the original, sparse ER models.

The black empty squares in Fig. 3 indicate a reference compression threshold $f^{*}_{0}=10/B$, i.e. grouping symbols observed less than 10 times. We consider that frequencies larger than $f_{0}$ are reliably estimated 666We estimate the standard deviation $\sigma$ of the osberved frequency based on a Binomial variable with probability $\hat{f}$ in a sample of size $B$: $\sigma=\sqrt{({\hat{f}}\,(1-{\hat{f}})/B}\simeq\sqrt{{\hat{f}}/B}$. Probabilities larger than $f_{0}$ correspond to signal-to-noise ratios ${\hat{f}}/\sigma>\sqrt{f^{*}_{0}\,B}>3$. For $f_{0}>f^{*}_{0}$ the compression procedure discards information on ’reliable’ colors, and a loss in performance is expected. For both inference procedures the increase in KL divergence in Fig. 3 becomes visible only at a large value of the cut frequency $f_{0}\gtrsim 0.001$, independently from the sampling depth. We indeed expect that the increase of the KL divergence due to the color compression procedure be a monotonic function of the cut frequency $f_{0}$, so it is irrilevant at small cut frequency. For high values of $f_{0}$ the increase of the KL divergence with respect to the uncompressed model is of course much more significant at high $B$.

Table 2 shows that the strong regularization ($\gamma_{J}\approx N/B$) used in the fully connected PLM inference is essential to reduce overfitting: a PLM model inferred with a weak regularization ($\gamma_{J}\approx 1/B$) gives indeed very large KL divergences, especially at low sampling depth. Color compression, acting as an additional regularization, helps in this latter case to reduce overfitting and to lower the KL divergence. As shown in Appendix VIII.5.1 and Fig.15 the KL divergences for weak regularized PLM model reach indeed a minimal value, at intermediates or large compression thresholds $f_{0}$ depending on sampling depth, compatible, but slightly larger, than the one obtained for the usual strongly regularized PLM of Table 2.

IV.2 Low-order Statistics

In this Section we discuss the generative properties of the inferred models, in particular its ability to reproduce the low order statistics of the original model : the site frequencies $f_{i}(a)=\displaystyle{\sum_{generated\>{\mathbf{a}}}\delta(a_{i},a)/B_{gen}}$ and covariances $cov_{ij}(a,b)=\displaystyle{\sum_{generated\>{\mathbf{a}}}\left[\delta(a_{i},a)\delta(a_{j},b)/B_{gen}\right]-f_{i}(a)f_{j}(b)}$. To benchmark the generative power of the inferred model as a function of the color compression two sets of 20000 configurations are generated by Markov Chain Monte Carlo, respectively with the real and with the inferred model for each $B$, $f_{0}$, and graph realization, and their low order statistics are compared.

Figure 4 shows, for the models inferred from the reference data set, the comparison of the frequencies and covariances computed from the configurations generated by the real model (test sequences) and by the models inferred with ACE and PLM, without color compression (top panels) and with $f_{0}=0.01$ (bottom panels). As shown in the figures PLM covariances are downscaled, due to the strong regularization, as happens for the couplings (Section IV.4). Moreover PLM assigns smaller frequencies to the unseen Potts states (left panels of Fig. 4 in log-log scale). This is probably due to the fact that the pseudocount used during decompression seems to be well fixed for the weak regularization used in ACE but not for the large regularization used in PLM. The insets in Fig 4 show that, contrary to spACE, zero covariances are set to non-zero values with PLM due to overfitting.

To have a more systematic comparison, we analyzed their absolute mean square error defined as:

[TABLE]

These quantities are then computed for different $B$ and $f_{0}$ and averaged over 10 graph realizations.

Figure 5 shows the Mean Square Errors (MSE) for the frequencies and covariances. spACE has again better performances than PLM 777Similar results are obtained when considering the Pearson correlations between the real and inferred frequencies or correlations rather than the absolute errors (not shown).. The MSE for the full Potts model (full lines) shows a little dependence on the compression frequency $f_{0}$; meanwhile the MSE restricted to explicitly modeled colors (dotted lines) generally increases at large compression performances. The only exception is that frequencies are better reproduced with PLM at strong color compression because the decompression procedure (Eq. 12), acting as an independent model, correctly assigns fields to grouped states. The above observations are simply explained by the fact that the contribution to the MSE of a color increases with its frequency. The MSEs on the full Potts model, after decompression (full lines in Fig. 5 ) are quite independent from $f_{0}$, as the averages in the MSE are dominated by the large amount of colors with small frequencies, see Fig. 4. On the contrary the MSE on explicitly modeled colors (dot lines in Fig. 5) shows a large increase with the compression threshold because the mean is restricted to more and more frequent colors.

IV.3 Interaction networks

In this section we focus on the reconstruction of the interaction network and the prediction of pairs of sites that are connected, or “in contact”, in the interaction graph. The original ER graph is sparse, with an average connectivity of about 2.5. We can predict contacting sites as those site pairs with large couplings, as traditionally done for protein structures sulkowska2012genomics ; nugent2012accurate ; hopf2012three . To this end, we compute the Frobenius norm of the (10$\times$10) inferred and decompressed coupling matrix between each pair of sites $i,j$, $F_{ij}=\sqrt{\sum_{a,b}J_{ij}(a,b)^{2}}\,.$

The heat map of the Frobenius norms of the couplings inferred by ACE and PLM at different color compression gives very similar results and allow us to identify the largely coupled sites by ACE and PLM (see Figs. 16 in Appendix). The only difference between ACE and PLM is that spACE infer a sparse network with zero Frobenius norms on the majority of links, while PLM with the $L_{2}$ norm regularization described in Eq. 5 infers a dense Frobenius norm matrix, There is therefore no straightforward separation between pairs of sites predicted to be in interaction or not.

To gain more insight into these predictions, as done for protein structure morcos2011direct ; ekeberg2014fast ; cocco2018inverse , we can sort links by decreasing Frobenius norm and follow the precision obtained progressively including the corresponding links in the so called Positive Predicted Value (PPV) curve (see Appendix  VIII.7).This is shown in Figs. 6 for a number of links up to the last inferred one for ACE or to $N_{0}=50$ for PLM, for the reference case.

The Frobenius norm is computed in the consensus gauge, which will be used to compare the couplings and fields in the next section and in the zero sum gauge used on protein structure prediction dunn2008mutual ; morcos2014coevolutionary ; cocco2018inverse . Performance can be further improved with the Average Product Correction (APC) (defined in Appendix  VIII.7) when using zero sum gauge. When inferring sparse networks, APC only corrects the ranking of the predictions in the PPV curve, but it does not change the overall number of site pairs predicted to be in interaction, nor the global precision, on the contrary APC on the zero sum gauge largely improve PLM precision, as expected. Figure 6 (bottom) shows the average PPV, over all the data realization and as a function of $f_{0}$, at the number of real links $PPV(N_{0}=50)$ for PLM and at the lesser between the predicted $N_{pred}$ or real $N_{0}$ number of links with ACE.

As expected, the plots show that the contact prediction improves with sampling, and that the APC significantly improves the results for PLM in zero sum gauge. PLM generally gives higher PPV, especially at high sampling depth $B=10^{4}$ and $B=10^{5}$ where the reconstruction error for spACE are due to the sparsity threshold. We have verified that, at large sample and when inferring fully connected netwoeks ACE has the same PPV as PLM. We see that, for both algorithms, the performance is usually stable against the introduction of color compression.

IV.4 Couplings and Fields

In Fig. 7 we compare the fields and the couplings of the real model (x-axis) and the inferred Potts model (y-axis) obtained for the reference data of the graph ER05 sampled at $B=1000$ without color compression (top panels) and with $f_{0}=0.01=10/B$ (bottom panels), respectively, with ACE and PLM. Different colors in Fig. 7 show Potts states (or Potts states pairs) occurring at different frequencies and therefore treated in the color compression procedure as explicitly modeled, grouped, or unseen in the configuration sample. For couplings, if at least one of the two Potts state is unseen, the pair is considered as unseen; if at least one site is grouped, the pair is considered as grouped; if both sites are explicitly modeled, the pair is considered as explicitly modeled. The comparison is performed in the consensus gauge.

Figure 7 shows that, as observed in Section  IV.3 the sparse procedure, spACE, misses some edges and the corresponding couplings are fixed to zero. PLM couplings are systematically smaller in amplitude than real ones, ending up in a tilted entry-by-entry comparison. This is due to the large regularization introduced to avoid overfitting.

Figure 8 shows the absolute errors on the fields and couplings defined through

[TABLE]

These errors measure the average distances from the diagonal of the points in the scatter plots of Fig. 7. We observe that the couplings and fields reconstruction performances are stable as a function of the color compression up to the reference value $f^{*}_{0}\simeq{10\over B}$, where they drop because the compression become too strong. Such drop is perfectly clear for the couplings parameters, especially within the small regularization used in spACE. Large compression threshold, as well as large coupling regularizations, degrades indeed the informations on the correlations of grouped, badly sampled, variables. The increase of field reconstruction error at all compressions is negligible for small sampling depth $B=500,B=1000$ and for both algorithm, showing that the decompression procedure introduced in Sec. III, correctly assign, in the limit of the available information, the large and negative fields for the grouped and unseen Potts symbols, as done by using prior information with the $L_{2}$ regularization and shown in Fig. 7. The dashed lines in Fig. 8, in agreement with Fig. 7, show that by restricting the coupling comparison to the explicitly modeled symbols, the better and better sampled ones at larger and larger $f_{0}$, the reconstruction indicators are better and better. In other words, even in the largely undersampled regime, parameters for well sampled colors are correctly inferred and are not affected by the poorly sampled states. This underlines the difference between sites and states in a Potts model. In the standard renormalization procedure cardy1996scaling when the number of sites are reduced in an effective “renormalized” Potts model the parameter values of the retained sites change. In contrast, in the space of Potts states, grouping some of them, and keeping the probabilities conserved, does not affect the others 888Similar results are obtained when considering the Pearson correlations between real and inferred parameters, rather than their absolute differences. ACE gives better results than PLM for parameter reconstruction, especially on couplings (not shown). This is due to the fact that spACE avoid overfitting of data and setting many non-zero couplings for non interacting sites in the real interaction graph..

V Application to sequence-based models of protein families

We now apply our inference approach to protein sequence data. Input samples are multiple sequence alignments of protein families, the nodes of the graph are the protein sites, and states are the 20 amino acids plus the insertion-deletion symbol ($q=21$). In this context, we aim at reconstructing the contact map morcos2011direct ; hopf2012three or the fitness landscape ferguson2013translating ; mann2014fitness ; figliuzzi2016coevolutionary ; hopf2017mutation . In particular, we would like to compare the change of energy corresponding to single point mutations with respect to a wild-type protein sequence to the experimentally measured changes of fitness of the protein. Differently from ER data, sequences are not uniformly sampled. A reweighting procedure, described in Appendix VIII.8, is usually introduced to reduce the initial number $B$ of sequences in the sample to an effective number of independent ones $B_{eff}.$

We here consider three protein families, whose fitnesses have been systematically assessed against single-point mutations. We specify for each of them the number of sites in the alignment $N$, the number of sequences $B$ and the number of effective sequences $B_{eff}$:

• •

WW ( $N=31$, $B=8251$, $B_{eff}=2590$) is a protein domain that mediates specific interactions with protein ligands. Here fitness has been measured in terms of the capability to bind a certain ligand araya2012fundamental .

• •

PDZ ($N=84$, $B=24795$, $B_{eff}=4240$) is a protein domain present in signaling proteins. Here fitness has been measured in terms of binding affinity mclaughlin2012spatial ;

• •

RRM ($N=82$, $B=70780$, $B_{eff}=18800$)is an RNA recognition motif; fitness was estimated through growth rate measurements in melamed2013deep .

Alignments and experimental fitness measures used in this section were taken from hopf2017mutation . The PLM procedure was applied with the same large regularization used for the artificial data, $\gamma_{J}=N/B$, $\gamma_{h}=0.1/B$. The SpACE procedure was applied with a threshold cutoff $N_{2}^{max}=3N$, and a smaller regularization $\gamma_{J}=10/B$, $\gamma_{h}=0.1/B$ ; however the relative error $\epsilon_{max}$ (Eq.15) was generally too large even in its local minima, indicating that the procedure had not converged. As shown in barton2016adaptive a Boltzmann Machine Learning (BML) procedure was further used, starting from the spACE inferred parameters as an initial guess, to better reproduce the low order statistics of the data and therefore the quality of the inferred model. Mutational cost of single mutations is predicted by their energetic cost, corresponding to the value of the field of the mutated amino acid, after having transformed the field parameters in the gauge of the wild type sequence through Eq. 7 cocco2018inverse .

Contrary to what happens for synthetic data, where the true model is known, the relationship between inferred energies and experimental fitness values may be nonlinear, so we use as a quality measure of the inference the Spearman correlation coefficient between them rather than the Pearson. The Spearman values of Fig. 9 for the full model (no compression) and PLM are in agreement to the ones previously obtained hopf2017mutation with PLM algorithm and a slightly different procedure, giving Spearman values of 0.6, 0.57, 0.5 for RRM1, WW and PDZ respectively. ACE + BML performances are compatible with PLM results, better for families with a large number of sequences such as RRM, and slightly worse in the PDZ case. The impact of color compression on the prediction is then investigated and shown in Fig. 9 and has to be compared to the quality of field reconstruction in Fig. 8, after decompression (full lines). For the WW and PDZ families, as for PLM with the ER data with $B=1000$ and $B=500$, the predictions are not affected by large compression thresholds. Very large compressions may even improve performances, as for the to WW family, where inference with almost binary Potts variables (at $f_{0}=0.1$ $<q_{i}>=2.5$ see Fig. 1) give better predictions than with a full $q=21$-state Potts model. A similar effect was observed also for HIV fitness predictions in mann2014fitness . The fact that large color-compressions can improve fitness predictions may be related to the experimental uncertainties and limited resolution on the fitness measures and also to the difficulty in estimating the effective number of independent sequences in the data. For the better sampled $RRM$ family, predictions worsen at large compression $f_{0}>10/B_{eff}$, as expected and observed in Fig. 8 for the well sampled artificial data (at $B=10^{4}$ and $B=10^{5}$).

VI Gain in computational time for synthetic and protein sequence data

In this section we study how the computational time scales with the sample size and the color compression frequency threshold for spACE and PLM on synthetic and protein sequence data. Times have been obtained on a processor Intel® Xeon(R) CPU E5-2690 v4 $@$ 2.60GHz x 56. The limiting factors for computational time with the ACE expansion are three: the number of colors $<q>$ in the inferred model, due to the computation of the partition function which grows as $<q>^{K}$ on clusters of size $K$, the overall number of clusters which enters in the construction rule to and the number of configurations sampled by Monte-Carlo to calculate the relative errors in the reconstruction of the statistics of the data.

For PLM the average over the sampled configurations for the pseudo-likelihood and moments calculation determines a linear dependence on the sample size $B$, and on the number of parameters. Fig. 10 and Table 4 show the computational times for the reference ER graph realization as a function of $f_{0}$ and B. On such data spACE (Fig. 10) weakly depends on $K$ and $\langle q\rangle$,which take small values, the limiting step is therefore the number of Monte Carlo sampled configurations. Having a large number of MC steps (here fixed to 500000) is important to correctly estimate the relative errors at very large sample size, given the small value of the sampling variances. The computational time show therefore a weak linear dependence on the compression, as shown in Fig. 13 (Appendix). The computational times for spACE and PLM on real proteins are compared in Tables  5 and in Fig. 11. As shown in Table 5 for real proteins a larger numbers and larger sizes of clusters are summed on in the ACE expansion, as compared to the artificial data in Table 4. The total number of processed clusters ($\hat{N}_{cl}$) and their average size $\langle q\rangle$ are clearly the limiting factor for the running time, making ACE slower than PLM. PLM further benefits from parallel computing, here results obtained on 1 and 20 are compared, while for spACE the computation has also been followed by time consuming MC-learning. Fig. 11 shows the large reduction in computational time within spACE by color compression, such compression makes often feasible the numerical computation, especially for non sparse ACE inference (see Fig. 13) which would take prohibitively long times without significant compression.

VII Conclusion

The present work reports an extensive numerical benchmarking of inferred color compressed Potts model with 2 algorithms (ACE and PLM ) on synthetic data, as well as on protein sequence data. Knowing the ground-truth model that generated the data, we can assess the inference performance at different compression strengths, by (1) computing the Kullback-Leibler divergence between the real and inferred models; (2) checking the reconstruction of low-order statistics; (3) testing the reconstruction of the structure of the interaction network, and of the couplings and field parameters.

We will in the following resume and discuss the comparison between PLM and ACE inference methods before discussing the results on color compression.

A first important advantage shown by ACE with respect to PLM inference is that ACE can be easily stopped at large value of the cluster inclusion threshold to reconstruct a good sparse model for the artificial data, which have been generated from sparse random graphs.Imposing sparsity regularization naturally reduces the number of coupling parameters in the inferred model, leaving the possibility to choose a small value for the $L_{2}$ regularization parameter on the left non-zero couplings. Such double regularization scheme allows a very good reconstruction of the overall model, the model parameters and the statistics. On the contrary for PLM , the inferred model is fully connected and therefore a large $L_{2}$ regularization is needed to avoid overfitting, in the optimal model reconstruction, with the consequences that all the amplitudes of coupling parameters are systematically underestimated, and the statistics of the data less well reproduced. Further theoretical investigations are needed to obtain the optimal regularization value found for the inference of a fully connected model, as a function of the sparsity of the original graph, the number of sites and color and the number of data and will be carried out on a forthcoming paper fantome2019 . It would be also interesting to introduce a double regularization scheme for PLM algorithm. The situation is different on protein sequence data because the spACE procedure with the simple stopping criterium implemented here does not converge to small statistical reconstruction errors, without additional MC learning algorithm. The resulting inferred network is no more sparse and no gain of performances of ACE with respect to PLM are achieved at larger computational costs. spACE procedure can be improved to better impose sparsity constraint in a more general case franz2019fast ; cocco2019 . A question which deserves further investigation is if the inference of sparse connectivity graph is more appropriate for parameter reconstruction in the large undersampling regime, even for data generated with models, such as protein sequence data, which are not necessary sparse. Preliminary results for fitness prediction on protein sequence data seems to indicate that sparse models gives very good performances cocco2019 . Next we will resume conclusion and discussion on color compression.

The central finding of the present work is that color compression does not degrade the studied performances in a very large range of frequency compression cut off, while largely reducing the dimensionality of the inferred model and its computational time. The reduction of computational time obtained thanks to color compression often become essential for solving the inverse problem in reasonable times. For example, ACE inference of models for many HIV proteins barton2015scaling ; barton2016relative ; barton2016entropy ; louie2018fitness has been possible thanks to color compression. The color compressed version of the PLM algorithm introduced adapted here from the routine of the group of Aurell and collaborators ekeberg2013improved ; ekeberg2014fast can be crucially important when dealing with large protein sequences or with whole genome inference gao2018correlation ; pensar2019high with a much larger number of variables than single protein sequence data. The color compression and decompression procedures introduced here are not restricted to pairwise graphical models. They could be used in other machine learning approaches, on protein sequence data, such as restricted Botzmann machines tubiana2019learning , or variational auto-encoders riesselman2018deep .

Acknowledgements. We thank Lorenzo Posani for useful discussions. This work benefited from the financial support of the grants PSL ProTheoMicS and RBMPro ANR-17-CE30-0021-01.

VIII Appendix

VIII.1 Reminder about Adaptive Cluster Expansion and the inclusion threshold

In the ACE inference procedure the cross entropy is expanded as the sum of cluster contributions. Defining a cluster as a sub-set of variables: $\Gamma=\{i_{1},\ldots,i_{k}\},k\leq N$, we can formally write the cross-entropy as the sum of cluster contributions:

[TABLE]

where the sum is over all nonempty clusters of the $N$ variables. The cluster cross-entropy contributions $\Delta S_{\Gamma}$ are recursively defined through

[TABLE]

Here $S_{\Gamma}$ denotes the minimum of the cross entropy (4) restricted only to the variables in $\Gamma$. Thus, $S_{\Gamma}$ depends only on the frequencies $p_{i}(a)$, $p_{ij}(a,b)$ with $i,j\in\Gamma$. Provided that the number of variables in $\Gamma$ is small (typically $\lesssim 10$ for $q=10$ Potts state as in the present work) numerical maximization of the likelihood restricted to $\Gamma$ is tractable. The definition of $\Delta S_{\Gamma}$ ensures that the sum over all clusters $\Gamma$ in (18) yields the cross entropy for the entire system of $N$ variables. As detailed in cocco2012adaptive ; barton2016adaptive , a recursive construction rule is used to avoid, before selection, the computation of all cluster entropies. Such rule consists in building up clusters of size $k$ by combining selected clusters selected of size $k-1$. The ACE expansion consists in truncating the expansion in Eq. (18) by fixing a cluster inclusion threshold $t$ and summing up in Eq.(19) only cluster contribution with $|\Delta S_{\Gamma}|>t$.

VIII.2 ACE and SpACE pseudocodes

The ACE algorithm, described in detail in barton2016adaptive , is based on the selection and summation of individual cluster contributions to the Cross Entropy. It is built by the recursion of the following routine. Given a list $L_{k}$ of clusters of size $k$, beginning with the list of all the $N\;(N-1)/2$ possible clusters of size $k=2$ :

For each cluster $\Gamma\in L_{k}$

1. (a)

Compute $S_{\Gamma}$ by numerical minimization of Eq. (4) restricted to $\Gamma$. 2. (b)

Record the parameters minimizing Eq. (4), called ${\bf J}_{\Gamma}$. 3. (c)

Compute $\Delta S_{\Gamma}$ using Eq. (19). 2. 2.

Add all clusters $\Gamma\in L_{k}$ with $|\Delta S_{\Gamma}|>t$ to a new list $L^{\prime}_{k}(t)$. 3. 3.

Construct a list $L_{k+1}$ of clusters of size $k+1$ from overlapping clusters in $L^{\prime}_{k}(t)$.

The rule used by default for constructing new clusters of size $k+1$ from selected clusters of size $k$ is the so called strict rule: a new cluster is added only if all of its $k+1$ subclusters of size $k$ belong to $L^{\prime}_{k}(t)$ . The above process is then repeated until no new clusters can be constructed (for ACE) or also until the maximal number of 2-site clusters $N_{2}^{max}$ is in the selected list $L^{\prime}_{k}(t)$ (for SpACE ). After the summation of clusters terminates, the approximate value of the parameters minimizing the cross-entropy, given the current value of the threshold, is computed by

[TABLE]

Then a Monte Carlo simulation is run to estimate the model one and two point correlations functions, which are compared to the empirical ones, taking into account their expected statistical fluctuations, by the computations of relative errors defined in Eq. (13, 14, 15) see also Fig. 2 in main text. When the maximal error is smaller than one the algorithm stops. Within the SpACE approximation even if a relative error of 1 is not reached, the smaller relative errors among the ones obtained for logarithmically spaced threshold intervals is chosen.

VIII.3 Cluster expansion and computational time as a function of the color compression $f_{0}$ for fully connected graphs.

Fig. 12 shows that the behavior of the maximal relative reconstruction error $\epsilon_{max}$ as a function of the cluster inclusion threshold $t$ when changing the color compression threshold $f_{0}$. The presence of two relative minima corresponding to a sparse and a fully connected models fitting the data is observed for all the values of $f_{0}$, see the two stars in Fig. 12. Moreover the threshold $t$ corresponding to the sparse inferred graph is largely independent of the level of color compression.

Fig. 10 in the main text shows a mild computational gain as a function of the color compression when inferring a sparse interaction network (large-threshold minima). Such gain is generally huge when the expansion converges only at low threshold values and sums up clusters of larger and larger sizes $K$. The numerical computation of the cross entropy requires indeed the sums over a number of $q^{K}$ configurations for $K$ Potts variables with $q$ states each. To illustrate this effect in Fig. 13 we show the reduction in computational time when the cluster expansion is stopped at the small threshold value corresponding to a fully connected inferred graph. One can reach a 1000-fold computational time reduction with large color compressions. As shown in Fig. 13 the expansion was stopped to maximal relative error of order 10 at small threshold $t$ after 11 days while it took 50 minutes to infer a good quality, fully connected model for the maximal color compression $f_{0}=0.1$. Note that the computational time to reach the sparse good model shown in Fig. 10 is smaller due to the reduced number of clusters. For large interconnected models color compression can therefore be essential to reach convergence in a reasonable amount of time and infer a model that reproduces the statistics of the data.

VIII.4 Kullback-Leibler Divergence from the ACE expansion

The computation is done in the Ising case for the simplicity of the notations, the generalization to the Potts case being straightforward. We denote $\boldsymbol{J}^{B}=\{J_{ij}^{B},h_{i}^{B}\}$ the inferred parameters at sample size $B$, and $\boldsymbol{J}^{true}=\{J_{ij}^{true},h_{i}^{true}\}$ the true underlying model parameters. The inferred cross-entropy at sampling $B$ writes

[TABLE]

where the sum is over all possible configurations $\boldsymbol{\sigma}=\{\sigma_{1},...,\sigma_{N}\}$. The inferred probability distribution at finite sampling $B$ is

[TABLE]

The Kullback-Leibler (KL) divergence between the true and the inferred distributions writes

[TABLE]

However, Eqs. (21) & (22) give

[TABLE]

The KL divergence between the true and the inferred distributions then writes

[TABLE]

Moreover, a reasonable approximation is

[TABLE]

because the true underlying parameters are recovered by the inference method in the perfect sampling case: $P_{\boldsymbol{J}^{B\rightarrow\infty}}(\boldsymbol{\sigma})\rightarrow P_{\boldsymbol{J}^{true}}(\boldsymbol{\sigma})$. Therefore,

[TABLE]

where $\left\langle\cdot\right\rangle^{B}=\sum_{\boldsymbol{\sigma}}\cdot P_{\boldsymbol{J}^{B}}(\boldsymbol{\sigma})$, and $\left\langle\cdot\right\rangle^{\infty}=\sum_{\sigma}\cdot P_{\boldsymbol{J}^{B\rightarrow\infty}}(\boldsymbol{\sigma})\approx\sum_{\sigma}\cdot P_{\boldsymbol{J}^{true}}(\boldsymbol{\sigma})$. It naturally generalizes to the $q$-states Potts case:

[TABLE]

The artificial data are in a compressed representation (cf. Section III). The complete inferred parameters are recovered as explained in Eq. (12).

VIII.5 Assignment of fields to zero-frequency states after inference

In section III we have discussed the decompression method used in the paper. In particular, we have seen that a pseudo-count is associated to the unseen states to assign them a field with respect to the reference of the grouped/compressed state or the least probable state and, in principle, this is different to what implicitly done when the model is inferred without color compression. To better understand the difference between the two approaches let us consider a simplified example of an independent-site model. Without color compression, the fields are obtained as the minimum of:

[TABLE]

which, for the unseen colors in the gauge of $Z=\sum_{a=1}^{q}\;e^{h_{i}(a)}=1$, gives:

[TABLE]

where $Lw(y)$ is the Lambert function, solution of $xe^{x}=y.$ On the other hand, the field which we assign to these symbols during color decompression is, in the same gauge,

[TABLE]

where we set, as in the rest of the paper, $\alpha=0.1$. In this independent-site approximation there is then a shift between the two procedures given by $\Delta h=h_{i}(a)^{\Gamma}-h_{i}(a)^{p}$, that depends from the pseudocount $\alpha$ and from the value of the regularization $\gamma_{h}$. In table 6 we give these shifts for the two values of $\gamma_{h}$ used respectively by ACE ($\gamma_{h}=0.01/B$) and PLM ($\gamma_{h}=0.1/B$).

If, in the approximation of independent-sites, the difference $\Delta h$ is the same in all gauges, the specific values of $h_{i}(a)^{\Gamma}$ and $h_{i}(a)^{p}$ of table 6 are specific of the $Z=1$ gauge. To have a comparison in the consensus gauge as done in the rest of the paper one has to subtract $h_{i}(a)^{\Gamma}$ and $h_{i}(a)^{p}$ the field of the most common color $c_{i}$ on the considered site. In the independent-site approximation this is just $h_{i}(c_{i})=log(p_{i}(c_{i}))$, and makes that unseen colors of different sites are found to have different fields. In Fig. 14 we plot the fields for the unseen symbols with a color compression $f_{0}=0.01$ (green diamonds) and $f_{0}=0$ (blue squares) versus the one for no color-compression for PLM and ACE (same fields as in Fig. 7 of the main paper), in the consensus gauge. We can see a systematic shift towards lower values (at least for PLM, to be checked for ACE). We can compare this shift with the theoretical shift obtained with the independent model as described above. Even if we neglect the terms due to the couplings we can well reproduce such shift as the difference between the field obtained with the pseudo-count with respect to the one obtained with the regularization, there is a good agreement between what observed and the theoretical shift for independent variables. In particular the shift is smaller for the regularization chosen by the ACE procedure.

VIII.5.1 KL divergence for PLM at low $L_{2}$ regularization as a function of the sampling depth $B$ and of the color compression threshold $f_{0}$

In this section we study what happens with PLM at lower $L_{2}$ regularization, e.g. $\gamma_{J}=1/B$ as the one used for ACE, when the threshold for color compression is varied. Without color compression, the performance obtained for $\gamma_{J}=1/B$ becomes significantly worse, see Table 3.

Figure 15 shows the average KL divergence between the true model and the inferred one for several color compression frequencies at low $L_{2}$ regularization ($\gamma_{J}=1/B$). For all the sample sizes, we observe the existence of an optimal value of $f_{0}$. Especially at small sampling depth $B\leq 1000$, large color compression leads to a very significant decrease in the KL divergence. However, in spite of the improvement due to color compression, the KL divergences do not reach the minimal values obtained with strong $L_{2}$ regularizations (dashed lines in Fig. 15), showing that a good choice of the regularization is always essential.

VIII.6 Graphical reconstruction: Contact maps and $Fscore$ for spACE graph reconstruction as a function of color compression

In Figs. 16 we compare the real contact maps with the heat map of the Frobenius norms of the couplings inferred by ACE and PLM with no color compression and $f_{0}=0.01$ for the reference data (ER05, $B=1000$). Results at different color compressions are very similar, such as the identification of largely coupled sites by ACE and PLM.

To encompass both the precision and the recall in a single measure it is possible to use the Fscore, which is the harmonic mean between the two and gives

[TABLE]

where $TP$ is the number of true predicted contacts, $N_{pred}$ is the number of predicted contacts and $N_{0}$ is the number real contacts. The Fscore for ACE inference are plotted in Fig.17.

VIII.7 Definition of Positive Predicted Value and Average Product Correction

The Positive Predicted Value (PPV) curve is defined as:

[TABLE]

Where $TP(n)$ is the number of true predicted edges in the top $n$ pairs. the Average Product Correction (APC) dunn2008mutual ; morcos2014coevolutionary ; cocco2018inverse ,

[TABLE]

Here the dot indicates the average over the corresponding variables, e.g. $F_{i,.}$ is the average of $F_{ij}$ over the second index $j$.

VIII.8 Reweighting procedure

To reduce sampling bias, we decrease the statistical weight of sequences having many similar ones in the MSA. More precisely, the weight of each sequence is defined as the inverse number of sequences within Hamming distance $d_{H}<xL$, with an arbitrary but fixed $x\in(0,1)$:

[TABLE]

for all $m=1,...,M$. The weight equals one for isolated sequences, and becomes smaller the denser the sampling around a sequence is. Note that $x=0$ would account to removing double counts from the MSA. The total weight

[TABLE]

can be interpreted as the effective number of independent sequences.