Protein Direct Coupling Analysis (DCA), which predicts residue-residue contacts based on covarying positions within a multiple sequence alignment, has been remarkably effective. This suggests that there is more to learn from sequence correlations than is generally assumed, and calls for deeper investigations into DCA and perhaps into other types of correlations. Here we describe an approach that enables such investigations by measuring, as an estimated p-value, the statistical significance of the association between residue-residue covariance and structural interactions, either internal or homodimeric. Its application to thirty protein superfamilies confirms that direct coupling (DC) scores correlate with 3D pairwise contacts with very high significance. This method also permits quantitative assessment of the relative performance of alternative DCA methods, and of the degree to which they detect direct versus indirect couplings. We illustrate its use to assess, for a given protein, the biological relevance of alternative conformational states, to investigate the possible mechanistic implications of differences between these states, and to characterize subtle aspects of direct couplings. Our analysis indicates that direct pairwise correlations may be largely distinct from correlated patterns associated with functional specialization, and that the joint analysis of both types of correlations can yield greater power. Data, programs, and source code are freely available at http: //evaldca. igs. umaryland. edu.
Contacts among residues largely determine a protein’s three-dimensional structure. For proteins sharing a common structure, such contacts generally produce correlated substitution patterns between residue pairs. Over evolutionary time substitutions at one residue position often result in compensating substitutions at other positions in order to maintain critical interactions. This allows the prediction of protein structural contacts based upon multiple sequence alignment (MSA) covariance analysis. Early approaches were only partially successful, with a major shortcoming the confounding effect of indirect correlations: When residues at positions i and j correlate, as do those at positions j and k, then residues at positions i and k may also correlate even though they fail to interact directly. Direct Coupling Analysis (DCA) and related methods [1–8] have overcome this problem by disentangling direct correlations from indirect coupling effects. As used here, the term DCA refers to all such approaches. DCA constitutes a major breakthrough in protein structure prediction and is currently being applied successfully on a large scale [9]. DCA programs employ a variety of algorithmic strategies, including sparse inverse covariance estimation (PSICOV) [4], pseudo-likelihood maximum entropy optimization (EVcouplings-PLM) [5,6] (CCMpred) [10] and multivariate Gaussian modeling (GaussDCA) [11]. DCA methods are evaluated by comparing those residue pairs with the highest direct coupling (DC) scores to residue-to-residue contacts within protein structures. Currently this involves using, for example, ROC curves [11], the Matthews correlation coefficient [12], F1 scores, or often the positive prediction value (PPV). Such measures are applied by labeling data points according to a binary classification scheme; for DCA, those residue pairs that are a specified distance apart within a benchmark structure (e. g. , ≤ 5 Å) are labeled as positives and other pairs as negatives. However, there are reasons to criticize such measures in particular circumstances [13]. In particular, it is not clear how to assess the significance of such measures when comparing different proteins or distinct structures. To standardize such comparisons, it is desirable to obtain a measure of statistical significance, which also provides insight into how surprised we should be with a given result. As illustrated here, one can use such a measure to determine whether it is better to base DC scores on an MSA of more closely related proteins rather than on an entire superfamily MSA. Given a set of structures for a protein superfamily, a significance measure can help identify those of greatest interest: Direct couplings between pairs of residues presumably are due to selective constraints maintaining functionally important structural interactions. Hence, those protein structures that exhibit the most biologically relevant interactions should achieve the highest level of significance. One could therefore use a significance measure to select among alternative structural models generated by homology or by ab initio structure prediction methods. One may also adapt such a measure to evaluate the degree to which high DC scores are associated with properties other than 3D structural contacts. As illustrated here, for example, one may determine whether those residues most distinctive of a particular protein family are overrepresented among the highest DC-scoring residue pairs. Here we describe a method to estimate, in various contexts, the statistical significance of the correspondence between DC scores and either protein structural contacts or other protein properties. Unlike the current practice of selecting for analysis an arbitrary number of the highest scoring pairs (e. g. , 1. 5 times the MSA length [14]), our approach determines the optimal number of such pairs automatically based on a statistical criterion, while adjusting automatically for the number of multiple hypotheses tested. Unlike binary classification schemes, our approach takes into account both the order of each residue pair based on DC scores and their ranks based on 3D pairwise distances; hence, it treats the structurally closest residue pairs having high DC scores as of higher biological relevance than such pairs having low DC scores. By providing a quantitative measure of significance, our approach can detect subtle yet important features of the data that qualitative measures would fail to distinguish from background noise. We illustrate this approach by investigating: the relative performance of alternative methods; the biological relevance of alternative structures; subtle structural changes associated with the transition state of Ran GTPase; the contribution of homo-oligomer interfaces to aggregate DC scores; DCA’s dependence on the sequences included in the input MSA; and the correspondence between DCA pairwise correlations and correlated patterns associated with protein functional specialization.
Abstractly, given an array of elements ordered by a primary criterion (e. g. , as used here, DC scores), we ask how well it agrees with a secondary criterion (e. g. , 3D pairwise distances) that distinguishes and ranks a subset of the elements. More specifically, we seek to identify an optimal initial cluster of elements of the array (defined by a cut), as measured by a relevant p-value. Our approach is based upon Initial Cluster Analysis (ICA) [15]; see Methods. For reference, Table 1 provides a summary of the variables used below. ICA answers the question: Given a random array of length L, containing D' 1' s (representing distinguished elements), and L—D' 0' s, what initial cluster, consisting of elements up to and including a cut point X, contains the most surprising number d of' 1' s, and what is its probability of occurring? (Below, we call the d' 1' s in an initial cluster “left-distinguished elements. ”) For L = 18 and D = 7, for example, one such array is “101101100000010001”, with optimal cut point X = 7 (underlined), yielding d = 5. Here we note that, in practice, to distinguish elements within our array, we frequently rank all the elements, and distinguish those with rank ≤ D. We then might denote our example array as “401603200000070005” with digits > 0 denoting the ranks of distinguished elements. ICA ignores these ranks when choosing the optimal X, whereas we would prefer the d distinguished elements to the left of X to have superior ranks (i. e. , lower numbers) than those to the right. To generalize ICA to exploit ranking information we incorporate a ball-in-urn model to calculate a ranking specific p-value Pb. For a specific cut point X that yields d left-distinguished elements, we imagine first coloring red, among all D distinguished elements, those d elements with superior ranks (e. g. , with the smallest pairwise distances); and then recording the number R that are red among the left-distinguished elements. Ideally, all the left-distinguished elements will outrank the remaining distinguished elements, yielding R = d, but more generally higher values of R are better; in the example of the previous paragraph, D = 7, d = 5 and R = 4. Given the null hypothesis that rankings are random, we may then use the cumulative hypergeometric distribution to calculate the probability Pb that ≥ R of the left-distinguished elements are red: Pb=[∑i=Rd (di) (D−dd−i) ]÷ (Dd). This corresponds to drawing d balls from an urn containing D balls, of which d are red; note that the number of balls drawn here equals the number colored red. A low value of Pb is reported for a cut with a surprising number, among its d left-distinguished elements, having the d smallest pairwise distances. Before it corrects for optimizing over all possible cuts, ICA can be understood as calculating a p-value Pa for finding d distinguished elements to the left of a cut point X. Because the calculation of Pa ignores ranking information, it will be independent of Pb, and these two p-values may therefore be combined to yield a joint p-value PJ [16–19] using the formula PJ=PaPb (1−lnPaPb). Low values of PJ may arise from low values of Pa, or Pb, or of both. PJ can provide a statistically stronger measure than Pa alone of the congruence of two orderings, here derived from DC scores and 3D distances. The p-values P we report in this paper correspond to PJ, after it has been corrected for optimization over the multiple cut points X considered, as described in [15]. One may wish to optimize as well over various values of D, but in the current application larger values of D are then almost always preferred, due to the indirect couplings considered below. We therefore choose a fixed D, based upon a maximum allowed 3D distance within a reference structure. To summarize, in order to apply the theory above to the question of how well DCA actually uncovers direct contacts within proteins, we proceed as follows. Given an MSA, a method to calculate DC scores for all column pairs, and a reference structure corresponding to one of the sequences in the MSA, we consider only those pairs of MSA columns separated by ≥ m intervening positions within the reference sequence, with m = 5 by default. Ordering these column pairs by descending DC score yields our array of elements, of length L. We then distinguish those D elements whose 3D distance per the reference structure is ≤ r Å, with r = 5 by default, and rank them by increasing distance. (This distance is defined as the minimum between sidechain atoms, including hydrogens, of the paired residues. For glycine, Cα and its attached hydrogen serve as the sidechain atoms.) ICA' s original Pa depends only upon the specification of these pairs as distinguished, whereas PJ takes account as well of their rankings, through the ball-in-urn derived Pb. As we will show below, for this application Pa is, in general, far smaller than Pb. However, we have found PJ to provide, in general, greater statistical power than Pa for analyzing protein sequence-structural relationships, and our focus in this paper is to illustrate its use. Currently DCA performance is often evaluated using the positive prediction value (PPV), defined as the percentage of observed reference-structure 3D contacts corresponding to a fixed number x (e. g. , x = 100 [11,20]) of the highest DC-scoring column pairs. In contrast, the cut point X is not fixed but chosen to optimized significance. Because the number of column pairs grows with increasing MSA length ℓ, x is often chosen, using a parameter F, as x = F × ℓ. Typical values of F range from 0. 5 to 1. 5 [10,14]. Since we propose the S-score as a replacement for PPV, we compare these two metrics below in several ways. To aid these comparisons, we define Px as the probability, based on the cumulative hypergeometric distribution, of d being at least the value observed for a constant value of x = F × ℓ. We define PF as the estimated joint p-value that combines Pb and Px: PF=PxPb (1−lnPxPb) =PF×lPb (1−lnPF×lPb). Implementation and availability. We implemented these algorithms and statistical models in C++ as the STARC (Statistical Tool for Analysis of Residue Couplings) program, which, along with the source code, is freely available at http: //evaldca. igs. umaryland. edu. Here and below, for an estimated or theoretical p-value P we define a corresponding s-score as S = −log10 P. Our theory should yield accurate p-values and s-scores for randomly generated, or shuffled arrays. However, in the present application many column pairs within an MSA are interrelated (e. g. , {i, j}, {j, k} and {i, k}), possibly affecting their DC scores as well as the corresponding distances derived from a structure. To test whether computed p-values remain valid given such interrelationships we generated, based on randomization of each of six MSAs with six corresponding structures, sets of random p-values, as described in Methods. We define Ŝ, as a function of S, to be -log10 of the proportion of observed (simulation-based) s-scores that are greater than or equal to S. If our p-value calculations are accurate, Ŝ should equal S to within stochastic error. In Fig 1 we plot, for S from 2 to 5, the Ŝ obtained from 100,000 p-values for MSAs in which the residues in each column and the order of the columns for each MSA were randomly permuted (termed column-permuted MSAs). This operation retains the distribution of column relative entropies observed in the original MSA. For comparison, we plot as well the Ŝ obtained from an equivalent number of shuffled DC arrays (see Methods). These arrays abolish interrelationships among DC-scores, and so conform better to theory. The straight, solid line represents the agreement of Ŝ with theory, and dashed curves represent error ranges of two standard deviations. As can be seen, within stochastic error, Ŝ agrees well with theory for the shuffled arrays. (Because we can generate p-values rapidly for shuffled arrays, we have confirmed the accuracy of Ŝ in this case for S ≤ 8.) For the MSA with the largest number of sequences (146,217) among the six used for column-permuted simulations, and corresponding to structure 1wznA, the values of Ŝ deviate consistently below the error bounds corresponding to slightly inflated s-scores. When we randomly removed all but 5,117 (3. 5%) of the sequences in this MSA, however, this effect was essentially eliminated (see 1wznA plot in Fig 1). Values of Ŝ for the other column-permuted MSA simulations exhibit less of a tendency to deviate outside of the error bounds. Based on these examples, it appears that for large alignments, s-scores may be slightly inflated, but it is not clear why this is so. We cannot be sure that values of Ŝ for randomized MSAs will conform to theory beyond the range tested. However, we may apply s-scores in a manner similar to that of Z-values. A Z-value is the distance between a raw score and the population mean in units of standard deviation. One may convert a Z-value into a p-value under the assumption (based on the Central Limit Theorem) that the variables are drawn from a normal distribution. Although this assumption is typically invalid for raw scores far away from the mean, Z-values still provide a useful metric for assessing significance. Extreme s-scores likewise provide a useful measure of statistical significance even though the true distribution may depart from the theoretical distribution used here. We ran the STARC program on the output from four DCA programs, EVcouplings (EVC) [5,6], GaussDCA with Frobenius norm ranking (GSF) [11], PSICOV (PCV) [4], and CCMpred (CCM) [10], each applied to thirty protein domain MSAs with reference 3D contacts ≤ 5 Å (Table 2). For a given MSA, better performing programs should typically generate more significant results, and thus generally higher s-scores. We have studied, through the score S, the correspondence between a multiple alignment' s DC scores and the pairwise distances implied by the structure for a particular sequence in the alignment. However, to calculate S, there are typically many structures to choose among, and these may differ in important particulars. Recent studies [27–32] have demonstrated that high DC scoring pairs that are distant in certain benchmark 3D structures may come into contact within alternative conformations or across homo-oligomer interfaces, and have thereby provided insight into protein biophysical and dynamic properties. Other studies [33,34] have combined DCA with correlation analyses involving larger groups of structurally and/or functionally correlated residues, thereby generating further insight. Here we illustrate the application of our method to these sorts of studies. To the degree to which DC scores capture the pairwise correlations imposed by the functional requirements common to a protein family, we expect the S yielded by a particular structure to reflect the degree to which that structure exhibits critical interactions characteristic of the family. In other words, S may measure the degree to which a specific structural conformation is biologically relevant. To investigate this, we consider three cases—human Ran GTPase, Gna1 N-acetyltransferase from C. elegans, and the bacterial (E. coli) clamp loader complex. Using available structures for each of these, we add hydrogen atoms using the Reduce program [21] to better discriminate among residue-to-residue contact distances. A previous DCA analysis [31] found that the heavy atom distance distribution for directly coupled residue pairs exhibited local maxima at 2. 8 Å and 3. 7 Å, which were interpreted as corresponding to the donor-acceptor distance of hydrogen bonds and to hydrophobic interactions, respectively. Here we choose to focus on hydrogen bond interactions. Since our analyses explicitly model hydrogen atoms, we calculate S using a maximum structural distance of 2. 6 Å, which, based on the sum of the van der Waals radii for hydrogen plus either nitrogen or oxygen [35], corresponds to an upper bound on the hydrogen-acceptor distance of hydrogen bonds.
STARC-computed s-scores quantify a DCA method’s ability to detect 3D residue-to-residue contacts. When used in combination with the Wilcoxon signed rank test, they yield a significance measure of the performance of one method versus another and can quantify, as well, a method’s tendency to detect indirect couplings. Larger domains tend to yield higher s-scores due to enhanced statistical power. However, this is not a confounding factor when s-scores are used to compare different DCA methods applied to the same MSA and 3D structure or to compare alternative structures based on the same DCA method and MSA. The following advantages of S-scores over SF-scores and PPV may be noted: S-scores are not biased toward any particular method, but rather correspond to the optimal value of F for each method and reference structure under consideration; this may reveal important features that would otherwise be overlooked. S-scores avoid ranking inconsistencies due to one’s choice of fixed values of F. S-scores can tap into additional information regarding each structure’s possible biological relevance, as illustrated here. S-scores take into consideration not only the number and arrangement of (false negative) contacting residue pairs to the right of a potential cut point, but also the ordering of those pairs based on their 3D distances within a reference structure. And, unlike PPV, S-scores provide a measure of statistical significance. Of course, researchers also have the option of computing SF, thereby obtaining both a measure of significance and an assessment of program performance as functions of F. We could further develop the STARC statistical model by considering the arrangement of the d distinguished pairs before X. A pair with a higher DC score should be more likely than one with a lower DC score to correspond to a 3D interaction. Ideally, the d pairs should thus be arranged in order of decreasing DC score. To measure how closely a DCA method’s output comes to achieving this configuration, we may first define a permutation π by ranking the d distinguished pairs based on 3D distance, with smaller distances receiving superior ranks (i. e. , lower numbers), and then define τ (π) =12∑i=1d (π[i]−i) 2. One may show that τ is an integer function that for random permutations is symmetrically distributed about its mean μ = (d3−d) /12, with standard deviation σ=μ/d−1. For d ≤ 16 one may compile exact p-values for τ by exhaustive enumeration, and for d > 16 estimate them using either a Gaussian approximation or Monte Carlo simulation. However, it is unclear whether these p-values are independent of Pa and Pb, and whether there is biological benefit to including this order in our statistical model. We plan to investigate these questions. An important potential application of our approach, which is beyond the scope of this study, is the evaluation of MSA accuracy without the need for benchmark alignments, which typically contain a relatively small number of sequences and whose accuracy may be uncertain [46]. Our proposed approach would proceed on the assumption that, given available structures, more accurate MSAs will yield higher values of S. We are developing this approach, which should benefit from the large amount of sequence and structural data becoming available. Our analysis of Ran, Gna1 and the DNA clamp loader complex suggests that S may be useful for evaluating the biological relevance of alternative structural conformations of the same protein and for characterizing the nature of conformation-specific interactions. Viewing direct couplings as functionally imposed constraints and proteins as molecular machines, S may measure the degree to which a particular crystal structure captures a protein in a mechanistically important state. If so, then analyzing in what ways various residue pairs contribution to S may provide mechanistic clues. Likewise, comparative analyses among MSAs corresponding to a protein’s subfamily, family and superfamily may provide mechanistic clues regarding functional specialization. Our analysis here also suggests that one may use STARC to search for the most biologically relevant among the many structures often available for a major protein superfamily.
For the thirty STARC analyses in Table 2, we obtained high quality crystal structures from the Protein Data Bank (PDB) (www. rcsb. org/pdb). The pdb and chain identifiers are given in column 1 of Table 2. Likewise, the coordinates for the Ran, Gna1 and DNA clamp loader analyses were obtained from the PDB; their pdb identifiers are given in Tables 4,5 and 6, respectively. For all analyses, hydrogen atoms were added using the Reduce program [31], except for the pdb coordinate file for 3f1lA in which hydrogens were already present. Hence, residue-to-residue distances are based on any two atoms, including hydrogens, albeit ignoring main chain to main chain interactions. This allows better discrimination among hydrogen bond interactions based on subtle differences in contact distances. EVcouplings (EVC) was run over the EVfold website (http: //evfold. org) using the pseudo-likelihood maximization (PLM) option with default settings. For each analysis, taking as input the sequence corresponding to the reference structure as the query, EVcouplings uses jackhmmer [45] to create an MSA, from which it then computes the DC scores. The score file and the corresponding PDB coordinates serve as the input to STARC. We also used the jackhmmer alignment as input to the other programs. The GaussDCA program was run with Frobenius norm ranking (with default parameters); this was done interactively under Julia (www. julialang. org). PSICOV version 2. 4 was run using the author recommended –p and –d 0. 03 options and the jackhmmer alignment reformatted by the fasta2aln program, which is included with the PSICOV package (http: //bioinf. cs. ucl. ac. uk/downloads/PSICOV). CCMpred version 0. 3. 2 (https: //travis-ci. org/soedinglab/CCMpred) was run with default settings again using the reformatted alignment. Note that the output from GaussDCA, CCMpred and PSICOV does not include the query sequence, which, along with the DC scores, were provided as input to STARC. We performed two types of simulations for each of the six MSAs in Fig 1, which are labeled by their corresponding pdb identifiers, 3fhkF, 1ijxA, 4cmlA, 1k30A, 1olzA, 1wznA, and which correspond to analyses in Table 2. These MSAs vary substantially in their numbers of aligned columns (127 to 481) and aligned sequences (1,042 to 146,217), and in the degree of shared sequence similarity. For the first type of simulation, we randomly shuffled the DC score array for each of 100,000 runs. This simulation corresponds to the theory behind the ICA algorithm, which is described in the next section. The second type of simulation corresponds more closely to a STARC analysis by computing a DC score array from a simulated MSA. For each MSA, we first randomly permuted the residues in each aligned column and then randomly permuted the order of the columns in the MSA (termed a column-permuted MSA). Next, using these simulated MSAs as input to the CCMpred program [10], DC scores were computed for each of 100,000 runs. Finally, for each run, STARC was applied using as input the DC scores and the corresponding protein structure. We describe Initial Cluster Analysis (ICA) in detail elsewhere [15], but summarize the approach briefly here. ICA seeks to determine whether a set of distinguished elements within a linear array is clustered significantly near the start of the array and, if so, what is the most significant initial cluster of these elements. Abstractly, given a linear array of length L containing D' 1' s (the distinguished elements) and L-D' 0' s, it considers a generative model in which the' 1' s occur with particular and differing probabilities before and after a cut point X in the array. For any particular X it is relatively easy to calculate a likelihood L (X) of the array of data, and one may optimize L (X) by simply evaluating it for all possible X. However, the values of L (X) for close values of X are highly correlated, dependent upon a calculable" density of independent trials" ρ (X). Because ρ (X) is not constant but rather grows approximately as the reciprocal of X' s distance from 0 or L, simply optimizing L (X) inherently favors, a priori, small or large values of X. Therefore, if one' s application suggests no such bias, choosing to optimize L (X) /ρ (X) rather than L (X) for a given array of' 0' s and' 1' s may be a better strategy. This is referred to in [15] as using" flattened priors" , and is the approach we take here. ICA estimates the effective total number of independent trials implicit in either optimization, which it uses in calculating a p-value for the optimal X from its corresponding Pa. This provides a mathematically principled way to define an optimal initial cluster of distinguished elements, balancing the claims of very short and dense clusters with those of longer but sparser clusters. We have extended ICA here by taking account not only of distinguished elements within an array, but of a ranking assigned to these elements as well. Thus we seek here initial clusters not only with a high density of distinguished elements, but clusters in which these elements have relatively better rankings. Our s-score may be understood as providing a measure of the congruence between two orderings, as well as, simultaneously, an assessment of statistical significance. We evaluated the performance of alternative DCA methods using the Wilcoxon signed-rank test [22], first dividing each S by the total number of residue pairs L. For CCMpred, EV-couplings and GaussDCA, these normalized s-scores then approximately follow a Gaussian distribution, as indicated by the Shapiro-Wilk test statistic [47] (p = 0. 52,0. 60, and 0. 09, respectively). For PSICOV the Shapiro-Wilk test score corresponded to p = 0. 04, which is slightly below the acceptance threshold of p > 0. 05. The STARC program uses a modified version of the Initial Cluster Analysis (ICA) algorithm [15] to find the optimal score S, as described above. Alternatively, as an option, it will calculate SF given a specified F. STARC converts PSICOV and GaussDCA formatted DC score files into EVcouplings format automatically; this requires as input the query sequence in fasta2aln format. We modified the CCMpred source code and recompiled the program to generate PSICOV-formatted DC score files. The source code for STARC is freely available at: http: //evaldca. igs. umaryland. edu/. BPPS [39,48] partitions the sequences in a superfamily MSA into families and subfamilies. It uses Markov chain Monte Carlo (MCMC) sampling to stochastically move sequences between subgroups, while modifying each subgroup’s characteristic pattern. BPPS also identifies and removes unrelated or aberrant sequences. We applied BPPS here to generate both family and subfamily MSAs for sequences of interest. Here we also use STARC to assess the correspondence between pairs of BPPS-defined pattern residues and high DC-scoring pairs.
The success of Direct Coupling Analysis (DCA) for protein structure prediction suggests that multiple sequence alignments implicitly contain more structural information than had previously been realized, and prompts deeper investigations of the sequence correlations uncovered by either DCA or other approaches. To aid such investigations and thereby broaden the utility of and improve DCA, we describe an approach that measures the statistical significance of the association between DCA and either 3D structure or correlated patterns associated with functional specialization. This approach can be used to evaluate the relative performance of DCA methods, their ability to distinguish direct from indirect couplings, and the potential biological relevance and mechanistic implications of alternative conformations and homodimeric interactions.
quantum chemistry
crystal structure
enzymes
split-decomposition method
dna clamps
condensed matter physics
enzymology
multiple alignment calculation
protein structure
dna
dna structure
crystallography
ground state
research and analysis methods
sequence analysis
solid state physics
sequence alignment
bioinformatics
proteins
structural proteins
guanosine triphosphatase
chemistry
molecular biology
physics
biochemistry
hydrolases
computational techniques
nucleic acids
database and informatics methods
genetics
biology and life sciences
physical sciences
macromolecular structure analysis
Statistical investigations of protein residue direct couplings