Improved protein structure prediction using predicted interresidue orientations

Predicting interresidue geometries from MSAs using a deep neural network.

Unlike most other approaches to contact/distance predictions from MSAs, in addition to C_β–C_β distances, we also sought to predict interresidue orientations (Fig. 1A). Orientations between residues 1 and 2 are represented by 3 dihedral (ω, θ₁₂, θ₂₁) and 2 planar angles (φ₁₂, φ₂₁), as shown in Fig. 1A. The ω dihedral measures rotation along the virtual axis connecting the C_β atoms of the 2 residues, and θ₁₂, φ₁₂ (θ₂₁, φ₂₁) angles specify the direction of the C_β atom of residue 2 (1) in a reference frame centered on residue 1 (2). Unlike d and ω, θ and φ coordinates are asymmetric and depend on the order of residues (1–2 and 2–1 pairs yield different coordinates, which is the reason why the θ and φ maps in SI Appendix, Fig. S1 are asymmetric). Together, the 6 parameters d, ω, θ₁₂, φ₁₂, θ₂₁, and φ₂₁ fully define the relative positions of the backbone atoms of 2 residues. All of the coordinates show characteristic patterns (SI Appendix, Fig. S1), and we hypothesized that a deep neural network could be trained to predict these.

The overall architecture of the network is similar to those recently described for distance and contact prediction (3, 4, 7, 12). Following RaptorX-Contact (4, 12) and AlphaFold (7), we learn probability distributions over distances and extend this to orientation features. The central part of the network is a stack of dilated residual-convolutional blocks that gradually transforms 1- and 2-site features derived from the MSA of the target to predict interresidue geometries for residue pairs (Fig. 1B) with C_β atoms closer than 20 Å. The distance range (2 to 20 Å) is binned into 36 equally spaced segments, 0.5 Å each, plus one bin indicating that residues are not in contact. After the last convolutional layer, the softmax function is applied to estimate the probability for each of these bins. Similarly, ω, θ dihedrals and φ angle are binned into 24, 24, and 12, respectively, with 15° segments (+ one no-contact bin) and are predicted by separate branches of the network. Branching takes place at the very top of the network, with each branch consisting of a single convolutional layer followed by softmax. The premise for such hard parameter sharing at the downstream layers of the networks is that correlations between the different objectives (i.e., orientations and distance) may be learned by the network, potentially yielding better predictions for the individual features. We used cross-entropy to measure the loss for all branches; the total loss is the sum over the 4 per-branch losses with equal weight. Previous work (4) implicitly captured some orientation information by predicting multiple interresidue distances (C_β–C_β, C_α–C_α, C_α–C_g, C_g–C_g, and N–O), but in contrast to our multitask-learning approach, a separate network was used for each of the objectives. Our network was trained on a nonredundant (at 30% sequence identity) dataset from PDB consisting of 15,051 proteins (structure release dates before 1 May 2018). The trained network is available for download at https://github.com/gjoni/trRosetta.

We couple the derivation of residue–residue couplings from MSAs by covariance matrix inversion to the network by making the former part of the computation graph in TensorFlow (13). Sequence reweighting, calculation of one-site amino acid frequencies, entropies, and coevolutionary couplings and related scores take place on the GPU, and the extracted features are passed into the convolutional layers of the network (most previous approaches have precomputed these terms). We took advantage of our recent observation (14) that with proper regularization, covariance matrix inversion yields interresidue couplings (Methods) with only minor decrease in accuracy compared to pseudolikelihood approaches like GREMLIN (15) (the latter are prohibitively slow for direct integration into the network). Since the MSA-processing steps are now cheap to compute (compared to the forward and backward passes through the network during parameter training), this coupled network architecture allows for data augmentation by MSA subsampling during training. At each training epoch, we use a randomly selected subset of sequences from each original MSA, so that each time the network operates on different inputs.

Structure modeling from predicted interresidue geometries.

Following AlphaFold, we generated 3D structures from the predicted distances and orientations using constrained minimization (Fig. 1C). Discrete probability distributions over the predicted orientation and distance bins were converted into interresidue interaction potentials by normalizing all of the probabilities by the corresponding probability at the last bin (Methods) and smoothing using the spline function in Rosetta. These distance- and orientation-dependent potentials were used as restraints, together with the Rosetta centroid level (coarse-grained) energy function (16), and folded structures satisfying the restraints were generated starting from conformations with randomly selected backbone dihedral angles by 3 rounds of quasi-Newton minimization within Rosetta. Only short-range (sequence separation <12) restraints were included in the first round; medium-range (sequence separation <24) restraints were added in the second round, and all were included in the third. A total of 150 coarse-grained models were generated using different sets of restraints obtained by selecting different probability thresholds for inclusion of the predicted distances and orientations in modeling.

The 50 lowest-energy backbone + centroid models were then subjected to Rosetta full-atom relaxation, including the distance and orientation restraints, to add in side chains and make the structures physically plausible. The lowest-energy full-atom model was then selected as the final model. The structure generation protocol is implemented in PyRosetta (17) and is available as a web server at https://yanglab.nankai.edu.cn/trRosetta/.

Accuracy of predicted interresidue geometries.

We tested the performance of our network on 31 free-modeling (FM) targets from CASP13. (None of these were included in the training set, which is based on a pre-CASP PDB set.) The precision of the derived contacts, defined as the fraction of top L/n (n = 1, 2, 5) predicted contacts realized in the native structure, is summarized in Table 1 and SI Appendix, Table S1. For the highest probability 7.5% of the distance/orientation predictions (Fig. 2C), there is a good correlation between modes of the predicted distance/orientation distributions and the observed values (Fig. 2C): Pearson r for distances is 0.72, and circular correlation r_c (18) for ω, θ, and φ are 0.62, 0.77, and 0.60, respectively. The predicted probability of the top L long- + medium-range contacts correlates well (r = 0.84) with their actual precision (Fig. 2B). This correlation between predicted probability and actual precision allows us to further improve the results by feeding a variety of MSAs generated with different e-value cutoffs or originating from searches against different databases, into the network and selecting the one that generates predictions with the highest predicted accuracy.

Comparison with baseline network.

We evaluated our extensions to previous approaches by generating a baseline model to predict distances only, with no MSA subsampling and selection; the contact prediction accuracy of this network is comparable to previously described models (3, 12, 19, 20). Incorporating MSA subsampling during training and extending the network to also predict interresidue orientations improve contact prediction accuracy by 1.7 and 2.2%, respectively. Subsequent alignment selection improves performance an additional 3.1% on the CASP13 FM set (Table 1, last row). The improvements described above, together with increasing the number of layers in the network, increase the accuracy of predicted contacts by 7.6% over the baseline network on the CASP13 FM set. Although we ensured that there is no overlap between the training and test sets by selecting pre-CASP PDBs only (before 1 May 2018), our model was trained at a later date when more sequences were available; we also included metagenomic sequence data. Hence, we may be overestimating the gap in performance between our method and those used by other groups in CASP13; future blind tests in CASP will be important in confirming these improvements. Nevertheless, the gain in performance with respect to the baseline model is independent of the possible variations in the training sets and sequence databases. All of the targets in the Continuous Automated Model Evaluation (CAMEO) validation set below are more recent than both structural and sequence data in the training set.

Accuracy of predicted structure models.

We tested our method on the CASP13 FM targets, with results shown in Fig. 3. The average TM-score (21) of our method is 0.625, which is 27.3% higher than that (0.491) by the top Server group Zhang-Server (Fig. 3A). Our method also outperforms the top Human group A7D by 6.5% (0.625 vs. 0.587; Fig. 3B). The relatively poor performance on T1021s3-D2 (the outlier in the upper triangle of Fig. 3B) reflects the MSA-generation procedure: the majority of sequence homologs in the full-length MSA for T1021S3 only covers the first of the 2 domains; performance is significantly improved (TM-score increased from 0.38 to 0.63; the TM-score of A7D model is 0.672) using a domain-specific MSA. An example of the improved performance of our method is shown in Fig. 3C for the CASP13 target T0950; the TM-score of this model is 0.716, while the highest values obtained during CASP13 are: RaptorX-DeepModeller (0.56), BAKER-ROSETTASERVER (0.46), Zhang-Server (0.44), and A7D (0.43).

Fig. 3A deconstructs the contributions to the improved performance of the different components of our approach. When modeling is only guided by the distance predictions from the baseline network (no orientations and no MSA subsampling and selection; “baseline” bar in Fig. 2A), the TM-score is 0.537, lower than A7D but significantly higher than Zhang-Server and RaptorX. When predicted distances from the complete network are used, the TM-score increases to 0.592, higher than that of A7D. When the orientation distributions are included, the TM-score is further increased to 0.625. The folding is driven by restraints; very similar models are generated without the Rosetta centroid terms, and very poor models are generated without the restraints. To compare our Rosetta minimization protocol (trRosetta) to CNS (8), we obtained predicted distance restraints and structure models for all CASP13 FM targets from the RaptorX-Contacts server, which uses CNS for structure modeling (4), and used the distance restraints to generate models with trRosetta. The average TM-score of the trRosetta models is 0.45 compared to 0.36 for the RaptorX CNS models; the improvement is likely due to both improved sampling and the supplementation of the distance information with the general Rosetta centroid energy function.

Comparison between distance and orientation-based folding.

Both predicted distance and orientation can guide folding alone. The average TM-score of coarse-grained models for the CASP13 FM targets is 0.57 when folding with predicted orientation alone and 0.55 when folding with predicted distance only (SI Appendix, Fig. S2A). Relaxation improved the TM-score to 0.58 and 0.59 for orientation and distance guided folding, respectively (SI Appendix, Fig. S2B). The differences in quality of models generated using either source of information alone suggest that the 2 are complementary, and indeed better models are generated using both distance and orientation information (SI Appendix, Fig. S2).

Validation on hard targets from the CAMEO experiments.

We further tested our method on 131 hard targets from the CAMEO experiments (22) over the 6 mo between 8 December 2018 and 1 June 2019. The results for contact prediction are summarized in Table 1 and Fig. 2A; as in the case of the CASP13 targets, our method improves over the baseline network. The results for structure modeling are shown in Fig. 3D. The contributions of different components to our method are presented in SI Appendix, Fig. S4. On these targets, the average TM-score of our method is 0.621, which is 8.9 and 24.7% higher than Robetta and HHpredB, respectively. We note that the definition of “hard” is looser than the CASP definition; a hard target from CAMEO can have close templates in PDB. Making the definition of “hard” more stringent by requiring the TM-score of the HHpredB server to be less than 0.5 reduces the number of targets to 66. On this harder set, the TM-score for our method is 0.534, 22% higher than the top server Robetta and 63.8% higher than the baseline server HHpredB. Fig. 3E shows an example of a CAMEO target where our method predicts very accurate models (5WB4_H). For this target, the TM-scores of the template-based models by HHpredB, IntFOLD5-TS, and RaptorX are about 0.4. In comparison, the TM-score of our predicted model is 0.921, which is also higher than the top server Robetta (0.879).

Accuracy estimation for predicted structure models.

We sought to predict the TM-score of the final structure model using the 131 hard targets from CAMEO. We found that, unlike direct coupling-based methods such as GREMLIN, the depth of the MSA did not have a good correlation with the accuracy of the derived contacts. Instead, a high correlation (Pearson r = 0.90) between the average probability of the top-predicted contacts and the actual precision was observed (SI Appendix, Fig. S3A). The average contact probability also correlates well with the TM-score of the final structure models (r = 0.71; SI Appendix, Fig. S3B). To obtain a structure-based accuracy metric, we rerelaxed the top 10 models without any restraints. The average pairwise TM-score between these 10 nonconstrained models also correlates with the TM-score of the final models (r = 0.65; SI Appendix, Fig. S3C). Linear regression against the average contact probability and the extent of structural displacement without the restraints gave a quite good correlation between predicted and actual TM-score (r = 0.84; SI Appendix, Fig. S3D). We used this method to provide an estimated model accuracy.

Refinement of predicted models.

As noted above, CASP13 showed that protein structure-refinement methods can consistently improve models for cases where the sampling problem is more tractable (smaller monomeric proteins). We first evaluated the iterative hybridization protocol (23) previously used to improve models generated using direct contacts predicted from GREMLIN on the entire set of CASP13 and CAMEO targets (SI Appendix, Fig. S5). Incorporating our network-derived distance predictions resulted in consistent improvement in model quality when the starting model’s TM-score was over 0.7, in a few cases by more than 10% in TM-score. We also tested the incorporation of the network-derived distance restraints into the more compute-intensive structure refinement protocol we used in CASP13 (10) on the CASP13 FM targets with an estimated starting TM-score >0.6 that were not heavily intertwined oligomers and not bigger than 250 residues. Consistent improvements were observed on a set of 6 such targets (SI Appendix, Fig. S6), with an average TM-score improvement of about 4%. The net improvement in prediction for these targets using the combination of our structure-generation method and refinement using the distance predictions is indicated by the red points in Fig. 3B.

Training set for the neural network.

To train the neural network for the prediction of distance and orientation distributions, a training set consisting of 15,051 protein chains was collected from the PDB. First, we collected 94,962 X-ray entries with resolution ≤ 2.5 Å (PDB snapshot as of May first 2018), then extracted all protein chains with at least 40 residues, and finally removed redundancy at 30% sequence identity cutoff, resulting in a set of 16,047 protein chains with the average length of 250 amino acids. All of the corresponding primary sequences were then used as queries to collect MSAs using the iterative procedure described below. Only chains with at least 100 sequence homologs in the MSA were selected for the final training set.

Independent test sets.

Two independent test sets are used to test our method. The first is the 31 FM domains (25 targets) from CASP13 (first target released on 1 May 2018). The second one is from the CAMEO experiment. We collected 131 CAMEO hard targets released between 8 December 2018 and 1 June 2019, along with all of the models submitted by public servers during this period. Note that for the CASP13 dataset, the full protein sequences rather than the domain sequences are used in all stages of our method to mimic the situation of the CASP experiments.

MSA generation and selection.

The precision of predicted distance and orientation distribution usually depends on the availability of an MSA with ‘good’ quality. A deep MSA is usually preferable but not always better than a shallow MSA (see the examples provided in ref. 3). In this work, 5 alternative alignments are generated for each target. The first 4 are generated independently by searching the Uniclust30 database (version 2018_08) with HHblits (version 3.0.3) (28) with default parameters at 4 different e-value cutoffs: 1e⁻⁴⁰, 1e⁻¹⁰, 1e⁻³, and 1. The last alignment was generated by several rounds of iterative HHblits searches with gradually relaxed e-value cutoffs (1e⁻⁸⁰, 1e⁻⁷⁰,…, 1e⁻¹⁰, 1e⁻⁸, 1e⁻⁶, and 1e⁻⁴), followed by the hmmsearch (version 3.1b2) (29) against the metagenome sequence database (20) in case not enough sequences were collected at previous steps. The metagenome database includes about 7 billion protein sequences from the following resources: 1) JGI Metagenomes (7,835 sets), Metatranscriptomes (2,623 sets), and Eukaryotes (891 genomes); 2) UniRef100; 3) NCBI TSA (2,616 sets); and 4) genomes manually collected from various genomic centers and online depositories (2,815 genomes). To avoid attracting distant homologs at early stages and making alignment unnecessarily deep, the search was stopped whenever either of the 2 criteria were met: at least 2,000 sequences with 75% coverage or 5,000 sequences with 50% coverage (both at 90% sequence identity cutoff) were collected. The final MSAs for the test datasets are available at https://yanglab.nankai.edu.cn/trRosetta/.

Protein structure representation.

In addition to the traditional interresidue distance matrices, we also make use of orientation information to make the representation locally informative. For a residue pair (i, j), we introduce ω dihedral between C_α, C_β of one residue and C_β, C_α of the other, as well as 2 sets of spherical coordinates centered at each of the residues and pointing to the C_β atom of the other residue. These 6 coordinates (d, ω, θ_ij, φ_ij, θ_ji, φ_ji) are sufficient to fully define the relative orientation of 2 residues with respect to one another. Additionally, as described below, any biasing energy term defined along these coordinates can be straightforwardly incorporated as restraints in Rosetta.

Input features.

All of the input features for the network are derived directly from the MSA and are calculated on-the-fly. The 1D features include: 1) one-hot-encoded amino acid sequence of the query protein (20 feature channels), 2) position-specific frequency matrix (21 features: 20 amino acids + 1 gap), and 3) positional entropy (1 feature). These 1D features are tiled horizontally and vertically and then stacked together to yield 2 × 42 = 84 2D feature maps.

Additionally, we extract pair statistics from the MSA. It is represented by couplings derived from the inverse of the shrunk covariance matrix constructed from the input MSA. First we compute 1-site and 2-site frequency countsfi(A)=1Meff∑m=1MwmδA,Ai,m and fi,j(A,B)=1Meff∑m=1MwmδA,Ai,mδB,Aj,m, where A and B denote amino acid identities (20 + gap), δ is the Kronecker delta, indices i, j run through columns in the alignment, and the summation is over all M sequence in the MSA; w_m is the inverse of the number of sequences in the MSA, which share at least 80% sequence identity with sequence m (including itself); Meff=∑m=1Mwm. We then calculate the sample covariance matrix

and find its inverse (also called the precision matrix) after shrinkage (i.e., regularization by putting additional constant weights on the diagonal):

(More details on tuning the regularization weight in Eq. 2 are provided in SI Appendix, Fig. S8). The 21 × 21 coupling matrices si,jA,B of the precision matrix (Eq. 2) are flattened, and the resulting L×L×441 feature matrix contributes to the input of the network. The above couplings (Eq. 2) are also converted into single values by computing their Frobenius norm for nongap entries:

followed by the average product correction (APC):

where s.,j*, si,.*, and s.,.∗ are row, column, and full averages of the si,j* matrix, respectively. The coefficient in Eq. 2 was manually tuned on a nonredundant set of 1,000 proteins to maximize accuracy of the top L predicted contacts. From our experience, the final results are quite stable to the particular choice of the regularization coefficient in Eq. 2. To summarize, the input tensor has 526 feature channels: 84 (transformed 1D features) + 441 (couplings; Eq. 2) + 1 (APC score; Eq. 4).

Network architecture.

The network takes the above L×L×526 tensor as the input and applies a sequence of 2D convolutions to simultaneously predict 4 objectives: 1 distance histogram (d coordinate) and 3 angle histograms (ω, θ and φ coordinates). After the first layer, which transforms the number of input features down to 64 (2D convolution with filter size 1), the stack of 61 basic residual blocks with dilations are applied. Dilations cycle through 1, 2, 4, 8, and 16 (12 full cycles in total). After the last residual block, the network branches out into 4 independent paths—one per objective—with each path consisting of a single 2D convolution followed by softmax activation. Since maps for d and ω coordinates are symmetric, we enforce symmetry in the network right before the corresponding 2 branches by adding transposed and untransposed feature maps from the previous layer. All convolution operations, except the first and the last, use 64 3 × 3 filters; ELU activations are applied throughout the network.

Training.

We use categorical cross-entropy to measure the loss for all 4 objectives. The total loss is the sum over the 4 individual losses with equal weight (= 1.0), assuming that all coordinates are equally important for structure modeling. During training, we randomly subsample the input MSAs, uniformly in the log scale of the alignment size. Big proteins of more than 300 amino acids long are randomly sliced to fit 300 residue limit. Each training epoch runs through the whole training set, and 100 epochs are performed in total. Adam optimizer with the learning rate 1e⁻⁴ is used. All trainable parameters are restrained by the L₂ penalty with the 1e⁻⁴ weight. Dropout keeping probability 85% is used. We train 5 networks with random 95/5% training/validation splits and use the average over the 5 networks as the final prediction. Training a single network takes ∼9 d on one NVIDIA Titan RTX GPU.

Converting distance and orientation distribution to energy potential.

The major steps for structure modeling from predicted distributions are shown in Fig. 1C. For each pair of residues, the predicted distributions are converted into energy potential following the idea of Dfire (30). For the distance distribution, the probability value for the last bin, i.e., [19.5, 20], is used as a reference state to convert the probability values into scores by the following equation:

where p_i is the probability for the ith distance bin, N is the total number of bins, α is a constant (= 1.57) for distance-based normalization, and d_i is the distance for the ith distance bin. For the orientation distributions, the conversion is similar but without normalization, i.e.,

All scores are then converted into smooth energy potential by the spline function in Rosetta and used as restraints to guide the energy minimization. The range for distances is [0, 20 Å] with a bin size of 0.5 Å, while for orientations, the ranges are [0, 360°] for θ and ω, and [0, 180°] for φ, all with a bin size of 15°; corresponding cubic spline curves are generated from the discrete scores defined by Eqs. 5 and 6. For the distance-based potential, the AtomPair restraint is applied. For the θ- and ω-based potential, the Dihedral restraint is applied. For the φ-based potential, the Angle restraint is applied.

Quasi-Newton–based energy minimization and full atom-based refinement.

To speed up the modeling, coarse-grained (centroid) models are first built with the quasi–Newton-based energy minimization (MinMover) in Rosetta. A centroid model is a reduced representation of protein structure, in which the backbone remains fully atomic but each side chain is represented by a single artificial atom (centroid). The optimization is based on the L-BFGS algorithm (lbfgs_armijo_nonmonotone). A maximum of 1,000 iterations is used, and the convergence cutoff is 0.0001. Besides the restraints introduced above, the following Rosetta energy terms are also used: ramachandran (rama), the omega and the steric repulsion van der Waals forces (vdw), and the centroid backbone hydrogen bonding (cen_hb). More details about these energy terms can be found in ref. 16. The weights for the AtomPair, Dihedral, and Angle restraints, rama, omega, vdw, and cen_hb, are 5, 4, 4, 1, 0.5, and 1, respectively. The final models are selected based on the total score which includes both Rosetta energy and restraints scores.

The MinMover algorithm is deterministic but can be easily trapped into local minima. It is sensitive to the initial structure and restraints. Two strategies are proposed to introduce randomization effect, and those models trapped into local minima can be discarded based on total energy. The first strategy is to use different starting structures with random backbone torsion angles (10 are tried). The second strategy consists of using different sets of restraints. For each residue pair, we only select a subset of restraints with probability higher than a specified threshold (from 0.05 to 0.5, with a step of 0.1).

For each starting structure, 3 different models are built by selecting different subsets of restraints based on sequence separation s: short range (1 ≤ s < 12), medium range (12 ≤ s < 24), and long range (s ≥ 24). The first one is progressively built with short-, medium-, and long-range restraints. The second one is built with short- + medium-range restraints and then with long-range restraints. The last one is built by using all restraints together.

In total, 150 (= 10 × 5 × 3) centroid models were generated. The top 10 models (ranked by total energy) at each of the probability cutoff are selected for full-atom relax by FastRelax in Rosetta. In this relax, the restraints at probability threshold 0.15 are used together with the ref2015 scoring function. The weights for the AtomPair, Dihedral, and Angle restraints are 4, 1, and 1, respectively.

Supplementary File