Ramon et al. BMC Bioinformatics (2019) 20:410 https://doi.org/10.1186/s12859-019-2991-2 RESEARCH ARTICLE Open Access HIV drug resistance prediction with weighted categorical kernel functions Elies Ramon1* , Lluís Belanche-Muñoz2 and Miguel Pérez-Enciso1,3 Abstract Background: Antiretroviral drugs are a very effective therapy against HIV infection. However, the high mutation rate of HIV permits the emergence of variants that can be resistant to the drug treatment. Predicting drug resistance to previously unobserved variants is therefore very important for an optimum medical treatment. In this paper, we propose the use of weighted categorical kernel functions to predict drug resistance from virus sequence data. These kernel functions are very simple to implement and are able to take into account HIV data particularities, such as allele mixtures, and to weigh the different importance of each protein residue, as it is known that not all positions contribute equally to the resistance. Results: We analyzed 21 drugs of four classes: protease inhibitors (PI), integrase inhibitors (INI), nucleoside reverse transcriptase inhibitors (NRTI) and non-nucleoside reverse transcriptase inhibitors (NNRTI). We compared two categorical kernel functions, Overlap and Jaccard, against two well-known noncategorical kernel functions (Linear and RBF) and Random Forest (RF). Weighted versions of these kernels were also considered, where the weights were obtained from the RF decrease in node impurity. The Jaccard kernel was the best method, either in its weighted or unweighted form, for 20 out of the 21 drugs. Conclusions: Results show that kernels that take into account both the categorical nature of the data and the presence of mixtures consistently result in the best prediction model. The advantage of including weights depended on the protein targeted by the drug. In the case of reverse transcriptase, weights based in the relative importance of each position clearly increased the prediction performance, while the improvement in the protease was much smaller. This seems to be related to the distribution of weights, as measured by the Gini index. All methods described, together with documentation and examples, are freely available at https://bitbucket.org/elies_ramon/catkern. Keywords: HIV, Drug resistance prediction, Categorical kernel, Weighted kernel, PI, NRTI, NNRTI, INI, Machine learning, Support vector machine, Random Forest, Kernel PCA Background greatly thanks to antiretroviral therapy. Antiretroviral drugs HIV is a retrovirus that infects human immune cells, lower the virus level in blood by targeting different stages causing a progressive weakening of the immune system. of the virus life cycle. The most important classes of anti- When untreated, the affected person develops acquired retroviral drugs are protease inhibitors (PIs), which target immunodeficiency syndrome (AIDS), which leads to a the protease, and nucleoside and non-nucleoside reverse rise of opportunistic infections and, finally, death. HIV transcriptase inhibitors (NRTIs and NNRTIs, respectively) has infected more than 35 million people worldwide and which target the reverse transcriptase. Other classes of anti- is considered a global pandemic [1]. Despite the efforts, retroviral drugs are the integrase inhibitors (INIs) and the to date there is no definitive cure that eradicates the fusion inhibitors. virus from the organism. However, the lifespan and quality Some of the main reasons why HIV is so difficult to of life of many people that live with HIV have expanded fight are its short life cycle (1–2 days), high replication rate (108–109 new virions each day), and high mutation * Correspondence: elies.ramon@cragenomica.es rate (10− 4–10− 5 mutations per nucleotide site per repli- 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB cation cycle) caused because reverse transcriptase lacks Consortium, Campus UAB, 08193 Bellaterra, Barcelona, Spain Full list of author information is available at the end of the article proofreading activity. This permits the fast emergence of © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Ramon et al. BMC Bioinformatics (2019) 20:410 Page 2 of 13 new HIV variants, some of which may be resistant to the of cross-resistance information in the form of ensemble drug treatment [2]. These variants can be transmitted, methods has also been reported to improve resistance pre- and some studies show that ~ 10% of patients who had diction [20–22]. never been on antiretroviral therapy carry at least one Nevertheless, HIV sequence data specificities pose sig- resistant HIV [3]. Cross-resistance (simultaneous resist- nificant challenges to resistance prediction. First, se- ance to two or more drugs, often of the same class) is quence data is categorical in nature. However, most also a common phenomenon. It is therefore advisable to machine learning algorithms are designed to cope with do a resistance test before the treatment to find the best numeric data (DT and RF being exceptions), thus obli- drug choice [2, 4], especially in developing countries, as ging to perform some kind of pre-processing. A typical recommended by the WHO and the International AIDS approach is to recode each position into m or m − 1 Society-USA Panel [3]. A resistance test can be per- “dummy variables”, which can take the values 0 or 1 [5]. formed in vitro, obtaining HIV samples from the patient Usually, m is the number of all possible alleles that can and using them to infect host cells cultured in presence be potentially found in a position (i.e., m = 20 in protein of increasing levels of drug concentration. The virus sus- sequences). However, some authors restrict the dummy ceptibility is then obtained empirically as the IC50 [4] variables to the drug associated mutations already and usually delivered as the relative IC50 (resistance of appearing in the literature [6, 10, 12]. A very different the virus variant compared to the wild type). Another approach is found in [14], where each amino acid was strategy is to infer the HIV variant resistance from its se- codified as an integer ranging 1–22 (the 20 canonical quence. This can be either gene sequence or the trans- amino acids plus two extra characters B and Z). Other lated protein sequence; this latter approach eliminates encodings have been used with HIV sequence data, like the noise of synonymous mutations. In any case, as gen- amino acid composition frequencies, reduced amino acid ome sequencing is cheaper, faster and more widely available alphabets or physicochemical properties [5, 16, 20]. than performing an in vitro drug susceptibility test, much Another challenge is the presence of mixtures of al- effort has been invested in developing algorithms that pre- leles (normally two, rarely three or four) in at least one dict the drug resistance from the virus sequence [5]. position of the viral sequence for most clinical samples. The first attempts of automatic prediction can be In the case of HIV, this event indicates that the patient traced back, at least, to the early 2000s [6]. These carries two or more virus variants [4]. It is well estab- approaches were rule-based: study the mutational lished that HIV tends to generate viral swarms of closely profile of the HIV variant to look for known major related viruses (quasispecies), as a consequence of its drug-associated resistance mutations (lists of these high mutation rate [2]. Mixtures introduce ambiguity in mutations are periodically updated and can be found the genotype-phenotype correlation [6] and a problem in reviews, e.g., [7]). The rule-based algorithms con- of technical nature: the vast majority of machine learn- tinue to be used to this day because of their inter- ing methods are not able to deal directly with these pretability. Some publicly available examples are the “multiallelic” codes. To our knowledge, algorithms so far Stanford HIVdb, Rega or ANRS softwares [5]. How- have handled allele mixtures with some sort of previous ever, the aforementioned high mutation rate of the pre-processing of the data, e.g., keeping only the most HIV, which favors the emergence of large numbers frequent amino acid of the mixture [19], replacing the of new resistance mutations and complex mutational positions by a missing value [17], excluding the affected patterns, makes the rule-based approach suboptimal. sequences [15] or expanding the data to obtain all the In this scenario machine learning methods can be possible sequences that could be generated with the ob- extremely helpful, especially in recent years with the served mixtures [11, 14, 18]. increasing size of available data. This second ap- In this paper, we propose the use of kernel functions proach is also very popular and there exists machine specifically adapted to the aforementioned HIV data in- learning software to predict resistance online [8, 9]. tricacies, and able to integrate the relevance of the major Different methods have been proposed, the most resistance associated protein residues. Kernels are math- common ones being Linear Regression [10, 11], Arti- ematical functions with interesting properties. They can ficial Neural Networks (ANN) [10, 12–14], Support be coupled to numerous machine learning algorithms, Vector Machines (SVMs) [10, 15, 16], Decision Trees the so-called kernel methods, and provide a framework (DT) [10, 17] and their ensemble counterpart, Ran- to deal with data of virtually any type (e.g. vectors, dom Forests (RF) [15, 16, 18, 19]. Some machine strings, graphs). They can also encode complementary learning studies have complemented the sequence data knowledge about a problem, as long as some mathemat- with structural information, e.g., [11, 15, 16, 18], or have ical conditions are satisfied [23]. Our aim using kernel benefited from the knowledge about major drug associ- functions that address the aforementioned HIV data par- ated mutations to perform feature selection. The inclusion ticularities was not only to improve prediction, but also Ramon et al. BMC Bioinformatics (2019) 20:410 Page 3 of 13 reduce pre-processing, thus preserving the data integrity Table 1 Final number of HIV isolates per drug and lowering the risk of inserting spurious patterns. Drug Data Drug Data Drug Data size Drug Data size size size Methods ATV 1058 SQV 1603 DDI 1511 BIC 84 Datasets and data pre-processing DRV 665 TPV 766 TDF 1236 CAB 0 The Genotype-Phenotype Stanford HIV Drug Resistance FPV 1559 3TC 1482 EFV 1511 DTG 209 Database [24] is a public dataset with sequences from IDV 1607 ABC 1513 ETR 502 EVG 577 HIV isolates and its relative susceptibility to several anti- LPV 1372 AZT 1502 NVP 1517 RAL 636 retroviral drugs. We retrieved the PhenoSense dataset from Stanford webpage (version date: 2019-2-20). The NFV 1654 D4T 1510 RPV 181 data is split in four databases (PI, NRTI, NNRTI and INI), which contain between 1,000–3,500 HIV isolates. given by the inner product of two vectors in input space, INI is a new addition to the Stanford database and in- x and y: cludes some of the drugs most recently approved for therapeutic use. The complete dataset contains eight k Lin ðx; y Þ ¼ xT y ð1Þ protease inhibitors: atazanavir (ATV), darunavir (DRV), In our case, x and y represent the protein sequence of fosamprenavir (FPV), indinavir (IDV), lopinavir (LPV), two HIV isolates, recoded as dummy variables [25]. We nelfinavir (NFV), saquinavir (SQV) and tipranavir (TPV); used this kernel as the linear method of reference. An al- five integrase inhibitors: bictegravir (BIC), cabotegravir ternative expression is: (CAB), dolutegravir (DTG), elvitegravir (EVG) and ralte- gravir (RAL); and two classes of reverse transcriptase in- X d hibitors: six NRTIs, lamivudine (3TC), abacavir (ABC), k Lin ðx; y Þ ¼ w i xi y i ð2Þ i¼1 zidovudine (AZT), stavudine (D4T), didanosine (DDI) and tenofovir (TDF); and four NNRTIs, efavirenz (EFV), etra- where d is the length of the sequence. This expression virine (ETR), nevirapine (NVP) and rilpivirine (RPV). Se- stresses the possibility of assigning a weight wi to each quence length is 99 amino acids in the case of PI database, protein position, as it is known that not all positions 288 in the case of INI database and 240 in the case of contribute equally to the virus resistance [2]. Weights NRTI and NNRTI databases. The dataset contains the are nonnegative and sum to one. We considered two op- strain virus resistance (relative IC50) to each drug, and the tions: the simplest one was to consider that all positions sequence of the protein targeted by this drug. We built have the same importance, i.e., assigning equal weight 1/ the regression models for each drug separately, taking d to all variables. The second one was including additional each polymorphic protein position as a predictor variable information into the kernels, using RF mean decrease in and the drug resistance value as the target variable. Since node impurity as a metric for position importance. the distributions of resistances are highly skewed we used the log-transformed values, as recommended in [5]. Re- RBF kernel dundant viruses obtained from the same patient were re- It is a nonlinear kernel, usually defined as: moved to minimize bias. We deleted all sequences 2 affected by events that changed protein length (protein k RBF ðx; y Þ ¼ e−γ jjx−yjj ð3Þ truncations, insertions and deletions). These events were Where ||x − y|| is the squared Euclidean distance be- 2 uncommon in the dataset and affected less than 5% of tween two vectors, and γ > 0 is a hyperparameter. As in HIV sequences. Also, we removed all isolates with one or the case of the linear kernel, the original data was more missing values. Missing values are present in the tar- recoded. We also introduced the possibility of weighting get variables as well as in the sequences, because not all the positions: HIV isolates have been tested for all drugs. The final num- Pd ber of data instances for each drug is shown in Table 1. 2 k RBF ðx; y Þ ¼ e−γ i¼1 wi ðxi −yi Þ ð4Þ To ensure a minimum of data rows for training/test parti- tions and cross-validation, we did not consider drugs with The RBF kernel is a widely accepted default method a sample size lower than 100. [23, 25], so we used it as a benchmark to compare with the categorical kernels. Methods We compared the performance of a nonlinear, nonkernel Overlap kernel method (RF) to a kernel method: SVMs. SVMs can be This is the most basic categorical kernel. This kernel as- either linear or nonlinear, depending on the kernel used. signs 1 if the two instances compared are equal and 0 The linear kernel is the simplest of all kernel functions, otherwise. Ramon et al. BMC Bioinformatics (2019) 20:410 Page 4 of 13  1 if xi ¼ yi prediction, we implemented the stacking algorithm de- k Ov ðxi ; yi Þ ¼ ð5Þ 0 if xi ≠yi scribed in [22] for continuous outcomes. This meta-learner approach consists of two principal steps. In the first step, where xi and yi represent the alleles of a given protein single drug models are built from the training data as usual. position i in two HIV sequences, x and y. In the second step, the fitted values (i.e. predictions of the training data) of all drugs obtained in step 1 are used as in- Jaccard kernel put to a new (stacked) model, being each drug a different The Jaccard index measures the similarity between two predictor. The method that integrates the single drug finite sets and is a valid kernel function [26]. We used it models in step 2 and delivers the definitive predictions is to handle allele mixtures, while in the rest of methods called a combiner algorithm. Data size largely varied be- we randomly sampled one allele of the mixture. Letting tween drugs (see Table 1), even within the same drug class, again i denote a given protein position (so that Xi and Yi so we chose Decision Trees (DT) as our combiner algo- are non-empty sets of alleles in the i-th position for iso- rithm, as they can easily handle missing data. We combined lates x and y) then: the drugs within the same database (PI, NRTI, NNRTI and INI) and applied this stacking methodology to our previ- j X i ∩Y i j k Jac ðX i ; Y i Þ ¼ ð6Þ ously proposed weighted kernels (Linear, RBF, Overlap and j X i ∪Y i j Jaccard). When ∣Xi ∣ = ∣ Yi ∣ = 1, i.e., none of the individuals have an allele mixture at that i-th position, Jaccard re- Experimental setup and model tuning duces to the Overlap kernel. Unlike Overlap, the Jaccard To assess the performance of the methods used, each data- kernel can deal simultaneously with allele mixtures and base was split at random in two partitions: training set categorical data. (60% of the database) and test set (40%). Hyperparameter optimization was done by a 10 × 10 cross-validation on the “RBF-like” categorical kernels training set. Once the optimum hyperparameter was found, For the whole protein sequences, we can aggregate all the final model was built using the whole training set. To single position Overlap and Jaccard evaluations as the assess the model performance, the NMSE (Normalized convex combination of kernels evaluations (Eq. 5 or 6) Mean Square Error) between the actual and the predicted and position weights. This results in a valid kernel func- drug resistances of the test set was computed: tion, since the product of a positive scalar and a kernel is a kernel, and the sum of kernels is also a kernel. To ensure that the only difference between categorical ker- NMSE ðobserved; predicted Þ P nels and RBF was the categorical part, we introduced an ðobserved−predicted Þ2 exponential factor and the hyperparameter γ, in a way ¼ ð8Þ ðN−1Þ  var ðobserved Þ analogous to (3) and (4): Pd k cat ðx; y Þ ¼ e−γ eγ i¼1 wi k ðxi ;yi Þ ð7Þ NMSE can be understood as the fraction of target variance not explained by the model. This is also a valid kernel function, since the exponen- We repeated the whole process 40 times, each time tial of a kernel gives another kernel, and where e−γ nor- with different 60/40 randomly split training/test parti- malizes the kernel matrix, keeping the evaluations tions, to obtain an error distribution. Kernel position between 0 and 1. The final versions of the Overlap and weights were calculated using the training set only. the Jaccard kernels are obtained replacing the k(xi, yi) Note that only the Jaccard kernel can directly handle term by (5) or (6), respectively. In our analyses, we com- allele mixtures; for the rest of kernels and the RF, we pared weighted and unweighted versions for all linear, generated 40 versions of the database randomly sam- RBF, Overlap and Jaccard kernels. Thus we can ensure a pling one allele at a time. Then, the 40 replicates fair comparison between the categorical and the nonca- were used to compute all the models except Jaccard, tegorical kernels. which could deal directly with the database without further preprocessing. This way we can ensure an Stacked models honest comparison between Jaccard and the rest of So far, we have built prediction models for each inhibitor kernels and methods. separately. As mentioned in the Introduction, it is reported All analyses were implemented in the R statistical that there exists some degree of relationship between the re- computing language [27]. A documented package imple- sistance of different drugs (e.g. in case of cross-resistance). menting these methods is available at https://bitbucket.org/ To check whether the use of this information can improve elies_ramon/catkern/. Ramon et al. BMC Bioinformatics (2019) 20:410 Page 5 of 13 Visualization Figure 1 shows the NMSE distribution boxplot for four Kernel PCA is a kernel method obtained by coupling ker- representative drugs: FPV (PI database, panel a), DDI nel functions to a Principal Components Analysis. We (NRTI database, panel b), NVP (NNRTI database, panel c) used the Jaccard kernel PCA to visually check whether se- and EVG (INI database, panel d). The remaining 17 box- quences that are considered more similar by the kernel plots can be found at Additional file 1: Figures S1-S17. function are also similar in their drug resistance. As this method is for visualization purposes only, we did not sep- Performance overview arate training and test sequences. Thus, we used the mean NMSE varied greatly across drugs and methods. The kernel weights of the 40 training sets to compute the best prediction was achieved for 3TC, with an average weighted Jaccard. NMSE ranging 0.07–0.16 depending on the method used To check whether the important protein positions (i.e. (Additional file 1: Figure S8). The drug with worst pre- kernel weights) detected by RF could have an structural diction error was DTG, with an average NMSE ranging relevance, we highlighted our top ranking positions on 0.65–0.75 (Additional file 1: Figure S16). This was also the tridimensional structure of the protein. Pictures of the second drug with lowest data size (Table 1). Not un- protein-drug complexes were generated with Molsoft expectedly, methods applied to drugs with low N had ICM-Browser v.3.7–2 using structural data obtained considerably worse performance overall (especially DTG, from RCSB Protein Data Bank. RPV, ETR and TPV, but also TDF and to some extent DRV). In the PI database, errors were fairly similar Performance comparison to other approaches across all drugs and around 0.12–0.20 on average (e.g. We compared our SVM plus weighted Jaccard with the Figure 1a), with the sole exception of TPV, with an aver- ANN approach described in [14], which to our knowledge age NMSE ranging 0.30–0.45. In turn, predictive perfor- achieves the best performance so far in this dataset. We mances for the integrase and reverse transcriptase used the R interface to keras to implement the ANN. First, inhibitors were far more variable across drugs. Overall, we followed the specifications described in [14] about the the best method was the SVM with the Jaccard kernel range of candidate architectures (1–3 hidden layers, with (either in its weighted or in its unweighted version), 2–10 nodes per layer, for all drugs), number of epochs and which achieved the best performance in 20 out of 21 early stopping. As our dataset version and data pre-process- drugs. ing differ from [14], we also evaluated a different range of hyperparameters: three fixed ANN architectures (one hid- Unweighted case den layer with 30 nodes, two hidden layers with 20 and 10 Nonlinear kernels performed much better than the lin- nodes respectively, and three hidden layers with 30, 20 and ear kernel in almost all drugs, with the only exception of 10 nodes) with the L2 regularization parameter λ. Both ap- ETR and D4T. Categorical kernels outperformed RBF, proaches (from now on referred to as ANN1 and ANN2) although RBF was close to Overlap (or even marginally were trained and tested as for the rest of methods (see: better) in some cases. Among categorical kernels, the Data and dataset pre-processing), with the previously de- Jaccard kernel performed better than Overlap in all in- scribed 40 replicates, allele mixture treatment, training/test hibitors, sometimes by a large margin, as in the cases of ratio and 10 × 10 cross-validation to choose the best num- SQV, 3TC, AZT, EFV, NVP, RAL or EVG (Fig. 1 c and ber of layers and nodes per layer (in the case of ANN1) or d). Predictive performances of unweighted kernels and λ (in the case of ANN2). We chose the best architecture of RF were markedly different in protease with respect obtained in training within ANN1 and ANN2 options for to integrase and transcriptase inhibitors. RF was consist- each drug. ently worse than kernel methods for the PI database (e.g. Figure 1a), whereas RF performance was comparable or Results better than those of kernel methods in both reverse tran- As expected, HIV protein sequences showed a large vari- scriptase and integrase inhibitors (e.g. Figure 1b, c and d). ability. As many as 93% of the protease positions were polymorphic and, among these, the number of different Weighted case observed alleles varied between 2 and 16. In the case of Figure 2 shows three representative examples of the reverse transcriptase, 89% of the positions were poly- weights obtained from RF. The remaining plots are morphic and the number of alleles per polymorphic pos- shown in Additional file 2: Figures S18-S35. We ascer- ition ranged between 2 and 14. Integrase was the least tained that RF detected most of the major resistance variable protein: 75% of the positions were polymorphic associated positions described in the literature (e.g. re- and, in these positions, the number of alleles ranged be- view in [2]). Overall, a higher percentage of relevant po- tween 2 and 8. Almost 60% of the sequences had at least sitions were identified in protease inhibitors than in one allele mixture. both reverse transcriptase and integrase inhibitors. To Ramon et al. BMC Bioinformatics (2019) 20:410 Page 6 of 13 a FPV (PI) b DDI (NRTI) 0.40 Methods 0.20 RF RBF Jaccard Linear Overlap 0.35 Unweighted Weighted 0.18 NMSE 0.30 0.16 0.14 0.25 0.12 0.20 RF LIN RBF wLIN wRBF RF LIN RBF wLIN wRBF OV JAC wOV wJAC OV JAC wOV wJAC c NVP (NNRTI) d EVG (INI) 0.40 0.30 0.35 0.25 0.30 NMSE 0.25 0.20 0.20 0.15 0.15 0.10 0.10 RF LIN RBF wLIN wRBF RF LIN RBF wLIN wRBF OV JAC wOV wJAC OV JAC wOV wJAC Fig. 1 NMSE distributions for a PI (FPV, panel a), an NRTI (DDI, panel b), an NNRTI (NVP, panel c) and an INI (EVG, panel d). Note that NMSE scale varies between panels evaluate this numerically, we computed the Gini index of the error in RAL and EVG drugs but not in DTG. In the RF importance distributions for each of the drugs. This summary, Jaccard was the best weighted kernel followed index is shown in Fig. 2 and Additional file 2. We also no- by Overlap, RBF and Linear. ticed differences regarding the location of the important positions in the tridimensional structures of protease Factors affecting prediction error (Fig. 3a) and reverse transcriptase (Fig. 3b). The most im- To investigate the relevance of each factor in prediction, portant protease positions according to RF are distributed we fitted the following linear model to NMSE obtained over the whole structure, whereas in the case of the reverse in each replicate across all kernels and drugs (40 repli- transcriptase they are located at the drug binding site. cates × 21 drugs × 8 kernels): As for predictive performance, weighting was more ef- fective in integrase and reverse transcriptase inhibitors NMSE  N þ K þ W þ GINI þ ε ð9Þ than in protease inhibitors. In NRTI and NNRTI data- bases, weighted kernels outperformed RF in all cases, where N is the drug data size (Table 1), K is a class vari- whereas their unweighted counterparts did not. This able with the kernel used (Linear, RBF, Overlap or was particularly the case for 3TC, DDI (Fig. 1b), EVG Jaccard), W = 0 or 1 depending on whether the kernel (Fig. 1d) and especially NVP (Fig. 1c), where weighting was unweighted or weighted, respectively, and GINI is decreased the Jaccard kernel error by around 50%. In the standardized Gini index of RF weights. Table 2 sum- contrast, the effect of weighting was less marked in the marizes the coefficients and their significance. We found PI database: similar errors were obtained (e.g. Figure 1a) that all factors are significant and behave additively (in- for all drugs but TPV, where the error actually increased. teractions were not significant; results not shown). As In the INI database, weighting decreased dramatically expected NMSE decreases with N but, interestingly, also Ramon et al. BMC Bioinformatics (2019) 20:410 Page 7 of 13 a FPV 0.4 RF weights 0.3 * Mean Gini index = 0.85 * ** * 0.2 * ** *** * 0.1 0.0 0 10 20 30 40 50 60 70 80 90 100 b NVP 0.4 0.3 * Mean Gini index = 0.93 RF weights * * 0.2 * **** 0.1 0.0 * ** 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 c EVG 0.4 Mean Gini index = 0.93 0.3 * RF weights 0.2 * * * * *** 0.1 0.0 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210 220 230 240 250 260 270 280 290 Protein position Fig. 2 RF relative importance of each protein position for three drugs: a protease inhibitor (a), a reverse transcriptase inhibitor (b) and an integrase inhibitor (c). Standard error across the 40 replicates is marked with error bars. Asterisks highlight the major drug related positions reported in the literature [2] with Gini index, i.e., prediction improves when there are Sample size is one of the most important factors in only a few positions of large effect. Categorical kernels any experimental design, and the main one influencing were consistently better than noncategorical ones and total cost. Figure 4 panels d, e and f show the residuals Jaccard was the best option in all cases. Weighting pro- of model NMSE ~ K + W + GINI vs. N. Although Table 2 tein positions significantly lowers the error, although only shows that the NMSE decreases with sample size for all in reverse transcriptase and integrase inhibitors (as also drugs and proteins, a clear trend appears only for re- observed in Fig. 1 and Additional file 1: Figures S1-S17). verse transcriptase inhibitors. In this case, a law of To visualize the impact of Gini index not ascribable to diminishing returns is observed, and adjusted NMSE de- the effects of data size (N) and the kernel used (K), we crease with N is very small for N > ~ 600. plotted the residuals of model NMSE ~ N + K + ε against GINI (Fig. 4 panels a, b and c). For protease inhibitors, the Kernel PCA Gini effect is confined to TPV drug (red dots in Fig. 4a). Even if weighting increases prediction accuracy overall, The effect is rather linear for reverse transcriptase inhibi- the effect was markedly different when we compare re- tors, although NMSE variability was larger than average verse transcriptase and integrase with protease (Table 2). for RPV (red dots), the drug with lowest N. In the case of In the latter protein, weighted kernels were not clearly integrase inhibitors, Gini takes values in a narrow range superior. To further investigate this issue, we performed and does not seem to have an impact on prediction. As in a PCA on the Jaccard kernel. Figure 5 shows the results the case of RPV, large variability in NMSE values is ob- of for FPV (a protease inhibitor, panels a and b) and served in DTG (blue dots), which is the drug with second NVP (a reverse transcriptase inhibitor, panel c and d), lowest sample size. both with unweighted and weighted Jaccard kernels. The Ramon et al. BMC Bioinformatics (2019) 20:410 Page 8 of 13 Fig. 3 a Wild type protease (in yellow and blue) with an inhibitor (NFV, in green) (PDB code: 3EKX). We highlight the ten most important positions according to RF: 10, 90, 54, 46, 71, 88, 84, 30, 20 and 82. These positions are scattered throughout the protein and only a few belong to the drug binding site (e.g. 30, 82 and 84). Mutations at the binding site reduce the affinity for the inhibitor, but can impair the protease catalytic activity as a collateral damage. Mutations in distant residues are typically concurrent with these binding site mutations and often have a compensatory role (e.g. stabilizing the protease structure or restoring the catalytic activity). Position 30 appears to be important only in the case of the NFV drug, while the other positions are found in all (or almost all) protease inhibitors. This agrees with the literature [2]. b Binding pocket of the reverse transcriptase (in yellow) with an NNRTI (NVP, in pink) (PDB code: 3V81). We highlight the five most important positions for NVP according to RF: 103, 181, 190, 188 and 101. All these positions reside in the NNRTI binding pocket of the enzyme, and also appear in the other NNRTIs analyzed. Thus, in EFV, we find 100 (but not 181) in the top 5; and in ETR, we have 179 instead of 188 (also highlighted). Positions 103 and 101 are located near the entry of the inhibitor binding pocket and, when mutated, interfere with the entrance of the inhibitor to the binding site. Y181 and Y188 have a crucial contribution the NVP binding via stacking interactions between its side chains and the inhibitor aromatic groups. G190 mutations lead to resistance through steric hindrance, because of the substitution by a more voluminous side chain. L100 effect is also related to steric hindrance [2] remaining figures can be found at (Additional file 4: Fig- Intriguingly, we found that the stacked versions of SVM ures S36-S54). Unweighted kernel PCA results, overall, with weighted kernels have similar performances to those in a good, spectrum-like separation between resistant of the individual models. This suggests that all the infor- and susceptible isolates for protease inhibitors, whereas mation of the sequence has been already extracted in the weighted kernels can improve dramatically the separ- first step, and so stacking the models was of no additional ation in the case of reverse transcriptase. The integrase value. inhibitors RAL and EVG behave similarly to reverse transcriptase inhibitors, while DTG (which has a very Performance comparison to other approaches small sample size) do not achieve a good separation ei- Figure 6 shows the performance comparison between ther in the weighted or the unweighted kernel PCAs. our best method (SVM with weighted Jaccard kernel) with the ANN1 and ANN2 (see “Performance compari- Stacked models son to other approaches” in Material and methods). We compared the performances of four methods (SVM ANN2 tends to have better performance than ANN1, plus weighted Linear, RBF, Overlap and Jaccard kernels) especially in drugs with small sample size, but also pre- with those of their stacked counterparts in Additional file 3: sents greater standard errors in some drugs. In the case Tables S1 (mean NMSE) and S2 (NMSE standard error). of protease inhibitors (panel a) both ANN1 and ANN2 Table 2 Linear model coefficient estimates and p-values All drugs PIs NRTIs + NNRTIs INIs −4 −5 −4 N increment −2.1·10 *** −2.0·10 *** −1.8·10 *** −1.3·10− 3 *** −2 −3 −2 Unweighted → Weighted −2.0·10 *** 2.0·10 −3.4·10 *** −3.1·10− 2 *** Gini Index increment −4.9·10− 3 *** −6.0·10− 2 *** −5.4·10− 2 *** 1.7·10−2 ** −2 −2 −2 Jaccard → Linear 4.5·10 *** 2.1·10 *** 5.0·10 *** 9.4·10− 2 *** Jaccard → RBF 2.5·10− 2 *** 1.3·10− 2 *** 2.8·10− 2 *** 4.7·10− 2 *** −2 −3 −2 Jaccard → Overlap 1.9·10 *** 9.3·10 *** 2.1·10 *** 3.6·10− 2 *** Legend: Significance codes are 0 *** 0.001 ** 0.01 * 0.05 Ramon et al. BMC Bioinformatics (2019) 20:410 Page 9 of 13 a Protease Inhibitors b Reverse Transcriptase Inhibitors c Integrase Inhibitors ATV FPV LPV SQV ABC D4T EFV NVP TDF Drugs DTG EVG RAL Drugs Drugs DRV IDV NFV TPV AZT DDI ETR RPV X3TC 0.90 0.96 0.94 0.85 0.93 Gini index Gini index Gini index 0.92 0.90 0.80 0.90 0.87 0.75 −0.1 0.0 0.1 0.2 −0.2 0.0 0.2 0.4 0.6 −0.2 0.0 0.2 0.4 NMSE ~ N + K Residuals NMSE ~ N + K Residuals NMSE ~ N + K Residuals d Protease Inhibitors e Reverse Transcriptase Inhibitors f Integrase Inhibitors ATV FPV LPV SQV ABC D4T EFV NVP TDF Drugs DTG EVG RAL Drugs Drugs DRV IDV NFV TPV AZT DDI ETR RPV X3TC 0.50 Model NMSE ~ K + W + GINI Residuals Model NMSE ~ K + W + GINI Residuals Model NMSE ~ K + W + GINI Residuals 0.6 0.1 0.25 0.3 0.00 0.0 0.0 −0.25 −0.1 800 1100 1400 1700 400 700 1000 1300 200 300 400 500 600 N N N Fig. 4 a, b and c NMSE residuals (observed − fitted values) of the linear model containing only data size (N) and kernel (K) vs. Gini index. Each color represents a different drug. Note different scale for the Gini index between panels. d, e and f Residuals (observed − fitted values) of the linear model containing K, W and GINI vs. data size (N). Each color represents a different drug are only marginally worse than the weighted Jaccard which to our knowledge achieved the best performance so SVM, with the exception of the FPV drug. In the case of far in this dataset [14]. We observed that weighted Jaccard reverse transcriptase and the integrase inhibitors (panels outperforms ANN in all drugs, and that the ANN predic- b, c and d), the difference between the performance of tion performances were worse than those originally re- weighted Jaccard and the ANN increases. The latter ported (which had R2 values ranging between of 0.85 and method presents higher NMSE and larger standard er- 0.99). It has to be stressed, however, that we used different rors, especially for 3TC, DDI, TDF, the NNRTIs, and the versions of the dataset (the version used by [14], for in- INIs. stance, did not contain information about the INIs) and that we followed very different strategies concerning pre- Discussion processing. In [14] a pre-processing with removal of out- Recent results on predicting HIV drug resistance as a re- liers and rare variant filtering is performed, which can result gression problem can be found in [14, 18]. Shen et al. in a loss of generalizability, as is acknowledged by the au- [18] used RF and computed the 5-fold cross-validation thors. Another reason for the discrepancy is probably the R2. Sheik Amamuddy et al. [14] used ANN and com- treatment of allele mixtures, as we discuss next. puted the R2 of the test set without replicates. The two In this work, we present a novel approach to predict approaches were based in a previous version of the drug resistance in HIV, using kernel functions that directly Stanford dataset (version date: 2014-9-28) and share a address the presence of allele mixtures and the categorical similar treatment of amino acid mixtures based on se- nature of the data. Previous work handled these two issues quence expansions. We did a comparison with the ANN, using several pre-processing strategies. Categorical data Ramon et al. BMC Bioinformatics (2019) 20:410 Page 10 of 13 a FPV unweighted Jaccard b FPV weighted Jaccard 10 5 Resist 2nd PC − 7.27 % 2nd PC − 10.5 % 2 5 1 0 0 0 −1 −5 −5 −10 −5 0 5 −5 0 5 10 1st PC − 12.7 % 1st PC − 34.6 % c NVP unweighted Jaccard d NVP weighted Jaccard 10 10 5 Resist 2nd PC − 4.31 % 2nd PC − 9.45 % 2 0 0 1 −10 0 −5 −1 −10 −20 −10 0 10 −10 0 10 20 30 1st PC − 10.6 % 1st PC − 27.1 % Fig. 5 The Jaccard kernel PCA in a protease inhibitor (FPV, panels a and b) and a reverse transcriptase inhibitor (NVP, panels c and d). Panels a and c correspond to unweighted Jaccard, and b and d to weighted Jaccard. Dot color represents the actual log-resistance value for each specific drug; in red the more resistant, and in green the least resistant. Sequences with missing resistance value are in gray are systematically recoded into numeric data, usually in widespread phenomenon: about 60% of the sequences the form of dummy data or, in [14], assigning an integer had at least one mixture. Mixtures introduce ambigu- to each category. Here, we have shown that addressing the ity in the genotype-phenotype correlation since it categorical nature of the data and the presence of allele makes impossible to know the actual sequences of mixtures does lower the test error in comparison to the strains. Also, the quasispecies distribution may have dummy variable approach (Table 2). In fact, even the sim- undergone undefined modifications during the in vitro plest categorical kernel (i.e. the Overlap kernel) improves assay [28]. Previous approaches to deal with this issue in- prediction upon the standard RBF kernel, although the ex- cluded keeping the most frequent amino acid of the mix- tent of the improvement depends on the specific drug. It ture [19] and sequence expansion [11, 14, 18]. The latter has to be stressed that recoding the categorical data into strategy consists on expanding the data to sequences with dummy variables increases the dimensionality of the prob- single amino acids at each mixture location until all pos- lem, thus increasing computation needs and leading to sible combinations have been exhausted. These “derived” sparse datasets. As this effect depends on the number of sequences share the resistance value, i.e., the resistance of different categories of the variables, categorical methods the original sequence. This approach dramatically enlarges may be more useful when data has more than few categor- data size (in the aforementioned works, minimum by a ies. Coding the different alleles as an integer does not in- 10x factor in the protease inhibitors and almost a 30x in crease the dimensionality either, but introduces an order the reverse transcriptase inhibitors). This might be one of without biological meaning among the amino acids. the main reasons for the discrepancy between the ANN The treatment of amino acid mixtures is more chal- performance computed in this work and in [14]. Without lenging. In the data analyzed we observed that it is a expansion, the data size ranges between 200 and 1500, but Ramon et al. BMC Bioinformatics (2019) 20:410 Page 11 of 13 a b 1.0 Method 1.0 0.9 ANN1 0.9 ANN2 0.8 0.8 wJAC 0.7 0.7 NMSE NMSE 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 FPV ATV IDV LPV NFV SQV TPV DRV EFV NVP ETR RPV c d 1.0 1.0 0.9 0.9 0.8 0.8 0.7 0.7 NMSE NMSE 0.6 0.6 0.5 0.5 0.4 0.4 0.3 0.3 0.2 0.2 0.1 0.1 0.0 0.0 3TC ABC AZT D4T DDI TDF RAL EVG DTG Fig. 6 Mean NMSE values and their corresponding standard errors for the SVM + weighted Jaccard kernels (red), ANN1 (light gray) and ANN2 (dark gray). PIs are shown in panel a, NRTIs in panel c, NNRTIs in panel b and INIs in panel d the number of (dummy) variables is almost 2000 in the shown that Jaccard is clearly the best among kernels PIs, and more than 4000 in the other drugs. The higher assessed and that also improves the RF results, in most number of variables compared to observations might have cases by a large margin. Both Overlap and Jaccard are adversely affected the ANN performance in comparison basic kernel functions, but our kernel definition (7) is to the original work and, also, in comparison to SVMs, as general enough to replace them for more sophisticated the latter are less prone to over-fitting. Furthermore, the categorical kernels, perhaps with improved prediction expansion potentially biases the dataset by over represent- performance. ing sequences with mixtures (especially those with a larger An additional theoretical proposal was to weigh kernel number of mixtures and/or alleles per mixture) and it can positions according to its inferred influence on drug re- generate HIV variants not found in the patient. Expansion sistance. Here we employed RF decrease in impurity as also increases the difficulty of the training/test splitting weights but numerous options are equally justified and because all expansions of the same sequence must be so additional research on this topic is warranted. Using placed either in the training set or in the test set; other- RF we were able to identify, from protein sequence wise, the independence of both sets is lost. In our work, alone, important positions for the drug resistance that we preferred keeping only one amino acid of the mixture, have a structural meaning (Fig. 3). We observed a dis- which is allegedly the most conservative pre-processing tinct effect of weighting in protease inhibitors and tran- choice. This differs from e.g. [19], because we keep one scriptase reverse inhibitors that correlates with the amino acid at random, while they pick the most frequent distribution of the importances. At least part of this one, which is sound if mixtures are considered a technical behavior might be due to differences in the mutational artifact. However, in case of HIV, this event mostly reflects pattern between the two enzymes in regards to drug re- the coexistence of actual HIV variants in the body of the sistance. In the reverse transcriptase, the major resist- patient [2, 4, 6, 28] and the ambiguity lies in the resistance ance mutations tend to be located in specific positions, value delivered via the in vitro test. In any case, part of the particularly at the drug binding sites of the N-terminal original information is lost by picking one of the allele of side, weakening the affinity between drug and enzyme. the mixture. This does not happen when using the Jaccard As early as 1998, it was noted that a single mutation of kernel, which naturally handles allele mixtures. We have the reverse transcriptase may confer high resistance to Ramon et al. BMC Bioinformatics (2019) 20:410 Page 12 of 13 drugs like 3TC and NVP [28], whereas the virus acquires DRV: Darunavir; DT: Decision Trees; DTG: Dolutegravir; EFV: Efavirenz; resistance to protease inhibitors by accumulating muta- ETR: Etravirine; EVG: Elvitegravir; FPV: Fosamprenavir; HIV: Human immunodeficiency virus; IC50: Half maximal inhibitory concentration; tions. First, primary resistance mutations arise at the ac- IDV: Indinavir; INI: Integrase inhibitor; LPV: Lopinavir; NFV: Nelfinavir; tive site pocket and the surrounding residues. But, as NMSE: Normalized Mean Square Error; NNRTI: Non-nucleoside reverse these mutations often cause conformational changes, transcriptase inhibitors; NRTI: Nucleoside reverse transcriptase inhibitors; NVP: Nevirapine; PCA: Principal Components Analysis; PI: Protease inhibitors; additional secondary mutations that compensate the im- RAL: Raltegravir; RF: Random Forests; RPV: Rilpivirine; SQV: Saquinavir; paired catalytic activity and stabilize the protease tend to SVM: Support Vector Machine; TDF: Tenofovir; TPV: Tipranavir; WHO: World be selected in turn [2]. There are at least 36 important Health Organization residues (out of a total of 99) involved in protease drug Acknowledgements resistance mutations and (unlike reverse transcriptase) A partial version of this work was presented at the 7th International Work- they are distributed along the whole sequence [2]. These Conference IWBBIO 2019 in Granada, Spain (May 8-10, 2019), and is available on the conference proceedings (LNBI Proceedings Part II). differences may explain why RF, and therefore the weighted categorical kernels, performed better at the Authors’ contributions NRTI and NNRTI databases. Further, the estimate of the LBM and MPE conceived and supervised research. ER developed code and variable importance is more reliable when few relevant performed research. ER and MPE wrote the manuscript with help from LBM. All authors read and approved the final manuscript. protein positions have a large impact on resistance. In contrast, the compensatory secondary mutations of the Funding protease probably introduce some degree of correlation This work was funded by project AGL2016–78709-R (Ministry of Economy between protein positions, which may explain why and Competitiveness, Spain) to MPE. ER has funding from a FI-AGAUR PhD studentship grant (Generalitat de Catalunya). We acknowledge financial sup- weighting in PI database does not result in a clear im- port from the Spanish Ministry of Economy and Competitiveness, through provement of performance. the “Severo Ochoa Programme for Centres of Excellence in R&D” 2016–2019 (SEV-2015-0533)”, and from the EU through the BFU2016–77236-P (MINECO/ AEI/FEDER, EU). Conclusions Machine learning is an effective approach to predict Availability of data and materials The datasets analyzed during the current study are available in the HIV drug resistance, and a straightforward alternative to Genotype-Phenotype Stanford HIV Drug Resistance Database repository, the much slower and expensive in vitro assay. Results https://hivdb.stanford.edu/pages/genopheno.dataset.html. Structural data can show that kernels that take into account both the cat- be found at https://www.rcsb.org/structure/3ekx and https://www.rcsb.org/ structure/3v81. Code used in this manuscript is available at https://bitbucket. egorical nature of the data and the presence of mixtures org/elies_ramon/catkern. consistently result in the best prediction model. As for the introduction of position weights, we found that the Ethics approval and consent to participate amount of improvement was a function of the number Not applicable. of positions with large effect on drug resistance, which Consent for publication may be related to the known different mutational pat- Not applicable. terns regarding drug resistance among the viral proteins. Using more sophisticated categorical kernels and/or ker- Competing interests MPE is a member of the editorial board (Associate Editor) of BMC nels able to take into account structural information Bioinformatics. may improve even more the resistance prediction. Author details 1 Centre for Research in Agricultural Genomics (CRAG), CSIC-IRTA-UAB-UB Additional files Consortium, Campus UAB, 08193 Bellaterra, Barcelona, Spain. 2Computer Science Department, Technical University of Catalonia, Carrer de Jordi Girona Additional file 1: Figures S1-S17. NMSE distribution for drugs ATV, 1-3, 08034 Barcelona, Spain. 3Institució Catalana de Recerca i Estudis DRV, IDV, LPV, NFV, TPV, SQV, 3TC, ABC, AZT, D4T, TDF, EFV, ETR, RPV, RAL Avançats (ICREA), Passeig de Lluís Companys 23, 08010 Barcelona, Spain. and DTG (PDF 671 kb) Received: 28 February 2019 Accepted: 11 July 2019 Additional file 2: Figures S18-S35. RF weights for drugs ATV, DRV, IDV, LPV, NFV, TPV, SQV, 3TC, ABC, AZT, D4T, DDI, TDF, EFV, ETR, RPV, DTG and RAL (PDF 530 kb) References Additional file 3: Table S1. Average NMSE of stacked methods. Table 1. UNAIDS. Facts sheet – world AIDS day 2018. UNAIDS. http://www. S2. NMSE Standard error of stacked methods (PDF 42 kb) unaids.org/sites/default/files/media_asset/UNAIDS_FactSheet_en.pdf. Additional file 4: Figures S36-S54. Kernel PCAs for drugs ATV, DRV, Accessed 02 Feb 2019. IDV, LPV, NFV, TPV, SQV, 3TC, ABC, AZT, D4T, DDI, TDF, EFV, ETR, RPV, DTG, 2. Iyidogan P, Anderson KS. Current perspectives on HIV-1 antiretroviral drug EVG and RAL (PDF 2075 kb) resistance. Viruses. 2014;6(10):4095–139. 3. German Advisory Committee Blood (Arbeitskreis Blut). Subgroup ‘assessment of pathogens transmissible by Blood'. Human Abbreviations immunodeficiency virus (HIV). Transfus Med Hemother. 2016;43(3):203–22. 3TC: Lamivudine; ABC: Abacavir; AIDS: Acquired immunodeficiency 4. Shafer RW, Dupnik K, Winters MA. Eshleman SH. A guide to HIV-1 reverse syndrome; ANN: Artificial Neural Networks; ATV: Atazanavir; AZT: Zidovudine; transcriptase and protease sequencing for drug resistance studies. HIV Seq BIC: Bictegravir; CAB: Cabotegravir; D4T: Stavudine; DDI: Didanosine; Compend. 2001;2001:1–51. Ramon et al. BMC Bioinformatics (2019) 20:410 Page 13 of 13 5. Bonet I. Machine learning for prediction of HIV drug resistance: a review. Curr Bioinforma. 2015;10(5):579–85. 6. Schmidt B, et al. Simple algorithm derived from a geno−/phenotypic database to predict HIV-1 protease inhibitor resistance. AIDS (London, England). 2000;14(12):1731–8. 7. Wensing AM, Calvez V, Gunthard HF, Johnson VA, Paredes R, et al. 2017 update of the drug resistance mutations in HIV-1. Top Antivir Med. 2017;24(4):132–3. 8. Beerenwinkel N, et al. Geno2pheno: estimating phenotypic drug resistance from HIV-1 genotypes. Nucleic Acids Res. 2003;31(13):3850–5. 9. Riemenschneider M, Hummel T. Heider D. SHIVA - a web application for drug resistance and tropism testing in HIV. BMC Bioinformatics. 2016;17(1):314. 10. Rhee SY, Taylor J, Wadhera G, Ben-Hur A, Brutlag DL, Shafer RW. Genotypic predictors of human immunodeficiency virus type 1 drug resistance. Proc Natl Acad Sci U S A. 2006;103(46):17355–60. 11. Yu X, Weber IT, Harrison RW. Prediction of HIV drug resistance from genotype with encoded three-dimensional protein structure. BMC Genomics. 2014;15(Suppl 5):S1. 12. Wang D, Larder B. Enhanced prediction of lopinavir resistance from genotype by use of artificial neural networks. J Infect Dis. 2003;188(5):653–60. 13. Pasomsub E, Sukasem C, Sungkanuparph S, Kijsirikul B, Chantratita W. The application of artificial neural networks for phenotypic drug resistance prediction: evaluation and comparison with other interpretation systems. Jpn J Infect Dis. 2010;63(2):87–94. 14. Sheik Amamuddy O, Bishop NT, Tastan Bishop Ö. Improving fold resistance prediction of HIV-1 against protease and reverse transcriptase inhibitors using artificial neural networks. BMC Bioinformatics. 2017;18(1):369. 15. Masso M, Vaisman II. Sequence and structure based models of HIV-1 protease and reverse transcriptase drug resistance. BMC Genomics. 2013;14(Suppl 4):S3. 16. Khalid Z, Sezerman OU. Prediction of HIV drug resistance by combining sequence and structural properties. IEEE/ACM Trans Comput Biol Bioinform. 2018;15(3):966–73. 17. Beerenwinkel, et al. Diversity and complexity of HIV-1 drug resistance: a bioinformatics approach to predicting phenotype from genotype. Proc Natl Acad Sci U S A. 2002;99(12):8271–6. 18. Shen C, Yu X, Harrison RW, Weber IT. Automated prediction of HIV drug resistance from genotype data. BMC Bioinformatics. 2016;17(Suppl 8):278. 19. Tarasova O, Biziukova N, Filimonov D, Poroikov V. A computational approach for the prediction of HIV resistance based on amino acids and nucleotide descriptors. Molecules. 2018;23(11);2751. 20. Heider D, Senge R, Cheng W, Hüllermeier E. Multilabel classification for exploiting cross-resistance information in HIV-1 drug resistance prediction. Bioinformatics. 2013;29(16):1946–52. 21. Riemenschneider M, Senge R, Neumann U, Hüllermeier E, Heider D. Exploiting HIV-1 protease and reverse transcriptase cross-resistance information for improved drug resistance prediction by means of multi- label classification. BioData Min. 2016;9:10. 22. Xing L, Lesperance M, Zhang X. Simultaneous prediction of multiple outcomes using revised stacking algorithms. Version 1. 29 Jan 2019. arXiv:1901.10153v1 [q-bio.QM]. 23. Scholkopf B, Vert JP, Tsuda K. Kernel methods in computational biology. 1st ed. Cambridge: The MIT Press/A Bradford Book; 2004. 24. Genotype-Phenotype Stanford University HIV Drug Resistance Database. https://hivdb.stanford.edu/pages/genopheno.dataset.html. Accessed 30 May 2019. 25. Belanche LA, Villegas MA. Kernel functions for categorical variables with application to problems in the life sciences. Front Artif Intell Appl. 2013;256:171–80. 26. Bouchard M, Jousselme AL. Doré PE. A proof for the positive definiteness of the Jaccard index matrix. Int J Approx Reason. 2013;54(5):615–26. 27. R Core Team. R: a language and environment for statistical computing. Vienna: R Foundation for Statistical Computing. https://www.R-project.org/. Accessed 15 Feb 2019 28. Perrin L, Telenti A. HIV treatment failure: testing for HIV resistance in clinical practice. Science. 1998;280(5371):1871–3. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.