Protein–protein interactions are fundamental to many biological processes. Experimental screens have identified tens of thousands of interactions, and structural biology has provided detailed functional insight for select 3D protein complexes. An alternative rich source of information about protein interactions is the evolutionary sequence record. Building on earlier work, we show that analysis of correlated evolutionary sequence changes across proteins identifies residues that are close in space with sufficient accuracy to determine the three-dimensional structure of the protein complexes. We evaluate prediction performance in blinded tests on 76 complexes of known 3D structure, predict protein–protein contacts in 32 complexes of unknown structure, and demonstrate how evolutionary couplings can be used to distinguish between interacting and non-interacting protein pairs in a large complex. With the current growth of sequences, we expect that the method can be generalized to genome-wide elucidation of protein–protein interaction networks and used for interaction predictions at residue resolution.https://doi.org/10.7554/eLife.03430.001
DNA is often referred to as the ‘blueprint of life’, as this molecule contains the instructions that are required to build a living organism from a single cell. But these instructions largely play out through the proteins that DNA encodes; and most proteins do not work alone. Instead they come together in different combinations, or complexes, and a single protein may participate in many complexes with different activities.
Proteins are so small that it is difficult to get clear information about what they look like. Visualizing protein complexes is even harder. Most protein–protein interactions remain poorly understood, even in the best-studied organisms such as humans, yeast, and bacteria.
Proteins are made from smaller molecules, called amino acids, strung together one after the other. The order in which different amino acids are arranged in a protein determines the protein’s shape and ultimately its function. Like DNA, protein sequences can change over time. Sometimes, the sequence of one protein changes in a way that prevents it binding to another protein. If these two proteins must work together for an organism to survive, the second protein will often develop a compensating change that allows the protein–protein complex to reform.
Identifying pairs of changes in the sequences of pairs of proteins suggests that the two proteins interact and gives some information about how the proteins fit together. Different species can have copies of the same proteins that have slightly different sequences. Since the DNA sequences from many different organisms are already known, there are now many opportunities to find sites in pairs of proteins that have evolved together, or co-evolved, over time.
To find sites that seem to have co-evolved, Hopf et al. used a computer program based on an approach from statistical physics to look at pairs of proteins that were already known to form complexes. Co-evolving sites were found in over 300 pairs of proteins; including 76 where the structure of the complex was already known. When sites that were predicted to be co-evolving were then mapped to these known complex structures, the co-evolving sites were remarkably close to the true protein–protein contacts. This indicates that the information from the co-evolved sequences is sufficient to show how two proteins fit together.
Hopf et al. then turned their attention to 82 pairs of proteins that were thought to interact, but where a structure was unavailable. For 32 of these pairs, structures of the entire complex could be predicted, showing how the two proteins might interact. Furthermore, when other researchers subsequently worked out the structure of one of these complexes, the prediction was a good match to the solved complex structure.
The machinery of life is largely made up of proteins, which must interact in ever-changing but precise ways. The new methods developed by Hopf et al. provide a new way to discover and investigate the details of these interactions.https://doi.org/10.7554/eLife.03430.002
A large part of biological research is concerned with the identity, dynamics, and specificity of protein interactions. There have been impressive advances in the three-dimensional (3D) structure determination of protein complexes which has been significantly extended by homology-inferred 3D models (Mosca et al., 2012; Webb et al., 2014) (Hart et al., 2006; Zhang et al., 2012). However, there is still little, or no, 3D information for ∼80% of the currently known protein interactions in bacteria, yeast, or human, amounting to at least ∼30,000/∼6000 incompletely characterized interactions in human and Escherichia coli, respectively (Mosca et al., 2012; Rajagopala et al., 2014). With the rapid rise in our knowledge of genetic variation at the sequence level, there is an increased interest in linking sequence changes to changes in molecular interactions, but current experimental methods cannot match the increase in the demand for residue-level information of these interactions. One way to address the knowledge gap of protein interactions has been the use of hybrid, computational–experimental approaches that typically combine 3D structural information at varying resolutions, homology models, and other methods (de Juan et al., 2013), with force field-based approaches such as RosettaDock, residue cross-linking, and data-driven approaches that incorporate various sources of biological information (Kortemme and Baker, 2002; Dominguez et al., 2003; Kortemme and Baker, 2004; Kortemme et al., 2004; Svensson et al., 2004; Chaudhury et al., 2011; Schneidman-Duhovny et al., 2012; Velazquez-Muriel et al., 2012; Karaca and Bonvin, 2013; Rodrigues et al., 2013; Webb et al., 2014). However, most of these approaches depend on the availability of prior knowledge and many biologically relevant systems remain out of reach, as additional experimental information is sparse (e.g., membrane proteins, transient interactions, and large complexes). One promising computational approach is to use evolutionary analysis of amino acid co-variation to identify close residue contacts across protein interactions, which was first used 20 years ago (Gobel et al., 1994; Pazos and Valencia, 2001), and subsequently used also to identify protein interactions (Pazos et al., 1997; Pazos and Valencia, 2002). Others have used some evolutionary information to improve a machine learning approach to developing docking potentials (Faure et al., 2012; Andreani et al., 2013; Andreani and Guerois, 2014). These previous approaches relied on a local model of co-evolution that is less likely to disentangle indirect and therefore incorrect correlations from the direct co-evolution, as has been described in work on residue–residue interactions in single proteins (Marks et al., 2012). More recently, reports using a global model have been successful in identifying residue interactions from evolutionary co-variation, for instance between histidine kinases and response regulators (Burger & van Nimwegen, 2008; Skerker et al., 2008; Weigt et al., 2009), and this approach has only recently been generalized and used to predict contacts between proteins in complexes of unknown structure, in an independent effort parallel to this work (Ovchinnikov et al., 2014). In principle, just a small number of key residue–residue contacts across a protein interface would allow computation of 3D models and provide a powerful, orthogonal approach to experiments.
Since the recent demonstration of the use of evolutionary couplings (ECs) between residues to determine the 3D structure of individual proteins (Marks et al., 2011; Morcos et al., 2011; Aurell and Ekeberg, 2012; Jones et al., 2012; Kamisetty et al., 2013), including integral membrane proteins (Hopf et al., 2012; Nugent and Jones, 2012), we reason that an evolutionary statistical approach such as EVcouplings (Marks et al., 2011) could be used to determine co-evolved residues between proteins. To assess this hypothesis, we built an evaluation set based on all known binary protein interactions in E. coli that have 3D structures of the complex as recently summarized (Rajagopala et al., 2014). We develop a score for every predicted inter-protein residue pair based on the overall inter-protein EC score distributions resulting in accurate predictions for the majority of top ranked inter-protein EC pairs (inter-ECs) and sufficient to calculate accurate 3D models of the complexes in the docked subset (Figure 1A). This approach was then used to predict evolutionary couplings for 32 complexes of unknown 3D structures that have a sufficient number of sequences. Predictions include the structurally unsolved interactions between the a-, b-, and c-subunits of ATP synthase, which are supported by previously published experimental results.
We first investigated whether co-evolving residues between proteins are close in three dimensions by assessing blinded predictions of residue co-evolution against experimentally determined 3D complex structures. We follow this evaluation by then predicting co-evolved residue pairs of interacting proteins that have no known complex structure.
To compute co-evolution across proteins, individual protein sequences must be paired up with each other that are presumed to interact, or being tested to see if they interact. Without this condition, proteins could be paired together that do not in fact interact with each other and therefore detection of co-evolution would be compromised. Given that the evolutionary couplings method depends on large numbers of diverse sequences (Hopf et al., 2012), some assumption must be made about which proteins interact with each other in homologous sequences in other species. Since it is challenging to know a priori whether particular interactions are conserved across many millions of years in thousands of different organisms, we use proximity of the two interacting partners on the genome as a proxy for this, with the goal of reducing incorrect pairings.
To assemble the broadest possible data sets to test the approach and make predictions, we take all known interacting proteins assembled in a published data set that contains ∼3500 high-confidence protein interactions in E. coli (Rajagopala et al., 2014). After removing redundancy and requiring close genome distance between the pairs of proteins this results in 326 interactions, see ‘Materials and methods’ (Figure 1B, Figure 1—figure supplement 1, Supplementary file 1 and 2),
The paired sequences are concatenated and statistical co-evolution analysis is performed using EVcouplings (Marks et al., 2011; Morcos et al., 2011; Aurell and Ekeberg, 2012), that applies a pseudolikelihood maximization (PLM) approximation to determine the interaction parameters in the underlying maximum entropy probability model (Balakrishnan et al., 2011; Ekeberg et al., 2013; Kamisetty et al., 2013), simultaneously generating both intra- and inter-EC scores for all pairs of residues within and across the protein pairs (Figure 1A). Evolutionary coupling calculations in previous work have indicated that this global probability model approach requires a minimum number of sequences in the alignment with at least 1 non-redundant sequence per residue (Marks et al., 2011; Morcos et al., 2011; Hopf et al., 2012; Jones et al., 2012; Kamisetty et al., 2013). Our current approach allows complexes with fewer available sequences to be assessed (minimum at 0.3 non-redundant sequences per residue) by using a new quality assessment score to assess the likelihood of the predicted contacts to be correct. The EVcomplex score is based on the knowledge that most pairs of residues are not coupled and true pair couplings are outliers in the high-scoring tail of the distribution (See ‘Materials and methods’, Figure 2A,B, Figure 2—figure supplement 1 and 2). The score can intuitively be understood as the distance from the noisy background of non-significant pair scores, normalized by the number of non-redundant sequences and the length of the protein (‘Materials and methods’, equations 1 and 2). If the number of sequences per residue is not controlled for, there is a large bias in the results, overestimating performance with low numbers of sequences (Figure 2B,C). The precise functional form of the correction for low numbers of sequences was chosen non-blindly after observing the dependencies in the test set.
Of the 329 interactions identified that are close on the E. coli genome, 76 have a sufficient number of alignable homologous sequences and known 3D structures either in E. coli or in other species. This set was used to test the inter-protein evolutionary coupling predictions (Supplementary file 1). The relationship between the EVcomplex score and the precision of the corresponding inter-protein ECs suggests that on average 74% (69%) of the predicted pairs with EVcomplex score greater than 0.8 will be accurate to within 10 Å (8 Å) of an experimental structure of the complex (Figure 2C). Most complexes have at least one inter-protein predicted contact above the selected score threshold of 0.8 (53/76 complexes). Three complexes have more than 20 predicted inter-protein residue contacts which are over 80% accurate, namely the histidine kinase and response regulator system (78 residue pairs), t-RNA synthetase (32 residue pairs), and the vitamin B importer complex (21 residue pairs), with precision over 80% (complex numbers 330, 019, 130 respectively, Figure 2D, Figure 2—figure supplements 3–8, Supplementary file 1).
We suggest that users of EVcomplex consider predicted contacts that lie below the threshold of 0.8 in the context of other biological knowledge, where available, or in comparison to other higher scoring contacts for the same complex. In this way additional true positive inter-residue contacts can be distinguished from false positives. For instance, the ethanolamine ammonia-lyase complex (complex 065) has only 3 predicted inter-protein residue pairs above the score threshold, but in fact has 5 additional correct pairs with EVcomplex scores slightly below the threshold of 0.8 which cluster with the 3 high-scoring contacts on the monomers, indicating that they are also correct.
Some of the high confidence inter-protein ECs in the test set are not close in 3D space when compared to their known 3D structures. These false positives may be a result of assumptions in the method that are not always correct. This includes (1) the assumption that the interaction between paired proteins is conserved across species and across paralogs, and (2) that truly co-evolved residues across proteins are indeed always close in 3D, which is not always the case. In addition, the complexes may also exist in alternative conformations that have not necessarily all been captured yet by crystal or NMR structures, for instance in the case of the large conformational changes of the BtuCDF complex (Hvorup et al., 2007).
To test whether the computed inter-protein ECs are sufficient for obtaining accurate 3D structures of the whole complex, we selected 15 diverse examples (with 5 or more inter-protein residue contacts) for docking (Table 1, Figure 3, Figure 3—figure supplement 1, Supplementary file 3) with HADDOCK (Dominguez et al., 2003; de Vries et al., 2007). The docking procedure is fast and generates 100 3D models of each complex using all residue pairs with EVcomplex scores above the selection threshold. We additionally dock negative controls to assess the amount of information added to the docking protocol by evolutionary couplings (500 models per run, no constraints other than center of mass, see ‘Materials and methods’). The best models for all 15 complexes docked with evolutionary couplings have interface RMSDs under 6 Å, 12/15 have the best scoring model under 4 Å and the top ranked models for 11/15 are under 5 Å backbone interface RMSD compared to a crystal or NMR structure interface. Over 70% of the generated models are close to the experimental structures of the complexes (<4 Å backbone iRMSD), compared to less than 0.5% in the controls (and these were not high–ranked) (Figure 3—figure supplement 1, Supplementary file 3, ‘Supplementary data’). Not surprisingly complexes that have the largest numbers of true positive predicted contacts perform the best when docking. For example, the ribosomal proteins RS3 and RS14 have 11 true positive inter-protein ECs and result in a top ranked model only 1.1 Å iRMSD from the reference structure. More surprisingly, other complexes with a lower proportion of true positive inter-protein contacts, such as Ubiquinol oxidase (6 out of 11) or the epsilon and gamma subunits of ATP synthase (8 out of 15) also produced accurate predicted complexes, with an iRMSD of 1.8 and 1.4 Å respectively. The docking experiments therefore demonstrate that inter-protein ECs, even in the presence of incorrect predictions, can be sufficient to give accurate 3D models of protein complexes, but more work will be needed to quantify the likelihood of successful docking from the predicted contacts.
The top 10 inter-EC pairs between MetI and MetN are accurate to within 8 Å in the MetNI complex (PDB: 3tui [Johnson et al., 2012]), resulting in an average 1.4 Å iRMSD from the crystal structure for all 100 computed 3D models (Table 1, Supplementary file 3 and ‘Supplementary data’). The top 3 inter-EC residue pairs (K136-E108, A128-L105, and E74-R124, MetI-MetN respectively) constitute a residue network coupling the ATP-binding pocket of MetN to the membrane transporter MetI. This network calculated from the sequence alignment corresponds to residues identified experimentally that couple ATP hydrolysis to the open and closed conformations of the MetI dimer (Johnson et al., 2012) (Figure 4A). The vitamin B12 transporter (BtuC) belongs to a different structural class of ABC transporters, but also uses ATP hydrolysis via an interacting ATPase (BtuD). The top 5 inter-ECs co-locate the L-loop of BtuC close to the Q-loop ATP-binding domain of the ATPase, hence coupling the transporter with the ATP hydrolysis state in an analogous way to MetI-MetN. The identification of these coupled residues across the different subunits suggests that EVcomplex identifies not only residues close in space, but also particular pairs that are constrained by the transporter function of these complexes (Kadaba et al., 2008; Johnson et al., 2012).
The ATP synthase ε and γ subunit complex provides a challenge to our approach, since the ε subunit can take different positions relative to the γ subunit, executing the auto-inhibition of the enzyme by dramatic conformational changes (Cingolani and Duncan, 2011). In a real-world scenario, where we might not know this a priori, there may be conflicting constraints in the evolutionary record corresponding to the different positions of the flexible portion of ε subunit. EVcomplex accurately predicts 6 of the top 10 inter-EC pairs (within 8 Å in the crystal structure 1fs0 (Rodgers and Wilce, 2000) or 3oaa (Cingolani and Duncan, 2011)), with the top 2 inter-ECs εA45-γL215 and εA40-γL207 providing contact between the subunits along an inter-protein beta sheet. The location of the C-terminal helices of the ε subunit is significantly different across 3 crystal structures (PDB IDs: 1fs0 [Rodgers and Wilce, 2000], 1aqt [Uhlin et al., 1997], 3oaa [Cingolani and Duncan, 2011]). The top ranked intra-ECs support the conformation seen in 1aqt, with the C-terminal helices packed in an antiparallel manner and tucked against the N-terminal beta barrel (Figure 4B, green circles) and do not contain a high ranked evolutionary trace for the extended helical contact to the γ subunit seen in 1fs0 or 3oaa (Figure 4B, gray box). Docking with the top inter-ECs results in models with 1.4 Å backbone iRMSDs to the crystal structure for the interface between the N-terminal domain of the ε subunit and the γ subunit (Table 1, Supplementary file 3). εD82 and γR222 connect the ε-subunit via a network of 3 high-scoring intra-ECs between the N- and C-terminal helices to the core of the F1 ATP synthase. In summary, these examples suggest that inter-protein evolutionary couplings can provide residue relationships across the proteins that could aid identification of functional coupling pathways, in addition to obtaining 3D models of the complex.
A total of 82 protein complexes with unknown 3D structure of the interaction that satisfy the conditions for the current approach, i.e., have sufficient sequences and are close in all genomes, were predicted using EVcomplex (all residue–residue inter-protein evolutionary couplings scores are available in ‘Supplementary data’). 32 of these have high EVcomplex scores with at least one predicted contact (Figure 5, Figure 5—figure supplement 1 and 2, and Supplementary file 4). Analysis of the inter-EC predictions for known 3D complex structures shows that protein pairs with more high-scoring ECs (EVcomplex score ≥0.8) have a higher proportion of true positives (Figure 2D). Hence, the protein complexes in the set of unknown structures with more high-scoring inter-ECs are most likely to have predicted ECs that indicate residue pairs close in 3D (column Q, Supplementary file 2, the exact pairs can be found in Supplementary file 4). Three examples of predictions with multiple high-scoring inter-ECs include MetQ-MetI, UmuD-UmuC, and DinJ–YafQ. The top 15 inter-ECs between MetQ and MetI are from one interface of MetQ to the MetI periplasmic loops, or the periplasmic end of the helices, consistent with the known binding of MetQ to MetI in the periplasm.
The UmuD and UmuC complex is induced in the stress/SOS response facilitating the cleavage of UmuD to UmuD’ (between C24 and G25) to form UmuD’2 which then interacts with UmuC (DNA polymerase V) in order to copy damaged DNA (Beuning et al., 2006). The truncated dimer form (UmuD’2) has at least two contrasting conformations where the N-terminal arm is placed on opposite sides of the dimer in one conformation or in close proximity in the alternative (Figure 5—figure supplement 3). For 6/7 ECs above the score threshold, residues in UmuD predicted to interface with UmuC are co-located on one face of the dimer. Two residues (Y33, I38) are located in the N-terminal arm of UmuD that, after cleavage of the 24 N-terminal amino acids, may become available for binding UmuC. Since UmuD switches functions after this cleavage and can then bind UmuC, these inter-ECs may identify the critical residues for translesion synthesis function (Beuning et al., 2006). Although the ECs from this UmuD arm to UmuC involve residues in two separate domains of UmuC (S 415 and Y 74), intra-monomer evolutionary couplings predict that these residues are close in UmuC (Figure 5—figure supplement 3A, black rectangles). The relative positions of the contacting residues within each monomer therefore support the plausibility of the accuracy of the interaction interface.
Whilst this manuscript was in review, the 3D structure of the previously unsolved biofilm toxin/antitoxin DinJ–YafQ complex was published (PDB: 3mlo (Liang et al., 2014)), showing the intertwining of subunits in a heterotetrameric complex. 17/19 predicted EC residue pairs are within 8 Å in this 3D structure (Supplementary file 4 and ‘Supplementary data’). In general, the agreement between our de novo predicted inter-protein ECs and available experimental data serves as a measure of confidence for the predicted residue pair interactions and suggests that EVcomplex can be used to reveal 3D structural details of yet unsolved protein complexes given sufficient evolutionary information.
To investigate whether the EVcomplex score can also distinguish between interacting and non-interacting pairs of proteins, we use the E. coli ATP synthase complex as a test case. The ATP synthase structure is of wide biological interest (reviewed in Walker, 2013) with a remarkable 3D structural arrangement, but completion of all aspects of the 3D structure has remained experimentally challenging (Baker et al., 2012) (Figure 6A). As a demonstration exercise, we calculated evolutionary couplings for all 28 possible pair combinations of different ATP synthase subunits (centered around the E. coli ATP synthase) and transformed the ECs into EVcomplex scores for all inter-protein residue pairs (experimentally determined stoichiometry: α3β3γδεab2c10, Supplementary file 5 and ‘Supplementary data’). Using the default EVcomplex score threshold of 0.8 to discriminate between interacting and non-interacting pairs of subunits, 24 of the 28 possible interactions between the subunits are correctly classified as interacting or non-interacting. The four incorrect predictions (namely: ε and c, γ and c, ε and β, b and β, for which there is some experimental evidence) are not identified as interacting using the 0.8 EVcomplex threshold. Choosing a threshold lower than 0.8 does identify 2 of these as interacting but also introduces new false positives. The ε and β interaction in the crystal structure 3oaa (Cingolani and Duncan, 2011) is a special case in that it involves a highly extended conformation of the last two helices of the ε subunit that reach up into the enzyme making contacts with the β subunit. The false negative EVcomplex score for this pair could be a result of the transience of their interaction or reflect a more general problem of lack of conservation of this interaction across the aligned proteins from different species. In total 80% of the interacting residue pairs in the known 3D structure parts of the synthase complex (7 pairs of subunits) are correctly predicted (threshold: 10 Å minimum atom distance between two residues). This exercise of predicting the presence or absence of interaction between any two proteins indicates the potential of the EVcomplex method in helping to elucidate protein–protein interaction networks from evolutionary sequence co-variation and identify interacting subunits of large macromolecular complexes.
While much of the 3D structure of ATP synthase is known (Walker, 2013), the details of interactions between the a-, b-, and c-subunits have not yet been determined by crystallography. We analyse the details in these interactions, as the EVcomplex scores between these subunits are substantial (Figure 6B). We are fortunately able to provide a missing piece for this analysis, the unknown structure of the membrane-integral penta-helical a-subunit, using our previously described method for de novo 3D structure prediction of alpha-helical transmembrane proteins (Hopf et al., 2012). To our knowledge, there are no experimentally determined atomic resolution structures of the a-subunit of ATP synthase. A 3D model of the a-subunit is from 1999 (1c17(Rastogi and Girvin, 1999)) and was computed using five helical–helical interactions that were inferred from second suppressor mutation experiments, and then imposed as distance restraints for TMH2-5, revealing a four helical bundle (with no information for TMH1). Later, cross-linking experiments (Schwem and Fillingame, 2006) identified contacting residues from all pairs of helical combinations of TM2-TM5 (6 pairs), supporting the earlier 4 helical bundle topology. 7 of the 8 cross-linked pairs are either exactly the same pair (L120-I246) or adjacent to many pairs in the top L intra a-subunit evolutionary couplings (ECs).
In fact, the helix packing arrangement in the predicted structure of the a-subunit is consistent with the topology suggested on the basis of crosslinking studies (Long et al., 2002; DeLeon-Rangel et al., 2003; Fillingame and Steed, 2014), including the lack of contacts for transmembrane helix 1 with the other 4 helices (‘Supplementary data’).
The top inter-protein EC pair between subunits a and b, aK74–bE34, coincides with experimental crosslinking evidence of the interaction of aK74 with the b-subunit and the position of E34 of the b subunit emerging from the membrane on the cytoplasmic side (Long et al., 2002; DeLeon-Rangel et al., 2003). Indeed, 6 of the 13 high score ECs are in the same region as the experimental crosslinks, for instance between the cytoplasmic loop between the first two helices of the a-subunit and the b-subunit helix as it emerges from the membrane bilayer (DeLeon-Rangel et al., 2013), or between a239V in TM helix 5 and bL16 (Figure 6C, Figure 6—figure supplement 1, Supplementary file 6). Additionally, the top EC between the a- and c-subunits (aG213 – cM65) lies close to the functionally critical aR210–cD61 interaction (Dmitriev et al., 1999) on the same helical faces of the respective subunits (Figure 6C). This prediction of missing aspects of subunit interactions may help in the design of targeted experiments to complete the understanding of the intricate molecular mechanism of the ATP synthase complex.
A primary limitation of our current approach is its dependence on the availability of a large number of evolutionarily related sequences. If a protein interaction is conserved across enough sequenced genomes, using a single pair per genome can give accurate predictions of the interacting residues. However, if the protein pair is present in limited taxonomic branches, there may be insufficient sequences at any given time to make confident predictions. A solution to this could be to include multiple paralogs of the interacting proteins from each genome, but this requires correct pairing of the interaction partners, which is in general hard to ascertain. In addition, details of interactions may have diverged for paralogous pairs. Hence, in this current version of the method, we have imposed a genome distance requirement across all genomes for all homolog pairs in order to be less sensitive to these complications.
As the need to use genome proximity to pair sequences becomes less important with the increasing availability of genome sequences, there will be a dramatic increase in the number of interactions that can be inferred from evolutionary couplings, including those unique to eukaryotes. With currently available sequences (May 2014 release of the UniProt database), EVcomplex is able to provide information for about 1/10th of the known 3000 protein interactions in the E. coli genome. Once there are ∼10,000 bacterial genome sequences of sufficient diversity, one would have enough information to test each potentially interacting pair of homologs for evidence of interaction and, given sufficiently strong evolutionary couplings, infer the 3D structure of each protein–protein pair, as well as of complexes with more than two proteins. For any set of species, e.g., vertebrates or mammals, one can imagine guiding sequencing efforts to optimize species diversity to facilitate the extraction of evolutionary couplings. This can open the doors for more comprehensive and more rapid determination of approximate 3D structures of proteins and protein complexes, as well as for the elucidation in molecular detail of the most strongly evolutionarily constrained interactions, pointing to functional interactions.
Determining the three-dimensional models of complexes from the predicted contacts was successful in many of the tested instances. Using minimal computing resources and a small number of inter-EC-derived contacts, low interface positional RMSDs relative to experimental structures can be achieved. However, a significant number of proteins exist as homomultimers within larger complexes. To determine models of these complexes one must deconvolute homomultimeric inter-ECs from the intra-protein signal, which is an important technical challenge for future work.
The analysis of subunit interactions in ATP synthase in this work is a ‘proof of principle’ study showing that methods such as EVcomplex can determine which proteins interact with each other at the same time as specific residue pair couplings across the proteins (as also shown in the work by the Baker lab on ribosomal protein interactions (Ovchinnikov et al., 2014)). Understanding the networks of protein interactions is of critical interest in eukaryotic systems, such as networks of protein kinases, GPCRs, or PDZ domain proteins. An understanding of the distributions of interaction specificities is of high interest to many fields. Although we do not know how well our evolutionary coupling approach will handle less obligate interactions, results on the two-component signaling system (histidine kinase/response regulator) both here and in other work (Skerker et al., 2008; Weigt et al., 2009) suggest optimism.
The approximately scale-free EVcomplex score is a heuristic based on the distribution of raw EC scores from the statistical model, their dependence on sequence alignment depth and the length of the concatenated sequences. The score provides a simple way of accounting for these dependencies such that a uniform threshold, say 0.8, can be used for any protein pair with the expectation of reasonably accurate predictions. Since cutoff thresholds can be useful but overly sharp, we recommend investigating predicted contacts below the threshold used in this work, especially where there is independent biological knowledge to validate the predictions.
The work presented here is in anticipation of a genome-wide exploration and, as a proof of principle, shows the accurate prediction of inter-protein contacts in many cases and their utility for the computation of 3D structures across diverse complex interfaces. As with single protein (intra-EC) predictions, evolutionarily conserved conformational flexibility and oligomerization can result in more than one set of contacts that must be de-convoluted. Can evolutionary information help to predict the details and extent for each complex? A key challenge will be the development of algorithms that can disentangle evolutionary signals caused by alternative conformations of single complexes, alternative conformations of homologous complexes, and effectively deal with false positive signals. Taken together, these issues highlight fruitful areas for future development of evolutionary coupling methods.
Despite conditions for the successful de novo calculation of co-evolved residues, the method described here may accelerate the exploration of the protein–protein interaction world and the determination of protein complexes on a genome-wide scale at residue level resolution. The use of co-evolutionary analysis in computational models to determine protein specificity and promiscuity, co-evolutionary dynamics and functional drift will open up exciting future research questions.
The candidate set of complexes for testing and de novo prediction was derived starting from a data set of binary protein–protein interactions in E. coli including yeast two-hybrid experiments, literature-curated interactions and 3D complex structures in the PDB (Rajagopala et al., 2014). Three complexes not contained in the list were added based on our analysis of other subunits in the same complex, namely BtuC/BtuF, MetI/MetQ, and the interaction between ATP synthase subunits a and b. Since our algorithm for concatenating multiple sequence pairs per species assumes the proximity of the interacting proteins on the respective genomes of each species (See below), we excluded any complex with a gene distance >20 from further analysis. The gene distance is calculated as the number of genes between the interacting partners based on an ordered list of genes in the E. coli genome obtained from the UniProt database. The resulting list of pairs (∼350) was then filtered for pseudo-homomultimeric complexes based on the identification of Pfam domains in the interacting proteins (330). All remaining complexes with a known 3D structure (as summarized in Rajagopala et al., 2014) or a homologous interacting 3D structure (93) (identified by intersecting the results of HMMER searches against the PDB for both monomers) were used for evaluating the method, while complexes without known structure (236) were assigned to the de novo prediction set (Figure 1—figure supplement 1). The set with protein complexes of known 3D structure was further filtered for structures that only cover fragments (<30 amino acids) of one or both of the monomers and structures with very low resolution (>5 Å), which led to the re-assignment of Ribonucleoside-diphosphate reductase 1 (complex_002), Type I restriction-modification enzyme EcoKI (complex_012), RpoC/RpoB (complex_041), RL11/Rl7 (complex_165), the ribosome with SecY (complex_226, complex_250, and complex_255), and RS3/RS (complex_254) to the set of unknown complexes. Large proteins were run with the specific interacting domains informed by the known 3D structure, when the full sequence was too large for the number of retrieved sequences (for domain annotation see ‘Supplementary data’).
This set could serve as a benchmark set for future development efforts in the community.
Each protein from all pairs in our data set was used to generate a multiple sequence alignment (MSA) using jackhmmer (Johnson et al., 2010) to search the UniProt database (UniProt Consortium, 2014) with 5 iterations. To obtain alignments of consistent evolutionary depths across all the proteins, a bit score threshold of 0.5 * monomer sequence length was chosen as homolog inclusion criterion (-incdomT parameter), rather than a fixed E-value threshold which selects for different degrees of evolutionary divergence based on the length of the input sequence.
In order to calculate co-evolved residues across different proteins, the interacting pairs of sequences in each species need to be matched. Here, we assume that proteins in close proximity on the genome, e.g., on the same operon, are more likely to interact, as in the methods used previously matching histidine kinase and response regulator interacting pairs (Skerker et al., 2008; Weigt et al., 2009) (‘Supplementary data’). We retrieved the genomic locations of proteins in the alignments and concatenated pairs following 2 rules: (i) the CDS of each concatenated protein pair must be located on the same genomic contig (using ENA (Pakseresht et al., 2014) for mapping) and (ii) each pair must be the closest to one another on the genome, when compared to all other possible pairings in the same species. The concatenated sequence pairs were filtered based on the distribution of genomic distances to exclude outlier pairs with high genomic distances of more than 10k nucleotides (‘Supplementary data’). Alignment members were clustered together and reweighted if 80% or more of their residues were identical (thus implicitly removing duplicate sequences from the alignment). Supplementary file 1 and 2 report the total number of concatenated sequences, the lengths, and the effective number of sequences remaining after down-weighting in the evaluation and de novo prediction set, respectively.
Inter- and intra-ECs were calculated on the alignment of concatenated sequences using a global probability model of sequence co-evolution, adapted from the method for single proteins (Marks et al., 2011; Morcos et al., 2011; Hopf et al., 2012) using a pseudo-likelihood maximization (PLM) (Balakrishnan et al., 2011; Ekeberg et al., 2013) rather than mean field approximation to calculate the coupling parameters. Columns in the alignment that contain more than 80% gaps were excluded and the weight of each sequence was adjusted to represent its cluster size in the alignment thus reducing the influence of identical or near-identical sequences in the calculation. For the evaluation set, we can then compare the predicted ECs for both within and between the proteins/domains to the crystal structures of the complexes (for contact maps and all EC scores, see ‘Supplementary data’).
In order to estimate the accuracy of the EC prediction, we evaluate the calculated inter-ECs based on the following observations: (1) most pairs of positions in an alignment are not coupled, i.e., have an EC score close to zero, and tend to be distant in the 3D structure; (2) the background distribution of EC scores between non-coupled positions is approximately symmetric around a zero mean; and (3) higher-scoring positive score outliers capture 3D proximity more accurately than lower-scoring outliers (See also Figure 2). The width of the (symmetric) background EC score distribution can be approximated using the absolute value of the minimal inter-EC score. The more a positive EC score exceeds the noise level of background coupling, the more likely it is to reflect true co-evolution between the coupled sites. For each inter-protein pair of sites i and j with pair coupling strength ECinter(i, j), we therefore calculate a raw reliability score ('pair coupling score ratio', Figure 2B) defined by
Since the accuracy of evolutionary couplings critically depends both on the number and diversity of sequences in the input alignment and the size of the statistical inference problem (Marks et al., 2011; Morcos et al., 2011; Jones et al., 2012), we incorporate a normalization factor to make the raw reliability score comparable across different protein pairs. The normalized EVcomplex score is defined as
where Neff is the effective number of sequences in the alignment after redundancy reduction, and L (total number of residues) is the length of the concatenated alignment. Previous work on single proteins has shown that the method requires a sufficient number of sequences in the alignment to be statistically meaningful. We thus filter for sequence sufficiency requiring Neff/L > 0.3 (Table 1, Supplementary files 1 and 2). Predictions of coupled residues in the evaluation set were evaluated against their residue distances in known structures of protein pairs (Rajagopala et al., 2014) (See Supplementary file 7) in order to determine the precision of the method.
To interpret the EVcomplex prediction of interaction between subunits a and b of the ATP synthase as well as UmuC and UmuD, individual monomer models were built de novo for the structurally unsolved a-subunit of ATP synthase and UmuC using the EVfold pipeline as previously published (Marks et al., 2011; Hopf et al., 2012). In both cases coupling parameters were calculated using PLM (Balakrishnan et al., 2011; Ekeberg et al., 2013) and sequences were clustered and weighted at 90% sequence identity (the resulting models are provided in ‘Supplementary data’).
Following this same protocol EVcomplex scores were calculated for all possible 28 combinations of the 8 E. coli ATP synthase F0 and F1 subunits. Since we want to compare the computational predictions to some ‘ground truth’, as with the complexes for the rest of the manuscript, we used known 3D structures of the ATP synthase complex to assign whether or not the subunits interact (PDB: 3oaa, 1fs0, 2a7u; Supplementary file 7). Since we are also determining whether the subunits interact, not necessarily knowing full atomic detail residue interactions, we included subunit interactions that have been inferred from cryo-EM, crosslinking or other experiments, but do not necessarily have a crystal structure. These are represented as solid blue boxes, if the interaction is well established (DeLeon-Rangel et al., 2013; Schulenberg et al., 1999; Brandt et al., 2013; McLachlin and Dunn, 2000), or crosshatched blue if there is a lack of consensus in the community (Figure 6B, left panel).
For each possible interaction the EVcomplex score of the highest ranked inter-EC was considered as a proxy for the likelihood of interaction. Pairs with scores above 0.8 are considered likely to interact, between 0.75 and 0.8 weakly predicted, while interactions with scores below 0.75 are rejected as possible complexes (blue boxes, blue crosshatched, and white respectively in right panel of Figure 6B, and ‘Supplementary data’).
A diverse set of 15 complexes was chosen from the 22 in the evaluation set that had at least 5 couplings above a complex score of 0.8 and was subsequently docked (Supplementary file 3). Proteins that have been crystallized together in a complex could bias the results of the docking, as they have complementary positions of the surface side chains. Therefore, where possible we used complexes that had a solved 3D structure of the unbound monomer, namely GcsH/GcsT, CyoA, FimC, DhaL, AtpE, PtqA/PtqB, RS10, and HK/RR, and in all other cases the side chains of the monomers were randomized either by using SCWRL4 (Krivov et al., 2009) or restrained minimization with Schrodinger Protein Preparation Wizard (Sastry et al., 2013) before docking. For ubiquinol oxidase (complex_054) the unbound structure of subunit 2 (CyoA) only covers the COX2 domain. In this case docking was performed using this unbound structure plus an additional run using the bound complex structure with perturbed side chains.
We used HADDOCK (Dominguez et al., 2003), a widely used docking program based on ARIA (Linge et al., 2003) and the CNS software (Brunger, 2007) (Crystallography and NMR System), to dock the monomers for each protein pair with all inter-ECs with an EVcomplex score of 0.8 or above implemented as distance restraints on the α-carbon atoms of the backbone.
Each docking calculation starts with a rigid-body energy minimization, followed by semi-flexible refinement in torsion angle space, and ends with further refinement of the models in explicit solvent (water). 500/100/100 models generated for each of the 3 steps, respectively. All other parameters were left as the default values in the HADDOCK protocol. Each protein complex was run using predicted ECs as unambiguous distance restraints on the Cα atoms (deff 5 Å, upper bound 2 Å, lower bound 2 Å; input files available in ‘Supplementary data’). As a negative control, each protein complex was also docked using center of mass restraints alone (ab initio docking mode of HADDOCK) (de Vries et al., 2007) and in this case generating 10,000/500/500 models.
Each of the generated models is scored using a weighted sum of electrostatic (Eelec) and van der Waals (Evdw) energies complemented by an empirical desolvation energy term (Edesolv) (Fernandez-Recio et al., 2004). The distance restraint energy term was explicitly removed from the equation in the last iteration (Edist3 = 0.0) to enable comparison of the scores between the runs that used a different number of ECs as distance restraints.
All computed models in the docked set were compared to the cognate crystal structures by the RMSD of all backbone atoms at the interface of the complex using ProFit v.3.1 (http://www.bioinf.org.uk/software/profit/). The interface is defined as the set of all residues that contain any atom <6 Å away from any atom of the complex partner. For the AtpE–AtpG complex, we excluded the 2 C-terminal helices of AtpE as these helices are mobile and take many different positions relative to other ATP synthase subunits (Cingolani and Duncan, 2011). Similarly, since the DHp domain of histidine kinases can take different positions relative to the CA domain, the HK-RR complex was compared over the interface between the DHp domain alone and the response regulator partner. In the case of the unbound ubiquinol oxidase docking results, only the interface between COX2 in subunit 2 and subunit 1 was considered. Accuracy of the computed models with EC restraints was compared with computed models with center of mass restraints alone (negative controls) (Figure 3—figure supplement 1, Supplementary file 3).
Supplementary data 1: Concatenated alignments for complexes predicted in this work
Supplementary data 2: Genome distance distribution of concatenated sequences per alignment
Supplementary data 3: EVcomplex predictions for evaluation and de novo set
Supplementary data 4: Docking input files and top 10 predicted models for evaluation set
Supplementary data 5: ATP synthase predictions, ATP synthase subunit a model
Supplementary data 6: UmuC model
Evolution of protein interactions: from interactomes to interfacesArchives of Biochemistry and Biophysics 554:65–75.https://doi.org/10.1016/j.abb.2014.05.010
Arrangement of subunits in intact mammalian mitochondrial ATP synthase determined by cryo-EMProceedings of the National Academy of Sciences of USA 109:11675–11680.https://doi.org/10.1073/pnas.1204935109
Learning generative models for protein fold familiesProteins 79:1061–1078.https://doi.org/10.1002/prot.22934
Individual interactions of the b subunits within the stator of the Escherichia coli ATP synthaseThe Journal of Biological Chemistry 288:24465–24479.https://doi.org/10.1074/jbc.M113.465633
Structure of the ATP synthase catalytic complex (F(1)) from Escherichia coli in an autoinhibited conformationNature Structural & Molecular Biology 18:701–707.https://doi.org/10.1038/nsmb.2058
The role of transmembrane span 2 in the structure and function of subunit a of the ATP synthase from Escherichia coliArchives of Biochemistry and Biophysics 418:55–62.https://doi.org/10.1016/S0003-9861(03)00391-6
Structure of the subunit c oligomer in the F1Fo ATP synthase: model derived from solution structure of the monomer and cross-linking in the native enzymeProceedings of the National Academy of Sciences of USA 96:7785–7790.https://doi.org/10.1073/pnas.96.14.7785
HADDOCK: a protein-protein docking approach based on biochemical or biophysical informationJournal of the American Chemical Society 125:1731–1737.https://doi.org/10.1021/ja026939x
Improved contact prediction in proteins: using pseudolikelihoods to infer Potts modelsPhysical Review. E, Statistical, Nonlinear, and Soft Matter Physics 87:012707.
InterEvol database: exploring the structure and evolution of protein complex interfacesNucleic Acids Research 40:D847–D856.https://doi.org/10.1093/nar/gkr845
Identification of protein-protein interaction sites from docking energy landscapesJournal of Molecular Biology 335:843–865.https://doi.org/10.1016/j.jmb.2003.10.069
Half channels mediating H transport and the mechanism of gating in the F sector of Escherichia coli FF ATP synthaseBiochimica Et Biophysica Acta, 10.1016/j.bbabio.2014.03.005.
Data from: Sequence co-evolution gives 3D contacts and structures of protein complexesDryad, 10.5061/dryad.6t7b8.
Assessing the utility of coevolution-based residue-residue contact predictions in a sequence- and structure-rich eraProceedings of the National Academy of Sciences of USA 110:15674–15679.https://doi.org/10.1073/pnas.1314045110
A simple physical model for binding energy hot spots in protein-protein complexesProceedings of the National Academy of Sciences of USA 99:14116–14121.https://doi.org/10.1073/pnas.202485799
Computational design of protein-protein interactionsCurrent Opinion in Chemical Biology 8:91–97.https://doi.org/10.1016/j.cbpa.2003.12.008
Computational redesign of protein-protein interaction specificityNature Structural & Molecular Biology 11:371–379.https://doi.org/10.1038/nsmb749
Structural and functional Characterization of Escherichia coli toxin-antitoxin complex DinJ-YafQThe Journal of Biological Chemistry 289:21191–21202.https://doi.org/10.1074/jbc.M114.559773
Characterization of the first cytoplasmic loop of subunit a of the Escherichia coli ATP synthase by surface labeling, cross-linking, and mutagenesisThe Journal of Biological Chemistry 277:27288–27293.https://doi.org/10.1074/jbc.M202118200
Protein structure prediction from sequence variationNature Biotechnology 30:1072–1080.https://doi.org/10.1038/nbt.2419
Direct-coupling analysis of residue coevolution captures native contacts across many protein familiesProceedings of the National Academy of Sciences of USA 108:E1293–E1301.https://doi.org/10.1073/pnas.1111471108
Interactome3D: adding structural details to protein networksNature Methods 10:47–53.https://doi.org/10.1038/nmeth.2289
Accurate de novo structure prediction of large transmembrane protein domains using fragment-assembly and correlated mutation analysisProceedings of the National Academy of Sciences of USA 109:E1540–E1547.https://doi.org/10.1073/pnas.1120036109
Assembly information services in the European Nucleotide ArchiveNucleic Acids Research 42:D38–D43.https://doi.org/10.1093/nar/gkt1082
Similarity of phylogenetic trees as indicator of protein-protein interactionProtein Engineering 14:609–614.https://doi.org/10.1093/protein/14.9.609
Correlated mutations contain information about protein-protein interactionJournal of Molecular Biology 271:511–523.https://doi.org/10.1006/jmbi.1997.1198
IPython: a system for Interactive Scientific computingComputing in Science and Engineering 9:21–29.https://doi.org/10.1109/MCSE.2007.53
The binary protein-protein interaction landscape of Escherichia coliNature Biotechnology 32:285–290.https://doi.org/10.1038/nbt.2831
Structure of the gamma-epsilon complex of ATP synthaseNat Struct Biol 7:1051–1054.https://doi.org/10.1038/80975
Protein and ligand preparation: parameters, protocols, and influence on virtual screening enrichmentsJournal of Computer-aided Molecular Design 27:221–234.https://doi.org/10.1007/s10822-013-9644-8
The gammaepsilon-c subunit interface in the ATP synthase of Escherichia coli. cross-linking of the epsilon subunit to the c subunit ring does not impair enzyme function, that of gamma to c subunits leads to uncouplingThe Journal of Biological Chemistry 274:34233–34237.https://doi.org/10.1074/jbc.274.48.34233
Cross-linking between helices within subunit a of Escherichia coli ATP synthase defines the transmembrane packing of a four-helix bundleThe Journal of Biological Chemistry 281:37861–37867.https://doi.org/10.1074/jbc.M607453200
Activities at the Universal protein resource (UniProt)Nucleic Acids Research 42:D191–D198.https://doi.org/10.1093/nar/gkt1140
Assembly of macromolecular complexes by satisfaction of spatial restraints from electron microscopy imagesProceedings of the National Academy of Sciences of USA 109:18821–18826.https://doi.org/10.1073/pnas.1216549109
The ATP synthase: the understood, the uncertain and the unknownBiochemical Society Transactions 41:1–16.https://doi.org/10.1042/BST20110773
Modeling of proteins and their assemblies with the Integrative Modeling PlatformMethods in molecular biology 1091:277–295.https://doi.org/10.1007/978-1-62703-691-7_20
Identification of direct residue contacts in protein-protein interaction by message passingProceedings of the National Academy of Sciences of USA 106:67–72.https://doi.org/10.1073/pnas.0805923106
John KuriyanReviewing Editor; Howard Hughes Medical Institute, University of California, Berkeley, United States
eLife posts the editorial decision letter and author response on a selection of the published articles (subject to the approval of the authors). An edited version of the letter sent to the authors after peer review is shown, indicating the substantive concerns or comments; minor concerns are not usually shown. Reviewers have the opportunity to discuss the decision before the letter is sent (see review process). Similarly, the author response typically shows only responses to the major concerns raised by the reviewers.
Thank you for sending your work entitled “Sequence co-evolution gives 3D contacts and structures of protein complexes”““ for consideration at eLife. Your article has been favorably evaluated by John Kuriyan (Senior editor), working with a member of our Board of Reviewing Editors, and 3 reviewers.
The editors and the other reviewers discussed their comments before we reached this decision, and the Reviewing editor has assembled the following comments to help you prepare a revised submission.
This manuscript describes an approach to predict 3D residue-residue contacts at protein interfaces from an analysis of multiple sequence alignments. The described computational approach (using sequence co-evolution to predict contacts between interacting proteins and to generate three-dimensional models of complex structures) appears to be quite successful and is an important contribution that is of interest to a broad audience. The analyses, results, and conclusions are similar to work conducted during the same period by Baker & coworkers, published recently in eLife. Nevertheless, it is the consensus opinion of the reviewers that the present manuscript is potentially suitable for publication in eLife as well, provided that all of the comments raised by the reviewers can be addressed satisfactorily.
We have appended the comments from the three reviewers below. Although many individual points are raised, you will see that the principal concern of the reviewers is that the manuscript falls short of providing the reader with sufficient analysis to judge the validity of the conclusions. It should be possible for you to revise the manuscript to address these concerns without much in the way of new calculations, so we hope that it should be relatively straightforward to deal with these issues.
In addition, since the Baker paper came out recently, it should be given appropriate credit in a revised version. We recommend removing the word “new” from the manuscript for this reason, and instead simply state what was done.
1) The Abstract claims that the authors' approach can discriminate between interacting and non-interacting proteins. However, as far as I can see, the rest of the paper concerns only calculating the interface residues of proteins which are known to interact. No evidence is presented to support the claim to be able to calculate non-interactors. The authors' claim should be corrected accordingly.
2) The authors claim that the co-evolutionary approach has not previously been applied to prediction of protein-protein interfaces. This is incorrect, as the authors seem to be unaware of the previous work of Raphael Guerois' group, which seems rather surprising. The authors should acknowledge the prior work of Faure et al (2012) and Andreani et al. (2013), and they should modify their claim to novelty accordingly. It would also be appropriate to cite properly the prior work of the Baker group (Ovchinnikov et al., 2014) in the same paragraph (the authors currently mention this work only as a “note added during submission”), and the recent review by Andreani and Guerois (2014).
3) In several places the authors say that their approach exploits the fact that interacting proteins are often coded close to each other in a genome. Why is this assumption necessary? According to the description given in Methods, any pair of sequences which are presumed to interact could be concatenated and then used in the authors method. Please clarify.
4) If I understand the authors’ method correctly, they first concatenate the multiple sequences of the presumed interactors, and they then use their previous approach for detecting interacting pairs of columns in the multiple alignment. In this case, when concatenating the alignments for protein “A” with those of protein “B”, the “intra-EC” pairs would appear as A-A and B-B pairs, while the “inter-EC” pairs would come from A-B or B-A interactions. My question is then, are there any observable differences in the inter and intra scores? One might expect that the conservation scores would be lower for inter interactions than intra interactions, as an interface might show more “plasticity” than the core of a domain. It would be very interesting if the authors could comment on this, preferably supported by numerical results.
5) It is incorrect for the authors to refer to their own data set as a “benchmark”. Suggestion: remove all references to the term “benchmark”. However, if the authors make all their data available in a convenient way on-line, it could be proposed as a new “benchmark”.
6) It is wrong to claim that “Experimental evidence... agrees with EVcomplex predictions” (!). Please rewrite this to say that your predictions agree with the experimental evidence. But, if this is your claim, please also support it in some way. What is the experimental evidence? How to you calculate and quantify “agrees with”?
7) Same point in the main text. “The resulting model is consistent with the topologies from cross-linking”. How do you quantify consistent? RMSD from the earlier models? Please provide supporting details.
My main comments concern the way the data in the current manuscript are analyzed and presented. While the main Figures / Table 1 currently illustrate useful examples, they should also show overall analyses of the results over all cases in the benchmark datasets. In several cases this is essential to support the stated conclusions. It does not seem that this would require a substantial amount of new work or lengthy simulations, as most of the relevant data are in the Supplementary Materials (most data are in raw table format or with a separate Figure for each example, which makes it hard for the reader to gauge overall performance – this could be presented in a more easily digestible format in the main text). Specifically, I feel the paper could be substantially improved if the following points were addressed:
1) Please support the statement “Benchmark calculations here ... indicate that the number of sequences in the alignment is critical...” with a Figure showing overall performance versus # of sequences (Figure 3 shows dependency of accuracy on relative rank for only two examples).
2) Show a main Figure comparing predicted contacts against residue distances in the 30 known crystal structures of complexes that had a sufficient number of sequences (Figure 2 only shows contact plots for two examples, and structural pictures for several more, but no overall quantification). This type of analysis could also help to support the statement “For the top ranked benchmark complexes, the majority of the top 5 ECs is correct to within 8A...” with a more quantitative analysis (Table 1 only shows 7 of the 30 complexes in the dataset; Supplementary file 3 shows all individual contact maps, but no overall quantification).
3) “... these benchmarks show ... and demonstrate the criteria needed for successful prediction of unknown interactions”. To support this statement, it would be useful to have a main Figure / Table for each important criterion, analyzed over an entire benchmark set (as for the number of sequences in point 1).
4) Please add metrics to support statements that ECs are “completely consistent” with crosslinking data, or that a coupling “coincides with experimental evidence”. The comparison Table in Supplementary Figure 6 could be presented in more quantitative terms in the main text (for example, what does “crosslinking neighborhood” mean?).
5) What are the possible reasons for false-positive co-variation that does not correspond to contact?
6) The manuscript deals entirely with pairwise co-evolution, and seems to neglect the possibility of higher-order complexity in the sequence pattern.
The manuscript presents an upbeat and optimistic view of the success of the method, with a relatively small amount of critical discussion and little if any investigation of what one can learn from cases in which it does not perform as expected. (For example, a list of possible reasons for false-positive co-variation that does not correspond to contact is given, but an analysis of actual instances is not.) Likewise, a set of methods is presented that is relatively particular, and the reason why these specific methods were chosen over other possibilities is not provided. Moreover, the manuscript deals entirely with pairwise co-evolution, and seems to neglect the possibility of higher-order complexity in the sequence pattern. The manuscript could be substantially strengthened by adding analyses and insights in these issues.https://doi.org/10.7554/eLife.03430.031
- Thomas A Hopf
- Debora S Marks
- Chris Sander
- Charlotta P I Schärfe
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
- John Kuriyan, Reviewing Editor, Howard Hughes Medical Institute, University of California, Berkeley, United States
© 2014, Hopf et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.