Mathematical epidemiology of HIV/AIDS in cuba during the period 1986-2000

1.   Introduction

Hormone-binding protein (HPB) is a kind of protein that selectively and non-covalently binds to hormone. HPB is a soluble outer region of the growth hormone receptor (HR), and is an important component of the growth hormone (GH)-insulin-like growth factor axis ^[1]. The abnormal expression of HBP can cause a variety of diseases ^[2]. Due to the complex in vivo effects of HBP, its biological function is still not fully understood ^[1]. Therefore, accurate identification of HBP will be helpful to understand the molecular mechanisms and regulatory pathways of HBP.

Traditional methods to identify HBP were wet biochemical experiments, such as immunoprecipitation, chromatography, crosslinking assays, etc ^[3,4,5,6]. However, the disadvantages of these methods, such as time-consuming and expensive, make them are unable to keep up with the rapid growth of protein sequences in the post-genomic era. Therefore, it is necessary to develop automatics machine learning methods to identify HBP. As a pioneer work, Tang et al. developed a support vector machine-based method to identify HBP in which proteins were encoded using the optimal features obtained by adopting optimized dipeptide composition ^[7]. Subsequently, Basith et al. developed a computational predictor named iGHBP, in which an optimal feature set was obtained based on combining dipeptide composition and amino acid index value by adopting two-step feature selection protocol ^[8]. However, the overall accuracy was still far from satisfactory. In order to improve the performance for the identification of HBP, it is necessary to apply new feature extraction and selection methods to select optimal features to represent HBP.

In this paper, by examining 5 feature encoding methods and 2 feature selection methods, we investigated the advantages and disadvantages of various models for identifying HBP and then established a predictor called HBPred2.0 based on the optimal model. Finally, a user-friendly webserver was established for HBPred2.0. The paper is organized based on the following aspects (Figure 1): (1) The construction of benchmark dataset, (2) feature extraction and selection, (3) machine learning method, and (4) performance evaluation.

2.   Materials and methods

2.1.   Benchmark dataset and independent dataset

This paper adopted the benchmark dataset built by Tang et al. ^[7]. In the database, there are 123 hormone-binding proteins (HBPs) and 123 none hormone-binding proteins (non-HBPs). To verify the portability and validity of the model, we built a high quality independent dataset by obeying following rules. Firstly, we selected the 357 manually annotated and reviewed HBP proteins from Universal Protein Resource (UniProt) ^[9] using 'hormone-binding' as keywords in molecular function item of Gene Ontology. Subsequently, we excluded the proteins with sequence identity > 60% by using CD-HIT ^[10]. Thirdly, sequences that appear in the training dataset were excluded. As a result, 46 HBPs were obtained as independent positive samples. Negative samples were randomly selected from UniProt while using 'hormone' and 'DNA damage binding' as keywords in molecular function item of Gene Ontology, respectively. The sequence identities of negative samples are also ≤ 60%. Finally, 46 non-HBPs (37 hormone proteins and 9 DNA damage binding proteins) were randomly obtained. It should be noted that there is no similar sequences between the training and testing data. All data could be downloaded from http://lin-group.cn/server/HBPred2.0/download.html.

2.2.   Feature extraction methods

2.2.1.   Natural vector method (NV)

Suppose a sample protein P with L residues, it can be expressed as below.

where R_i represents the i-th amino acid residue of the sample protein P; i = (1, 2, …L). The Natural Vector Method (NV) method is briefly described as follows ^[11]:

For each of the 20 amino acid k, define:

where w_k(R_i) = 1, if R_i = k. otherwise, w_k(R_i) = 0.

Let n_k be the number of amino acid k in the protein sequence P, which can be calculated as:

Let s_(k)(i) be the distance from the first amino acid (regarded as origin) to the i-th amino acid k in the protein sequence. Let T_k be the total distance of each set of the 20 amino acids. Let μ_k be the mean position of the amino acid k. And they can be calculated as:

Let D₂^k be the second-order normalized central moments, which can be calculated as:

Thus, a sample protein P can be formulated as:

where the symbol T is the transposition of the vector.

2.2.2.   Composition transition distribution (CTD)

The CTD was first proposed for protein folding class prediction by Dubchak et al. in 1995 ^[12]. It's a global composition feature extraction method includes hydrophobicity, polarity, normalized van der Waals volume, polarizability, predicted secondary structure, solvent accessibility and so on. In this method, 20 amino acids were divided into 3 different groups: polar, neutral, and hydrophobic. For each of the amino acids attributes, three descriptors (C, T, D) were calculated. 'C' stands for 'Composition', which represents the composition percentage of each group in the peptide sequence, and thus can yield 3 features. 'T' stands for 'Transition', which represents the transition probability between two neighboring amino acids belonging to two different groups, and thus can yield 3 features. 'D' stands for 'Distribution', which represents the position (the first, 25%, 50%, 75%, or 100%) of amino acids in each group in the protein sequence, and thus can yield 5 features for each group (total 15 features).

In this paper, the sequence description of a sample protein P in term of hydrophobicity consists of 3 + 3 + 15 = 21 features.

2.2.3.   G-gap dipeptide composition (g-gap)

Adjacent dipeptide composition can only express the correlation between two adjacent amino acid residues. In fact, the amino acids with g-gap residues may be adjacent in three-dimensional space ^[13]. To find important correlations in protein sequences, we used the g-gap dipeptide composition that extends from adjacent dipeptides. A protein P can be formulated as below by using this method.

where the symbol T is the transposition of the vector; the v_i^g is the frequency of the i-th (i = 1, 2, …, 400) g-gap dipeptide and can be formulated as:

where n_i^g is the number of the i-th g-gap dipeptide; L is the length of the protein P; g is the number of amino acid residues separated by two amino acid residues.

In this paper, we studied the cases of g ranging from 1 to 9 because the case of g = 0 has been studied in reference ^[7].

2.2.4.   Pseudo amino acid composition (PseAAC)

The PseAAC method can not only include amino acid composition, but also the correlation of physicochemical properties between two residues ^[14,15]. In this paper, we adopted the type Ⅱ PseAAC, in which a sample protein P can be formulated as below.

where '9' is the number of amino acid physicochemical properties considered, namely, hydrophobicity, hydrophilicity, mass, pK1, pK2, pI, rigidity, flexibility and irreplaceability; 'λ' is the rank of correlation; 'x' is the frequencies for each element and is formulated as:

where ω is the weight factor for the sequence order effect; f_u is the frequency of the 400 dipeptides; τ_j is the correlation factor of the physicochemical properties between residues. More detailed information about the formula derivation process can be found in the reference ^[16].

In this paper, the parameter λ is from 1 to 95 with the step of 1, the parameter ω is from 0.1 to 1 with the step of 0.1. Therefore, 95×10 = 950 feature subsets based on PseAAC will be obtained.

2.2.5.   Tripeptide composition (TPC)

Tripeptide is composed of three adjacent amino acids in a protein sequence, which is a biosignaling with minimal functionality. By adopting TPC, a sample protein P can be formulated by:

where the symbol T is the transposition of the vector; the t_i is the frequency of the i-th (i = 1, 2, …, 8000) tripeptide and can be formulated as:

where n_i is the number of the i-th tripeptide; L is the length of the protein P.

2.3.   Feature selection methods

2.3.1.   Analysis of variance (ANOVA)

Feature selection is important to improve the classification performance. It can filter the noisy features ^{[17,18,19,20]}. We adopted the ANOVA method to select optimal features from g-gap dipeptide compositions and PseAAC. The ANOVA method calculated the ratio of the variance among groups and the variance within groups for each attribute ^[21,22]. The formula expressions can be described as follows:

where F(i) is the score of the i-th feature, a high F(i)-value means a high ability to identify the sample; S_w²(i) is the variance within groups; S_b²(i) is the variance among groups; and they can be calculated as follows:

where SS_b(i) is the sum of the squares between the groups; SS_w(i) is the sum of squares within the groups; K is the total number of classes; N is the total number of samples.

2.3.2.   Binomial distribution (BD)

We adopted the BD method to select optimal features from tripeptide composition ^[21]. In this algorithm, the confidence level (CL) of each feature can be calculated by:

where CL_ij is the confidence level for the i-th tripeptide in the j-th type; j denotes the type of samples (positive sample or negative sample); N_i is the total number of the i-th tripeptide in the dataset; the probability q_j is the relative frequency of type j in the dataset;

According to the formula as defined in Eq. (15), a high CL-value means a high ability to identify the sample. The BD method can extract the over-represented motifs, which is an excellent statistical method widely used in bioinformatics ^[23,24].

2.3.3.   Incremental feature selection (IFS) process

In general, if a model was built on a low-dimensional feature subset, it will not provide enough information. On the contrary, if a model was built on a high-dimensional feature subset, it can lead to information redundancy and overfitting problems. Therefore, the ANOVA and BD method with the IFS process and 5-fold cross-validation was applied to investigate the optimal feature set with the maximum accuracy ^[7,25,26,27] (Figure 2). We ranked all features according to the F(i)-values or CL-values and obtained new feature vectors, which are shown below.

The first feature subset contains the feature with the highest F(i)-value or CL-value, P^' = [g₁^']^T; By adding the second highest F(i)-value or CL-value to the first subset, the second feature subset P^' = [g₁^', g₂^']^T is formed. The procedure was repeated until all features were considered.

2.4.   Support vector machine (SVM)

The support vector machine (SVM) is a supervised machine learning method and has been widely used in bioinformatics ^{[28,29,30,31,32,33]}. Its main idea is to map the input features from low-dimensional space to a high-dimensional space through nonlinear transformation and find the optimal linear classification surface. For convenience, SVM software packages LibSVM can be download from https://www.csie.ntu.edu.tw/~cjlin/libsvm/. In the current study, the LibSVM-3.22 package was adopted to investigate the performance for identifying HBP. Besides, the radical basis function kernel was selected to perform predictions. The grid search spaces are [2^-5, 2¹⁵] with step of 2 for penalty parameter C and [2³, 2^-15] with step of 2^-1 for kernel parameter g.

2.5.   Performance evaluation

Three cross-validation methods, namely, the independent dataset test, the sub-sampling test, and the jackknife test, are widely used to investigate the performance of a predictor in practical application ^{[30,34,35,36,37,38,39,40,41]}. In order to save computing time, the 5-fold cross-validation test was adopted to calculate the optimal parameter C and g of SVM in this paper.

Five evaluation indexes were adopted to evaluate the models ^{[42,43,44,45,46,47,48,49]}. Sensitivity (S_n) is used to evaluate the model's ability to correctly predict positive samples. Specificity (S_p) is used to evaluate the model's ability to correctly predict negative samples. Overall Accuracy (Acc) reflects the proportion of the entire benchmark dataset that can be correctly predicted. The Matthew correlation coefficient (Mcc) is used to evaluate the reliability of the algorithm. Area under the ROC curve (AUC) reflects model's classification ability across decision values. They can be calculated as follows:

where TP, TN, FP, and FN represent the number of the correctly recognized positive samples, the number of the correctly recognized negative samples, the number of negative samples recognized as positive samples, and the number of positive samples recognized as negative samples, respectively.

3.   Results and discussion

3.1.   Performances of different features

In this study, we examined the performance of 5 feature extraction methods and their combinations. Based on CTD, NV, CTD+NV methods, protein samples can be expressed as 21-D (dimensional), 60-D and 81-D vector, respectively. The Accs of 60.16%, 70.33% and 67.07% were obtained by using SVM in the 5-fold cross-validation, respectively (as shown in Table 1). It was found that the prediction performances were far from satisfactory.

Based on the g-gap method, a protein sample can be expressed as a 400-D vector. By changing the value of g from 1 to 9, we obtained 9 feature subsets. Firstly, we investigated the performances of these 400-D features subsets based on SVM. The results were reported in Figure 3A. Subsequently, the ANOVA method with the IFS process was applied to investigate the optimal feature set, and the results were recorded in Figure 3B. One may notice that while g = 1, a maximum Acc of 80.89% was obtained when the top 144 features were used. Obviously, Accs were significantly increased by adopting ANOVA method. However, prediction performances still needed to improve.

Based on the PseAAC method, we obtained 95×10 = 950 (95 kinds of λ and 10 kinds of ω) feature subsets. Firstly, we investigated the performances of these 950 models by using SVM in the 5-fold cross-validation test and reported the results in Figure 4A. It was found that the maximum Acc of 76.83% was achieved when λ = 18 and ω = 0.1. In order to improve Acc, the ANOVA method was adopted to rank the 400 + 18 × 9 = 572 features. By adopting SVM with IFS, a maximum Acc of 84.15% was obtained when the top 194 features were used (Figure 4B). Although the result was encouraging, the Acc still has room to rise.

Based on the TPC method, 8000 features were extracted for each protein sequence. Considering that it would lead to overfitting problem, the BD method was adopted as the feature selection method. By adopting SVM with IFS process in the 5-fold cross-validation test, a maximum Acc of 97.15% was obtained when the top 1169 features were used (Figure 5). In this case, the S_n, S_p and Mcc are 96.75%, 97.56%, and 0.943, respectively The AUC reached 0.994, this result indicates that the performance of the model based on the optimal TPC is smart and reliable for identifying HBP.

3.2.   Comparison with other methods

In order to show the superiority of SVM to identify HBP, we compared its performance with those of other machine learning algorithms based on the same feature subset (i.e. 1169 optimal features). From Table 2, we can find that the SVM classifier could produce the best performance among these algorithms. Thus, the final model was constructed based on SVM.

It is also necessary to compare the methods proposed in this paper with existing methods. Table 3 shows the detailed results of different methods for identifying HBP. Based on the same benchmark dataset, Tang et al. achieved an Acc of 84.9% by using a SVM-based method, in which proteins sequences were encoded using the optimal 0-gap dipeptide composition features obtained by the ANOVA feature selection technique ^[7]. Basith et al. obtained an Acc of 84.96% in cross-validation test by training an extremely randomized tree with optimal features obtained from dipeptide composition and amino acid index values based on two-step feature selection ^[8]. Our proposed method could produce an Acc of 97.15% which is superior to the two published results, demonstrating that our method is more powerful for identifying HBP.

3.3.   Performance evaluation based on the independent dataset

For further comparing the performance of these methods, an independent dataset was used. The results were recorded in Table 4. One may observe that the HBPred2.0 predictor achieved the best performance among the three predictors, suggesting that HBPre2.0 has better generalization ability.

Specificity could reflect the discriminated capability of model on negative samples. From the Table 4, a higher specificity of the HBPred2.0 indicates that the model could produce less false positives.

4.   Conclusion

In this paper, we systematically investigated the performances of various features and classifiers on HBP prediction. By a great number of experiments, we obtained the best model by combining SVM with optimal tripeptide composition. This model could produce the overall accuracy of 84.78% on the independent data. Finally, Due to published database ^{[50,51,52,53]} and webserver ^{[54,55,56,57,58,59,60,61,62,63]} could provide more convenience for scientific community, we established a free webserver for the proposed method, called HBPred2.0, which can be free accessed form http://lin-group.cn/server/HBPred2.0/. We expect that the tool will help scholars to study the mechanism of HBP's function, and promote the development of related drug research.

Acknowledgments

This work was supported by the National Nature Scientific Foundation of China (61772119, 31771471, 61702430), Natural Science Foundation for Distinguished Young Scholar of Hebei Province (No. C2017209244), the Central Public Interest Scientific Institution Basal Research Fund (No. 2018GJM06).

Conflicts of interest

The authors declare that there is no conflict of interest.