HIV-1 tropism prediction by the XGboost and HMM methods

Human Immunodeficiency Virus 1 (HIV-1) co-receptor usage, called tropism, is associated with disease progression towards AIDS. Furthermore, the recently developed and developing drugs against co-receptors CCR5 or CXCR4 open a new thought for HIV-1 therapy. Thus, knowledge about tropism is critical for illness diagnosis and regimen prescription. To improve tropism prediction accuracy, we developed two novel methods, the extreme gradient boosting based XGBpred and the hidden Markov model based HMMpred. Both XGBpred and HMMpred achieved higher specificities (72.56% and 72.09%) than the state-of-the-art methods Geno2pheno (61.6%) and G2p_str (68.60%) in a 10-fold cross validation test at the same sensitivity of 93.73%. Moreover, XGBpred had more outstanding performances (with AUCs 0.9483, 0.9464) than HMMpred (0.8829, 0.8774) on the Hivcopred and Newdb (created in this work) datasets containing larger proportions of hard-to-predict dual tropic samples in the X4-using tropic samples. Therefore, we recommend the use of our novel method XGBpred to predict tropism. The two methods and datasets are available via http://spg.med.tsinghua.edu.cn:23334/XGBpred/. In addition, our models identified that positions 5, 11, 13, 18, 22, 24, and 25 were correlated with HIV-1 tropism.

while CM constructs gapped PSSMs and takes the 11/25 rule and net charge into consideration. Recently, many genotypic methods based on machine learning have also been published. The method Geno2pheno 4 combines two machine learning approaches, support vector machine (SVM) and decision trees, and uses clinical information such as viral loads and CD4-cell counts if available. Another method from the same laboratory, G2p_str 23 , combines SVM and Lasso regression and uses the amino acid structure feature. Hivcopred 24 is based on SVM light with the split amino acid composition feature. T-CUP2 25 employs random forests (RFs) with the structural information of hydrophobicity and electrostatic potential. Currently, Geno2pheno is the most widely used method and the only genotypic method recommended for usage in clinical routines by the European Consensus Group 5,8 .
Genotypic methods can predict R5 viruses (~90%) accurately, but are inaccurate in the prediction of X4-using viruses (~50-70%) 26 . Thus, more accurate tropism prediction methods are required. Here, we present two methods, XGBpred and HMMpred. We analyzed the HIV-1 tropism prediction ability of our methods and compared them with the Geno2pheno, G2p_str, Hivcopred, CM and WebPSSM methods. The results show that XGBpred is robust with the hard-to-predict dual tropic sequences.

Methods
Datasets. To construct the Newdb dataset, we extracted 6790 R5 tropic, 590 X4 tropic and 1125 dual tropic sequences from the Los Alamos HIV sequence database (http://www.hiv.lanl.gov/, last update: 10 Sep 2017). The tropisms of the sequences from the Los Alamos HIV sequence database have been phenotypically determined, none of them have been inferred from sequences. Then we removed sequences containing non-canonical residues, reserved sequences with lengths between 31 and 39, and dislodged duplicated sequences to guarantee the high quality of genotype/phenotype pairs. This process finally generated 2335 R5 and 663 X4-using (245 X4 and 418 dual) tropic sequences. The distribution of the six major subtypes in the Newdb dataset is shown in Table 1. To compare our methods with the Geno2pheno, G2p_str, Hivcopred, CM and WebPSSM methods, we used the datasets constructed in these studies, respectively. These datasets can be accessed in Supplementary Spreadsheet S1. The distributions of tropisms in different datasets are shown in Table 2.
Machine learning method: XGBpred. Extreme gradient boosting (XGboost), like RFs used by T-CUP2 25 , is an ensemble algorithm of decision trees 27 . The ensemble works by combining a set of weaker machine learning algorithms to get an improved machine learning algorithm in overall. The main difference between XGboost and RFs is the way of sampling. RFs are based on uniform sampling with return. Instead, XGboost gives higher weights to the wrongly predicted samples in the current weaker leaner, and then these samples will be paid more attention when training the next weaker leaner. In addition, XGboost adds regularization to avoid overfitting. Therefore, XGboost is a more complicated algorithm than RFs, and thus always outperforms.
Because XGboost is designed for vectors, it is necessary to convert V3 loop string sequences of different lengths to numerical vectors. For this task, we used many kinds of features to describe the characteristics of protein sequences, such as split amino acid composition 24 , dipeptide composition 28 , and net charge or hydropath 29 . We also proposed an additional set of features: the alignment score. The 35-dimensional alignment scores were generated by scoring alignments using the block substitution matrices BLOSUM62, BLOSUM90 or BLOSUM100 30 , and the alignments were generated by aligning sequences to the consensus sequence with  and linear sequences. Just as the PSSM profile, HMM also can be used to describe protein families. The HMM profile described by state-transition and symbol-emission probabilities performs better than PSSM in terms of sequence alignment and homology recognition because it can deal with gaps in protein families better by hidden state chains 32 .
HMM profile construction. We used the maximum likelihood estimation method to establish R5 and X4-using specific HMM profiles from R5 and X4-using tropic multiple sequence alignments generated by ClustalO 33 , respectively. In addition, we simply assigned columns that had more than half gap characters as insertion states.
The structure of HMM that we used was no transition allowed from D j to I j or from I j to D j+1 (This kind of structure performed better than the full structure, as shown in Supplementary Table S1). M, D, and I denote match, deletion and insertion states, respectively.
Where in, k and l are indices over states M, D, or I; a is an amino acid symbol or gap; â kl means the estimated probability of transiting from state k to state l, ê a ( ) k means the estimated probability of emitting residue a at state k, and A kl and E k (a) are the corresponding frequencies. In order to avoid the zero probability which represents it cannot happen in the future, we applied the Laplace's pseudo-count rule that added one to each frequency.
Sequence-profile alignment. We employed Viterbi algorithm 34 , a dynamic programing algorithm, to get two alignment scores S R5 and S non-R5 . Those alignment scores represent the optimal state pathway scores from the R5 and X4-using HMM profiles, respectively. the final score was defined as: Then the given sequence would be classified as R5 tropic if the final score S is higher than a threshold, otherwise it would be classified as X4-using tropic.

Ten-fold cross validation.
The widely-used 10-fold cross validation was used to evaluate the performance of our methods in this study, where the sequences were divided into 10 subsets randomly, one subset was used as the testing set, and the others were used as the training set. After ten repetitions, the final performance was average of the performances of those ten subsets.
Evaluation parameters. For evaluation, we used sensitivity, specificity, accuracy and Matthew's correlation coefficient (MCC). In particular, MCC is robust even when the size of classes varies widely 35 . An MCC value '0' corresponds to a completely random prediction, while '1' corresponds to a perfect perdition. These parameters were calculated using the following equations: where TP is the number of true positives, FP false positives, TN true negatives and FN false negatives. We regarded R5 tropic samples as positives in this study.
In contrast to the four threshold-dependent parameters, the receiver operating characteristic (ROC) curve, a threshold-independent parameter, illustrates the trade-off between sensitivity and specificity at various threshold settings. In this study, we used the area under the curve (AUC) to measure a predictive power, where 0.5 means a random method, and 1 means a perfect method 36 . www.nature.com/scientificreports www.nature.com/scientificreports/

Results
Performance on the Newdb dataset. The feature set and the model that gave the strongest predictive power for the XGBpred and HMMpred methods were found, respectively (Supplementary Tables S1 and S2). The performances of the two methods on the Newdb dataset in a same 10-fold cross validation test are shown in Fig. 1A and Table 3. XGBpred had a higher specificity, accuracy, MCC and AUC than HMMpred when having the same sensitivity. Furthermore, the specificity of XGBpred was higher than 80% (84.62%) at the sensitivity of 91.78%. Results from the two methods were highly consistent: they predicted same tropisms for 87.96% of total samples, and achieved 96.70% sensitivity, 83.39% specificity and 93.93% accuracy.
Considering the poorer performance of HMMpred, the score distributions of the two methods were plotted (Fig. 1B). As depicted, the scores of dual tropic sequences mostly placed in the middle of the scores of X4 and R5 tropic sequences. Furthermore, HMMpred generated higher scores for a considerable number of dual tropic samples than XGBpred. This phenomenon illustrates that it is hard for dual tropic sequences to be correctly classified, especially by HMMpred.
The performances of the two methods for the six major subtypes (subtypes B, C, D, 01_AE, A and 02_AG) in the Newdb dataset were analyzed due to the sequence divergence among different subtypes and the different number of sequences in each subtype (Fig. 1C). HMMpred for subtypes B and D showed much lower AUCs (0.8942 and 0.8839) than for subtypes C and 01_AE (0.9369 and 0.9486). The reason was that subtypes B and D contained more hard-to-predict dual tropic sequences (Table 1). This also resulted in a low AUC (0.5887) for subtype 02_AG, and a higher AUC (0.9029) for subtype A than for subtype D (0.8839) generated by HMMpred. In contrast, the performance of XGBpred was not so deeply influenced by dual tropic sequences. XGBpred had higher AUCs for the top four most common subtypes (subtypes B, C, D and 01_AE) than for subtypes A and www.nature.com/scientificreports www.nature.com/scientificreports/ 02_AG. In addition, The V3 loops of subtypes 01_AE and 02_AG come from subtypes E and A, respectively 37,38 . This can also further lead to the weaker predictive power for subtypes A and 02_AG as it is a trickier task to determine tropism for subtype A than the other subtypes 39,40 . Besides, the large biases existed between the number of R5 and X4-using samples for subtypes C and A (Table 1). Therefore, we also reported the mean average precision (mAP) of the two classes (Fig. 1C). AP is the areas under the precision-recall curve for a certain class. The bigger the mAP is, the better the method preforms. The mAPs and AUCs demonstrated the same tendency for the predictive power of our methods. Among all subtypes, just as AUCs, XGBpred and HMMpred showed the highest mAPs (0.9646, 0.9491) for subtype 01_AE. Moreover, for both XGBpred and HMMpred, the divergences between AUCs and mAPs for subtypes C and A were biggest. This may arise from the large biases between the amount of R5 and X4-using samples.
Comparison with other methods. In this section, to evaluate our methods, we compared with the previously published methods Geno2pheno, G2p_str 23 , Hivcopred 24 , CM 22 , and WebPSSM 21 by implementing our methods in a 10-fold cross validation test on the datasets used in these published methods, respectively. The exception was WebPSSM 21 where we used the training set from WebPSSM to model our methods in a 10-fold cross validation test and used the validation set from WebPSSM to test (Table 3).
First when comparing with the Geno2pheno and G2p_str methods 23 , XGBpred and HMMpred achieved AUCs of 0.8952 and 0.9002, respectively. Our methods had higher AUCs than Geno2pheno (0.860) and G2p_ str (0.892). In addition, XGBpred and HMMpred achieved specificities of 72.56% and 72.09% at the sensitivity of 93.73%. The specificities were obviously higher than the specificities of Geno2pheno (61.6%) and G2p_str (68.6%) at the same sensitivity. Second, when comparing with the Hivcopred method 24 , XGBpred had a higher AUC (0.9483) than Hivcopred (0.904), but HMMpred had a low AUC (0.8829) as on the Newdb dataset. Third, when comparing with the CM method 22 . Our methods were as accurate as the CM method on the CM dataset which only contains a small amount of hard-to-predict dual tropic samples (Table 2). Finally, when comparing with the WebPSSM method 21 , although the WebPSSM dataset is small, XGBpred had a higher AUC (0.9043) than WebPSSM (0.881), and HMMpred presented a similar AUC (0.8678) with WebPSSM.
Feature importance analysis. Given the high performance of XGBpred presented in the previous subsections, we discussed which features XGBpred provided with its predictive power (Fig. 2). We did not analyze the feature importance on the WebPSSM dataset as it contains few training samples (Table 2). In the XGBpred method, the feature alignment score in the 5 th position of the V3 loop appeared in the top three most important features on all datasets. Interestingly, amino acid Tyr in position 5 appeared more frequently in X4-using tropic than in R5 tropic sequences (Supplementary Fig. S1). Currently, X4-using tropism can be predicted by the 11/25 rule 19 . However, since position 5 was as same important as positions 11 and 25, the pragmatic 11/25/5 rule was proposed to predict a virus as X4-using tropic by the presence of a positively charged amino acid in positions 11 or 25, or by the presence of amino acid Tyr in position 5 of its V3 loop. Compared with the 11/25 rule, the 11/25/5 rule reduced sensitivities by 1.29%, 1.03%, 1.14% and 1.19% on the Newdb, G2p_str, Hivcopred and CM datasets while increasing specificities by 7.39%, 5.11%, 6.34% and 10.77%, respectively. The 11/25/5 rule also had higher accuracies and MCCs than the 11/25 rule on the four datasets, which indicates the influence of amino acid Tyr in position 5 with regard to viral tropism (Supplementary Table S3 (3), 302 (7), 306 (11), 308 (13), 315 (18), 317 (20), 319 (22), 321 (24), 322 (25) and 328 (32) are important for tropism. Furthermore, the feature importance distribution generated by XGBpred is a feasible method to judge whether a new discovered association pattern is of importance to co-receptor usage or not.

Discussion
In this study, we present two methods, XGBpred and HMMpred, for HIV-1 co-receptor usage prediction. XGBpred is based on machine learning, and HMMpred is based on statistics. XGBpred performed best on the Hivcopred and Newdb datasets containing larger proportions of hard-to-predict dual tropic samples in the X4-using samples, while HMMpred performed worst. In contrast, the predictive powers of the two methods were similar on the smaller G2p_str and CM datasets containing fewer dual tropic samples ( Table 3). The poor ability of HMMpred to predict tropism stemmed from the high probability that HMMpred incorrectly predicted dual tropic samples as R5 tropic (Fig. 1B and Supplementary Fig. S2). The profiles used in HMMpred may not be meticulous enough. Several reasons may account for this phenomenon. Firstly, the two sequence families are highly similar since even one amino acid substitution may change their tropisms 42,43 . Secondly, the characteristics of dual tropic sequences may be overwhelmed by R5 and X4 tropic sequences. Finally, the unavailability of X4-using tropic samples makes it uncertain to learn its accurate HMM profile. Moreover, as the number of samples increased, the gap of predictive powers between XGBpred and HMMpred became large (Tables 2 and 3). This corresponds to the fact that the machine learning based Geno2pheno method is more widely used than the statistics based 11/25 rule and WebPSSM. As a result, a machine learning based method, in particular XGBpred, is recommended to predict co-receptor usage as the number of samples continues to expand.
In an effort to further increase the predictive power, we also generated three meta methods by the means of stacking 44 . The scores generated by XGBpred, Hivcopred (SVM light ) and HMMpred were added as additional features to the new stacking based XGBpred models. Compared with the original XGBpred method, the new stacking-based XGBpred methods had slightly higher AUCs on the G2p_str dataset but lower AUCs on the other datasets (Supplementary Table S4). The poor performances of the meta methods may due to the poorer predictive abilities of Hivcopred and HMMpred than the original XGBpred method, and/or the dependence of the results generated by XGBpred, Hivcopred and HMMpred (Supplementary Table S5). This may stem from the fact that the V3 loop is not the sole determinant of viral tropism. Moreover, V1, V2, C4 and the bridge sheet regions of gp120 also have an impact on co-receptor usage 45,46 . To predict tropism, several methods gain a higher accuracy by employing other information in addition to the V3 loop, such as clinical information 47 , V2 loop sequences 48 and structure information 23,25,41 . Therefore, the stacking based method can be constructed to improve its predictive power by combining methods with different kinds of information.

G2p_str
Hivcopred CM Feature important score Figure 2. Distribution of feature importance scores. The top 30 most important features indicated by XGBpred on the Newdb, G2p_str, Hivcopred and CM datasets. S# means the alignment score in position #, and R1R2 represents a dipeptide.