TAFPred: Torsion Angle Fluctuations Prediction from Protein Sequences

1. Introduction

Protein is an organic molecule of carbon, hydrogen, nitrogen, oxygen, and sulfur [1]. The core carbon atom is coupled to a side chain group, an amine group, a carbonyl group, and a hydrogen atom [2] to form a protein molecule. Protein molecules are essential and comprise many structures and functions within the cell. They also play an important role in the cell, forming many structures and performing numerous functions [3]. Protein molecules, such as actin and tubulin, can be structural and help shape the cell. They can also be functional, such as enzymes that perform critical metabolic reactions. A protein’s tertiary structure refers to how it folds in space in three dimensions. The polypeptide chain may require the assistance of chaperone proteins after it is synthesized at the ribosome [4,5]. The chaperone proteins also form temporary hydrogen bonds with the polypeptide chain, encouraging it to take the proper shape. This enables proper folding and, as a result, proper protein function. Some protein molecules do not fold to their native state and remain in a flexible state. Protein structural flexibility enables several types of dynamic and functional motions which are essential for interactions between proteins and peptides, DNA, RNA, or carbohydrate.

Protein structure can be illustrated by backbone torsion angles (Figure 1): rotational angles about the N-Cα bond (φ) and the Cα-C bond (ψ) or the angle between Cαi-1-Cαi-Cαi + 1 (θ) and the rotational angle about the Cαi-Cαi + 1 bond (τ) [6]. Prediction of Cα atom-based angle has demonstrated their potential usefulness in model quality assessment and structure prediction [7,8].

Several methods have been developed to predict backbone torsion angles. Angle predictions have been shown to be useful in fold recognition [9, 10] and fragment-based [11] or fragment-free structure prediction [12]. ANGLOR [13] utilizes support vector machines and neural networks for predicting the value of φ and ψ separately. TANGLE [14] uses a support vector regression method to predict backbone torsion angles (φ, ψ). Li et al. [15] predicted protein torsion angles by using four deep learning architectures, consisting of a deep neural network (DNN), a deep restricted Boltzmann machine (DRBN), and a deep recurrent neural network (DRNN), and a deep recurrent restricted Boltzmann machine (DReRBM). In addition, Heffernan et al. [7] captured the nonlocal interactions and yielded the highest reported accuracy in angle prediction by using long short-term memory bidirectional recurrent neural networks. A good prediction of angle probability may provide significant information on structural flexibility and intrinsic protein disorder in the extreme scenario [14]. In recent times, there have been notable advancements in the field of protein structure prediction using deep learning techniques. Specifically, AlphaFold [16], OmegaFold [17], and ESMFold [18] have demonstrated remarkable performance in predicting the three-dimensional (3D) structure of structured proteins. However, it should be noted that these methods excel primarily in predicting structured proteins [16]. On the contrary, the prediction of Phi and Psi angle fluctuations holds the potential for aiding in the prediction of unstructured or disordered protein structures.

However, to the best of our knowledge, only one research [19] presents work on backbone torsion angle fluctuation which is derived from the variation of backbone torsion angles. Because most proteins lack a known structure, the need for locating flexible (potentially functional) regions of a protein is the driving force behind the sequence-based prediction of torsion angle fluctuation. Moreover, using predicted torsion angles and flexibility as restraints can aid in protein structure and disordered region predictions. So, there is a dire need to improve the existing method to predict torsion angle fluctuations from protein sequences. The only method we found is developed by Zhang et al. [19]. They represented a neural network method for backbone torsion angle fluctuation based on sequence information only. Their model achieved ten-fold cross-validated correlation coefficients of 0.59 and 0.60 and the mean absolute errors (MAEs) of 22.7° and 24.3° for the angle fluctuation of φ and ψ, respectively.

In this work, we developed a machine learning method [20], TAFPred, to predict the backbone torsion angle fluctuation. Various features are directly extracted from protein sequences. A sliding window is used to include information from the neighbor residues. Furthermore, in TAFPred, we utilize a genetic algorithm (GA)-based feature selection method to extract several relevant features from the protein sequence. Finally, we train an optimized Light Gradient Boosting Machine to predict the backbone torsion angle fluctuation. We believe this is the second work that presents a sequence-based prediction method for backbone torsion angle fluctuation. We hope it will help further the protein structure and protein disorder prediction.

2. Materials and Methods

This section describes the dataset, feature extraction method, performance evaluation metrics, feature window selection, and finally, we describe the selected method for training the model. Figure 2 shows the workflow of the proposed TAFPred method.

2.1. Dataset

We collected 1268 protein chains from the author [19]. The protein chains are determined by the Nuclear Magnetic Resonance (NMR) method from the precompiled CulledPDB lists by PISCES using a sequence identity threshold of 25%. 997 protein chains are selected [19] by removing the chains with less than 5 NMR models, smaller than 25 amino acids, and consisting of nonstandard amino acid types. Finally, 936 protein chains are obtained by removing chains for which feature could not be obtained (referred to as NMR936) [21]. The backbone torsion angle fluctuation is derived from the variation of backbone torsion angles from different NMR models.

2.2. Feature extraction

We extracted several relevant profiles from the protein sequences, i.e., the Residue profile, Conservation profile, Physiochemical profile, Structural profile, and Flexibility profile. Here we briefly describe each of the profiles.

Residue profile. Twenty different numerical values are used to represent 20 standard amino acids (AA) types, yielding one feature per amino acid. The importance of this feature in solving bioinformatics problems has been shown in previous studies [22,23,24].

Physiochemical profile. In this work, five highly compact numeric patterns reflecting polarity, secondary structure, molecular volume, codon diversity, and electrostatic charge are extracted from [25] and used as features to represent the respective properties of each amino acid.

Conservation profile. The conservation profile of the protein sequence is obtained in terms of a normalized position-specific scoring matrix (PSSM) from the DisPredict2 program [24]. The PSSM is a matrix of L × 20 dimensions. The high scores suggest highly conserved locations, while scores around zero or negative indicate a less conserved position. The PSSM score was used to calculate monogram (MG) and bi-gram (BG) features. In terms of transition probabilities from one amino acid to another, the MG and BG properties can be used to characterize the portion of a protein sequence that can be conserved within a fold. From the DisPredict2 tool, we collect 1-D MG and 20-D BG characteristics.

Structural profile. Numerous biological problems have been solved using local structural features such as predicted secondary structure (SS) and accessible surface area (ASA) of amino acids. Here, the predicted ASA and SS probabilities for helix (H), coil (C), and beta-sheet (E) at the residue level is obtained from the DisPredict2 program. Moreover, we collect a separate set of SS probabilities for H, C, and E at the residue level from the BalancedSSP [26] program, as it provides a balanced prediction of these SS types. Thus, we extracted seven total structural properties (one ASA per amino acid and six predicted SS probabilities) as a structural profile of protein sequences.

Flexibility profile. Previous studies have demonstrated that an intrinsically disordered region (IDR) contains PTM sites, sorting signals and plays an important role in regulating protein structures and functions [27,28,29]. In this study, we used a disorder predictor named DisPredict2 [24] to accurately predict the protein’s disordered regions and obtain the disorder probability as a feature. To further improve the feature quality, we obtained two predicted backbone angle fluctuations, dphi (ΔΦ) and dpsi (ΔΨ), the DAVAR program [19].

Energy profile by Iqbal and Hoque [24], proposed a novel method that uses contact energy and predicted relative solvent accessibility (RSA) to estimate position-specific estimated energy (PSEE) of amino acid residues from sequence information alone. They showed that the PSEE could be used to distinguish between a protein’s structured and unstructured or intrinsically disordered regions. We utilized the PSEE score per amino acid as a feature in our study since it has been empirically demonstrated to have the ability to address a number of biological issues.

2.3. Machine Learning Methods

We analyzed the performance of eight individual regression methods: i) Light Gradient Boosting Machine Regressor (LightGBM) [30]; ii) Extreme Gradient Boosting Regressor (XGB) [31]; iii) Extra Tree Regressor (ET) [32]; iv) Decision Tree Regressor [33]; v) K-Nearest Neighbors Regressor [34,35]; vi) Convolutional Neural Network (CNN) [36]; and Long Short-Term Memory (LSTM) [37]; and Deep Neural Network (TabNet) [38]. The Light Gradient Boosting Machine Regressor (LightGBM) performs better, as shown in Table 2 and Table 3 under the Results Section.

2.4. Feature Selection Using Genetic Algorithm (GA)

We collected a feature vector of 179 dimensions (Figure 3) from different tools during the feature extraction process, which is significantly large. Therefore, to reduce the feature space and select the relevant features that could help improve the classification accuracy, we used a Genetic Algorithm (GA), a class of evolutionary algorithms, for feature selection. A detailed description of the feature selection approaches is provided below.

A GA is a population-based stochastic search technique that mimics the natural process of evolution. It contains a population of chromosomes, each representing a possible solution to the problem under consideration. In general, a GA operates by initializing the population randomly and iteratively updating the population through various operators, including elitism, crossover, and mutation, to discover, prioritize, and recombine good building blocks in parent chromosomes and finally obtain fitter ones [39,40,41].

Encoding the solution of the problem under consideration in the form of chromosomes and computing the fitness of the chromosomes are two important steps in setting up the GA. The length of the chromosome space is equal to the length of the feature space. Moreover, to compute the fitness of the chromosome, we use the LightGBM algorithm [33]. LightGBM was chosen because of its fast execution time and reasonable performance compared to other machine learning classifiers. During feature selection, the values of LightGBM parameters: max_depth, eta, silent, objective, num_class, n_estimators, min_child_weight, subsample, scale_pos_weight, tree_method, and max_bin were set to 6, 0.1, 1, ‘multi:softprob’, 2, 100, 5, 0.9, 3, ‘hist’ and 500, respectively and the rest of the parameters were set to their default value. The values of the LightGBM parameters mentioned above were identified through the hit-and-trial approach. In our implementation, the objective fitness is defined as:

$$ob{j}_{fit}=1-MAE+PCC$$

(3)

2.5. Performance evaluation

The performances of all the machine learning methods have been examined using a 10-fold cross-validation approach with the evaluation metric shown in Table 1. We measure the performance of torsion angle fluctuation predictions by calculating the Pearson Correlation Coefficient (PCC) and Mean Absolute Error (MAE) with the following equations:

3. Results

In this section, we first show the performance of different machine learning methods. Then, we present the performance of the best model with optimized hyperparameters. Next, we present the applied sliding window technique results to find the optimum window size. Finally, we compared the proposed method with the state-of-the-art method.

3.1. Comparison between different methods

We experimented with eight machine learning methods. The performance comparison of the individual regressors on the training dataset for phi angle fluctuation is shown in Table 2. Most of the methods perform better than the state-of-the-art method [19] except Decision Tree Regressor. Table 2 further shows that the LightGBM is the best-performing regressor among eight regressors implemented in our study in terms of mean absolute value (MAE) and Pearson correlation coefficient (PCC). Moreover, LightGBM improves by 6.59% and 24.50% in terms of MAE and PCC, respectively, compared to the existing method.

Table 2. Results from different machine learning methods (Phi Angle).

Methods / Metric	MAE	PCC	MAE (% imp.)	PCC (% imp.)	Average (% imp.)
State-of-the-art method [19]	0.126	0.598	-	-	-
Extra Trees Regressor	0.122	0.741	3.57%	23.88%	13.73%
XGB Regressor	0.123	0.727	2.67%	21.57%	12.12%
KNN Regressor	0.129	0.681	-2.30%	13.89%	5.79%
Decision Tree Regressor	0.167	0.527	-24.38%	-11.84%	-18.11%
LSTM	0.125	0.678	1.13%	13.35%	7.24%
CNN	0.166	0.608	-24.21%	1.68%	-11.27%
Tabnet	0.117	0.736	7.26%	23.09%	15.18%
LightGBM Regressor	0.118	0.745	6.59%	24.50%	15.54%

Best score values are boldfaced. Here, ‘imp.’ stands for improvement. The ‘% imp.’ represents the improvement in percentage achieved by TAFPred compared to the state-of-the-art method. Likewise, the ‘Average (% imp.)’ represents the average percentage improvement achieved by TAFPred for both MAE and PCC. Additionally, ‘(-)’ denotes that the % imp. or (Average % imp.) cannot be calculated.

Table 2. Results from different machine learning methods (Phi Angle).

Methods / Metric	MAE	PCC	MAE (% imp.)	PCC (% imp.)	Average (% imp.)
State-of-the-art method [19]	0.126	0.598	-	-	-
Extra Trees Regressor	0.122	0.741	3.57%	23.88%	13.73%
XGB Regressor	0.123	0.727	2.67%	21.57%	12.12%
KNN Regressor	0.129	0.681	-2.30%	13.89%	5.79%
Decision Tree Regressor	0.167	0.527	-24.38%	-11.84%	-18.11%
LSTM	0.125	0.678	1.13%	13.35%	7.24%
CNN	0.166	0.608	-24.21%	1.68%	-11.27%
Tabnet	0.117	0.736	7.26%	23.09%	15.18%
LightGBM Regressor	0.118	0.745	6.59%	24.50%	15.54%

Best score values are boldfaced. Here, ‘imp.’ stands for improvement. The ‘% imp.’ represents the improvement in percentage achieved by TAFPred compared to the state-of-the-art method. Likewise, the ‘Average (% imp.)’ represents the average percentage improvement achieved by TAFPred for both MAE and PCC. Additionally, ‘(-)’ denotes that the % imp. or (Average % imp.) cannot be calculated.

Table 3 compares the individual regressors’ performance for psi angle fluctuations. We also found that the LightGBM regressor performs best than other methods. LightGBM attains an MAE of 0.127 and PCC of 0.733. Moreover, The LightGBM Regressor improves by 6.59% and 24.50% in terms of MAE and PCC, respectively, compared to the state-of-the-art method.

Table 3. Results from different machine learning methods (Psi Angle).

Methods / Metric	MAE	PCC	MAE (% imp.)	PCC (% imp.)	Average (% imp.)
State-of-the-art method [19]	0.135	0.602	-	-	-
Extra Trees Regressor	0.131	0.729	2.77%	21.10%	11.94%
XGB Regressor	0.132	0.715	2.22%	18.73%	10.48%
KNN Regressor	0.139	0.670	-2.63%	11.24%	4.31%
Decision Tree Regressor	0.179	0.511	-24.65%	-15.11%	-19.88%
LSTM	0.132	0.665	2.29%	10.48%	6.38%
CNN	0.144	0.702	-6.46%	16.61%	5.07%
Tabnet	0.126	0.724	7.24%	20.28%	13.76%
LightGBM Regressor	0.127	0.733	6.09%	21.84%	13.96%

Best score values are boldfaced. Here, ‘imp.’ stands for improvement. The ‘% imp.’ represents the improvement in percentage achieved by TAFPred compared to the state-of-the-art method. Likewise, the ‘Average (% imp.)’ represents the average percentage improvement achieved by TAFPred for both MAE and PCC. Additionally, ‘(-)’ denotes that the % imp. or (Average % imp.) cannot be calculated.

Table 3. Results from different machine learning methods (Psi Angle).

Methods / Metric	MAE	PCC	MAE (% imp.)	PCC (% imp.)	Average (% imp.)
State-of-the-art method [19]	0.135	0.602	-	-	-
Extra Trees Regressor	0.131	0.729	2.77%	21.10%	11.94%
XGB Regressor	0.132	0.715	2.22%	18.73%	10.48%
KNN Regressor	0.139	0.670	-2.63%	11.24%	4.31%
Decision Tree Regressor	0.179	0.511	-24.65%	-15.11%	-19.88%
LSTM	0.132	0.665	2.29%	10.48%	6.38%
CNN	0.144	0.702	-6.46%	16.61%	5.07%
Tabnet	0.126	0.724	7.24%	20.28%	13.76%
LightGBM Regressor	0.127	0.733	6.09%	21.84%	13.96%

Best score values are boldfaced. Here, ‘imp.’ stands for improvement. The ‘% imp.’ represents the improvement in percentage achieved by TAFPred compared to the state-of-the-art method. Likewise, the ‘Average (% imp.)’ represents the average percentage improvement achieved by TAFPred for both MAE and PCC. Additionally, ‘(-)’ denotes that the % imp. or (Average % imp.) cannot be calculated.

3.2. Hyperparameters optimization

We optimized the LightGBM Regressor parameters: learning_rate, estimators, max_depth, num_leaves, max_bin, feature_fraction, etc., to achieve the best 10-fold cross-validation performance and for sampling hyperparameters and pruning efficiently unpromising trials. We have used the custom objective function as PCC+(1-MAE) for optimization. The best values of the parameters: learning_rate, estimators, max_depth, num_leaves, max_bin, and feature_fraction were found to be 0.014, 2561, 19, 380, 138, and 0.52, respectively.

3.3. Feature window selection

Here, we applied a widely used feature windowing technique to include the neighboring residue features. We examined a suitable sliding window size that determines the appropriate number of residues around a target residue that helps the model attain improved performance. We designed several models with different window sizes (ws) (1, 3, 5, and so on). We used the following custom metric as the objective function to measure the performance.

$$ob{j}_{fit}=1-MAE+PCC$$

(4)

Figure 4 shows the performance of the optimized LightGBM regressor for different window sizes for the Phi angle. The LightGBM regressor slightly improves window size 3, and the performance gradually decreases.

Figure 5 shows the performance of the optimized LightGBM regressor for different widow sizes for psi angle. The LightGBM regressor performance improves for window size 3, and then the performance gradually decreases. So, we select window size 3 to train the final model.

3.4. Comparison with the State-of-the-art method

Here, we compare the performance of the proposed method, TAFPred, with an existing state-of-the-art method [19] proposed by Zhang et al.. Table 4 shows that our proposed method improves by 10.08% in MAE, and 24.83% in PCC in the phi angle compared to the state-of-the-art method [19].

Table 5 shows that our proposed method improves by 9.93% in MAE and 22.37% in PCC in psi angle compared to the state-of-the-art method. Our proposed method significantly outperforms the existing state-of-the-art method and can more accurately predict the protein’s backbone torsion angle fluctuations.

4. Discussion

In this section, we explored diverse characteristics associated with the distribution of torsion-angle fluctuation. We examined the correlation between Δφ and Δψ, as well as the connection between torsion-angle fluctuation and disordered regions, utilizing our newly generated dataset.

The distribution of torsion-angle fluctuation

Figure 1 displays the distribution of torsion-angle fluctuation, with the dataset divided into 10 bins. The distributions are non-uniform, and the majority of residues exhibit angle fluctuations below 0.2. This observation indicates that stable protein structures are characterized by a limited presence of flexible residues.

Figure 6. The torsion-angle fluctuation is depicted in its distribution, with the data points divided into 10 bins. The fluctuations of the phi and psi angles are visually represented using red and green colors, respectively.

Figure 6. The torsion-angle fluctuation is depicted in its distribution, with the data points divided into 10 bins. The fluctuations of the phi and psi angles are visually represented using red and green colors, respectively.

Relationship between Δφ and Δψ

We further examined the relationship between the Δφ and Δψ angles, which represent the fluctuation of neighboring rotational angles in the protein backbone for the same residue. A chemical bond linkage results in a correlation between these angles, as it is impossible to alter one torsion angle without affecting the other. Consequently, a notable and significant correlation was observed between them. Consistently, the majority of residues exhibited minimal fluctuations below 0.2, as anticipated.

Figure 7. The relationship between Δφ and Δψ is shown in the figure, revealing that the majority of residues exhibit small fluctuations below 0.2.

Figure 7. The relationship between Δφ and Δψ is shown in the figure, revealing that the majority of residues exhibit small fluctuations below 0.2.

Relationship between torsion-angle fluctuation and disordered regions

We thoroughly investigated the connection between torsion-angle fluctuation and disordered regions. To gather disordered probability data, we utilized the SPOT-Disordered2 method. The figures provide clear evidence of the close relationship between phi and psi angle fluctuations and the presence of disordered regions. In the majority of samples, regions with low fluctuations exhibit a low disordered probability, while regions with higher fluctuations display relatively higher disordered probability, as illustrated in Figure 8 and Figure 9.

Figure 8. Relationship between torsion-angle fluctuation in the Phi angle and disordered regions. The disordered probability was obtained from the SPOT-Disordered2 tool. The figure illustrates that regions with low disordered probability exhibit correspondingly low fluctuations in the Phi angle, and conversely, regions with high disordered probability show higher fluctuations in the Phi angle.

Figure 8. Relationship between torsion-angle fluctuation in the Phi angle and disordered regions. The disordered probability was obtained from the SPOT-Disordered2 tool. The figure illustrates that regions with low disordered probability exhibit correspondingly low fluctuations in the Phi angle, and conversely, regions with high disordered probability show higher fluctuations in the Phi angle.

Figure 9. Correlation between torsion-angle fluctuation in the Psi angle and the presence of disordered regions. The disordered probability was obtained through the utilization of the SPOT-Disordered2 tool. The figure clearly illustrates that regions with low disordered probability exhibit lower fluctuations in the Psi angle, while regions with high disordered probability tend to have higher fluctuations in the Psi angle.

Figure 9. Correlation between torsion-angle fluctuation in the Psi angle and the presence of disordered regions. The disordered probability was obtained through the utilization of the SPOT-Disordered2 tool. The figure clearly illustrates that regions with low disordered probability exhibit lower fluctuations in the Psi angle, while regions with high disordered probability tend to have higher fluctuations in the Psi angle.

5. Conclusions

In this study, we explored eight machine learning methods, including a recently published Deep Neural Network (TabNet) [1], and found that the Light Gradient Boosting Machine Regressor (LightGBM) performs best based on MAE and PCC. To optimize LightGBM Regressor, we used state-of-the-art sampling and pruning algorithms for hyperparameter tuning. Moreover, a custom objective function is used for optimization, and a sliding window technique is used to extract more information from the neighbor residues for improved performance. Our proposed method, TAFPred, shows an average improvement of 15.54% and 13.96% in both metrics (MAE and PCC) on phi and psi angles, respectively, compared to the state-of-the-art method [19]. In the future, we also plan to investigate the effect of torsion angle fluctuation in disorder proteins. We strongly believe the developed method will be helpful to the researcher in protein structure prediction and disordered prediction.