Sparse Group Selection and Analysis of Function-Related Residue for Protein-State Recognition

2.1. Data Collection Using Molecular Dynamics Simulation

In this study, we employed the second PDZ domain (PDZ2) in the human PTP1E protein as the data source to develop feature selection and prediction models. The PDZ2 protein is a typical dynamics-driven allosteric protein upon binding with its allosteric effectors. This type of protein has extensively experimental and computational investigations in protein research. We adopted the PDZ2 protein as a test case to develop and evaluate the proposed structured sparse learning based feature selection and classification model to identify key features that potentially drive the allostery of the PDZ2 protein. As the developed machine learning method is a general approach, it can be conveniently used as an effective machine learning tool for other protein research studies.

The initial structures for PDZ2 protein are 3LNX and 3LNY for unbound and bound states, respectively. Both unbound and bound states are immersed in explicit water boxes using TIP3P model (ref), and sodium cations and chloride anions were added to the simulation boxes to neutralize the simulation boxes. Total of 13 independent sets simulations each with length 34ns were carried out for both state and subjected to machine learning model analysis. The canonical ensemble (NVT) Langevin MD simulations were used for the production run. For all simulations, 2 femtosecond (fs) step size was used and bond for hydrogen atoms were constrained. Frames were saved every 10 picoseconds. Periodic Boundary Condition (PBC) was applied in the simulations. All simulations were carried out using CHARMM simulation package version 40b1 and the CHARMM22 force field. For all the 34ns simulation, the initial 4 nanoseconds were treated as equilibrium. Therefore, each simulation set subjected to analysis is 30ns with 3000 frames extracted. Among 13 simulations of each state, ten simulations were randomly selected as training sets, and remaining three simulations were used as testing sets.

In previous work⁴⁵, there are 60,000 training observations and 18,000 testing observations. The training and testing observations were combined together as a single dataset for better evaluation with cross validation (as introduced later). The total number of features is 4743: the 1st − 4371th features are the inter-distance among the 94 residues, the 4372nd−4464th features are the sine values of Phi angles along the protein backbone, the 4465th − 4557th features are the cosine values of Phi angles along the protein backbone, the 4558th − 4650th features are the sine values of the Psi angles along the backbone, and the 4651st − 4743th features are the cosine values of the Psi angles along the backbone. Detail can be found in Table 1. To summarize, each residue can be assigned to one of the 19 groups, as shown in Table 2. PDZ2 structure with these 19 groups is illustrated in Fig. 1. In this study, the features are grouped into feature groups. The details of the 189 inter-residual groups can be found in Table 3.

2.2. Machine Learning Methods

2.3. Bayesian Hyperparameter Optimization

Bayesian Hyperparameter Optimization is used to tune the hyperparameters of sparse group lasso model and SVM. Comparing with the random or grid search, Bayesian hyperparameter optimization can efficiently conduct a search of a global optimization problem at finding the hyperparameters.

In the sparse group lasso model, the parameter λ₁ control the within-group sparsity of model vector, the parameter λ₂ control the between-group sparsity of model vector. It means that the sparsity of feature groups to be chosen can be controlled by adjusting this number. The parameters λ₁ and λ₂ are adjusted by the Bayesian hyperparameter optimization method with the training dataset using a five-fold cross validation to provide a realistic estimation of prediction errors and to prevent over-fitting. By adjusting the parameters(λ₁, λ₂), the number of features to be chosen in each group and the number of groups to be chosen can be controlled respectively while maintaining similar classification accuracy to the most optimal result. After running all possible combinations of parameters, the combination with the highest accuracy would be chosen. It can be a powerful application for chemists to design and consider different new design regarding protein residues.

2.4. Experimental Design

As shown in Fig. 3, we conduct standard 5-folds cross-validation on the entire dataset. In each time, the 4 folds were used as training set which is used for feature selection conducted by the sparse group lasso method. The training set can be further split into 5 folds to adjust the parameter λ₁ and the λ₂ by Bayesian hyperparameter optimization method. Based on the selected features by the sparse group lasso method, the SVM classification model are trained on the training set to make the final classification. The remaining one fold of the data set was used as testing set to calculate the classification accuracy, the sensitivity and the specificity with the trained SVM model.

2.5. Comparison with Baseline Methods

To illustrate the advantages of the proposed method, the method in the previous work⁴⁵ was compared. Feature selection using an extra-trees classifier in Scikit-Learn package⁶⁹ of Python was first applied on all 4743 features. After that, two other baseline prediction models were built and evaluated.

• Decision Tree: Often referred to as CART or Classification and Regression Trees, decision tree is a non-parametric machine learning algorithm. It uses a decision tree (as a predictive model) to go from observations about an item (represented in the branches) to conclusions about the item’s target value (represented in the leaves).

• Artificial neural network (ANN): An ANN is based on a collection of connected units or nodes called artificial neurons which loosely model the neurons in a biological brain. The word ”artificial” was added to differentiate the human neural network from the neural network for machine learning. In essence, ANN can be considered as a black-box learning algorithm. When appropriately tuned, ANN can outperform other machine learning algorithms if vast amount of data is available.

More details of the baseline methods can be found in the previous work⁴⁵.

3.1. Evaluation Criterion

There are four evaluation indexes used in this regard. 1) Accuracy (Acc) refers to the ratio between the number of correctly classified number of protein states and total number of protein states; 2) Sensitivity (Sen) stands for the proportion of positives that are correctly identified; 3) Specificity (Spe) refers to the proportion of negatives that are correctly classified number of protein states; and 4) Density of model vector β (or feature vector) is defined as the ratio between the number of non-zero elements of the feature vector and that of the length. Model is deemed to be good if accuracy, sensitivity and specificity are high and density of feature vector is low (but not too low that the necessary features are not included in the selection result).

Sen, Spe and Acc are calculated according to the following formulate:

where TP refers to the number of true positives, TN stands for the number of true negatives, FP refers to the number of false positives, FN stands the number of false negatives.

For an example of understanding what overall density and group density of model vector is, it is assumed that there are 10 features in 2 groups for differentiating an apple from an orange, and the selection result (in the form of model vector β) is as follows:

Since all elements in feature vector of group 1 (β₁) are all zero, feature group 1 is deemed to be irrelevant in differentiating an apple from an orange. On the other hand, feature group 2 deserves more attention as two elements out of the total five in the feature vector of group 2 (β₂) are not zero. Feature density of group 1 is 0, whereas feature density of group 2 is 2/5 = 0.4 (2 out of all 5 elements are non-zero). Last but not least, overall density is 2/10 = 0.2. The same is applied to the problem of protein-state recognition.

3.2. Performance of different features

Using the proposed method, we use the Bayesian optimization search method to adjust the parameter λ₁ and λ₂ in different range to control the sparsity and thus obtain various number of selected features. Given a certain range of the parameter λ₁ and λ₂, we can obtain the number of feature density and the corresponding accuracy, sensitivity and specificity. The results are summarized in Table 4. According to Table 4, as the number of features go up, the accuracy is higher. It can be observed that 28 selected features achieve 81.04% for classification accuracy, 79.90% for sensitivity and 82.18% for specificity. Furthermore, the classification accuracy achieves 91.50% when the number of selected features is 1936, and the corresponding sensitivity is 91.68%.

3.3. Performance comparison with Baseline methods

We compare our proposed method against the previous work⁴⁵. As shown in Table 5, we can see that the proposed sparse group lasso method can achieve better classification accuracy with less selected features comparing with baseline methods. When the number of selected features is 28, the accuracy of our method is improvement of approximately 1.04%. In addition, the accuracy of our method improves 11.5% when the number of selected features is 1936.

3.4. Feature analysis

In these studies, feature selection is more or less an implicit procedure to ensure the prediction quality and accuracy. Due to the large size of proteins with much more atoms in small molecules, including all the atomic distances as features for machine learning models is not feasible. Although in some special cases, some predefined reaction coordinates could be constructed manually as features. In general, a robust feature selection procedure would be desirable for many machine learning prediction models for protein simulation. However, more systematic and thorough investigation for feature selection of high dimensional system presented in the current study is still necessary.

There are 94 residues which can be categorized in 19 groups (Table 2). This work focuses on the feature groups which contribute the most to the differentiation of protein states. Therefore, if the density of a model vector for a particular feature group is significantly higher than zero, that particular feature group deserves more attention for further interpretation.

Using SGL method, total of 28 features were selected to achieve 81.4% of accuracy. 28 features are inter-residue distances, as shown in Table 6. The inter-residue distances are illustrated in Fig. 4. The distance features are located between Groups 1 (Loop 1) and 3 (Loop 2), Groups 1 (Loop 1) and 4 (Beta-strand 2), Groups 1 (Loop 1) and 8 (Beta-strand 3), Groups 1 (Loop 1) and 9 (Loop 5), Groups 1 (Loop 1) and 10 (Alpha-helix 2), Groups 1 (Loop 1) and 11 (Loop 6), Groups 2 (Beta-strand 1) and 14 (Beta-strand 5), Groups 2 (Beta-strand 1) and 16 (Alpha-helix 3), Groups 10 (Alpha-helix 2) and 15 (Loop 8), Groups 10 (Alpha-helix 2) and 16 (Alpha-helix 3), Groups 10 (Alpha-helix 2) and 17 (Loop 9), Groups 10 (Alpha-helix 2) and 18 (Beta-strand 6), Groups 10 (Alpha-helix 2) and 19 (Loop 10). In Table 6, it also shows the corresponding protein residues. Feature density of each group is calculated as the number of non-zero elements over the length of the model vector of the particular group. A higher feature density indicates that the importance of the group is more high. The final coefficient is multiple of coefficients from the 25 fold, which indicates the importance of each selected feature. For example, if α_i is the model vector from fold i, then the final coefficient element-wise products from all 25 folds ($\mathit{\alpha }={\mathit{\alpha }}_1\odot {\mathit{\alpha }}_2\odot {\mathit{\alpha }}_3\odot {\mathit{\alpha }}_4\odot \dots \odot {\mathit{\alpha }}_{25}$)

In Fig. 5, different colors of ideogram represents different groups. Each group includes different number of protein residues. It shows not only the inter-residual feature density (pairwise distance) but also the inter-group feature density for various combination of parameters with accuracy higher than 81%. Fig. 5(a), Fig. 5(b), Fig. 5(c), Fig. 5(d), Fig. 5(e), Fig. 5(f) indicate the visualization of inter-residual feature density and inter-group feature density when the number of selected features is 28, 35, 94, 158, 373 and 1936 respectively. The Inter-Group Feature Density illustrated in Fig. 5 provides straightforward representations for the important groups among all the groups. In addition, Inter-Residual Feature Density illustrated in Fig. 5 could provide much finer presentation about important feature distribution of the system.

On the whole, the inter-residual features between alpha-helix 2 and alpha-helix 3, alpha-helix2 and beta-strand 6 exist in most of the plots, indicating its importance in discriminating the protein states. The protein residues P1, K2, I6, N16, I20, G24, K38, P42, Q43, A45, A46, G50, S65, L66, E67, Q73, A74, V75, E76, T77, L78, R79, N80, V84, V85, H86, L87, L88, L89 from the above mentioned groups are critical for allostery of PDZ2 in various experimental studies. Such observation is indeed in agreement with the result of inter-group feature selection as in Table 6. With the increase of accuracy, the more inter-residual features are chosen.

3.5. Discussion

One of the innovations in the current study is grouping and labeling large number of atomic distances from protein (4371 for PDZ2 as the model system in this study) with relatively small number of secondary structures (19 for PDZ2, Fig. 1). This further labeling of the general features significantly simplified the data processing and interpretation of the machine learning model developed in this study. In the current study, Groups 1, 10, 16, and 18 are identified with the most contributions to the top important features (Table 6). Some of the residues from these Groups have been reported as critical for allostery of PDZ2 in various experimental studies, including Ala45 in Group 10,⁷⁰ and Q73, Ala74, Val75, E76, Thr77, Leu78, R79 in Group 16^70–73. In addition, there are more residues associated with high importance revealed in the current study that are also identified as key allosteric residues in the experimental characterization of PDZ2 allostery: N16, I20, G24, Q43, L66, N80^71,73.

The key residues identified through sparse group lasso also large overlap with findings from several computational studies. The most prominent observation is that the Alpha-helix 3 as Group 16 consisting residues Q73 through R79 highlighted in this study is also identified as containing many hot residues through Perturbation Response Scanning study²³, and part of allosteric path through Protein Structure Network (PSN) and Elastic Network Model (ENM)-based strategy.⁷⁴ Similarly, part of beta-strand 2 as group 4 (residues 20 and 24) and beta-stand 6 as group 18 (residues V84, V85, L87 and L88 ) are associated with important features in this study, and are also identified as hot residues,²³ or part of allosteric path⁷⁴. Numerous other residues associated with high importance from this study are also highlighted by these computational studies, including P1, K2, I6, K38.

These agreements between our refined feature selections and experimental as well as other computational studies about PDZ2 provide necessary and reassuring validation for this study, rendering our strategy as a general approach applicable for other allosteric proteins. These presentations based on the rigorous feature selection procedure presented in this study, provide comprehensive and innovative ways to build reliable and informative machine learning models for complicated biological systems.