Author response: Sequence co-evolution gives 3D contacts and structures of protein complexes (2014)

Full text Figures and data Side by side Abstract eLife digest Introduction Results Discussion Materials and methods Data availability References Decision letter Author response Article and author information Metrics Abstract 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 eLife digest 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 Introduction 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. Figure 1 with 1 supplement see all Download asset Open asset Co-evolution of residues across protein complexes from the evolutionary sequence record. (A) Evolutionary pressure to maintain protein–protein interactions leads to the co-evolution of residues between interacting proteins in a complex. By analyzing patterns of amino acid co-variation in an alignment of putatively interacting homologous proteins (left), evolutionary couplings between co-evolving inter-protein residue pairs can be identified (middle). By defining distance restraints on these pairs, the 3D structure of the protein complex can be inferred using docking software (right). (B) Distribution of E. coli protein complexes of known and unknown 3D structure where both subunits are close on the bacterial genome (left), allowing sequence pair matching by genomic distance. For a subset of these complexes, sufficient sequence information is available for evolutionary couplings analysis (dark blue bars). As more genomic information is created through on-going sequencing efforts, larger fractions of the E. coli interactome become accessible for EVcomplex (right). A detailed version of the workflow used to calculate all E. coli complexes currently for which there is currently enough sequence information is shown in Figure1—figure supplement 1. https://doi.org/10.7554/eLife.03430.003 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. Extension of the evolutionary couplings method to protein complexes 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. Figure 2 with 8 supplements see all Download asset Open asset Evolutionary couplings capture interacting residues in protein complexes. (A) Inter- and Intra-EC pairs with high coupling scores largely correspond to proximal pairs in 3D, but only if they lie above the background level of the coupling score distribution. To estimate this background noise a symmetric range around 0 is considered with the width being defined by the minimum inter-EC score. For the protein complexes in the evaluation set, this distribution is compared to the distance in the known 3D structure of the complex that is shown here for the methionine transporter complex, MetNI. (Plots for all complexes in the evaluation set are shown in Figure 2—figure supplement 1 and 2.) (B) A larger distance from the background noise (ratio of EC score over background noise line) gives more accurate contacts. Additionally, the higher the number of sequences in the alignment the more reliable the inferred coupling pairs are which then reduces the required distance from noise (different shades of blue). Residue pairs with an 8 Å minimum atom distance between the residues are defined as true positive contacts, and precision = TP/(TP + FP). The plot is limited to range (0,3) which excludes the histidine kinase—response regulator complex (HK–RR)—a single outlier with extremely high number of sequences. (C) To allow the comparison across protein complexes and to estimate the average inter-EC precision for a given score threshold independent of sequence numbers, the raw couplings score is normalized for the number of sequences in the alignment, resulting in the EVcomplex score. In this work, inter-ECs with an EVcomplex score ≥0.8 are used. Note: the shown plot is cut off at a score of 2 in order to zoom in on the phase change region and the high sequence coverage outlier HK-RR is excluded. (D) For complexes in the benchmark set, inter-EC pairs with EVcomplex score ≥0.8 give predictions of interacting residue pairs between the complex subunits to varying accuracy (8 Å TP distance cutoff). All predicted interacting residues for complexes in the benchmark set that had at least one inter-EC above 0.8 are shown as contact maps in Figure 2—figure supplement 3–8. https://doi.org/10.7554/eLife.03430.005 Blinded prediction of known complexes Evolutionary covariation reveals inter-protein contacts 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). Docking is accurate with few pairs of predicted contacts 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. Table 1 EVcomplex predictions and docking results for 15 protein complexes https://doi.org/10.7554/eLife.03430.014 EVcomplex contactsDocking quality (iRMSD)Complex nameSubunitsSeqs†ECs‡TP rate§Top ranked model#Best model¶Carbamoyl-phosphate synthaseCarB:CarA2.3170.881.91.9Aminomethyltransferase/Glycine cleavage system H proteinGcsH:GcsT2.950.25.45.4Histidine kinase/response regulatorKdpD:CheY (T. maritima)95.4780.722.12.0Ubiquinol oxidaseCyoB:CyoA1.0110.551.81.2Outer membrane usher protein/Chaperone proteinFimD:FimC3.660.833.23.0Molybdopterin synthaseMoaD:MoaE3.681.04.44.1Methionine transporter complexMetN:MetI1.9140.861.51.2Dihydroxyacetone kinaseDhaL:DhaK1.4120.426.72.4Vitamin B12 uptake B12 uptake synthase and complex of the of non-redundant sequences in concatenated alignment normalized by alignment with EVcomplex score for inter-ECs above score threshold. iRMSD of model from known structure, for docked model with best HADDOCK score. iRMSD across all Figure 3 with 1 supplement see all Download asset Open asset Blinded prediction of evolutionary couplings between complex subunits with known 3D structure. with EVcomplex score ≥0.8 on a selection of benchmark complexes subunits in and inter-ECs in pairs than 8 Å by The predicted inter-ECs for these complexes were then used to 3D models of the complex using protein–protein docking. For the complexes for which 3D structures were predicted using are shown in Figure 3—figure supplement 1. residue networks provide of functional constraints The top 10 inter-EC pairs between and are accurate to within 8 Å in the complex et al., 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 and a residue coupling the of to the membrane transporter This from the sequence alignment to residues identified experimentally that ATP to the and conformations of the et al., (Figure The vitamin B12 transporter to a different structural of but also ATP an interacting The top 5 inter-ECs the of close to the of the coupling the transporter with the ATP in an way to The of these coupled residues across the different subunits suggests that EVcomplex identifies not only residues close in but also particular pairs that are by the transporter of these complexes et al., 2008; et al., 2012). Figure 4 Download asset Open asset Evolutionary couplings give accurate 3D structures of complexes. EVcomplex predictions and comparison to crystal structure for (A) the transporter from E. coli and (B) the interaction of E. coli ATP synthase complex contact predicted inter-ECs with EVcomplex score ≥0.8 and to the chosen inter-EC and blue top and lower to close pairs in the complex crystal for minimum atom distance of Å for contacts and for Å within the with an EVcomplex score ≥0.8 are also on the subunits of the complex on and blue couplings than 8 Å in lower of the top ranked model from 3D docking the complex crystal structure and of the interface region with coupled residues The ATP synthase and complex a to our the can take different to the the of the by conformational changes and In a where we might not know this a there may be constraints in the evolutionary corresponding to the different of the of EVcomplex 6 of the top 10 inter-EC pairs 8 Å in the crystal structure and or and with the top 2 inter-ECs and contact between the subunits an inter-protein The of the of the is significantly different across 3 crystal structures and et al., and The top ranked the in with the in an and against the (Figure and do not a high ranked evolutionary for the extended contact to the in or (Figure Docking with the top inter-ECs results in models with 1.4 Å backbone to the crystal structure for the interface between the of the and the (Table 1, Supplementary file and the a of 3 high-scoring between the and to the of the ATP In these examples suggest that inter-protein evolutionary couplings can provide residue across the proteins that could of functional coupling in to obtaining 3D models of the complex. prediction of unknown complexes of interactions for 32 protein pairs with high-scoring evolutionary couplings A of 82 protein complexes with unknown 3D structure of the interaction that the for the current have sufficient sequences and are close in all were predicted using EVcomplex 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 Figure supplement 1 and and Supplementary file of the inter-EC predictions for known 3D complex structures that protein pairs with more high-scoring ECs score have a higher proportion of true positives (Figure the protein complexes in the set of unknown structures with more high-scoring inter-ECs are most likely to have predicted ECs that residue pairs close in 3D Supplementary file the pairs can be found in Supplementary file Three examples of predictions with high-scoring inter-ECs include and The top 15 inter-ECs between and are from one interface of to the or the of the with the known binding of to in the Figure 5 with 3 supplements see