Accurate Simulation and Detection of Coevolution Signals in Multiple Sequence Alignments

Background While the conserved positions of a multiple sequence alignment (MSA) are clearly of interest, non-conserved positions can also be important because, for example, destabilizing effects at one position can be compensated by stabilizing effects at another position. Different methods have been developed to recognize the evolutionary relationship between amino acid sites, and to disentangle functional/structural dependencies from historical/phylogenetic ones. Methodology/Principal Findings We have used two complementary approaches to test the efficacy of these methods. In the first approach, we have used a new program, MSAvolve, for the in silico evolution of MSAs, which records a detailed history of all covarying positions, and builds a global coevolution matrix as the accumulated sum of individual matrices for the positions forced to co-vary, the recombinant coevolution, and the stochastic coevolution. We have simulated over 1600 MSAs for 8 protein families, which reflect sequences of different sizes and proteins with widely different functions. The calculated coevolution matrices were compared with the coevolution matrices obtained for the same evolved MSAs with different coevolution detection methods. In a second approach we have evaluated the capacity of the different methods to predict close contacts in the representative X-ray structures of an additional 150 protein families using only experimental MSAs. Conclusions/Significance Methods based on the identification of global correlations between pairs were found to be generally superior to methods based only on local correlations in their capacity to identify coevolving residues using either simulated or experimental MSAs. However, the significant variability in the performance of different methods with different proteins suggests that the simulation of MSAs that replicate the statistical properties of the experimental MSA can be a valuable tool to identify the coevolution detection method that is most effective in each case.

only some parts -a branch or a specific group of sequences -of the evolutionary tree.
In particular, the automatic branch assignment of MSAvolve can satisfy a user's request to identify clusters in the experimental MSA using more than 2 principal components (PCs) of the covariance or other type of distance matrix ( Figure S9).
MSAvolve generates MSAs that match very well the overall similarity score (OSS, see below) of the experimental MSA, regardless of whether the ancestor is created using the background frequencies of amino acids in the family or the HMM emission probabilities. In addition, MSAvolve builds coevolution matrices of each simulated MSA by counting the mutations that segregate at specific times during the evolution of a protein. This function can be equated to that of an external observer that witnesses the historical process of evolution. There are two main observers in MSAvolve: one is a 'low power' observer, who counts all the co-segregating mutations regardless of whether they originate from chance, recombination, or functional/structural demands (the latter are the mutations under constraint of coevolution). There is also a 'high power' observer, who counts selectively the true co-segregating mutations of functional/structural significance.
Although it is possible to build ancestor sequences of arbitrary length, in most cases, the MSA simulation starts with reading the experimental MSA of a real protein family and converting it to a numeric format (each aa is represented by a number between 1 and 20, gaps are represented as number 21). The total length (including gaps) of each aligned sequence in the real MSA defines the length of each sequence in the simulated MSA. First, a Hidden Markov Model (HMM) of the real MSA is built: the HMM provides the emission probabilities and background frequencies of the 20 amino acids at each position in the MSA. A non-parametric kernel smoothing distribution with a smoothing bandwidth of 0.01 is fitted to the vector of emission probabilities for each position of the MSA. An object representing the fitted distribution is stored for each MSA position, and an ancestor sequence is built by random drawing (rounded to the nearest integer) from the probability distribution of each position. Alternatively, MSAvolve can be prompted to produce an ancestor also from the equilibrium frequencies.
Next, a phylogenetic tree built from the matrix of pairwise distances is used to determine the number of main branches in the real MSA, which will be replicated in the simulated MSA (an arbitrary number of branches can however be set manually). Several options are available to calculate pairwise distances and trees, as offered in the Matlab Bioinformatics Toolbox (which is required for MSAvolve to function). Alternatively, branches can be determined from a spectral analysis of the covariance matrix of the sequences in the MSA, as originally described by Halabi et al. [1]. The covariance matrix provides also a convenient way to define an overall similarity score (OSS) between all the sequences in the MSA using the definition of multidimensional mutual information, I, developed by Wang and Shen [2] to express the similarity between the images in a set:  increases C stochastic with respect to C structure/function . On the other end, cycles cannot be made arbitrarily small, because it would imply that no mutations are ever deleterious (any mutation is fixed), and in the limit, from the practical standpoint of how the programs works, all mutations would appear to be the product of functional or structural coevolution. This effect can be exploited by the user to test the 'strength' of a coevolution detection method. Weaker methods are expected to detect progressively fewer true covarions than stronger methods as the mutational rates increases, and the number of cycles decreases. Furthermore, since the mean entropy of the columns in the MSA grows monotonically with the overall mutation rate, to simulate data that are analogous to a given real MSA, it is useful to find a value of the product "mutation_rate x mutation_cycles" that gives an entropy similar to that of the real data. However, as previously discussed, while the final MSAs obtained maintaining this product constant will have the same column entropy, the corresponding coevolution matrices will have significantly different levels of C stochastic .
MSAvolve alternates rounds of mutations (consisting of multiple mutation cycles) with rounds of recombination. In the beginning of our example recombination occurs between the three branches of the MSA (Figure S11;  [3][4][5][6] can be used to predict which fragments of homologous proteins can be recombined without disturbing the integrity of the structure [7].

8
Once the cross-over point are defined, MSAvolve carries out recombination in two steps: 1. It selects at random the number and identity of the sequences that will receive the donated fragment.
2. It selects at random which fragment will be donated and which sequence will donate the fragment to the other sequences.  [5,6,7,8,9,10] may have spread from sequence 3 to sequence 1 and 2. If this spread produces a change at positions [5,6] in sequence 1 and positions [5,6,9] in sequence 2, the recombinant coevolution matrix will report an increase by 2 at indices [5,6] and an increase by 1 at indices [5,9] and [6,9].