Large-scale determination of previously unsolved protein structures using evolutionary information

Abstract

The prediction of the structures of proteins without detectable sequence similarity to any protein of known structure remains an outstanding scientific challenge. Here we report significant progress in this area. We first describe de novo blind structure predictions of unprecendented accuracy we made for two proteins in large families in the recent CASP11 blind test of protein structure prediction methods by incorporating residue-residue co-evolution information in the Rosetta structure prediction program. We then describe the use of this method to generate structure models for 58 of the 121 large protein families in prokaryotes for which three-dimensional structures are not available. These models, which are posted online for public access, provide structural information for the over 400,000 proteins belonging to the 58 families and suggest hypotheses about mechanism for the subset for which the function is known, and hypotheses about function for the remainder.

Keywords: B. subtilis; E. coli; H. salinarum; S. solfataricus; biophysics; co-evolution; evolutionary biology; genomics; protein family; protein fold; structural biology; structure prediction.

Figure 1.. Accurate blind structure prediction of CASP11 targets T0806 and T0824.

(A) Location of the most strongly co-evolving residue pairs. Lines connect residue pairs with normalized coupling strength greater than 1.0; yellow, distance less than 5 Å; orange, distance less than 10 Å and red, greater than 10 Å in the models. (B) CASP11 submitted models, colored from N to C terminus (blue to red). (C) X-ray crystal structures. For T0806, the Cα RMSD over the full-length protein is 3.6 Å and 2.9 Å over 223 aligned residues. For T0824; the Cα RMSD over the full-length protein is 4.2 Å and 2.7 Å over 77 aligned residues. For statistics on all five models submitted during CASP, see Figure 1—source data 1. DOI: http://dx.doi.org/10.7554/eLife.09248.003

Figure 2.. Conserved residues tend to cluster in the predicted structures.

Residue conservations from multiple sequence alignments were mapped to predicted structures using Al2Co (Pei and Grishin, 2001) and are colored in rainbow from blue (variable) to red (conserved). The most conserved residues (red or orange), displayed as spheres to highlight their positions, tend to line interaction surfaces and indicate potential functional sites. DOI: http://dx.doi.org/10.7554/eLife.09248.008

Figure 3.. Predicted structure of the Cytochrome bd oxidase complex.

(A) Location of the top co-evolving residue pairs in our model. For clarity, the monomers have been pulled apart slightly. (B) Location of conserved and experimentally characterized residues (Borisov et al., 2011) on structure model. (C) Residues that coordinate heme in CydA are in the same location as histidines (yellow spheres) in Cytochrome b561 (PDB: 4O6Y). H126 (red sphere) overlaps one of these histidines and is the proposed as a heme d coordination site. For clarity, both the model of CydA and the structure of Cyt b562 (4O6Y) are trimmed to highlight the four helix bundle(s). DOI: http://dx.doi.org/10.7554/eLife.09248.009

Figure 4.. Predicted structure of the tartrate dehydratase heterotetramer composed of two copies each of ttdA and ttdB.

(A) Co-evolving residue pairs. The monomers have been pulled apart to reveal the contacts. (B) The ttdA subunit (rainbow) contains a 4Fe-4S cluster (white spheres) that is near the interface with ttdB (green). DOI: http://dx.doi.org/10.7554/eLife.09248.010

Figure 5.. Succinate-acetate/proton symporter SatP (YaaH).

(A) Co-evolving residue pairs in homo-oligomer model. (B) Co-evolving residue pairs in SatP monomer model. (C) SatP co-evolution-based model places known acetate selective residues (magenta spheres) lining the channel. Conserved residues (gray spheres) line the periplasmic surface. The 6TMH channels are formed by threefold pseudo-symmetric TMH hairpins. (D) Proton-gated UreI channel protomer. C-alpha positions at the periplasmic constriction site (magenta spheres) and the cytoplasmic constriction site (pink spheres) are highlighted. The SatP model has the same fold as UreI (C vs D). DOI: http://dx.doi.org/10.7554/eLife.09248.011

Figure 6.. Lipid II flippase (FtsW) in complex with the transmembrane domain of Peptidoglycan synthase (FtsI).

FtsW is an essential cell division protein that transports lipids across the cytoplasmic membrane and is required for localization of FtsI. (A) Location of the top co-evolving residue pairs. (B) White spheres indicate conserved positions in FtsW that when mutated to alanine result in loss of flippase activity. (C, D) The last seven transmembrane (TM) helices of FtsW (TMH4-TMH10) adopting a similar topology as TMH4-TMH7 and TMH10-12 of the TM domain of STT3 (PDB: 3WAK). Both the model of FtsW and 3WAK was trimmed over the aligned helices for clarity. (C) Two residues from the corresponding TMH4 of FtsW (R145 and K153) are essential for flippase activity. (D) The side chain of R145 overlaps the residues that coordinate the divalent metal in the conserved DxH motif at the N-terminus of STT3 TMH4. (E) The model of RfaL adopts a similar fold as ftsW. DOI: http://dx.doi.org/10.7554/eLife.09248.012

Figure 7.. Prolipoprotein diacylglyceryl transferase (LGT).

(A) Predicted contacts indicated on model, (B) model with conserved positions at which alanine mutations result in loss in activity indicated in white spheres; five of these are clustered at the periplasmic end of the model. DOI: http://dx.doi.org/10.7554/eLife.09248.013

Figure 8.. UppP catalyzes the dephosphorylation of undecaprenyl diphophate (UPP).

(A) Location of the top co-evolving residue pairs. (B) Spheres in white indicate conserved residues experimentally shown to decrease activity to <1% (Chang et al., 2014); all these residues are in the core in the model. YfcA, a protein of unknown function is a very distant sequence homologue of UppP (they are in different PFAM families); (C) the predicted structure of YfcA has the same fold as UppP with prominent broken helices (highlighted in blue and yellow). DOI: http://dx.doi.org/10.7554/eLife.09248.014

Figure 9.. PrsW is an intramembrane protease that is crucial in the resistance to antimicrobial peptides.

(A) Location of the top co-evolving residue pairs. Our model of PrsW (B) contains a 4TMH substructure ([TMHs 3–6], (C) which is very similar to a substructure of type II CAAX prenyl protease [TMHs 4–7; D]). The predicted active site residues in PrsW motifs EExxK (TMH3; blue spheres, C) and HxxxD (TMH6; red spheres, C) are in positions similar to those of conserved residues in motifs EExxxR (TMH4; blue spheres, D) and HxxxN (TMH7; red spheres, D) of type II CAAX prenyl protease. Another conserved histidine in the fifth TMH of PrsW (but absent in type II CAAX prenyl protease) is also located in the predicted active site of PrsW (black sphere). DOI: http://dx.doi.org/10.7554/eLife.09248.015

Figure 10.. Our model of the inner membrane protein YeiH (A, B) is structurally similar to the structure of the antiporter NapA (C).

Lower panels: TM helices of core domains are highlighted in white and magenta: while these helices cross over each other in NapA (right), the core of our model of YeiH (left) is formed by two pairs of broken helices (TM5, 6 and TM8, 9) that exit the membrane on the same side. DOI: http://dx.doi.org/10.7554/eLife.09248.016

Figure 11.. Our model of RarD has a similar architecture to EmrE but different fold.

(A, B) Full-length RarD and EmrE homodimer. (C) RarD internal repeat and EmrE monomer. The N-terminus helix (blue) is swapped in RarD relative to EmrE due to helix insertion (gray). DOI: http://dx.doi.org/10.7554/eLife.09248.017

Figure 12.. Predicted structures of YqfA and YhhN have topologies similar to G protein-coupled receptors (GPCR-like).

The core seven TM helix (TMH) fold exhibited by members of the GPCR superfamily is colored in rainbow from the N- to the C-terminus. (A) Bacteriorhodopsin binds retinal (pink spheres) in a pocket formed by TMH3-7 [PDB ID: 1m0k]. (B) The agonist (pink spheres) binding site of P2Y12 receptor is formed by the same set of helices [PDB ID: 4pxz]. (C) A co-evolution-based structure model for YhhN has the GPCR topology with an N-terminal TMH extension. Conserved residues that might form an active site (magenta spheres) cluster in a similar place as the YqfA catalytic residues. (D) Our co-evolution-based structure model for YqfA has a GPCR like topology and clusters residues that may form an active site (magenta spheres mark the Calpha position) in a region that corresponds to the GPCR ligand-binding pocket. (E) Side and top view of the TMalign superposition of YqfA model (in rainbow) over the recently released 3wxw (in white) human ortholog. The N- and C-terminal loops were trimmed for clarity. The TMalign score between the model and the homolog is 0.8. DOI: http://dx.doi.org/10.7554/eLife.09248.018

Figure 13.. YfiP predicted structure has methyltransferase-like fold with knot.

Left: the top co-evolving residues pairs. Middle: conserved residues (magenta) surround the AdoMet-binding site and a conserved Cys could bind a Fe4S cluster. Right: 3nk7 methyltransferase bound to AdoMet. DOI: http://dx.doi.org/10.7554/eLife.09248.019

Figure 14.. Bacillus subtilis YitE model.

(A) The top co-evolving residues. YitE (B) has architecture similar to aquaporin (PDB:2B6P) (C), including the internal pseudo-symmetry (blue vs green), but completely different connectivity. DOI: http://dx.doi.org/10.7554/eLife.09248.020

Figure 15.. Escherichia coli protein YgdD.

Our model of the E. coli protein YgdD trimer (C) is based on predicted contacts satisfied within the monomer (B) and between monomers in the homo-trimer (A). Structural similarity to heme copper oxidase (D) along with a weak HHpred sequence match over part of the protein suggests that YgdD is evolutionarily related to heme–copper oxidases. DOI: http://dx.doi.org/10.7554/eLife.09248.021

Figure 16.. Dependence of the accuracy of predicted contacts on the normalized GREMLIN score (sco), the effective number of sequences (seq), the length (len), and the sequence separation (sep).

Contacts are defined based on amino acid specific Cβ-Cβ distance cutoffs as described in SI Table 3 in Kamisetty et al. (2013). (A) Observed vs predicted accuracies over a large data set of proteins of known structure with deep alignments (Supplementary file 3), sub sampled to different extents (seq/√(len) = 4 (red), 8 (green), 15 (purple), 32 (cyan), and 96 (orange)). Circles represent observed contact prediction accuracies, solid lines, a fit to a sigmoid function of the normalized coupling value, the number of sequences, the length, and the sequence separation (see Figure 16—figure supplement 1 and Figure 16—figure supplement 2). (B) Observed vs predicted accuracies in an independent data set of variable length alignments for 7047 pdb chains (Supplementary file 3), using maximum number of sequences obtained with HHblits as opposed to subsampling a large alignment. Circles again represent observed contact prediction accuracies; solid lines, the predicted accuracy using the model obtained by fitting to the data in (A). The contact prediction accuracy is correctly modeled for the independent data set, justifying its use on the unknown cases described in this article. The Equation use to calculate P(contact|sco,seq,len,sep) is $\text{P}\left(\text{contact}|\text{sco},\text{seq},\text{len},\text{sep}\right)\approx \frac{0.89\left(1–\text{P}\left(\text{contact}|\text{sep}\right)\right)}{1+\text{exp}\left(-0.58{\left(\frac{\text{seq}}{\sqrt{\text{len}}}\right)}^{0.50}\left(\text{sco}-5.46{\left(\frac{\text{seq}}{\sqrt{\text{len}}}\right)}^{-0.53}\right)\right)}+\text{P}\left(\text{contact}|\text{sep}\right).$ DOI: http://dx.doi.org/10.7554/eLife.09248.022

Figure 16—figure supplement 1.. Contact prediction accuracy is better correlated with (#sequences/sqrt(length)) than with (#sequences/length).

Accuracy is computed for the top 3L/2 GREMLIN predictions, with sequence separation ≥3, based on Cβ-Cβ amino acid specific distance as described in SI Table 3 in Kamisetty et al. (2013). The number of sequences after reducing the redundancy to 80% is shown. A set of 7047 pdb chains (see Supplemental file 3) was divided into two groups by length (less than 150 and greater than 400). (A) Larger proteins with similar number of sequence were less accurate then the smaller proteins. (B) #Sequences/length as often used does not accurately account for length dependence. There is a clear separation between the blue and green distributions. (C) #Sequences/√length better accounts for the length dependency. The blue and green distributions overlap. DOI: http://dx.doi.org/10.7554/eLife.09248.023

Figure 17.. The Rc metric used to assess fit of predicted contacts to a model.

The expected total GREMLIN score if the structure was native was estimated by summing sco*P(contact|sco, seq, len, sep) over all contacts with sep ≥6. To evaluate the fit of a particular model to a predicted contact set, we take the ratio of the actual total GREMLIN score of the model to the expected total score computed as above; we refer to this ratio of observed and expected contact scores as ‘Rc’. Blue line: the distribution of Rc in native structures with 4L–10L sequences; Red line: distribution of Rc after randomly reassigning contact predictions to structures. Rc values less than 0.7 are very infrequently observed for native structures; we use this value as a cutoff to evaluate the fit of a predicted contact set to a model. DOI: http://dx.doi.org/10.7554/eLife.09248.005