pathCHEMO, a generalizable computational framework uncovers molecular pathways of chemoresistance in lung adenocarcinoma

Abstract

Despite recent advances in discovering a wide array of novel chemotherapy agents, identification of patients with poor and favorable chemotherapy response prior to treatment administration remains a major challenge in clinical oncology. To tackle this challenge, we present a generalizable genome-wide computational framework pathCHEMO that uncovers interplay between transcriptomic and epigenomic mechanisms altered in biological pathways that govern chemotherapy response in cancer patients. Our approach is tested on patients with lung adenocarcinoma who received adjuvant standard-of-care doublet chemotherapy (i.e., carboplatin-paclitaxel), identifying seven molecular pathway markers of primary treatment response and demonstrating their ability to predict patients at risk of carboplatin-paclitaxel resistance in an independent patient cohort (log-rank p-value = 0.008, HR = 10). Furthermore, we extend our method to additional chemotherapy-regimens and cancer types to demonstrate its accuracy and generalizability. We propose that our model can be utilized to prioritize patients for specific chemotherapy-regimens as a part of treatment planning.

Keywords: Computational models; Lung cancer; Predictive medicine.

Fig. 1

Schematic representation of pathway altered on both transcriptomic and epigenomic levels. Pathway genes affected on transcriptomic and epigenomic levels in G alpha (s) signalling events pathway are represented by ovals, where their colors correspond to either over-expression (red), under-expression (blue) or no differential expression (white). Small satellite circles represent over-methylation (red) or under-methylation (blue)

Fig. 2

Integrative genome-wide transcriptomic and epigenomic analysis identifies candidate molecular pathways of chemotherapy response. a Schematic representation of the integrative transcriptomic and epigenomic analysis: first, patients are defined by their response to chemotherapy; second, our method integrates patients’ transcriptomic and epigenomic profiles; third, candidate pathways affected on both transcriptomic and epigenomic levels are identified; and finally, our method employs multi-modal validation of candidate pathways. b Box and whisker plot depicting p-value cutoff for query carboplatin–paclitaxel response composite methylation pathway signature (x-axis) and NESs from the corresponding GSEA comparison between composite methylation and expression pathways signatures (y-axis), based on analysis in TCGA-LUAD patient cohort. Arrow indicated optimal p-value threshold, which results in the strongest GSEA enrichment. c GSEA comparing carboplatin–paclitaxel response composite expression pathway signature (reference) and carboplatin–paclitaxel response composite methylation pathway signature (query, NES p ≤ 0.001), based on analysis in TCGA-LUAD patient cohort. Horizontal red bar indicates leading edge pathways altered on both transcriptomic and epigenomic levels. NES and p-value were estimated using 1000 pathway permutations. d ROC analysis comparing ability of the 7 candidate pathways to predict carboplatin–paclitaxel where their activity is defined based on their expression values (green) or methylation values (blue). AUROC is indicated

Fig. 3

Transcriptomic and epigenomic alterations in candidate pathways of carboplatin–paclitaxel response. a Representative molecular pathways altered on both transcriptomic and epigenomic levels, visualized through circlize R package. Genes from the leading edge in each pathway are represented as differentially expressed (pink), methylated (gray) and both differentially expressed and methylated (yellow). Width of each connecting line is proportional to the extent of differential expression and differential methylation. Pathways are depicting as follows: chemokine receptors bind chemokines pathway (19 differentially expressed genes, four differentially methylated genes, and eight differentially expressed and methylated genes); mRNA splicing pathway (21 differentially expressed genes, 39 differentially methylated genes, and 28 differentially expressed and methylated genes); and G alpha (s) signaling events pathway (37 differentially expressed genes, eight differentially methylated genes, and four differentially expressed and methylated genes). b In the seven candidate pathway network representation, nodes correspond to the genes, which are connected to central pathway-membership circles (i.e., indicating pathway membership). Gene colors describe differential expression (pink), differential methylation (grey) and both differential expression and methylation (yellow). Network was constructed with ggnetwork R package

Fig. 4

Candidate molecular pathways stratify patients based on response to carboplatin–taxane in an independent cohort. a Validation strategy, as follows: first, employment of molecular transcriptomic and epigenomic profiling of patients; second, predicting patients’ risk of developing chemoresistance; and finally, informed clinical decision making based on patients personalized risks. b t-SNE clustering of lung adenocarcinoma patients treated with carboplatin–taxane (e.g., paclitaxel) from the Tang et al. validation cohort (n = 39 biologically independent patient samples), based on activity levels of seven candidate pathways. Among two groups green group (n = 21 biologically independent patient samples) corresponds to patients with low composite activity levels of candidate pathways and orange group (n = 18 biologically independent patient samples) corresponds to patients with high composite activity levels of candidate pathways. c Kaplan–Meier survival analysis to estimate difference in response to carboplatin–taxane (e.g., paclitaxel) between two patient groups is identified in b. Log-rank p-value and number of patients in each group are indicated. d Two random models indicate non-random predictive ability of our model in the Tang et al. validation cohort: random model 1 (steel-blue) is defined based on to seven pathways selected at random, and random model 2 (goldenrod) is defined based on to equally sized patient groups selected at random

Fig. 5

Comparative performance analysis confirms robust predictive ability of pathCHEMO. a, b Comparison of pathCHEMO (turquoise) to other commonly utilized methods, including Panja et al. Epi2GenR (yellow), Zhong et al. SVM (light blue), Yu et al. PRES random forest (dark blue) using (a) ROC analysis (with AUROC indicated) and (b) Kaplan–Meier and Cox proportional hazards model (with log-rank p-value and hazard ratio indicated) in Tang et al. validation cohort. c Multivariable Cox proportional hazards analysis demonstrating adjustment of seven candidate pathways for common covariates (i.e., age, gender and stage at diagnosis). Hazard p-value is indicated. d Multivariable Cox proportional hazards analysis demonstrating adjustment of seven candidate pathways for signatures of lung cancer aggressiveness, including Larsen et al. (54 lung adenocarcinoma markers), Beer et al. (50 lung adenocarcinoma markers), and Tang et al. (12 non-small cell lung cancer markers). Hazard p-value is indicated

Fig. 6

pathCHEMO accurately identifies pathways of treatment resistance across chemo-regimens and cancer types. Treatment related Kaplan–Meier survival analysis in a cisplatin–vinorelbine treated lung adenocarcinoma (LUAD) patients in the Zhu et al. patient cohort (n = 39 biologically independent patient samples), b cisplatin–vinorelbine treated lung squamous cell carcinoma (LUSC) patients in the Zhu et al. patient cohort (n = 26 biologically independent patient samples), and c FOLFOX (folinic acid, fluorouracil, and oxaliplatin) treated colorectal adenocarcinoma (COAD) patients in the Marisa et al. patient cohort (n = 23 biologically independent patient samples), demonstrating ability of identified candidate pathways (for each analysis) to predict treatment response. Log rank p-value and number of patients in each group are indicated