The gene expression data and drug response of 1,001 cell lines screened with 265 drugs in the GDSC dataset (Iorio et al., 2016) were used as training sets, and ten sets of gene expression data and the response of PDX/patients treated with chemotherapies and targeted therapies were used for testing. The test sets consisted of three clinical trial datasets for docetaxel, erlotinib, and sorafenib (Geeleher et al., 2014), four sets from PDX Encyclopedia datasets (Gao et al., 2015), and three sets of TCGA patients (Weinstein et al., 2013; Ding et al., 2016). The sources for both training and testing datasets are detailed in Data Availability Statement. All datasets are publicly available.

Gene expression profiles of cell lines, which were RMA-normalized, log-transformed, and aggregated to the level of genes, were downloaded from the GDSC database. The gene expression profiles of the first three test datasets were preprocessed by Geeleher et al. (2014), and those of the remaining seven test sets (from MOLI) were converted to TPM and log-transformed by Sharifi-Noghabi et al. (2019).

The GED of cell lines from the GDSC dataset was first standardized by the mean and standard error (s.e.) of each gene. Next, the GED of each cell line was normalized by the house-keeping gene GAPDH across cell lines and homogenized with the GED of test sets by the ComBat() function from the sva library in R (Geeleher et al., 2014). The details of each dataset, such as the drug name, the number of samples, and the number of genes in the training and test sets, are provided in Table 1. Similar to existing methods (Geeleher et al., 2014; Sharifi-Noghabi et al., 2019), we included only genes present in both training and test sets for the subsequent analysis; the four columns from the right-hand side of Table 1 show the number of overlapping genes.

For a given drug d, let X d , Y d consist of X d ∈ R n d × p expression profiles of p genes of n d cell lines (PDXs or human tumors) and drug response values Y d ∈ 0 , 1 n d for drug d. Here, we dichotomized the drug response of cell lines (IC₅₀) into 0 (resistant) or 1 (sensitive) if a drug response is greater or less than or equal to its maximum drug concentration (given in the GDSC website), respectively (Iorio et al., 2016).

The proposed feature selection consisted of three procedures: 1) supervised DA (Mourragui et al., 2019), 2) differential expression between sensitive and resistant cells, and 3) the ratio of “between-group to within-group sums of squares” (the BW ratio) (Dudoit et al., 2002), where the two groups refer to sensitive cell lines and resistant cell lines for each drug. As cell lines (the training sets) are different from patients and mouse models in the test sets, we applied DA to sift genes whose conditional distributions given the label Y across domains were not significantly different. The intuition for feature selection procedures 2) and 3) is stated in Supplementary methods. In a pilot study, we also studied Logit (KNN) trained by genes and sifted by unsupervised DA in combination with the latter two proposed feature selection procedures. However, the predictors employing unsupervised DA performed in a manner considerably inferior to those employing supervised DA, given that the remaining procedures were kept the same. Thus, we used supervised DA in this study.

Specifically, for each drug and gene X, the Kolmogorov–Smirnov (KS) test for equality of the conditional distribution of selected genes X given the label Y in the source and target domains, F S X | Y and F T X | Y , was conducted at P ≥ α , where α = 0.6 , 0.7 , 0.8 , and 0.9. It should be noted that DA used all the information about the responses of the samples in both domains. In this study, we used stringent cutoffs P ≥ α and α = 0.6 , 0.7 , 0.8, and 0.9, as the source (cell lines) and target (PDX/patients) domains are quite different. Let X i denote g e n e   i . It should be noted that the aforementioned KS test is equivalent to using the distance measure sup F S X i | Y = 1 − F T X i | Y = 1 and sup F S X i | Y = 0 − F T X i | Y = 0 , where F S X i | Y   a n d   F T X i | Y denote the conditional distribution of X i given Y in the source and target domains, respectively.

Given the features X = X 1 , X 2 , … , X p , we define P S Y = 1 and P T Y = 1 to be the population proportion of responders in the source and target domains, respectively. Furthermore, let P ^ S Y = 1 | X and P ^ T Y = 1 | X be the estimated proportions of responders given X in the source and target domains, respectively, and r = P S Y = 1 / P S Y = 0 / P T Y = 1 / P T Y = 0 be the odds ratio of these two domains. Proposition 1 states that if the selected features satisfy the required DA condition and the odds ratio equals to 1, then the features are invariant across the source and target domains and vice versa.

Proposition 1

assumed that the features X = X 1 , X 2 , … , X p satisfy the DA condition and marginal conditional distributions of X i | Y are independent for i =1, …, p. Then, P S Y = 1 | X = P T Y = 1 | X if and only if r = 1.

The proof is given in Supplementary information.

Genes that passed the DA selection were then prioritized by their differential expression among all overlapping genes in sensitive cell lines versus resistant cell lines. To reveal explainable classifiers, we kept at most the top-ranked 1,000 genes with the smallest false discovery rate (FDR) values obtained from the two-sample t-test and sorted these genes by the BW ratio. The BW ratio for a gene j of the cell line i in group k is defined as follows: B W j = ∑ i ∑ k I y i = k x ¯ k j − x ¯ . j 2 ∑ i ∑ k I y i = k x i j − x ¯ k j 2 , where x ¯ . j and x ¯ k j denote the average expression level of gene j across all cell lines and the cell lines of group k only, respectively.

For fixed top-ranked p genes, where p ranged from 50 with step size 10 to 200 (denoted by 50 (10)200), 200 (20)400, and 400 (100)1,000) genes of the cell lines, we trained the hyperparameter λ (the penalty constant of logit regression) and p using 5-fold CV with ten repeats. We used grid-search to tune the hyper-parameter as follows. First, let λ = 10 a 0 and a 0 ∈ − 3 , 0 with step size 10^0.3, i.e., we ran 5-fold CV of LogitDA with grid points 10^–3, 10^–2.7, …, and 10⁰ and found the grid point whose associated CV score was the maximum, which was termed 10 a 1 , e.g., a 1 = − 2.7 . Second, we further evaluated LogitDA with grid points in 10 a 1 − 0.05 , 10 a 1 + 0.05 and step size 10^0.01. The grid point 10 a max , whose corresponding CV score is the maximum, determines the tuned hyperparameter λ = 10 a max . The logit model with the highest averaged CV AUC determined p and λ , which yielded one LogitDA. Then, we applied the LogitDA to a test set to determine its prediction AUC. The aforementioned procedures were repeated ten times (with different seeds for CV) to obtain the mean and s.e. of the prediction AUC.

For the classifier KNN, we used the distance measure 1 − rho, where rho is Spearman’s rho between any two cell lines with selected p genes because the default Euclidean distance did not work well for the GED of any two cell lines in our pilot study. For each drug and fixed top-ranked p gene, where p = 50 (10)200, 200 (20) 400, and 400 (100)1,000 genes of the cell lines, we first trained the hyperparameter K of KNN, via 5-fold CV with ten repeats and using the GED and drug response of cell lines, and computed the AUC in the cross-validation experiments. The hyperparameter K was determined by the experiment with the highest averaged CV score. We then fitted all data into this KNN classifier with each top-ranked p gene. Of all the top-p-ranked KNN classifiers trained, the one with the highest averaged CV score determined the value of p, which was one trained KNNDA predictor. We repeated the aforementioned procedures ten times to yield the mean and s.e. of the prediction AUC of KNNDA.

The following lemma and Proposition 2 established the theoretical foundation for adjusting the prediction probability cutoff when the drug response rates between cell lines and tumors differ.

Lemma. When r > 1 and P S Y = 1 | X > P T Y = 1 | X , the prediction probability is overestimated. Similarly, when r < 1 and P S Y = 1 | X < P T Y = 1 | X , the prediction probability is underestimated.

The proof is given in Supplementary information.

Proposition 2

assumed that the predictors X = X 1 , X 2 , … , X p satisfy the DA condition and marginal conditional distributions of X i | Y are independent for i =1, …, p. When the odds ratio between the source and target domains is r ≠ 1 , the cutoff of the prediction probability P S Y = 1 | X T should be adjusted to r / r + 1 .

The proof is given in Supplementary information.

As Proposition 2 suggests, the prediction probability cutoff will deviate from 0.5 when the population proportion of responders in the target and source domains differ. Due to a lack of information on the ratio of responders in the target domain (patients), we estimated it using test data directly in this study. However, a better estimate of the ratio can be obtained once more responses to these drugs are released. We note that Proposition 2 considers a continuous prediction probability function. It does not apply to the KNN classifier that makes a prediction on a test sample based on the majority voting for its K-nearest neighbors.

Finally, we applied the trained predictors LogitDA with α = 0.7 and KNNDA with α = 0.7 to the ten external test sets. To compare with the results of the baseline ridge regression (Geeleher et al., 2014) and MOLI complete (Sharifi-Noghabi et al., 2019), we repeated the experiments ten times for estimating the s.e. of the predictors for each drug.

In this study, we aimed to investigate the following questions: Do logistic ridge regression and KNN with adequately selected features outperform a deep learning-based predictor, MOLI complete (Sharifi-Noghabi et al., 2019), in terms of prediction AUC on external test sets (PDX and patient data)? Do the proposed predictors, LogitDA and KNNDA, work well for targeted therapies and/or chemotherapies? Information about the training and test datasets of the ten studied drugs is provided in Table 1.

After the features (namely, genes in this study) were selected by the proposed procedures (see Methods for details), we trained logistic ridge regression (KNN) with 5-fold CV using the GED of the prioritized features of GDSC cell lines screened with seven drugs, which included docetaxel, erlotinib, sorafenib, cetuximab, gemcitabine, paclitaxel, and cisplatin, in a total of ten sets. These drugs were chosen because we planned to compare LogitDA and KNNDA to the baseline logistic ridge regression (Geeleher et al., 2014) and MOLI complete (Sharifi-Noghabi et al., 2019).

As the training set (GDSC cell lines) is quite different from the test sets (PDX and patient data) (Sharifi-Noghabi et al., 2019), the cutoff for domain adaptation should be strict. Nevertheless, this threshold should allow sufficient features to pass so that a classifier can be adequately trained; see Methods for details. Therefore, we trained LogitDA with features ( X i ) that passed the KS test for equality of F S X i | Y and F T X i | Y with P ≥ α , where α = 0.6, 0.7, 0.8, and 0.9, and we denote the resulting predictor as LogitDA_{_α}, where α = 0.6, 0.7,0.8, and 0.9.

Supplementary Table S1 shows that the CV scores of LogitDA_{_0.6} were equivalent to those of LogitDA_{_0.7} for the ten drugs. However, to satisfy the DA condition required by Proposition 1, namely, the marginal conditional distribution of selected genes X given the label Y in both domains is equal, the value of α should be large, so we chose LogitDA_{_0.7}. Proposition 1 shows that features that satisfy the conditions will be domain-invariant. That is, if the features perform well in the source domain, they will also perform well in the target domain.

Moreover, we observed that each averaged CV score of LogitDA_{_0.7} was higher than that of LogitDA_{_0.8} and LogitDA_{_0.9}, except that they had the same CV score for cetuximab (PDX). Thus, among the LogitDA predictors, we suggest using LogitDA_{_0.7} for the prediction of the drug response of patients/PDX and denote it by LogitDA, henceforth, for simplicity; for details, we refer to Table 2.

Taking the training result of LigitDA into account, we trained the non-linear KNNDA with α = 0.7, 0.8, and 0.9 for the KS test and summarized the 5-fold CV result in Table 3. KNNDA_{_0.7} performed better than KNNDA_{_0.8}, as the former had higher (lower) averaged CV scores for five (two) drugs than the latter; the differences ranged from 1% to 4%. Moreover, KNND_{_0.8} outperformed KNNDA_{_0.9} in terms of higher averaged CV scores for nine of the ten drugs. Thus, we suggest using KNNDA_{_0.7} among these non-linear predictors for the test sets. For simplicity, we denote KNNDA_{_0.7} by KNNDA henceforth.

Next, Table 4 and Figure 2 report the prediction AUC of LogitDA and KNNDA for the ten test sets. The predictor LogitDA achieved a prediction AUC >0.8 for five drugs and predicted AUCs of 0.71 and 0.70 for docetaxel and sorafenib, respectively. In particular, LogitDA using the top-ranked 50, 130, 50, 100, and 650 genes resulted in prediction AUCs of 0.94, 0.93, 1.00, 0.83, and 0.81 for erlotinib, cetuximab (PDX), erlotinib (PDX), gemcitabine (PDX), and docetaxel (TCGA), respectively. This result shows that LogitDA may be useful for precision oncology, especially for the targeted therapies erlotinib and cetuximab.

Of the ten drugs, the predictor KNNDA achieved a prediction AUC >0.8 for four drugs. Specifically, KNNDA using the top-ranked 110, 220, 110, and 60 genes resulted in prediction AUCs of 0.87, 0.90, 0.95, and 1.00 for docetaxel, erlotinib, cetuximab (PDX), and erlotinib (PDX), respectively. This result shows that KNNDA may also be useful for precision oncology.

We further compared LogitDA to KNNDA. Of the ten drugs, LogitDA had a significantly higher (21% higher) (11% lower) prediction AUC compared to KNNDA for gemcitabine (PDX) (docetaxel) and performed equivalent to KNNDA for the remaining eight drugs. Thus, these predictors performed equivalently; we refer to Table 4 for details.

As shown in Table 4 and Figure 2, the predictors LogitDA and KNNDA performed significantly better (16%–35% and 13%–37% higher prediction AUC) than the baseline logistic ridge regression model (Geeleher et al., 2014) for nine and ten out of the ten drugs, respectively.

Next, we compared LogitDA and KNNDA to the deep neural network (DNN)-based method MOLI (Sharifi-Noghabi et al., 2019), which outperformed DNNs with early integration, with 5-fold CV and 10 repeats. Of the ten drugs in Table 4, LogitDA and KNNDA outperformed MOLI complete (expression data) for seven and eight drugs, respectively. In particular, LogitDA and KNNDA had 31%–61% and 44%–61% higher prediction AUCs compared to MOLI complete for docetaxel, cetuximab (PDX), erlotinib (PDX), and gemcitabine (PDX). Furthermore, LogitDA and KNNDA also had significantly higher (18%–21% and 14%–17%) prediction AUCs compared to MOLI complete for erlotinib and docetaxel (TCGA). LogitDA and KNNDA only performed significantly worse than MOLI complete (13% and 8% lower prediction AUC) for cisplatin (TCGA).

The prediction AUC of LogitDA for both cisplatin and gemcitabine (TCGA) was only 62%, which may be because the ratio of sensitive (responders) versus resistant (non-responders) samples is reversed from the training to the test sets (from about 1:2 to 10:1 for cisplatin); in other words, our assumption that the ratio of sensitive to resistant samples in both domains is equal was not met.

EGFR expression has been used as a biomarker to treat colorectal cancer (CRC) patients with wild-type KRAS in the US (patients with metastatic CRC and HNSCC in the EU). However, EGFR expression does not predict a response to cetuximab (Messersmith and Ahnen, 2008). The high prediction AUC of LogitDA for cetuximab (PDX) suggests that the fitted 130 genes may be promising for selecting KRAS wild-type patients with CRC for cetuximab, provided more test sets are validated.

Furthermore, we compared LogitDA (KNNDA) to MOLI complete (multi-omics data). LogitDA (KNNDA) has a significantly higher prediction AUC 19%–40% (19%–42%) in the test set of cetuximab (PDX), erlotinib (PDX), gemcitabine (PDX, LogitDA only), and docetaxel (TCGA) and performed equivalent (<6% differences of test AUC) to MOLI complete (multi-omics data) for the remaining drugs, except that KNNDA had 9% less test AUC for paclitaxel; please refer to Table 4 for details.

For targeted therapies such as erlotinib and cetuximab, Sharifi-Noghabi and colleagues further trained MOLI on multi-omics data of five drugs targeting the EGFR pathway (MOLI complete (pan-drug)), which consisted of >3,000 samples. It is interesting that LogitDA (KNNDA) (using merely hundreds of samples) outperformed MOLI (pan-drug) for erlotinib (PDX) and cetuximab (PDX) with 28% and 13% (28% and 15%) higher prediction AUCs, respectively.

In addition to the high prediction AUC for the aforementioned drugs, our approach also has the advantages of being interpretable and using much fewer (50–650) genes that are interpretable in comparison with the baseline logistic ridge regression and MOLI, which used more than 12,000 genes, except for docetaxel, for which ∼6,370 were used, as shown in Table 1. The use of much fewer genes (features) and hyperparameters may prevent LogitDA and KNNDA from overfitting problems.

For some of the ten drugs whose odds ratios of the source and target domains r = P S Y = 1 / P S Y = 0 / P T Y = 1 / P T Y = 0 deviate much from 1 (Supplementary Table S2A), Proposition 2 shows that their cutoff of the predicted probability should be adjusted to r/(r+1) to account for the differences of the ratios across domains. Therefore, we adjusted the cutoffs accordingly and obtained the prediction accuracy of the ten drugs in Table 5. Notably, for seven of the ten drugs, the resulting prediction accuracy is greater than or equal to 0.70. In particular, for 25 tumors treated with erlotinib, LogitDA achieved a prediction accuracy of 0.76, and its prediction accuracy increased to 0.85 if we focused on the 20 EGFR and KRAS wild-type patients with NSCLC; LogitDA correctly predicted all 12 resistant tumors and five of eight tumors sensitive to erlotinib(Supplementary Table S2B). To the best of our knowledge, to date, there is no efficient biomarker to predict the response to erlotinib of such patients (Geeleher et al., 2014), who were estimated to represent ∼30% of Caucasian patients with lung adenocarcinoma (Wang M. et al., 2021).

In a pilot study, we found that DA considerably improved the prediction AUC of logistic ridge regression (KNN) combined with feature selection using DE genes and BW ratio. Thus, it was of interest to evaluate the contribution of DA. We trained logistic ridge regression and KNN with the features selected by the two aforementioned feature selections (denoted as LogitDA-DA and KNNDA-DA, respectively) for the ten drugs. Supplementary Table S3 shows that LogitDA-DA (KNNDA-DA) used a few hundred genes to achieve equivalent test AUCs as the baseline logit model trained by more than 5,100–13,400 features for the ten drugs (Figure 2). Table 6 shows that DA increases the averaged prediction AUC of LogitDA (KNNDA) from 0.55 to 0.79 (0.57–0.78) over LogitDA-DA (KNNDA-DA), where the averaged prediction AUC was averaged over the ten drugs; the improvements are quite significant.

As LogitDA and KNNDA perform well in the prediction AUC for erlotinib and cetuximab (PDX), it is of interest to find the pathways in which the fitted genes of these predictors are involved. Thus, we first submitted the top-ranked 220 genes of LogitDA and KNNDA for erlotinib into the database Ingenuity Pathway Analysis (IPA; http://www.ingenuity.com) and uncovered the relevant pathways in Supplementary Table S4. Interestingly, several important metabolic pathways were discovered, e.g., Purine Nucleotides de novo Biosynthesis (Ali et al., 2020; Taha-Mehlitz et al., 2021) and Histidine Degradation VI (Tominaga et al., 2019; Han et al., 2020). Furthermore, pathways for epigenetic regulation (Yu et al., 2018; Du et al., 2019) and DNA repair (nucleotide excision repair enhanced pathway) (Dong et al., 2019; Wang T. et al., 2021; Sato et al., 2021) were also uncovered. The aforementioned pathways play essential roles in tumor malignancy and response to anti-cancer therapies.

The overlap of the aforementioned fitted genes and the uncovered pathways (in the molecules column of Supplementary Table S5) has been reported to contribute to tumor progression (cell proliferation, survival, invasion, and metastasis) and drug resistance. Specifically, LIG1 is an attractive target for personalization of ovarian cancer therapy (Ali et al., 2021), and decreased eEF2 phosphorylation, mediated by increased PP2A activity, contributes to resistance to HER2 inhibition (McDermott et al., 2014). ADSL has been suggested as a predictive biomarker of response to 6-mercaptopurine (under the brand name Purinethol) in a pre-clinical setting (Taha-Mehlitz et al., 2021).

Moreover, MTA3 downregulates SOX2OT, and the MTA3/SOX2-OT/SOX2 axis has been reported as a potential cancer stratification marker in human esophageal squamous cell carcinomas (Du et al., 2019). Finally, CAD, a key enzyme of de novo pyrimidine biosynthesis essential for cell proliferation, has been found to directly interact with the second generation of EGFR-TKI Afatinib, which also targets EGFR in the same pathway as erlotinib (Tu et al., 2021).

Similarly, we submitted the top-ranked 130 genes of LogitDA for cetuximab (PDX) and uncovered DNA repair, metabolic processes, and lysosome-associated pathways; the details of the pathways are listed in Supplementary Table S5. The overlap of the fitted genes and the uncovered pathways includes CDK7, IGF-1R, and others. In particular, CDK7 is a key regulator of transcription and cell-cycle control, and its deregulation in cancer has been linked to a worse prognosis (Jagomast et al., 2022). Inhibition of CDK7/12 promotes resistance emergence in response to targeted therapy in lung cancer cells (Rusan et al., 2018; Terai et al., 2018). Moreover, cetuximab therapeutically blocks EGFR, and this might concurrently induce the activation of IGF-1R, which could activate EGFR-downstream Akt signaling, thus mediating cetuximab resistance in gastric cancer cells (Li et al., 2015).