Francesca Rizzato

Alice Coucke

Eleonora De Leonardis

John P Barton

Jérôme Tubiana

Rémi Monasson

Simona Cocco

Inference of compressed Potts graphical models

We consider the problem of inferring a graphical Potts model on a population of 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 Potts states (colors) available to each variable is reduced, and interaction networks are made sparse. To achieve the color compression, 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, Adaptive Cluster Expansion (ACE) and PseudoLikelihood Maximization (PLM) on synthetic data obtained by sampling disordered Potts models on an Erdős-Rényi 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.

I. INTRODUCTION

.

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

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 ≤ q depends on the site i.

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 = N i q i + N i<j 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 [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF] or to their mutual information [START_REF] Haldane | Coevolutionary landscape of kinase family proteins: sequence probabilities and functional motifs[END_REF] or compression to a fixed number of colors [START_REF] Anton | Radi (reduced alphabet direct information): Improving execution time for direct-coupling analysis[END_REF], 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 [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF] 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

A. 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 a = a 1 , . . . a N ,

E(a|J) = - N i=1 h i (a i ) - N -1 i=1 N j=i+1 J ij (a i , a j ) (1)
and, consequently, a probability

P (a|J) = exp (-E(a|J)) Z(J) , (2)
where Z(J) = a exp (-E(a|J)) is the partition function and ensures that all probabilities sum to one. For simplicity, here we label the set of fields and couplings as 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 [START_REF] Jaynes | On the rationale of maximum-entropy methods[END_REF], 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.,

where δ(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.

B. 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

S(J|f ) = log Z(J) - N i=1 q a=1 h i (a)f s i (a) - N -1 i=1 N j=i+1 q a=1 q b=1 J ij (a, b)f s ij (a, b) , (4)
where, for simplicity, we indicate the set of single and pairwise frequencies as f and the set of fields and couplings as J.

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

∆S = γ h N i=1 q a=1 h i (a) 2 + γ J N -1 i=1 N j=i+1 q a=1 q b=1 J ij (a, b) 2 .
(

) 5
The regularization parameters γ J and γ h are related to the prior variances of fields (σ 2 h ) and couplings (σ 2 J ) through γ h = 1/(Bσ 2 h ), and γ J = 1/(Bσ 2 J ), where B is the number of configurations in the sample. When the penalties are relatively weak (γ ∼ O(1/B)), this regularization can be thought of as a weakly informative prior [START_REF] Gelman | A weakly informative default prior distribution for logistic and other regression models[END_REF] whose main purpose is to prevent pathologies in the inference.

C. Gauge invariance

The N • q frequencies f i (a) and 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:

h i (a) → h i (a) + H i + j( =i) K ij (a) J ij (a, b) → J ij (a, b) -K ij (a) -K ji (b) + κ ij
for any K ij (a), H i and κ ij . In the so called lattice-gas gauge, this freedom is used to define a gauge state c i at each site such that

J ij (a, c j ) = J ij (c i , b) = h i (c i ) = 0 , (6)
for all states a, b and sites i, j. The couplings and fields are transformed as follows:

h i (a) → h i (a) -h i (c i ) + j =i (J ij (a, c j ) -J ij (c i , c j )) , J ij (a, b) → J ij (a, b) -J ij (c i , b) -J ij (a, c j ) + J ij (c i , c j ) . (7)
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

q c=1 J ij (a, c) = q c=1 J ij (c, a) = q c=1 h i (c) = 0 , (8)
for all states a and all variables i, j. Fields and couplings can also be put in the so-called zero-sum gauge through

h i (a) → h i (a) -h i (•) + j( =i) [J ij (a, •) -J ij (•, •)] , J ij (a, b) → J ij (a, b) -J ij (•, b) -J ij (a, •) + J ij (•, •) , (9)
where g(•) 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 D in the descriptions of ACE and PLM.

D. Algorithms

The presence of the partition function Z in Eq. 4 precludes direct numerical minimization of the crossentropy when the system size is large, since this requires summing over all N i=1 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).

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 A and VIII B). The adaptive cluster expansion [START_REF] Cocco | Adaptive cluster expansion for the inverse ising problem: convergence, algorithm and tests[END_REF][START_REF] Cocco | Adaptive cluster expansion for inferring boltzmann machines with noisy data[END_REF][START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF] 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 crossentropy 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 γ J = 1/B and γ 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.

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 [START_REF] Ekeberg | Improved contact prediction in proteins: using pseudolikelihoods to infer potts models[END_REF][START_REF] Ekeberg | Fast pseudolikelihood maximization for direct-coupling analysis of protein structure from many homologous amino-acid sequences[END_REF] 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 [START_REF] Ekeberg | Improved contact prediction in proteins: using pseudolikelihoods to infer potts models[END_REF][START_REF] Hopf | Mutation effects predicted from sequence co-variation[END_REF][START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF] and in theoretical analyses [START_REF] Barton | Large pseudocounts and l 2-norm penalties are necessary for the mean-field inference of ising and potts models[END_REF], 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. III, and fixed γ J = N/B, γ 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:

γ J qj b=1 J ij (a, b) = γ h h i (a) , γ J qi a=1 J ij (a, b) = γ h h j (b) , qi a=1 h i (a) = 0 . ( 10
)
As for ACE, the PLM fields and couplings are subsequently transformed to the consensus gauge for comparison.

III. REGULARIZATIONS

A. Removing variable states

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

f i (a) > f 0 , (11)
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) ≡ q a =qi+1 f i (a ).

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. [START_REF] De Juan | Emerging methods in protein co-evolution[END_REF]. The field and couplings parametes on interacting sites are selected from Gaussian distributions of mean µ = 0 and standard deviations σ 2 J = 1 and σ 2 h = 5. Therefore if i and j interact J ij is a 10 × 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 × N , can be compared to the number of parameters to be inferred, M tot = qN + q 2 N (N -1)/2 1.2 10 5 . We have, for the previously listed values of B, B × 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.

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 = N i q i + N i<j 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 q i = q Potts model to the Ising model with only two colors q i = 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.

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 , the fields and the couplings are estimated through

h i (a ) = h i (q i + 1) + log f i (a ) f i (q i + 1) , J ij (a , b) = J ij (q i + 1, b) , (12)
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 ) 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") = α/B to the unseen symbols (a"). Herafter, we have fixed α = 0.1, in the expected range 0 < α < 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 E.

Gauge used in the ACE inference

ACE inference has been implemented in the lattice-gas gauge (Eq. 6) [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF]. 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.

B. 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 D 1 plays a similar role, while not affecting the amplitude of non-zero couplings.

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 δ(a i , a) and δ(a i , a)δ(b i , b) 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 [START_REF] Cocco | Adaptive cluster expansion for the inverse ising problem: convergence, algorithm and tests[END_REF][START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF], 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,

∆f i (a) = δ(a i , a) -f i (a) and ∆f ij (a, b) = δ(a i , a)δ(b i , b) -f ij (a, b)
, and the expected statistical fluctuations due to finite sampling,

σ i (a) = f i (a)(1 -f i (a))/B and σ ij (a, b) = f ij (a, b)(1 -f ij (a, b))/B. The relative error on frequencies is p = 1 N q i,a ∆f i (a) σ i (a) 2 . ( 13
)
The relative error on connected correlations,

c ij (a, b) = δ(a i , a)δ(b i , b) -δ(a i , a) δ(b i , b) , is c = 2 N (N -1)q 2 i<j,a,b ∆c i,j (a, b) σ c i,j (a, b) 2 , ( 14
)
where we estimate the standard deviation in the connected correlations as

σ c i,j (a, b) = σ ij (a, b) + f j (b)σ i (a) + f i (a)σ j (b). Finally, the maximum relative error is max = max {i,j,a,b} 1 2 log (M ) |∆f i (a)| σ i (a) , |∆f ij (a, b)| σ ij (ab) , (15)
where M = N q + (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 ( max < 5) at multiple values (marked by asterisks) of the expansion threshold t, see Table I. 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 I for the threshold corresponding to the minimal relative errors max . In particular the two relative minima better reproducing the data correspond to 2 different inferred networks. The minimum with 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 max = 1 is at low threshold t = 6.4 × 10 -5 where the expansion includes the maximal number of 2 site clusters N 2 = N × (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 A, the fully connected solution is overfitting the data.

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 max 2 of 2-site clusters, see Appendix VIII B. 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 max 2 = N × 2.5 ≈ 125. In practice to find the best sparse graph with a number of edges smaller than the prescribed value N max 2 we spanned the threshold range in regular intervals [START_REF]The first interval is I/τ < t < 1, the second is for 1/τ 2 < t < 1/τ and so on[END_REF], at which the local minima of the absolute error and the corresponding t parameters are recorded (see Table I) The parameters corresponding to the global minimum (under the sparsity conditions) are then chosen (t = 3.6 × 10 -2 for Table I). The spACE procedure is stable when changing N max 2 around its prior fixed average value. In the following we have stopped the algorithm = 200 (shown on Fig. 2), after checking that the value N max 2 = 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 Kmax , with their maximal size K max (see Sec. VI and Appendix VIII C): as shown in Fig. 2 K max = 3 at N max 2 = 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 I. 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 A 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 • 10 -3 , 10 -2 , 3 • 10 -2 , 10 -1 ]. Note that all frequency thresholds 0 < f 0 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.

A. 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:

D(P J real ||P J inf ) = a P J real (a) log P J real (a) P J inf (a) = log(Z inf ) -log(Z real ) + E inf (a) -E real (a) a in real ,
where a = {a 1 , . . . a N } is a configuration and • 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 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) [START_REF] Salakhutdinov | Learning and evaluating boltzmann machines[END_REF], 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) + α) -log(f i (cons i )) where cons i is the most common state at site i and α = 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 I), see Appendix VIII A. The KL divergences obtained directly from the ACE expansion and the ones obtained with importance sampling are very similar, as shown in Table II, 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.

KL Divergence for ACE models at different inclusion thresholds t

Table II displays the KL divergences for the reference data set and different cluster inclusion thresholds of Fig. 2 and in table I 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 I 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 II are significantly different from the ones obtained by the AIS method. 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 [40] 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 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 II shows that the strong regularization (γ J ≈ N/B) used in the fully connected PLM inference is essential to reduce overfitting: a PLM model inferred with a weak regularization (γ J ≈ 1/B) gives indeed very large KL divergences, especially at low sampling depth. Color compression, acting as an additional III: KL divergences between inferred and empirical distributions on the reference data sets at different B and with no color compression. The value obtained by weakly regularized PLM γ J = 1 B is compared to the ones obtained by the standards strongly regularized PLM inference and weakly regularized spACE inference, used in Fig. 3. regularization, helps in this latter case to reduce overfitting and to lower the KL divergence. As shown in Appendix VIII E 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 II.

PLM PLM ACE B (γJ = 1 B ,γ h = 0.002 B ) (γJ = 50 B ,γ h = 0.1 B ) (γJ = 1 B ,γ h = 0.

B. 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) = 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 D). 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:

∆f = i a (f gen i (a) -f test i (a)) 2 i q i , ∆corr = ij ab (cov test ij (a, b) -cov gen ij (a, b)) 2 ij q i • q j . ( 16
)
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 [41]. 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.

C. 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 [START_REF] Sułkowska | Genomics-aided structure prediction[END_REF][START_REF] Hopf | Three-dimensional structures of membrane proteins from genomic sequencing[END_REF][START_REF] Nugent | Accurate de novo structure prediction of large transmembrane protein domains using fragment-assembly and correlated mutation analysis[END_REF]. To this end, we compute the Frobenius norm of the (10×10) inferred and decompressed coupling matrix between each pair of sites i, j,

F ij = 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 [START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF][START_REF] Ekeberg | Fast pseudolikelihood maximization for direct-coupling analysis of protein structure from many homologous amino-acid sequences[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF], 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 G).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 [START_REF] Morcos | Coevolutionary information, protein folding landscapes, and the thermodynamics of natural selection[END_REF][START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF][START_REF] Dunn | Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction[END_REF]. Performance can be further improved with the Average Product Correction (APC) (defined in Appendix VIII G) 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 P P V (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. and the absolute error (∆, as in Eq. 17) are indicated.

D. 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 C 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

∆h = i a h inf i (a) -h real i (a) 2 i q i , ∆J = ij ab (J real ij (a, b) -J inf ij (a, b)) 2 ij q i • q j ; (17)
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 10 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 [START_REF] Cardy | Scaling and renormalization in statistical physics[END_REF] 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 [45].

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 [START_REF] Hopf | Three-dimensional structures of membrane proteins from genomic sequencing[END_REF][START_REF] Morcos | Direct-coupling analysis of residue coevolution captures native contacts across many protein families[END_REF] or the fitness landscape [START_REF] Ferguson | Translating hiv sequences into quantitative fitness landscapes predicts viral vulnerabilities for rational immunogen design[END_REF][START_REF] Mann | The fitness landscape of hiv-1 gag: advanced modeling approaches and validation of model predictions by in vitro testing[END_REF][START_REF] Figliuzzi | Coevolutionary landscape inference and the context-dependence of mutations in beta-lactamase tem-1[END_REF][START_REF] Hopf | Mutation effects predicted from sequence co-variation[END_REF]. 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 H, is usually introduced to reduce the initial number B of sequences in the sample to an effective number of independent ones B ef f . We here consider three protein families, whose fitnesses have been systematically assessed against singlepoint 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 ef f :

• WW ( N = 31, B = 8251, B ef f = 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 [START_REF] Araya | A fundamental protein property, thermodynamic stability, revealed solely from large-scale measurements of protein function[END_REF].

• PDZ (N = 84, B = 24795, B ef f = 4240) is a protein domain present in signaling proteins. Here fitness has been measured in terms of binding affinity [START_REF] Mclaughlin | The spatial architecture of protein function and adaptation[END_REF];

• RRM (N = 82, B = 70780, B ef f = 18800)is an RNA recognition motif; fitness was estimated through growth rate measurements in [START_REF] Melamed | Deep mutational scanning of an rrm domain of the saccharomyces cerevisiae poly (a)-binding protein[END_REF].

Alignments and experimental fitness measures used in this section were taken from [START_REF] Hopf | Mutation effects predicted from sequence co-variation[END_REF]. The PLM procedure was applied with the same large regularization used for the artificial data, γ J = N/B, γ h = 0.1/B. The SpACE procedure was applied with a threshold cutoff N max 2 = 3N , and a smaller regularization γ J = 10/B, γ h = 0.1/B ; however the relative error max (Eq.15) was generally too large even in its local minima, indicating that the procedure had not converged. As shown in [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF] 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 [START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF].

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 [START_REF] Hopf | Mutation effects predicted from sequence co-variation[END_REF] 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 [START_REF] Mann | The fitness landscape of hiv-1 gag: advanced modeling approaches and validation of model predictions by in vitro testing[END_REF]. 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 ef f , as expected and observed in Fig. 8 for the well sampled artificial data (at B = 10 4 and B = 10 5 ). 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 R 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 = 4N ) and PLM on 1 and 20 cores, for the data from the reference graph ER05, at 2 different sampling depth ER 1000 (N = 50B = 1000), ER100000 (N=50,B=100000), for different color compression thresholds f 0 . The number of clusters processed by the algorithm Ncl and the average number of colors explicitly modeled are also indicated q kept .

VI. GAIN IN COMPUTATIONAL TIME FOR SYNTHETIC AND PROTEIN SEQUENCE DATA

B = 10 3 B = 10 5 f 0 tACE Ncl kmax < qi >

WW PDZ RRM1

f 0 tACE Ncl kmax < qi > t = 93, 252, 246 (WW, PDZ, RRM1); the number of clusters processed by the algorithm Ncl and the average number of colors explicitly modeled are also indicated < q kept >. PLM times are given for both on 1 and 20 cores.

. 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 IV 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 q ,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 V and in Fig. 11. As shown in Table V 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 IV. The total number of processed clusters ( Ncl ) and their average size q 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 [START_REF] Fanthomme | Optimal regularizations for multivariate gaussian distributions[END_REF] . 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 [START_REF] Franz | A fast and accurate algorithm for inferring sparse ising models via parameters activation to maximize the pseudo-likelihood[END_REF][START_REF] Cocco | Functional couplings from sequence and mutational data[END_REF]. 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 [START_REF] Cocco | Functional couplings from sequence and mutational data[END_REF]. 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 [START_REF] Barton | Scaling laws describe memories of host-pathogen riposte in the hiv population[END_REF][START_REF] Barton | Relative rate and location of intra-host hiv evolution to evade cellular immunity are predictable[END_REF][START_REF] Barton | On the entropy of protein families[END_REF][START_REF] Louie | Fitness landscape of the human immunodeficiency virus envelope protein that is targeted by antibodies[END_REF] 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 [START_REF] Ekeberg | Improved contact prediction in proteins: using pseudolikelihoods to infer potts models[END_REF][START_REF] Ekeberg | Fast pseudolikelihood maximization for direct-coupling analysis of protein structure from many homologous amino-acid sequences[END_REF] can be crucially important when dealing with large protein sequences or with whole genome inference [START_REF] Pensar | High-dimensional structure learning of binary pairwise markov networks: A comparative numerical study[END_REF][START_REF] Gao | Correlation-compressed direct-coupling analysis[END_REF] 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 [START_REF] Tubiana | Learning protein constitutive motifs from sequence data[END_REF], or variational auto-encoders [START_REF] Riesselman | Deep generative models of genetic variation capture the effects of mutations[END_REF] .

where the sum is over all nonempty clusters of the N variables. The cluster cross-entropy contributions ∆S Γ are recursively defined through

∆S Γ = S Γ - Γ ⊂Γ ∆S Γ . (19)
Here S Γ denotes the minimum of the cross entropy (4) restricted only to the variables in Γ. Thus, S Γ depends only on the frequencies p i (a), p ij (a, b) with i, j ∈ Γ. Provided that the number of variables in Γ is small (typically 10 for q = 10 Potts state as in the present work) numerical maximization of the likelihood restricted to Γ is tractable. The definition of ∆S Γ ensures that the sum over all clusters Γ in (18) yields the cross entropy for the entire system of N variables. As detailed in [START_REF] Cocco | Adaptive cluster expansion for the inverse ising problem: convergence, algorithm and tests[END_REF][START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF], 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 |∆S Γ | > t.

B. ACE and SpACE pseudocodes

The ACE algorithm, described in detail in [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF], 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 :

1. For each cluster Γ ∈ L k (a) Compute S Γ by numerical minimization of Eq. ( 4) restricted to Γ.

(b) Record the parameters minimizing Eq. ( 4), called J Γ .

(c) Compute ∆S Γ using Eq. ( 19). 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 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 max 2 is in the selected list L 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

J(t) = k Γ∈L k (t) ∆J Γ , ∆J Γ = J Γ - Γ ⊂Γ ∆J Γ . ( 20
)
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. [START_REF] Sułkowska | Genomics-aided structure prediction[END_REF][START_REF] Hopf | Three-dimensional structures of membrane proteins from genomic sequencing[END_REF][START_REF] Nugent | Accurate de novo structure prediction of large transmembrane protein domains using fragment-assembly and correlated mutation analysis[END_REF] 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.

C.

Cluster expansion and computational time as a function of the color compression f0 for fully connected graphs.

Fig. 12 shows that the behavior of the maximal relative reconstruction error 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 FIG. 13: Reduction in computational time due to the color compression for fully connected ACE inference on 2 data sets obtained by sampling B=1000 configurations from Erdős-Rényi random graph models. Left: Computational time at the fully connected, low-threshold minimum, as a function of the color compression threshold f 0 . Right: Computational time relative to the one with no color compression as a function of the number of parameters. 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.

D. 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 J B = {J B ij , h B i } the inferred parameters at sample size B, and J true = {J true ij , h true i } the true underlying model parameters. The inferred cross-entropy at sampling B writes

S B = - σ P J B (σ) log P J B (σ) , (21)
where the sum is over all possible configurations σ = {σ 1 , ..., σ N }. The inferred probability distribution at finite sampling B is

P J B (σ) = exp N i=1 h B i σ i + N k<l J B kl σ k σ l Z B . (22)
The Kullback-Leibler (KL) divergence between the true and the inferred distributions writes

D(P J true ||P J B ) = σ P J true (σ) log P J true (σ) P J B (σ) = -S true - σ P J true (σ) i h B i σ i + k<l J B kl σ k σ l -log Z B = -S true + log Z B - σ P J true (σ) i h B i σ i + k<l J B kl σ k σ l .
However, Eqs. ( 21) & [START_REF] Hopf | Mutation effects predicted from sequence co-variation[END_REF] give

log Z B = S B + σ P J B (σ) i h B i σ i + k<l J B kl σ k σ l .
The KL divergence between the true and the inferred distributions then writes

D(P J true ||P J B ) =(S B -S true ) - σ P J true (σ) i h B i σ i + k<l J B kl σ k σ l + σ P J B (σ) i h B i σ i + k<l J B kl σ k σ l .
Moreover, a reasonable approximation is

S true = - σ P J true (σ) log P J true (σ) ≈ S B→∞ = - σ P J B→∞ (σ) log P J B→∞ (σ) , (23)
because the true underlying parameters are recovered by the inference method in the perfect sampling case: P J B→∞ (σ) → P J true (σ). Therefore,

D(P J true ||P J B ) = (S B -S ∞ ) + i h B i ( σ i B -σ i ∞ ) + k<l J B kl ( σ k σ l B -σ k σ l ∞ ) , (24)
where • B = σ •P J B (σ), and • ∞ = σ •P J B→∞ (σ) ≈ σ •P J true (σ). It naturally generalizes to the q-states Potts case:

D(P J true ||P J B ) =(S B -S ∞ ) + N i=1 q a=1 h B i (a)( σ ia B -σ ia ∞ ) + N k,l=1 k<l q c,d=1 J B kl (c, d)( σ kc σ ld B -σ kc σ ld ∞ ) . (25)
The artificial data are in a compressed representation (cf. Section III). The complete inferred parameters are recovered as explained in Eq. ( 12).

E. 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:

S ind = log q a=1 e hi(a) - q a=1 h i (a) p i (a) + γ h • q a=1 h i (a) 2 (26)
which, for the unseen colors in the gauge of Z = q a=1 e hi(a) = 1, gives:

h Γ u = -Lw( 1 2γ h ) (27)
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,

h p u = log( α B ) (28)
where we set, as in the rest of the paper, α = 0.1. In this independent-site approximation there is then a shift between the two procedures given by ∆h = h i (a) Γ -h i (a) p , that depends from the pseudocount α and from the value of the regularization γ h . In table VI we give these shifts for the two values of γ h used respectively by ACE (γ h = 0.01/B) and PLM (γ h = 0.1/B).

If, in the approximation of independent-sites, the difference ∆h is the same in all gauges, the specific values of h i (a) Γ and h i (a) p of table VI 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) Γ 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. In this section we study what happens with PLM at lower L 2 regularization, e.g. γ J = 1/B as the one used for ACE, when the threshold for color compression is varied. Without color compression, the performance obtained for γ J = 1/B becomes significantly worse, see Table III.

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 (γ 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 ≤ 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. In Figs. [START_REF] De Juan | Emerging methods in protein co-evolution[END_REF] 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 Fscore = 2 T P (N pred )

N pred + N 0 (29)
where T P 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.

G. Definition of Positive Predicted Value and Average Product Correction

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

P P V (n) = T P (n) n . (30)
Where T P (n) is the number of true predicted edges in the top n pairs. the Average Product Correction (APC) [START_REF] Morcos | Coevolutionary information, protein folding landscapes, and the thermodynamics of natural selection[END_REF][START_REF] Cocco | Inverse statistical physics of protein sequences: a key issues review[END_REF][START_REF] Dunn | Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction[END_REF],

F AP C ij = F ij - F i. F .j F .. . (31)
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.

H. 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 ∈ (0, 1):

w m = 1 ||{n|1 ≤ n ≤ M ; d H [(a n
1 , ..., a n L ), (a m 1 , ..., a m L )] ≤ xL}|| [START_REF] Barton | Ace: adaptive cluster expansion for maximum entropy graphical model inference[END_REF] 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

M ef f = M m=1 w m (33)
can be interpreted as the effective number of independent sequences.

  a1,...a N δ(a i , a)P (a 1 . . . , a N |J) = f i (a) , a1,...a N δ(a i , a)δ(a j , b)P (a 1 . . . , a N |J) = f ij (a, b) ,

FIG. 1 :

 1 FIG. 1: Mean number of colors per site <q i > (top) and relative reduction in the number of parameters to be inferred M cc /M tot (bottom), as a function of the f 0 : no compression, 0 (only unseen symbols are removed from inference) and 1/B < f 0 < 0.1. Data from a single ER realization (left), and sequence data for 3 protein families WW (N = 31), PDZ (N = 84), RRM (N = 82) (right) at different sample depths B.

2 FIG. 2 :

 22 FIG.2: Cluster expansion as function of the cluster inclusion threshold t for one data set, obtained by sampling one realization of an Erdős-Rényi random graph with B = 1000 configurations and applying ACE to compressed data with f 0 = 0.01. From top to bottom: i) relative reconstruction errors versus t on frequencies p , connected correlations c and maximal relative error max . Stars indicate possible solutions of the inverse models, reproducing well the empirical 1-and 2-point correlations ( max < 5), which are obtained at different expansion thresholds and corresponding to sparse graphs at large thresholds, or fully connected reconstructed graphs at small thresholds, see TableIii) Regularized cross entropy vs. t, iii) number of total clusters N cl and number of 2-site clusters N 2 included in the expansion vs. t, iv) maximal cluster size versus t. The spACE expansion stops at N 2 = 200.

FIG. 3 :

 3 FIG. 3: Kullback-Leibler divergence between real and inferred probability distributions averaged over 10 realizations plotted as function of the compression parameter f 0 for different sample sizes. The left plot is for ACE, the right plot for PLM. Error bars are standard deviations over the 10 realizations. Black empty squares correspond to the reference cutoff frequency f * 0 = 10/B.

FIG. 4 :

 4 FIG. 4: Reconstruction of average frequencies and covariances for the reference case (ER05, B=1000). Comparison between generated data and test set for no color compression (top) and f 0 = 0.01 for ACE ( left,a) and PLM ( right, b). The Pearson correlation coefficient (R) and the absolute error (∆, Eq. 16) are marked on top of the plots both for the full model and the for the reduced one (only explicitly modeled states).

  generated a δ(a i , a)/B gen and covariances cov ij (a, b) = generated a[δ(a i , a)δ(a j , b)/B gen ] -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.

5 FIG. 5 :

 55 FIG.5: Absolute error (Eq. 16) of frequencies (∆f , top panels) and covariances (∆corr, bottom panels) averaged over 10 ER realizations, as a function of the color compression for several sample sizes. Dashed lines: error on explicitly modeled Potts states only. Full lines: error on all parameters. Error-bars are standard deviations computed over the 10 realizations. Inference is performed respectively by ACE (left) and PLM (right).

  FIG.6: ACE and PLM interaction graph reconstruction as a function of the color compression f 0 . Positive Predicted value P P V (n = min {N 0 , N pred }) for ACE and PLM inference. Points and error bars are averages and standard deviations obtained on 10 ER realizations. For ACE all quantities are computed in the consensus gauge with (dotted line) and without (full line) APC. (Left and middle plot). Right Plot P P V (N 0 ) with PLM either in the consensus gauge without APC (full line), or in the zero-sum gauge with APC (dotted line).

FIG. 7 :

 7 FIG. 7: Comparison of inferred and real fields and couplings with ACE and PLM for one realization of ER graph with no color compression (top) and f 0 = 0.01 (bottom), for B=1000 sampled configurations. Parameters on explicitly modeled (kept), grouped, and unseen Potts states are colored differently. Left: field comparison. Right: coupling comparison. On top of each plot, the Pearson correlation coefficient (R) and the absolute error (∆, as in Eq. 17) are indicated.

5 FIG. 8 :

 58 FIG.8: Absolute errors (Eq. 17) on fields (∆h, top panels) and couplings (∆J, bottom panels) averaged over 10 ER realizations, as a function of the color compression for several sample sizes. Dashed lines: error on parameters related to explicitly modeled Potts states. Full lines: error on all parameters, after decompression of unseen and grouped Potts states (see Sec. III). Error-bars are standard deviations computed over the 10 realizations. Inference is performed respectively by ACE (left) and PLM (right).

FIG. 9 :

 9 FIG.9: Spearman correlation coefficient between the fitness predictions and the experimental measures found in literature, as a function of the color compression for ACE and PLM. The protein families used here are: WW (PF00397, in orange), PDZ (PF00595, in green) and RRM1 (PF00076 in light blue)

FIG. 10 :FIG. 11 :

 1011 FIG.10: Computational times averaged on all ER data for ACE and PLM for different B and color compression. Left: computational time, in seconds, using only 1 CPU, for ACE and PLM as a function B; a power law/ linear fits are added to ACE /PLM (dashed line). Right: computational time ratio between the compressed and decompressed inference for ACE (1 CPU) and PLM (25 parallel CPU) as a function of fraction of parameters to infer; dotted line: x=y; full line: linear PLM fit.

2 .

 2 Add all clusters Γ ∈ L k with |∆S Γ | > t to a new list L k (t).

3 .

 3 Construct a list L k+1 of clusters of size k + 1 from overlapping clusters in L k (t).

1 FIG. 12 :

 112 FIG. 12: Maximum relative error as a function of the expansion threshold for a particular graph realization (same used in Fig 5: ER05, sampled with B=1000),for different color compressionf 0 .

1 TABLEFIG. 14 :

 114 FIG. 14: Fields for the unseen Potts symbols for ER05 and B = 1000 in the non compressed model and the reference model with f 0 = 0.01 of Fig. 7. Dotted line x=y. Red Line: the shift given in the independent model for γ h = 1/10B as the one used with PLM.

F.

  Graphical reconstruction: Contact maps and F score for spACE graph reconstruction as a function of color compression

FIG. 15 :

 15 FIG.15: KL divergence between the true model and the one inferred with γ J = 1/B averaged over 10 realizations for several sample sizes B at different color compression thresholds f 0 (full line). Error bars are standard deviations over the 10 realizations. Horizontal dashed lines are there for comparison and correspond to the KL values obtained with γ J = 50/B without color compression.

  FIG. 16: ACE and PLM interaction graph reconstruction and PPV curve for one realization of ER graph. Left: contact maps. Upper triangular: contact map with real contacts (black empty squares), inferred true positive (full circles in empty squares), inferred false positive (full circles) in consensus gauge without Average Product Correction (APC). Lower triangular: Frobenius norm of the inferred parameters in consensus gauge with color-scale on the right. Right: Positive Predicted Value (PPV) curve in consensus and zero-sum gauge with and without APC. Top: no color compression. Bottom: f 0 = 0.01.

FIG. 17 :

 17 FIG. 17: Fscore for ACE inference. Points and error bars are averages and standard deviations obtained on 10 ER realizations. All quantities are computed in the consensus gauge with (dotted line) and without (full line) APC.

  Table I ii) Regularized cross entropy vs. t, iii) number of total clusters N cl and number of 2-site clusters N 2 included in the expansion vs. t, iv) maximal cluster size versus t. The spACE expansion stops at N 2 = 200.

	t	max N2 S λ c	Sc
	1	17	0	56 56
	9.6 × 10 -2 4.9

*

29 51.7 51.2 3.6 × 10 -2 3.9 * 55 50.6 49.9 2.2 × 10 -2 34 90 49.7 48.8 3.4 × 10 -5 1 * 1225 42.3 39

TABLE I :

 I 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. The Table gives the cluster inclusion threshold t, the number N 2 of 2-site clusters, the maximal relative error max , the regularized cross entropy S λ c and the cross entropy S c obtained with the cluster expansion at different thresholds t. The entropy of the model having generated the data is S = 50.3. Possible solutions ( max < 5) are indicated by asterisks. The optimal threshold determined by the spACE procedure is 3.6 × 10 -2 .

	for N max 2

TABLE II :

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

TABLE

  t 1 P LM t 20 P LM tACE Ncl kmax < qi > t 1

								P LM	t 20 P LM
	no 10 3 2511 4	10	39	8	9 10 2 1306 4	10 1.8 10 4 1.5 10 3
	0	10 3 2320 4	8.6	32	8	9 10 2 1306 4	9.9 1.8 10 4 1.5 10 3
	0.01 10 3 2065 4	5.8	22	7	9 10 2 1305 4	5.8	10 4	10 3
	0.1 8 10 2 1350 4	2.2	11	7	6 10 2 1314 4	2.3	2 10 3 3 10 2

TABLE IV :

 IV Running times, in seconds, for spACE (N max 2

TABLE V :

 V

Running times (in seconds) for ACE and PLM for the three studied protein families: WW (N=31 B=8251), PDZ (N=84,B=24795), RRM1 (N=82,B=70780) for different color compression thresholds f 0 . spACE has been stopped at N max 2

VIII. APPENDIX A. Reminder about Adaptive Cluster Expansion and the inclusion thresholdIn the ACE inference procedure the cross entropy is expanded as the sum of cluster contributions. Defining a cluster as a sub-set of variables: Γ = {i 1 , . . . , i k }, k ≤ N , we can formally write the cross-entropy as the sum of cluster contributions:S(J|f ) = Γ ∆S Γ ,(18)

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.