adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis

adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis Lorenzo Rosset , Roberto Netti , Anna Paola Muntoni , View ORCID Profile Martin Weigt , Francesco Zamponi doi: https://doi.org/10.1101/2025.01.31.635874 Lorenzo Rosset a Sorbonne Université, CNRS, Laboratory of Computational and Quantitative Biology , 75005 Paris, France b Laboratoire de Physique Théorique, Éole Normale Supérieure , 75231 Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Roberto Netti a Sorbonne Université, CNRS, Laboratory of Computational and Quantitative Biology , 75005 Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Anna Paola Muntoni c DISAT, Politecnico di Torino , 10129 Torino, Italy Find this author on Google Scholar Find this author on PubMed Search for this author on this site Martin Weigt a Sorbonne Université, CNRS, Laboratory of Computational and Quantitative Biology , 75005 Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Martin Weigt Francesco Zamponi d Dipartimento di Fisica, Sapienza Università di Roma , 00185 Rome, Italy b Laboratoire de Physique Théorique, Éole Normale Supérieure , 75231 Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: francesco.zamponi{at}uniroma1.it Abstract Full Text Info/History Metrics Preview PDF Abstract In this methods article, we provide a flexible but easy-to-use implementation of Direct Coupling Analysis (DCA) based on Boltzmann machine learning, together with a tutorial on how to use it. The package adabmDCA 2.0 is available in different programming languages (C++, Julia, Python) usable on different architectures (single-core and multi-core CPU, GPU) using a common front-end interface. In addition to several learning protocols for dense and sparse generative DCA models, it allows to directly address common downstream tasks like residue-residue contact prediction, mutational-effect prediction, scoring of sequence libraries and generation of artificial sequences for sequence design. It is readily applicable to protein and RNA sequence data. 1 Introduction Over the past decade, the increase in sequence data availability has significantly propelled the application of machine learning in biomolecular studies. Generative models, in particular, have emerged as highly effective tools. Even though the field of bioinformatics is moving towards the adoption of larger and larger deep neural network models, notably with the advent of Large Language Models (LLMs), Direct Coupling Analysis (DCA) methods [ 3 ] provide a simple, user-friendly, energy-efficient and interpretable tool that can be successfully adopted for generating functional protein sequences [ 14 ], predicting contacts in the ternary structure [ 11 ], predicting mutational effects in proteins [ 9 ] and unveiling coevolutionary signals in MSAs of homologous protein or RNA families [ 4 ]. This tutorial presents a new version of adabmDCA [ 12 ]. The package comes in three different languages: C++ (single-core CPU), Julia (multi-core CPU), and Python (GPU-oriented). They share the same front-end interface from the terminal allowing the user to install and use one of the three equivalent versions based on hardware or software constraints. We provide three different training routines: bmDCA : Trains a fully-connected DCA model [ 8 ]; eaDCA : Trains a DCA model on a sparse coupling network by progressively adding couplings during the training [ 2 ]; edDCA : Starts from a trained bmDCA model and iteratively removes the less informative couplings until the target sparsity is reached [ 1 ]. Additionally, we provide several routines for sampling and analyzing the generated sequences once a DCA model is trained, for constructing and evaluating - according to a DCA model - a single mutant library from a given wild type, and finally, for computing the pairwise contact scores, in terms of average-product corrected Frobenius norms of the DCA couplings [ 6 ]. This tutorial is organized as follows. In the first section, we present the theoretical framework of DCA models that is useful for understanding the parameters of the algorithm and the possible issues that one can encounter during the training. We also describe the typical preprocessing pipeline for preparing the data before using them to train the DCA model. In the second section, we describe the different routines and the most important hyperparameters that need to be set. Finally, in the last section, we summarize the main terminal commands for the different routines and we summarize all the input parameters. The reader who is already familiar with DCA modeling of biological sequences can skip the first part of the tutorial and go directly to the third section. 2 Boltzmann learning of biological models 2.1 Structure of the model and Boltzmann learning DCA models ( Figure 1 ) are probabilistic generative models that infer a probability distribution over sequence space. Their objective is to assign high probability values to sequences that are statistically similar to the natural sequences used to train the architecture while assigning low probabilities to those that significantly diverge. The input dataset is a Multiple Sequence Alignment (MSA), where each sequence is represented as a L -dimensional categorical vector a = ( a 1 , … , a L ) with a i ∈ {1, … , q }, each number representing one of the possible amino acids or nucleotides, or the alignment gap. To simplify the exposition, from here on, we will assume them to be amino acids. The following equation then gives the DCA probability distribution: Download figure Open in new tab Figure 1: Schematic representation of a DCA model. In this expression, Z is the normalization constant and E is the DCA energy function . The interaction graph 𝒢 = (𝒱, ℰ ) is represented in a one-hot encoding format: the vertices 𝒱 of the graph are the L × q combinations of all possible symbols on all possible sites, labeled by ( i, a ) ∈ {1, …, L } × {1, …, q }, and the edges E connect two vertices ( i, a ) and ( j, b ). The bias (or field) h i ( a ) corresponding to the amino acid a on the site i is activated by the Kronecker if and only if a i = a , and it encodes the conservation signal of the MSA. The coupling matrix J ij ( a, b ) represents the coevolution (or epistatic) signal between pairs of amino acids at different sites, and is activated by the term if a i = a and a j = b . Note that J ii ( a, b ) = 0 for all a, b to avoid a redundancy with the h i ( a ) terms, and that J ij ( a, b ) = J ji ( b, a ) is imposed by the symmetry structure of the model. The interaction graph 𝒢 can be chosen fully connected as in the bmDCA model, or it can be a sparse graph as in eaDCA and edDCA. Training the model The training consists of adjusting the biases, the coupling matrix, and the interaction graph to maximize the log-likelihood of the model for a given MSA, which can be written as where w ( m ) is the weight of the data sequence m , with , and are the empirical single-site and two-site frequencies computed from the data. Roughly speaking, f i ( a ) tells us what is the empirical probability of finding the amino acid a in the position i of a training sequence, whereas f ij ( a, b ) tells us how likely it is to find together the amino acids a and b at positions i and j , respectively, in a training sequence. For a fixed graph 𝒢, we can maximize the log-likelihood by iteratively updating the parameters of the model in the direction of the gradient of the log-likelihood, meaning where γ is a small rescaling parameter called learning rate . By differentiating the log-likelihood (2), we find the update rule for Boltzmann learning: where and are the one-site and two-site marginals of the model (1). Notice that the convergence of the algorithm is reached when p i ( a ) = f i ( a ) and p ij ( a, b ) = f ij ( a, b ). Monte Carlo estimation of the gradient The difficult part of the algorithm consists of estimating p i ( a ) and p ij ( a, b ), because computing the normalization Z of Eq. (1) is computationally intractable, preventing us from directly computing the probability of any sequence. To tackle this issue, we estimate the first two moments of the distribution through a Monte Carlo simulation. This consists of sampling a certain number of fairly independent sequences from the probability distribution (1) and using them to estimate p i ( a ) and p ij ( a, b ) at each learning epoch. There exist several equivalent strategies to deal with it. Samples from the model (1) can be obtained via Markov Chain Monte Carlo (MCMC) simulations either at equilibrium or out-of-equilibrium, where we start from N c configurations (we refer to them as chains ), chosen uniformly at random, from the data or the last configurations of the previous learning epoch, and update them using Gibbs or Metropolis-Hastings sampling steps up to a certain number of MCMC sweeps. It has been shown in [ 12 ] for Boltzmann machines and, in general, for energy-based models [ 5 ] that under certain conditions, training the model estimating the marginals from an out-of-equilibrium sampling leads to good enough generative models. Depending on the training set, these models can even be fairly close to what is achieved with perfectly equilibrated training. For this reason, in adabmDCA 2.0 , we implement a persistent contrastive divergence (PCD) scheme, in which the computation of the gradient is done with chains that may be slightly out of equilibrium. In particular, chains are persistent , i.e. they are initialized at each learning epoch using the last sampled configurations of the previous epoch. This is done because sampling from configurations that are already close to the stationary state of the model at the current training epoch is much more convenient. Furthermore, the number of sweeps to be performed should be large enough to ensure that the updated chains are close enough to an equilibrium sample of the probability (1). In practice, this is done by fixing the number of sweeps to a convenient value, k , that trades off between a reasonable training time and a fair mixing of the chains. Convergence criterium To decide when to terminate the training, we monitor the two-site connected correlation functions of the data and of the model, which are defined as When the Pearson correlation coefficient between the two reaches a target value, set by default at 0.95, the training stops. Once the model is trained, we can generate new sequences by sampling from the probability distribution (1), infer contacts on the tertiary structure by analyzing the coupling matrix, or propose and assess mutations through the DCA energy function E (see Section 4 ). While during the training we keep the number of Monte Carlo sweeps fixed, we underline that whenever we want to sample new sequences from a model (i.e. after convergence of training) we compute the mixing time as explained in 4.1 and we ensure that sequences are sampled in equilibrium. 2.2 Training sparse models What we have described so far is true for all the DCA models considered in this work, but we have not yet discussed how to adjust the topology of the interaction graph 𝒢. In the most basic implementation, bmDCA, the graph is assumed to be fully connected (every amino acid of any residue is connected with all the rest) and the learning will only tweak the strength of the connections. This results in a coupling matrix with L ( L −1) q 2 / 2 independent parameters, where L is the aligned sequence length, q = 21 for amino acid and q = 5 for nucleotide sequences. However, it is well known from the literature that the interaction network of protein families tends to be relatively sparse , suggesting that only a few connections should be necessary for reproducing the statistics of biological sequence data. This observation brings us to devising a training routine that produces a DCA model capable of reproducing the one and two-body statistics of the data with a minimal amount of couplings. To achieve this, we implemented two different routines: eaDCA that promotes initially inactive (i.e. zero) coupling parameters to active ones starting from a profile model, and edDCA, which iteratively prunes active but negligible coupling parameters starting from a dense fully-connected DCA model. Element activation DCA In eaDCA ( Figure 2-B ), we start from an empty interaction graph ℰ= ⊘, meaning that no connection is present. Each training step is divided into two different moments: first, we update the graph, and then we bring the model to convergence once the graph is fixed (a similar pipeline has been proposed in [ 2 ]). Download figure Open in new tab Figure 2: Schematic representation of the sparse model training. A) edDCA, the sparsification is obtained by progressively pruning contacts from an initial fully connected model. B) eaDCA, the couplings are progressively added during the training. To update the graph, we first estimate p ij ( a, b ) for the current model. Then, we decide how many couplings we want to target as the fraction factivate of the number of currently inactive couplings and, from all the possible quadruplets of indices ( i, j, a, b ), we select those for which p ij ( a, b ) is “the most distant” (according to the criterion introduced in [ 2 ]) from the target statistics f ij ( a, b ). We then activate the couplings corresponding to these quadru-plets, obtaining a new graph ℰ ′ ⊇ℰ. Notice that some of the selected couplings might be already active, so the model partially auto-regulates the number of new couplings that have to be activated at each step. Element decimation DCA In edDCA ( Figure 2-A ), we start from a previously trained bmDCA model and its fully connected graph 𝒢. We then apply the decimation algorithm, in which we prune connections from the edges ℰ until a target density of the graph is reached, where the density is defined as the ratio between the number of active couplings and the number of couplings of the fully connected model. Similarly to eaDCA, each iteration consists of two separate moments: graph updating and active parameter updating. To update the graph, we remove the fraction drate of active couplings that, once removed, produce the smallest perturbation on the probability distribution at the current epoch. In particular, for each active coupling, one computes the symmetric Kullback-Leibler distances between the current model and a perturbed one, without that target element. One then removes the drate elements which exhibit the smallest distances (see [ 1 ] for further details). Parameter updates in between decimations/activations In both procedures, to bring the model to conver-gence on the graph, we perform a certain number of parameter updates in between each step of edge activation or decimation, using the update rule (5). Between two subsequent parameter updates, k sweeps are performed to update the Markov chains. In the case of element activation we perform a fixed number of parameter updates, specified by the input parameter gsteps . Instead, when pruning the graph, we keep updating the parameters with the rule (5) until the Pearson correlation coefficient reaches a target value. 2.3 Input data and pre-processing adabmDCA 2.0 takes as input a multiple sequence alignment (MSA) of aligned amino acid or nucleotide sequences, usually forming a protein or RNA family. DCA implementations require the data to be saved in FASTA format [ 13 ]. adabmDCA 2.0 implements the three default alphabets shown in table 1 , but the user can specify an ad-hoc alphabet as far as it is compatible with the input MSA. View this table: View inline View popup Download powerpoint Table 1: Default alphabets implemented in adabmDCA 2.0 , where “-” indicates the alignment gap. View this table: View inline View popup Download powerpoint Table 2: Command line arguments for training bmDCA, eaDCA or edDCA. View this table: View inline View popup Download powerpoint Table 3: Command line arguments for sampling from a DCA model. View this table: View inline View popup Download powerpoint Table 4: Command line arguments for computing the DCA energies of an MSA. View this table: View inline View popup Download powerpoint Table 5: Command line arguments for generating the deep mutational scan of a provided wild type sequence. View this table: View inline View popup Download powerpoint Table 6: Command line arguments for computing the Frobenius contact matrix of a DCA model. An example of a FASTA file format is shown in Figure 3 . In particular, adabmDCA 2.0 correctly handles FASTA files in which line breaks within a sequence are present. Download figure Open in new tab Figure 3: Example of RNA sequences formatted in FASTA format. 2.3.1 Preprocessing Preprocessing pipeline The adabmDCA 2.0 code applies the following preprocessing pipeline to the input MSA: p.1 Remove the sequences having some tokens not included in the chosen alphabet; p.2 If needed, compute the importance weights for the sequences in the MSA; p.3 Apply a pseudocount to compute the MSA statistics. Their precise implementation is described in the following. Computing the importance weights The sequence weights are computed to mitigate as much as possible the systematic biases in the data, such as correlations due to the phylogeny or over-representation of some regions of the sequence space because of a sequencing bias. Given an MSA of M sequences, to compute the importance weight of each sequence a ( m ) , m = 1, … , M , we consider N ( m ) as the number of sequences in the dataset having sequence identity from a ( m ) greater or equal to 0.8 L (this threshold can be tuned by the user). Then, the importance weight of a ( m ) will be This reweighting allows us to give less importance to sequences found in very densely populated regions of the sequence space while enhancing the importance of isolated sequences. Alternatively, the user can decide to provide an external file containing the weights, or to ask the program to assign uniform weights across the sequences in the MSA. Pseudocount and reweighted statistics DCA models are trained to reproduce the one and two-site frequencies of the empirical data. To compute these, we introduce in the computation of the empirical statistics a small parameter α , called pseudocount, that allows us to deal with unobserved (pairs of) symbols in one (or two) column(s) of the MSA. The one and two-site frequencies are given by where and are computed from the MSA as in Eq. (3). 3 Implementation All the software implementations that we propose (Python, Julia, and C++) offer the same interface from the terminal through the adabmDCA command. The complete list of training options can be listed through the command Download figure Open in new tab The standard command for starting the training of a DCA model is Download figure Open in new tab where ∈ {bmDCA, eaDCA, edDCA} selects the training routine. By default, the fully connected bmDCA algorithm is used. edDCA can follow two different routines: either it decimates a pre-trained bmDCA model, or it first trains a bmDCA model and then decimates it. The corresponding commands are shown below (see section 3.2 ); is the FASTA file, with the complete path, containing the training MSA; is the path to a (existing or not) folder where to store the output files; is an optional argument. If provided, it will label the output files. This is helpful when running the algorithm multiple times in the same output folder. Once started, the training will continue until the Pearson correlation coefficient between the two-point connected correlations of the model and the empirical ones obtained from the data reaches a modifiable target value (set by default at target = 0.95 ). In figure 4 we show in a log-log plot the evolution of (one minus) the Pearson correlation coefficient as a function of the training time. The long-time evolution of the training is well approximated by a straight line, meaning that the Pearson growth follows a power law. As such, one should keep in mind that the training procedure converges very rapidly at the beginning (say, it reaches Pearson values of ∼0.9 in about 100 iterations in figure 4 ), while reaching higher values of the Pearson requires comparatively more algorithmic time. Therefore, if needed, it is possible to fit a coarse model in a very short amount of time by just lowering the target value somewhere around 0.9. Download figure Open in new tab Figure 4: One minus the Pearson correlation coefficient between the data and chains correlation matrices as a function of the training time for the protein family PF00072. The curve is well approximated by a power law decay. Output files By default the training algorithm outputs three kinds of text files: params.dat : file containing the parameters of the model saved in this format: – Lines starting with J represent entries of the coupling matrix, followed by the two interacting positions in the sequence and the two amino acids or nucleotides involved. – Lines starting with h represent the bias, followed by a number and a letter indicating the position and the amino acid or nucleotide subject to the bias. Note that inactive, i.e. zero couplings are not included in the file. chains.fasta : FASTA file containing the sequences corresponding to the last state of the Markov chains used during the learning; adabmDCA.log : .log file containing information collected during the training procedure. During the training, the output files containing the parameters and the chains are overwritten, and the log file is updated, every 50 gradient updates (epochs) for bmDCA , and every 10 graph updates for eaDCA and edDCA . Restore an interrupted training It is possible to start the training by initializing the parameters of the model and the chains at a given checkpoint. To do so, two arguments specifying the path of the parameters and the chains are needed: Download figure Open in new tab Importance weights It is possible to provide the algorithm with a pre-computed list of importance weights to be assigned to the sequences (see section 2.3.1 ) by giving the path to the text file to the argument -w . If this argument is not provided, the algorithm will automatically compute the weights using Eq. (7) and it will store them into the folder as weights.dat . The default sequence identity threshold used for computing the weights is set to 0.8, but this value can be modified using the argument --clustering seqid . It is also possible to ask the routine to use uniform weights for all sequences by means of the flag --no reweighting . Choosing the alphabet By default, the algorithm will assume that the input MSA belongs to a protein family, and it will use the preset alphabet defined in Table 1 (by default: --alphabet protein ). If the input data comes from RNA or DNA sequences, it has to be specified by passing respectively rna or dna to the --alphabet argument. There is also the possibility of passing a user-defined alphabet, provided that all the tokens match with those that are found in the input MSA. This can be useful if one wants to use a different order than the default one for the tokens, or in the eventuality that one wants to handle additional symbols present in the alignment. This is done using --alphabet ABCD- if for example the alphabet contains the symbols A, B, C, D, - . 3.1 eaDCA To train an eaDCA model, we just have to specify --model eaDCA . Two more hyperparameters can be changed: --factivate : The fraction of inactive couplings that are selected for the activation at each update of the graph. By default, it is set to 0.001. --gsteps : The number of parameter updates to be performed on the given graph. By default, it is set to 10. For this routine, the number of sweeps for updating the chains can be typically reduced to 5, since only a fraction of all the possible couplings have to be updated at each iteration. 3.2 edDCA To launch a decimation with default hyperparameters, use the command: Download figure Open in new tab where and are, respectively, the file names of the parameters and the chains (including the path) of a previously trained bmDCA model. The edDCA can perform two routines as described above. In the first routine, it uses a pre-trained bmDCA model and its associated chains provided through the parameters and . The routine makes sure everything has converged before starting the decimation of couplings. It repeats this process for up to 10,000 iterations if needed. If these parameters are not supplied, the second routine initializes the model and chains randomly, trains the bmDCA model to meet convergence criteria, and then starts the decimation process as described. Some important parameters that can be changed are: --gsteps : The number of parameter updates to be performed at each step of the convergence process on the given graph. By default, it is set to 10. --drate : Fraction of active couplings to be pruned at each graph update. By default, it is set to 0.01. --density : Density of the graph that has to be reached after the decimation. The density is defined as the ratio between the number of active couplings and the number of couplings of the fully connected graph. By default, it is set to 0.02. --target : The Pearson correlation coefficient to be reached to assume the model to have converged. By default, it is set to 0.95. 3.3 How to choose the hyperparameters? The default values for the hyperparameters are chosen to be a good compromise between having a relatively short training time and a good quality of the learned model for most of the typical input MSA, where for typical we mean a clean MSA with only a few gaps for each sequence and a not too structured dataset (not too clustered in subfamilies) that could create some ergodicity problems during the training. It may happen, though, that for some datasets some adjustments of the hyperparameters are needed to get a properly trained model. The most important ones are: Learning rate By default, the learning rate is set to 0.05, which is a reasonable value in most cases. If the resampling of the model is bad (very long thermalization time or mode collapse), one may try to decrease the learning rate through the argument --lr to some smaller value (e.g. 0.01 or 0.005). Number of Markov Chains By default, the number of Markov chains is set to 10000, which works well in most cases. Using fewer chains reduces the memory required to train the model, but it may also lead to a longer algorithm convergence time. To change the number of chains, we can use the argument --nchains . Number of Monte Carlo steps The argument --nsweeps defines the number of Monte Carlo chain updates (sweeps) between one gradient update and the following. A single sweep is obtained once we propose a mutation for all the residues of the sequence. By default, this parameter is set to 10, which is a good choice for easily tractable MSAs. The higher this number is chosen, the better the quality of the training will be, because in this way we allow the Markov chains to decorrelate more to the previous configuration. However, this parameter heavily impacts the training time, so we recommend choosing it in the interval 10 - 50. Regularization Another parameter that can be adjusted if the model does not resample correctly is the pseu-docount, α , which can be changed using the key --pseudocount . The pseudocount is a regularization term that introduces a flat prior on the frequency profiles, modifying the frequencies as in equations (8) and (9). If α = 1 we impose an equal probability to all the amino acids to be found in the residue at position i , while if α = 0 we just use the empirical frequencies of the data. By default, the pseudocount is set as the inverse of the effective number of sequences, 1 /M eff , where 4 Applications 4.1 Generate sequences Once we have a trained model, we can use it to generate new sequences. This can be done using the command: Download figure Open in new tab where output folder is the directory where to save the data and num gen is the number of sequences to be generated. The routine will first compute the mixing time ( t mix ) of the model by running a simulation starting from the sequences of the input MSA. After that, it will randomly initialize num gen Markov chains and run for nmix × t mix sweeps to ensure that the model equilibrates. It will save in the output directory a FASTA file containing the sequences sampled with the model and a text file containing the records used to determine the convergence of the algorithm. Figure 6 -Right shows the comparison between the entries of the covariance matrices obtained from the data and from the generated sequences. The Pearson correlation coefficient is the same used as a target for the training and the slope is close to 1, meaning that the model is able to correctly recover the two-sites statistics of the data MSA. Download figure Open in new tab Figure 5: Format for saving the coupling matrix J and biases h for an RNA domain. Download figure Open in new tab Figure 6: Analysis of a bmDCA model. Left : measuring the mixing time of the model using 10 4 chains. The curves represent the average overlap among randomly initialized samples (dark blue) and the one among the same sequences between times t and t/ 2 (light blue). Shaded areas represent the error of the mean. When the two curves merge, we can assume that the chains at time t forgot the memory of the chains at time t/ 2. This point gives us an estimate of the mixing time of the model, t mix . Notice that the times start from 1, so the starting conditions are not shown. Right : Scatter plot of the entries of the covariance matrix of the data versus that of the generated samples. Convergence criterion To determine the convergence of the Monte Carlo simulation, the following strategy is used. We extract nmeasure = N sequences from the data MSA according to their statistical weight and we make a copy of them, which are used to initialize a set of N chains. We then compare two sets of chains. The first set represents the chains simulated up to time t , and we denote them as A ( t ) = { a 1 ( t ), … , a N ( t )}, while the sequences of the second set are the chains simulated until time t/ 2, and we call them A ( t/ 2) = { a 1 ( t/ 2), … , a N ( t/ 2)}. With some abuse of notation, we define the intrachain correlation and autocorrelation as where σ ( i ) is a random permutation of the index i and is the normalized sequence identity (or overlap) between the sequences a and b . By construction, at the initialization we have SeqID( t, t/ 2) = 1 and SeqID( t ) somewhat close to the average sequence identity of the MSA. The convergence is obtained when chains are mixed , meaning that the system has completely forgotten the initial configuration. This requirement is satisfied when the statistics of a set of independent chains is the same as the one between the initialization and the evolved chains, meaning SeqID( t ) ∼= SeqID( t, t/ 2). The point at which the two curves merge is called mixing time , and we denote it as t mix . After reaching the mixing time of the model, the algorithm will initialize ngen chains at random. It will run a sampling for other nmix × t mix steps to guarantee complete thermalization, with nmix=2 by default. Together with the generated sequences, the script will output a text file containing the records of SeqID( t ) and SeqID( t, t/ 2) and their standard deviations ( figure 6 -Left). 4.2 Contact prediction One of the principal applications of the DCA models has been that of predicting a tertiary structure of a protein or RNA domain. In particular, with each pair of sites i and j in the MSA, adabmDCA 2.0 computes a contact score that quantifies how likely the two associated positions in the chains are in contact in the three-dimensional structure. Formally, it corresponds to the average-product corrected (APC) Frobenius norms of the coupling matrices [ 6 ], i.e. Note that the coupling parameters are usually transformed in a zero-sum gauge before computing the scores, and the gap symbol should be neglected while computing the sum in Eq. 13 [ 7 ]. The scores for all site pairs are provided in the output folder in a separate file called frobenius.txt . The first two columns indicate the site indices and the third one contains the associated APC Frobenius norm. The command for computing the matrix of Frobenius norms is Download figure Open in new tab 4.3 Scoring a sequence set At convergence, users can score a set of input sequences according to a trained DCA model by using the command line Download figure Open in new tab adabmDCA 2.0 will produce a new FASTA file in the output folder where the input sequences in have an additional field in the name that account for the DCA energy function computed according to the model in . Note that low energies correspond to good sequences. 4.4 Single mutant library Another possible application exploits the sequence-fitness score computable according to the energy function E in Eq. (1). This routine provides a single-mutant library for a given wild type to possibly guide Deep Mutational Scanning (DMS) experiments. In particular, adabmDCA 2.0 allows one to predict the fitness reduction (increase) in terms of Δ E = E (mutant) − E (wildtype) for positive (negative) value of Δ E , respectively [ 10 ]. To produce a FASTA file containing all weighted single-mutants one has to run Download figure Open in new tab where is the name of the FASTA file containing the wild type sequence, is a model file in a format compatible with adabmDCA 2.0 output, and corresponds to the folder that will contain the library file. The sequences in the output FASTA file are named after the introduced mutation and the corresponding Δ E ; for instance, >G27A | DCAscore: -0.6 denotes that position 27 has been changed from G to A and Δ E = −0.6. 5 Source Code The code and the documentation for the adabmDCA 2.0 package can be found at https://github.com/spqb/adabmDCA . We also prepared a Colab notebook with a simple tutorial for training and sampling a bmDCA model. 6 Quicklist 6.1 Commands Train a bmDCA model with default arguments: Download figure Open in new tab Restore the training of a bmDCA model: Download figure Open in new tab Train an eaDCA model with default arguments: Download figure Open in new tab Restore the training of an eaDCA model: Download figure Open in new tab Decimate a bmDCA model at 2% of density: Download figure Open in new tab Train and decimate a bmDCA model at 2% of density: Download figure Open in new tab Sample from a previously trained DCA model: Download figure Open in new tab Scoring a sequence set: Download figure Open in new tab Generating a single mutant library starting from a wild type: Download figure Open in new tab Computing the matrix of Frobenius norms for the contact prediction: Download figure Open in new tab 6.2 Script arguments Footnotes * These arguments are specific to the code implementation. References [1]. ↵ Pierre Barrat-Charlaix et al. “ Sparse generative modeling via parameter reduction of Boltzmann machines: Application to protein-sequence families ”. In: Physical Review E 104 . 2 ( Aug . 6, 2021 ), p. 024407 . doi: 10.1103/PhysRevE.104.024407 . OpenUrl CrossRef PubMed [2]. ↵ Francesco Calvanese et al. “ Towards parsimonious generative modeling of RNA families ”. In: Nucleic Acids Research 52 . 10 ( June 10, 2024 ), pp. 5465 – 5477 . issn: 0305-1048 . doi: 10.1093/nar/gkae289 . OpenUrl CrossRef PubMed [3]. ↵ Simona Cocco et al. “ Inverse statistical physics of protein sequences: a key issues review ”. In: Reports on Progress in Physics 81 . 3 ( Jan . 2018 ), p. 032601 . issn: 0034-4885 . doi: 10.1088/1361-6633/aa9965 . OpenUrl CrossRef PubMed [4]. ↵ Francesca Cuturello , Guido Tiana , and Giovanni Bussi . “ Assessing the accuracy of direct-coupling analysis for RNA contact prediction ”. In: RNA 26 . 5 ( 2020 ), pp. 637 – 647 . OpenUrl CrossRef [5]. ↵ Aurélien Decelle , Cyril Furtlehner , and Beatriz Seoane . “ Equilibrium and non-Equilibrium regimes in the learning of Restricted Boltzmann Machines ”. In: Advances in Neural Information Processing Systems . Vol. 34 . Curran Associates, Inc ., 2021 , pp. 5345 – 5359 . OpenUrl [6]. ↵ Magnus Ekeberg et al. “ Improved contact prediction in proteins: Using pseudolikelihoods to infer Potts models ”. In: Physical Review E 87 . 1 ( Jan . 11, 2013 ), p. 012707 . doi: 10.1103/PhysRevE.87.012707 . OpenUrl CrossRef PubMed [7]. ↵ Christoph Feinauer et al. “ Improving Contact Prediction along Three Dimensions ”. In: PLOS Computational Biology 10 . 10 ( Oct . 9, 2014 ), e1003847 . issn: 1553-7358 . doi: 10.1371/journal.pcbi.1003847 . OpenUrl CrossRef PubMed [8]. ↵ Matteo Figliuzzi , Pierre Barrat-Charlaix , and Martin Weigt . “ How Pairwise Coevolutionary Models Capture the Collective Residue Variability in Proteins? ” In: Molecular Biology and Evolution 35 . 4 ( Apr . 1, 2018 ), pp. 1018 – 1027 . issn: 0737-4038 . doi: 10.1093/molbev/msy007 . OpenUrl CrossRef [9]. ↵ Matteo Figliuzzi et al. “ Coevolutionary landscape inference and the context-dependence of mutations in betalactamase TEM-1 ”. In: Molecular biology and evolution 33 . 1 ( 2016 ), pp. 268 – 280 . OpenUrl CrossRef [10]. ↵ Thomas A. Hopf et al. “ Mutation effects predicted from sequence co-variation ”. In: Nature Biotechnology 35 . 2 ( Feb . 2017 ), pp. 128 – 135 . issn: 1546-1696 . doi: 10.1038/nbt.3769 . OpenUrl CrossRef PubMed [11]. ↵ Faruck Morcos et al. “ Direct-coupling analysis of residue coevolution captures native contacts across many protein families ”. In: Proceedings of the National Academy of Sciences 108 . 49 ( 2011 ), E1293 – E1301 . OpenUrl [12]. ↵ Anna Paola Muntoni et al. “ adabmDCA: adaptive Boltzmann machine learning for biological sequences ”. In: BMC Bioinformatics 22 . 1 ( Oct . 29, 2021 ), p. 528 . issn: 1471-2105 . doi: 10.1186/s12859-021-04441-9 . OpenUrl CrossRef PubMed [13]. ↵ WR Pearson and DJ Lipman . “ Improved tools for biological sequence comparison .” In: Proceedings of the National Academy of Sciences of the United States of America 85 . 8 ( Apr . 1988 ), pp. 2444 – 2448 . issn: 0027-8424 . OpenUrl [14]. ↵ William P Russ et al. “ An evolution-based model for designing chorismate mutase enzymes ”. In: Science 369 . 6502 ( 2020 ), pp. 440 – 445 . OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted February 05, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis Lorenzo Rosset , Roberto Netti , Anna Paola Muntoni , Martin Weigt , Francesco Zamponi bioRxiv 2025.01.31.635874; doi: https://doi.org/10.1101/2025.01.31.635874 Share This Article: Copy Citation Tools adabmDCA 2.0 – a flexible but easy-to-use package for Direct Coupling Analysis Lorenzo Rosset , Roberto Netti , Anna Paola Muntoni , Martin Weigt , Francesco Zamponi bioRxiv 2025.01.31.635874; doi: https://doi.org/10.1101/2025.01.31.635874 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Bioinformatics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41968) Biophysics (21460) Cancer Biology (18603) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19909) Epidemiology (2067) Evolutionary Biology (24326) Genetics (15614) Genomics (22513) Immunology (17741) Microbiology (40423) Molecular Biology (17193) Neuroscience (88645) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4825) Physiology (7647) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)