Establishment and molecular characteristics of lung cancer PDCs

We established a PDC collection of 102 samples (National Cancer Center; see Additional file 2, Table S1) from patients with advanced lung cancer enrolled in this study (Fig. 1A). PDCs were primarily collected from pleural effusions (92.2%), and secondarily from pericardial effusions (4.9%), ascites (2.0%), and tissues (1.0%; Additional file 2, Table S1). The pathologic type was adenocarcinoma (ADC; 84.3%), SCLC (5.9%), and miscellaneous types (squamous cell carcinoma, 4.9%; sarcomatoid carcinoma, 3.9%; and not otherwise specified, 1%). The drug panel used for response screening included 48 anti-cancer compounds of seven classes targeting angiogenesis (n = 3), the cell cycle (n = 8), DNA damage (n = 6), MAPK (n = 3), PI3K/AKT/mTOR (n = 3), protein tyrosine kinase (n = 7), and others (n = 18; see Additional file 2, Table S2). Among the 48 drugs, 16 were screened in all PDCs and 32 were screened in 38 PDCs.

Fig. 1

Genomic profiles of refractory lung cancer patient-derived cells (PDCs) for pharmacogenomic analysis. A Graphical workflow of pharmacogenomic analysis of lung cancer PDCs. B Gene expression and mutation profile comparison of our PDCs with TCGA samples. The scatterplot presents the principal component analysis for TCGA tumors, normal tissue samples, and PDCs. The bar plot indicates the mutation frequency of TCGA samples and PDCs. TCGA-T (orange) indicates tumor samples and TCGA-N (blue) indicates the adjacent normal tissue. C Genomic landscape heatmap of PDCs to harbor top-ranked somatic variants

To investigate the genomic characteristics, we identified somatic mutations and CNVs using target-seq (n = 98) and classified four molecular subtypes (C1–C4) from gene expression profiles (n = 102; Additional file 2, Table S1). Next, we checked whether the PDCs effectively recapitulate the mutations and expression profiles of the solid tumors. Owing to the dominance of LUAD in our cohort (> 80%), we compared PDC gene expression profiles with TCGA-LUAD dataset (n = 230). As expected, the PDC gene expression profile resembled with tumor samples and was separated from normal adjacent tissues (Fig. 1B) [10]. Constitutive somatic gene mutations were similar in PDCs and TCGA samples. The recurrence of TP53, RB1, and BRAF mutations was highly preserved in both PDC and TCGA samples. The EGFR mutation frequency was higher in PDCs, whereas the recurrence of KRAS, KEAP1, and STK11 mutations was lower than that of TCGA samples not shown in Fig. 1B. Thus, somatic mutations in TP53 (47%), EGFR (29%), and RB1 (8%) were frequently observed in PDC models (Fig. 1C). Moreover, MET (10%), CDK4 (6%), and MDM2 (6%) variants, as well as EML4-ALK (4%) and CD74-ROS (2%) fusion genes, were detected.

Before in-depth analysis, the most frequent TP53 and EGFR mutations were additionally categorized according to selectivity or functionality [22, 47]. TP53 mutations were divided into hotspot (known as gain-of-function; 5.9%; R175, G245, R248, R249, R273, R282) and non-hotspot (unknown or loss-of-function; 41.2%) mutations. EGFR mutations were categorized into (1) TKI-targetable single variant (TARGET: L858R and exon 19 deletion = 19.6%), (2) T790M acquisition (T790Maq: TARGET and T790M = 2.9%), (3) non-other-specified (NOS) single variant except TARGET (4.9%), and (4) NOS acquisition (NOSaq; TARGET and NOS = 2%; Table 1) [22, 47].

Table 1 Clinical characteristics of four RNA-subtypesMolecular subtype classification and targeted drug candidate identification

PDC samples could be classified into four molecular subtypes using gene expression profile by NMF clustering. To extrapolate the clinical characteristics for each subtype, we interrogated patient clinical profile encompassing histologic type, survival, smoking, and EGFR-TKI therapy record (Fig. 2 and Table 1). To uncover regulatory program for each molecular subtype, variant enrichment, and pathway regulation scores were assessed from multi-omics profile. In brief, subtype C1 was associated with a good outcome shown in Fig. 2A, and was activated in inflammatory and IL6-JAK-STAT3 signaling pathways; subtype C2 was associated with a modest outcome, dominance of TP53/EGFR wild-type, upregulation of epithelial-to-mesenchymal transition (EMT), and enrichment of osimertinib-resistant group (POST3) PDCs; subtype C3 was associated with the worst outcome, long-term smoking males, SCLC, TP53 hotspot mutation, MYC activation, and fasten G2M checkpoint; by contrast, subtype C4 exhibited the best survival, frequent TP53 non-hotspot mutation, the dominance of EGFR-TKI TARGET mutation, and NOTCH signaling activation (Fig. 2A, B). We also demonstrated the survival significance corresponding to C1–C4 subtype gene sets using additional transcriptome datasets (n = 1587; see Fig. S1A-B in Additional file 3 and Table S3 in Additional file 2). Upregulation of the C3 gene set concurrently exhibited the worst outcomes (P < 0.001 and hazard ratio [HR] = 2.6). When additionally demonstrating from known up-regulated genes of embryonic stem cell, we could observe the activation of human embryonic stem cell genes, and target genes’ upregulation of transcription factor MYC and SOX2 (Fig. S2 in Additional file 3). The activation of C1 and C2 subtype genes significantly showed better prognosis (P < 0.05 and HR < 0.79). Our subtype classification sustained a global prognostic signature for lung cancer. Finally, we could summarize the subtypes based on the following regulatory pathways: C1, inflammatory; C2, EMT-like; C3, stemness; and C4, EGFR-dominant.

Fig. 2

PDC molecular subtypes and their drug responses. A Overall survival curve plot of the four molecular subtypes. P-values were acquired using the log-rank test. Hazard ratio (HRs) 95% confidence intervals (CIs) were determined using the Cox model. B Heatmap indicating average scores of activated pathways for each subtype identified using the limma test (P < 0.001). C Volcano plot to test drug sensitivity of all possible pairs of subtypes and drugs. The x-axis is the log-scale fold change (log2 FC) to subtype’s average AUC divided by the rest samples’ average AUC. The y-axis is the P-value performed by Wilcoxon rank-sum test. Drug–subtype pairs are labeled. D Heatmap of drug–subtype pairs determined from response sensitivity tests for each subtype (P ≤ 0.05 and log2 FC ≤ 0.2). Only 30 drugs to pass P-value cutoff for each subtype were denoted. The size of each circle showed log-scale P-values, and circle was colorized by log2 FC

Sensitive drug candidates exhibited remarkable concordance with the previously identified regulatory pathways for each molecular subtype (P < 0.05 and |log2 fold change (FC)| < 0.2; FC was assessed to compare average drug AUC values between corresponding subtype and another group; Fig. 2C). The C1 inflammatory subtype was sensitive to only ruxolitinib (JAK1/2) that targets the JAT/STAT pathway (Fig. 2B-D). The C2 EMT-like subtype showed MAPK class drug resistance and sensitivity to dasatinib (angiogenesis and SRC inhibitor), cabozantinib (VEGFA inhibitor), miscellaneous class repotrectinib (ROS inhibitor), and XAV939 (WNT-TNKS-β-catenin inhibitor). The C3 stemness subtype was sensitive to the top five-ranked cell cycle inhibitors. The C4 EGFR-dominant subtype exhibited the strongest resistance to most drug classes except MAPK inhibitors (selumetinib and trametinib) and EGFR-TKIs (gefitinib and afatinib). Notably, targets of predicted drugs belonged to pathways that were found to be activated in each subtype (Fig. 2B).

Predicting drugs for variants and dissecting the molecular subtype of EGFR-mutated PDCs

To interrogate drug candidates for variants, we tested the difference in drug responses of single and co-occurrent mutation cases (P < 0.005 and |log2 FC| < 0.2; Fig. 3A). Co-mutated cases were also explored using the Fisher’s exact test (P < 0.25; see Additional file 3, Fig. S3). Unexpectedly, EGFR-TARGET mutated PDCs showed a relatively modest response to three EGFR-TKIs (afatinib P = 0.17, gefitinib P = 0.19, osimertinib P = 0.08). To uncover EGFR-mutated cells’ molecular features interfering with the EGFR-TKI response, we dissected TKI TARGET mutations according to our four molecular subtypes (Fig. 3B). The C4 EGFR-dominant subtype was the most sensitive to all EGFR-TKIs, and the C1 inflammatory subtype also showed an especially good response to afatinib and osimertinib. Mutated cases (n = 2) of the C3 stemness subtype was insufficient for statistical test. Finally, mutated cases (n = 6) of the C2 EMT-like subtype did not respond to any EGFR-TKIs. Additionally, both T790Maq and NOSaq (EGFR R776G and I744M with L858R) PDCs were TKI-sensitive. In the EGFR NOS type mutations, G719A was observed, and it comprised 11.5% among the NOS mutations in previous NSCLC study, and patient-derived xenografts demonstrated that the mutation was resistant to osimertinib [48]. The remaining NOS mutations excluded exon 18–21 and had a low possibility of finding a structure-based therapeutic target [48]. Therefore, we concluded that the remaining NOS group mutations were not targetable by EGFR-TKIs. Especially, EGFR TARGET mutated PDCs classified to EMT-like subtype exhibited low response to EGFR-TKIs. Thus, our molecular subtype of EGFR mutations revealed that these PDC models can be used to verify the heterogeneous tumor environment affecting drug responses.

Fig. 3

Drug candidates according to mutated cases. A Volcano plot to test drug sensitivity for single or double mutations (P ≤ 0.005 and log2 FC ≤ 0.2). B Heatmap to summarize the three EGFR-TKIs’ sensitivity of EGFR-target (L858R and exon 19 deletion; n = 20) mutated PDCs classified by four molecular subtypes. P-values (rectangle size) were acquired by Wilcoxon rank-sum test, and fold changes (FCs; color-coded) were calculated to compare target mutation group with the wild-type. The number of mutation cases for each subtype is denoted on the y-axis. C Waterfall plot of the alpelisib response. RB1, TP53, and RB1/TP53 mutated PDCs are denoted. SCLC (S) and molecular subtype are labeled. The drug response comparison between mutation and wild-type investigated in both PDCs and cell lines, summarized in a table. D Heatmap of the alpelisib gene signature (n = 12) extracted from machine-learning and GSEA (pathway “MITOTIC G1 PHASE AND G1 S TRANSITION”, P < 0.1). PDCs were divided into drug-sensitive and -resistant groups, and the normalized AUC values are presented as a dot plot

When interrogating drug candidates for mutations, BRAF variants were sensitive to trametinib (Fig. 2A). The TP53 non-hotspot exhibited resistance to olaparib (PARP inhibitor), repotrectinib, and dasatinib. RB1 variants responded well to alpelisib. Although TP53 mutation was not associated with any drug response (P < 0.28), RB1/TP53 co-mutated samples exhibited certain sensitivity that was stronger than RB1 (Fig. 3A). RB1/TP53 mutations were mostly enriched in SCLC (n = 4; P < 0.001) but were also observed in NSCLC (n = 3; Additional file 3, Fig. S3). We additionally investigated alpelisib sensitivity comparing PDCs with cell lines (n = 70; SCLC n = 7). Notably, although alpelisib inhibited SCLC (Fig. 3C), RB1/TP53 variants (P < 0.001) were more strongly inhibited than SCLC (P = 0.02). However, in contrast to RB1/TP53 (P < 0.01), we failed to identify SCLC sensitivity to alpelisib (P = 0.47; Fig. 3C) using PDCs. Meanwhile, another PIK3-AKT target drug capivasertib (AKT1) was sensitive in RB1/EGFR-NOS and EGFR-NOS cases. When dissecting details, two of the five (RB1/EGFR-NOS) PDCs was classified to SCLC type additionally harboring TP53 mutations (Fig. S3B in Additional file 3) [48]. Therefore, we could infer that RB1/EGFR-NOS PDC cases accompanied the dependency to SCLC type. Therefore, PI3K-AKT class drugs detected from SCLC or RB1/TP53 cells also affect EGFR-NOS mutations. Collectively, mutation-based drug detection result emphasize that the genomic characteristics could predict both the cancer drug response and histological lung cancer type SCLC (Additional file 2, Table S1). Furthermore, our PDC model predicted alpelisib as an appropriate treatment for SCLCs, which was not identified in cell lines (Fig. 3C).

To delineate the drug susceptibility of alpelisib from transcriptome, we extracted the transcriptome signature for each drug using machine-learning approaches. The response to alpelisib was modulated by a transcriptome signature enriched “MITOTIC G1 PHASE AND G1 S TRANSITION” pathway gene set (Fig. 3D, and Additional file 2, Table S4). Upregulation of MYC, E2F1, and CCNB1 expression led to resistance to alpelisib. Previously, we uncovered that stem-like C3 subtype exhibited both MYC up-regulation, and SCLC enrichment (Additional file 3, Fig. S2). MYC facilitates SCLC tumorigenesis [49]. Moreover, our SCLC PDCs predicted the response to alpelisib better than cell lines (Fig. 3C). Therefore, transcriptome gene signature extracted by machine-learning determined alpelisib response for SCLC like RB1/TP53.

FOXM1 over-expression in SCLC type and sensitivity to cell cycle inhibitors

For further in-depth exploration of the candidate drugs for treating different lung cancer types, we investigated the drug responses of ADC, SCLC, and miscellaneous cancers. SCLC was sensitive to drugs targeting cell cycle pathways (adavosertib, barasertib, berzosertib, and AZD7762) and DNA damage (vorinostat; P < 0.1 and |log2FC| < 0.3; Fig. 4A). The SCLC sensitivity to alpelisib was relatively lower than that to cell cycle drugs according to the FC value. The machine-learning signatures of three drugs (adavosertib, AZD7762, and berzosertib) consistently contained cell cycle pathways (Fig. 4B, and Additional file 2, Table S4). The transcriptome profile revealed the activation of cell cycle genes in SCLC. Moreover, the DNA damage response and RAS signal were downregulated in NSCLC (Additional file 3, Fig. S4). Collectively, these comprehensive pharmacogenomics result suggest the role of cell cycle inhibition in treating SCLC types.

Fig. 4

Sensitivity of lung cancer types to drug candidates. A Volcano plot for each lung cancer type. BP-value bar plot of gene set enrichment analysis of transcriptome signatures for three drugs acquired from machine learning. C Five genes were associated with the cell cycle pathway. PDC gene expression heatmap for five genes detected in machine-learning feature importance and “CYCLIN A/B1/B2 ASSOCIATED EVENTS” pathway, and normalized AZD7762 AUCs. D Scatter plots for AZD7762 AUC and FOXM1 expression for PDCs and cell lines. EFOXM1 expression comparison among four lung cancer cohorts between SCLC and LUAD. P-values were determined using the Wilcoxon rank-sum test between SCLC and LUAD. F Protein levels of FOXM1 and β-actin in siRNA-transfected H69 cells analyzed using western blotting. Ctrl represents the control, and siFOXM1-#6 and siFOXM1-#7 are FOXM1 siRNAs. G Cell cycle phase assessment according to control and FOXM1 siRNA conditions of H69 cells. After incubation, the cells were analyzed using flow cytometry to evaluate the DNA content. Representative DNA content profiles from three independent experiments are shown. Bar plots present the proportions of cells in cell cycle phases (upper panel) and cell counts (y-axis) based on propidium iodide (PI) staining (x-axis; lower panel). H Drug response curves for testing AZD7762 sensitivity in control siRNA (Ctrl) and siRNA #6- and #7-transfected H69 cells. The half-maximal inhibitory concentration (IC50) of AZD7762 siRNA-transfected cells is shown in the panel. All experiments were performed in quadruplicate. The values represent the mean ± SEM (Student’s t-test, *P < 0.05; ***P < 0.001)

PDCs and cell lines showing FOXM1 upregulation exhibited AZD7762 sensitivity (Fig. 4C, D; PDC R = –0.77; cell line R = –0.51). FOXM1 was clearly over-expressed in SCLC compared to that in LUAD (Fig. 4E) [10, 50,51,52]. To evaluate whether FOXM1 is a plausible target to AZD7762 for SCLC, we transfected H69 cells with small interfering RNAs (siRNAs) against FOXM1 (#6 and #7, siFOXM1) and control siRNA and measured the FOXM1 protein level using western blotting. The FOXM1 expression level was decreased in siFOXM1 (Fig. 4F). siFOXM1 cells also exhibited decreased cell proliferation compared to siControl-transfected cells. However, G2-M phase arrest was detected in siFOXM1 H69 cells (Fig. 4G). siFOXM1 became more resistant to AZD7762 (AUC FC > 1.1) than siControl-transfected cells (Fig. 4H). Another cell lines H209 also presented the same result (see Additional file 3, Fig. S5). Therefore, we concluded that FOXM1 plays essential role to SCLC participating in the G2-M phase, and can be targeted by cell cycle inhibitors.

Pharmacogenomic analysis according to EGFR-TKI treatment status

Among our PDCs, we identified 27 EGFR-TKI-treated cases with corresponding NGS results (independently obtained in the clinic using tissue samples). As described in methods, we these patients into four groups (Fig. 5A). Primarily, POST3 samples (86%) were enriched in the C2 EMT-like subtype (Table 1). EGFR mutation-calling failure was observed in the clinical tumor biopsy NGS record of some patients from which BASELINE (n = 2) and POST2 (n = 1) PDCs were derived. The mutation-calling failure likely originated from differences in the platform between the clinical biopsy and laboratory PDC testing. Moreover, TP53 non-hotspot mutations were present in all of the POST2 PDCs and were markedly absent in POST3 PDCs. EGFR T790Maq was identified in only POST2 PDCs (Fig. 5B red asterisk); interestingly, BRAF mutations highly co-occurred with EGFR T790Maq mutations (Fig. 5B and Additional file 3, Fig. S3). Thus, our PDC models significantly reflected EGFR mutation, acquisition, and extinction status according to EGFR-TKI therapies.

Fig. 5

Genomic profile and drug sensitivity of PDCs, and functional validation using cell lines. A Therapeutic groups were categorized into two arms. The first arm was BASELINE (yellow) to POST1 (light green). The second was BASELINE to POST2 (dark green) to POST3 (orange). B Heatmap presenting the mutation profile, RNA subtype, and tumor mutation burden (TMB) of each group. In the POST1 group, the red asterisk indicates the EGFR double mutation including T790M. C Heatmap of the DEGs for the four groups. DEGs and pathways (GSEA P < 0.1) for each group are denoted on the right. D Two scatter plots of scores extracted from EGFR-TKI resistance pathway gene pairs. ERBB2 and MET were upregulated in the POST2 group (dotted circle), and YAP/TAZ-AXL activated in the POST3 group (dotted circle). R indicates Pearson’s correlation coefficient, and the P value was acquired using the Wilcoxon rank-sum test from the AUC values. E Volcano plot to test drug sensitivity for each group. The x-axis is log2 fold change (FC) and the y-axis is P-value. F Drug response curves based on the viability of H1975, H1975_OR3, and OR4 cells treated with osimertinib. G Bar plots of YAP1 and its target gene expressions. Fold increases (x-axis) in mRNA expression in each gene (y-axis) were assessed using RT-PCR in H1975, H1975_OR3, and OR4 cells. H The protein level expression assessment in three cells using western blot analysis. I Cell migration across the transmembrane. Cell morphology images are shown from three cases. The bar plot indicates the average number of migrating cells (y-axis) for each condition (x-axis). J EpCAM expression (x-axis) determined using flow cytometry. Graphs were analyzed using the FlowJo program. K Drug response curves treated with XAV939. The AUC of each condition for XAV939 is shown in the panel. All experiments were performed in quadruplicate. The values represent the mean ± SEM (Student’s t-test, *P < 0.05; **P < 0.01; ***P < 0.001)

The four therapeutic groups exhibited distinct regulatory programs. To identify the resistance pathway, we applied two different methods. First, activated pathways were inferred from upregulated DEGs for each group (see Additional file 2, Table S5). The activation of ATM signaling and the DNA damage response according to CDK5 upregulation were enriched in the BASELINE PDCs. POST1 PDCs predisposed to RAS, IL-6, WNT, and PI3K-AKT-mTOR signals were enriched with upregulated genes (BRAP, KRAS, and NF1). POST2 PDCs were enriched in NOTCH and STAT3 signaling, including DLL4, EP300, and STAT3 over-expression. Finally, POST3 PDCs exhibited activated TGF-β signaling regulated by PDK1, SMAD2, TGFBR3, and ZEB1. The over-expression of TGF-β signaling and ZEB1 demonstrated that osimertinib-resistant cancer exhibited EMT pathway activation (Fig. 5C).

We next assessed the activities of 12 knowledge-based therapeutic resistance pathways for each sample according to the four EGFR-TKI therapeutic groups (Additional file 2, Table S6). The SERPINE1 signature was only activated in the BASELINE PDCs (Additional file 3, Fig. S6). The WNT, ERBB2, and MET pathway scores were increased in POST2 PDCs but were decreased in POST3 PDCs. The YAP/TAZ and PI3K-AKT pathways were elevated in POST3 PDCs. MET and ERBB2 gene expression also showed positive correlations with activation in the POST2 PDCs (R = 0.62), whereas the POST3 PDCs exhibited inactivation in both pathways (Fig. 5D). The EMT-like POST3 PDCs exhibited the over-expression of both YAP/TAZ and AXL (R = 0.86; Fig. 5D).

Among the four EGFR-TKI groups, we explored the sensitivity of POST3 PDCs to etoposide and XAV939 (Fig. 5E). To predict drugs for 27 PDCs, we prepared an additional drug screening dataset by including the extended PDC set without available sequencing data (n = 70; Additional file 3, Fig. S7A). To make up for the lack of genomic profiles, we collected EGFR-TKI therapy and clinical NGS test medical records. We could assign these samples to the four groups (BASELINE, n = 12; POST1, n = 23; POST2, n = 17; POST3, n = 18). The extended PDC dataset showed remarkable concordance in patient categorization obtained using NGS results and medical records. Finally, POST3 PDCs (PDC P < 0.2, extended PDC P < 0.083; Fig. 5E) were sensitive to etoposide and XAV939. We found that POST3 PDCs exhibited higher sensitivity to XAV939 than etoposide. The machine-learning XAV939 gene signature revealed that the response to XAV939 also involved EMT-like molecular features (i.e., collagen formation and extracellular matrix organization; Additional file 2, Table S4).

To evaluate the previous finding about XAV939 selectivity and pathways, we created two cell lines—EGFR-T790Maq (H1975) cells using POST2 PDCs and XAV939-resistant (H1975_OR3, H1975_OR4) cells using POST3 PDCs—with osimertinib exposure. Both cell lines acquired osimertinib resistance (FC > 1.7; Fig. 5F). TCGA pan-cancer and reverse transcription-quantitative polymerase chain reaction analyses (RT-PCR) revealed that osimertinib-resistant cells concurrently exhibited over-expression of YAP/TAZ target genes [46]. Ten target genes, including YAP1 and AXL, were considerably upregulated in the osimertinib-resistant cell lines compared to that in parent H1975 cells (Fig. 5G). In particular, AXL was elevated with the highest FC (5.1–9.0), and YAP1 and AXL proteins showed consistent over-expression in osimertinib-resistant cell lines (Fig. 5H).

Next, we identified osimertinib resistance-related molecular subtypes and cellular alterations in POST3 PDCs. Since POST3 PDCs were classified as an EMT-like subtype (Fig. 5C), we assessed the EMT features of osimertinib-resistant PDCs. Two osimertinib-resistant PDCs showed decreased E-cadherin and EpCAM (epithelial markers) and increased N-cadherin and Vimentin (mesenchymal markers; Fig. 5H) expression levels. Subsequently, using flow cytometry, we established that cell surface EpCAM expression was considerably decreased in these cell lines (Fig. 5I). We further conducted a migration assay to verify these EMT-associated molecular changes (Fig. 5J). As expected, the two osimertinib-resistant PDCs had greater migration ability than H1975 cells. Collectively, these findings demonstrated that the two successfully generated osimertinib-resistant cell lines were representative of POST3 PDCs. Next, we checked if XAV939 had similar inhibitory effects on POST3 PDCs and H1975, H1975_OR3, and H1975_OR4 cells. XAV939 was not cytotoxic to H1975 cells (AUC = 105.5); however, H1975_OR3 (93.2) and OR4 (87.8) cells were clearly more sensitive to XAV939 (Fig. 5K). Overall, we validated that osimertinib resistance facilitates YAP/TAZ and AXL activation, and an EMT-like phenotype. Thus, XAV939 may be used to treat osimertinib-resistant tumors.