The lack of knowledge about the relationship between tumor genotypes and therapeutic responses remains one of the most critical gaps in enabling the effective use of cancer therapies. Here, we couple a multiplexed and quantitative experimental platform with robust statistical methods to enable pharmacogenomic mapping of lung cancer treatment responses in vivo. The complex map of genotype-specific treatment responses uncovered that over 20% of possible interactions show significant resistance or sensitivity. Known and novel interactions were identified, and one of these interactions, the resistance of KEAP1-mutant lung tumors to platinum therapy, was validated using a large patient response data set. These results highlight the broad impact of tumor suppressor genotype on treatment responses and define a strategy to identify the determinants of precision therapies.

Significance:

An experimental and analytical framework to generate in vivo pharmacogenomic maps that relate tumor genotypes to therapeutic responses reveals a surprisingly complex map of genotype-specific resistance and sensitivity.

Efforts over the past decade have generated many novel cancer therapies (1, 2). However, patient responses are heterogeneous, with some patients responding well and others showing limited or no response (3, 4). Although it is believed that the genetic complexity of cancer underlies a significant portion of the variation in therapeutic response, the map of such pharmacogenomic interactions is currently lacking (5–7). Despite widespread tumor genotyping, only a few driver mutations currently inform clinical treatment decisions and clinical trial designs (8–10). This is driven by the fact that we do not yet know which tumor suppressor alterations influence sensitivity or resistance to specific therapies. The very premise that tumor suppressor genotype substantially affects therapeutic responses remains largely untested.

The pharmacogenomic landscape of cancer drug responses has been investigated using cell lines, patient-derived xenografts (PDX), and patient treatment outcome data (5, 11–14). However, such genotype–treatment interactions are notoriously difficult to measure using these systems for four major reasons: the large numbers of driver and passenger mutations, the observational instead of manipulative nature of the experiments, lack of the appropriate autochthonous in vivo environment, and the high stochasticity of tumor growth. Specifically, cell lines grown in vitro lack the appropriate in vivo environment, do not represent all cancer subtypes, and often carry additional alterations that arise during passaging (15). PDXs and human cell line transplantation models recapitulate some aspects of in vivo growth, but growth factor/receptor incompatibility, growth in nonorthotopic sites, and the obligate absence of the adaptive immune system compromise these approaches (16–18). Furthermore, human tumor–derived systems almost invariably have large numbers of mutations and genomic alterations. Thus, even large-scale analyses often lack the statistical power to glean cause-and-effect relationships between individual genomic alterations and therapeutic responses (5, 14). The same logic applies to patient treatment response data, which are generally too limited in scale to provide sufficient statistical power to confidently associate tumor suppressor genotypes with metrics of clinical response (19). Such data are particularly sparse for unapproved therapies (limited to clinical trial results) and are nonexistent for preclinical therapeutic candidates.

A cost-effective system that introduces defined genomic alterations, measures the response of a large number of isogenic tumors, and recapitulates the in vivo physiologic context could be valuable for uncovering genotype–treatment relationships. Here we present such a system based on tumor barcoding in genetically engineered mouse models. Genetically engineered mouse models of human cancer are important preclinical models, as they recapitulate the physiologic, tissue, and immunologic context of tumor growth (20, 21). These models uniquely enable the introduction of defined genomic alterations into adult somatic cells, which leads to the generation of autochthonous tumors (20). These tumors can recapitulate the genomic alterations, gene-expression state, histopathology, and therapy-refractive nature of corresponding human cancers (11, 22). Despite the potential value of these models in preclinical translation studies, the breadth of their utility has been limited in practice by the fact that they are neither readily scalable nor sufficiently quantitative (23–27).

To increase the scope and precision of in vivo cancer modeling and to assess tumor suppressor gene function in a multiplexed manner, we previously developed a system that couples tumor barcoding with high-throughput barcode sequencing (Tuba-seq; ref. 26). This method integrates CRISPR/Cas9-based somatic genome engineering and molecular barcoding into well-established Cre/Lox-based genetically engineered mouse models of oncogenic Kras-driven lung cancer (28). The initiation of lung tumors with pools of barcoded Lenti-sgRNA/Cre viral vectors enables the generation of many tumors of different genotypes in parallel. All neoplastic cells within each clonal tumor have the same two-component barcode, in which an sgID region identifies the sgRNA and a random barcode (BC) is unique to each tumor. Thus, high-throughput sequencing of the sgID-BC region from bulk tumor-bearing lungs can quantify the number of neoplastic cells in each tumor of each genotype (28). Previous Tuba-seq studies quantify tumor suppressor effects and their interaction with other tumor suppressor genes, focusing only on comparisons within mice (28–30). Comparisons of tumor distributions across mice are more challenging and required improvements in accuracy as well as new analytical methods.

Here, we optimize multiple key aspects of the Tuba-seq approach. The greatly improved accuracy in tumor calling enabled us to compare tumor size distributions between groups of mice, i.e., treated and untreated groups, and to generate a large-scale map that relates tumor genotype to therapeutic responses in vivo. We developed a new analytical and computational framework, pharmacogenomic tumor barcoding with high-throughput barcode sequencing (PGx-Tuba-seq). We quantify the treatment responses of tens of thousands of oncogenic KRAS-driven lung tumors of eleven different tumor suppressor genotypes to a diverse panel of therapies, and uncover a surprisingly complex pharmacogenomic map of resistance and sensitivity. PGx-Tuba-seq represents a more tractable method to uncover the therapeutic response of different tumor genotypes than previous in vitro and in vivo screening approaches.

### Mice, tumor initiation, and drug treatment

All animal experiments have been approved by Institutional Animal Care at Stanford University with protocol number 26696. Lung tumors were initiated by intratracheal delivery of the same lentiviral pools (26). 1.1 × 105 and 2.2 × 104 infectious unit/mouse were administered to each KrasLSL-G12D (K); R26LSL-Tomato (T) (hereafter KT), and KT;H11LSL-Cas9 mouse (31–33), respectively. Drug treatments were started 15 weeks after tumor initiation. For the main pharmacogenomic mapping experiment, mice were assigned to eight treatment arms or were left untreated for 3 weeks (Fig. 1A; Table 1).

Figure 1.

Optimization of tumor-barcoding coupled with Tuba-seq for the analysis of GSTRs in vivo. A, Overview of Tuba-seq pipeline to uncover GSTRs. The Lenti-TSPool/Cre viral pool contains barcoded vectors with sgRNAs targeting 11 putative tumor suppressors that are frequently mutated in human lung adenocarcinoma. Tumors are initiated in either KrasLSL-G12D/+;R26LSL-Tom (KT) or KrasLSL-G12D/+;R26LSL-Tom;H11LSL-Cas9 (KT;H11LSL-Cas9) mice. Following tumor development, mice are treated with therapies, and barcode sequencing libraries are prepared from each tumor-bearing lung. Multiple technical advances in the pipeline involve viral production, library preparation, sequencing, and analysis pipeline have been made, boosting the accuracy of our pipeline to enable many further applications. B, Stringent filtering effectively eliminated spurious tumors. Analysis of the barcodes associated with the sgID specific for the spike-in control cells (three cell lines with a defined sgID-BC added at 5 × 105 cell/sample as the benchmark) enables identification of recurrent barcode reads generated from sequencing and other errors (spurious tumors). Data are from a typical lane of 22 multiplexed Tuba-seq libraries from KT;H11LSL-Cas9 mice with Lenti-TSPool/Cre-initiated tumors. C, The relative size of tumors of each genotype in KT;H11LSL-Cas9 mice 18 weeks after tumor initiation with Lenti-sgTSPool/Cre. The relative sizes of tumors at the indicated percentiles were calculated from the tumor size distribution of all tumors in 5 mice. Error bars, 95% confidence intervals.

Figure 1.

Optimization of tumor-barcoding coupled with Tuba-seq for the analysis of GSTRs in vivo. A, Overview of Tuba-seq pipeline to uncover GSTRs. The Lenti-TSPool/Cre viral pool contains barcoded vectors with sgRNAs targeting 11 putative tumor suppressors that are frequently mutated in human lung adenocarcinoma. Tumors are initiated in either KrasLSL-G12D/+;R26LSL-Tom (KT) or KrasLSL-G12D/+;R26LSL-Tom;H11LSL-Cas9 (KT;H11LSL-Cas9) mice. Following tumor development, mice are treated with therapies, and barcode sequencing libraries are prepared from each tumor-bearing lung. Multiple technical advances in the pipeline involve viral production, library preparation, sequencing, and analysis pipeline have been made, boosting the accuracy of our pipeline to enable many further applications. B, Stringent filtering effectively eliminated spurious tumors. Analysis of the barcodes associated with the sgID specific for the spike-in control cells (three cell lines with a defined sgID-BC added at 5 × 105 cell/sample as the benchmark) enables identification of recurrent barcode reads generated from sequencing and other errors (spurious tumors). Data are from a typical lane of 22 multiplexed Tuba-seq libraries from KT;H11LSL-Cas9 mice with Lenti-TSPool/Cre-initiated tumors. C, The relative size of tumors of each genotype in KT;H11LSL-Cas9 mice 18 weeks after tumor initiation with Lenti-sgTSPool/Cre. The relative sizes of tumors at the indicated percentiles were calculated from the tumor size distribution of all tumors in 5 mice. Error bars, 95% confidence intervals.

Close modal
Table 1.

Treatments tested using the PGx-Tuba-seq platform.

Monotherapy Palbociclib 100 mg/kg Daily Oral gavage
Monotherapy Everolimus 10 mg/kg Daily Oral gavage
Monotherapy Phenformin 100 mg/kg Daily Oral gavage
Monotherapy Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
Monotherapy Trametinib 0.3 mg/kg Daily Oral gavage
Combination Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
+ Trametinib 0.3 mg/kg Daily Oral gavage
Combination Carboplatin 50 mg/kg Every five days Intraperitoneal injection
+ Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
Combination Carboplatin 50 mg/kg Every five days Intraperitoneal injection
+ Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
+ Trametinib 0.3 mg/kg Daily Oral gavage
Monotherapy Palbociclib 100 mg/kg Daily Oral gavage
Monotherapy Everolimus 10 mg/kg Daily Oral gavage
Monotherapy Phenformin 100 mg/kg Daily Oral gavage
Monotherapy Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
Monotherapy Trametinib 0.3 mg/kg Daily Oral gavage
Combination Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
+ Trametinib 0.3 mg/kg Daily Oral gavage
Combination Carboplatin 50 mg/kg Every five days Intraperitoneal injection
+ Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
Combination Carboplatin 50 mg/kg Every five days Intraperitoneal injection
+ Paclitaxel 20 mg/kg Every other day Intraperitoneal injection
+ Trametinib 0.3 mg/kg Daily Oral gavage

### Tuba-seq library generation

Genomic DNA was isolated from bulk tumor-bearing lung tissue from each mouse (26). Three benchmark control cell lines (∼5 × 105 cells/cell line) were added to each mouse lung sample prior to lysis to enable the calculation of the absolute neoplastic cell number within each tumor (28). To reduce the errors of the Tuba-seq pipeline by orders of magnitude, we implemented multiple critical changes to the library preparation, sequencing, and analysis (Tables 2 and 3). Q5 High-Fidelity 2x Master Mix (NEB, M0494X) was used to amplify the sgID-BC region from 32 μg of genomic DNA (34). To improve sequencing quality, we used unique dual-indexed primers and added 6–9 random nucleotides (Ns) to the flanking ends of both index primers before the sequence-specific primer regions (35). The libraries were pooled based on lung weight to ensure even read depth and sequenced on an Illumina HiSeq 2500 platform (Admera Health) with paired-end 150 bp reads.

Table 2.

Overview of our optimized Tuba-seq analysis pipeline for calling sgID-BCs from sequencing data and determining the number of neoplastic cells in each tumor.

Table 3.

Comparison of our current pipeline with our previous Tuba-seq pipeline.

ModulePrevious (Rogers et al., 2017)CurrentPurpose
Viral production Pooled Each viral vector was prepared separately Eliminate lentiviral template switching
Library preparation Taq polymerase Q5 polymerase Reduce PCR errors
Library preparation Single indexing Dual unique indexing Eliminate the impact of index hopping during sequencing on tumor calling
Sequencing Single-end Paired-end Reduce “spurious tumors” created by sequencing errors
Read processing and tumor calling DADA2 clustering Stringent filtering on reads Eliminate “spurious tumors” created by PCR and sequencing errors
Remove spurious tumors recursively based on hamming distance
Read processing and tumor calling No restriction on BC length Require exact length match Eliminate “spurious tumors” created by PCR and sequencing errors
ModulePrevious (Rogers et al., 2017)CurrentPurpose
Viral production Pooled Each viral vector was prepared separately Eliminate lentiviral template switching
Library preparation Taq polymerase Q5 polymerase Reduce PCR errors
Library preparation Single indexing Dual unique indexing Eliminate the impact of index hopping during sequencing on tumor calling
Sequencing Single-end Paired-end Reduce “spurious tumors” created by sequencing errors
Read processing and tumor calling DADA2 clustering Stringent filtering on reads Eliminate “spurious tumors” created by PCR and sequencing errors
Remove spurious tumors recursively based on hamming distance
Read processing and tumor calling No restriction on BC length Require exact length match Eliminate “spurious tumors” created by PCR and sequencing errors

### Processing reads to identify the sgID and barcode and removal of “spurious tumor” generated by read errors

We required both the forward and reverse sequencing reads to match perfectly within the BC region. FASTQ files were processed to identify the sgID and BC counts for each tumor. The sgID region identified the targeted tumor suppressor gene. The number of reads with each unique sgID-BC in each sample was summed to calculate each putative tumor's size. PCR and sequencing errors within the random barcode regions generate spurious tumors. We used stringent criteria to reduce and even eliminate the effects of PCR and sequencing errors on tumor calls, greatly reducing spurious tumors (Fig. 1B; Supplementary Fig. S1A) when quantifying relative tumor sizes (Fig. 1C; Supplementary Fig. S1B–S1D), showing larger effect sizes (Supplementary Fig. S2A–S2D).

### Developing unbiased procedures for detecting genotype-specific drug effects

Previous Tuba-seq analyses focused on comparing the sizes of tumors of different genotypes within individual mice (28, 30). Such analyses are largely robust to multiple sources of variation among mice (Supplementary Fig. S3A–S3D). Here we needed to compare tumor sizes between the untreated and treated groups when analyzing genotype-specific drug responses. We used the same viral pool to initiate tumors in all mice; therefore, the relative representation of transduced epithelial cells containing each Lenti-sgRNA/Cre is constant and does not vary across mice.

### Null model of tumor responses with no genotype specificity

We assume that the therapy affects all tumors proportionally to their sizes such that the size of each tumor changes from X to X1 = X × S after the treatment, where S is the proportion of remaining cancer cells. Under the null model (H0) of no genotype-specific drug responses, S is constant and does not depend on tumor genotype. Under the alternative model H1, S varies depending on the genotype: SsgID,j = SInert × (1 + Gj), with Gj representing the genotype-specific therapeutic response (GSTR) of tumors of specific genotypes to the drug j. If Gj > 0, the inactivation of the tumor suppressor gene confers relative resistance; if Gj < 0, the inactivation of the tumor suppressor gene confers relative sensitivity.

The most extreme treatment reduced tumor sizes by ∼87%. Although the depth of sequencing varied across mice and treatments, we wanted to reliably identify tumors in each treated and untreated mouse. Thus, we chose to use the cutoff of L = 1,000 cells in the untreated mice, allowing reliable detection and accurate size estimates of tumors in each mouse.

### Calculation of proportional size reduction as the drug effect

To estimate the value of the tumor reduction factor S that leads to the best match between the distributions of inert tumors (those with inert sgRNAs) in the treated and untreated groups, we calculated the value of S such that the median number of shrunk tumors across all the untreated mice was closest to the median of the number of observed tumors with the size above or equal to 1,000 cells across all the mice in the treated group.

### Using relative tumor number (ScoreRTN) to estimate GSTR

Our first approach defines response as the number of tumors that exceed a minimum size threshold (Fig. 2A and B). The null hypothesis for each genotype is that the number of tumors above the cutoff L in the untreated mice should match the number above the new cutoff L×S in the treated mice. If a GSTR exists, the tumors with a specific sgID are more resistant to the drug than the inert tumors, and more of such tumors should remain above the adjusted cutoff of L×S than expected, whereas if they are more sensitive, fewer such tumors should remain above the adjusted cutoff of L×S. We first calculate the ratio of the number of tumors above the cutoff L in the untreated mice of a particular sgID to that of the inert tumors (RTNi,j,L),

Although ScoreRTN and ScoreRGM may have an emphasis on different aspects of GSTR on tumor size distribution, it is helpful to have a single combined score. We calculated a combined score of GSTR (⁠|\hat{G}$|⁠) by taking the inverse variance weighted average of ScoreRTN and ScoreRGM, then converting it to the linear scale (Fig. 3). Figure 3. Tuba-seq quantifies GSTR to multiple therapies. A, Timeline of the experiment. Tumors were initiated in KT;H11LSL-Cas9 mice with the barcoded Lenti-sgTSPool/Cre. Three weeks of treatment was initiated after 15 weeks of tumor growth. The number of mice used for each treatment arm is shown. B–D, The estimated genotype-specific treatment response (Ĝ) calculated from the inverse variance weighted average of ScoreRTN and ScoreRGM for the pharmacogenomic mapping experiment (B), negative control experiment in KT mice (C), and palbociclib repeat experiment (D). Stars represent significant effects. Figure 3. Tuba-seq quantifies GSTR to multiple therapies. A, Timeline of the experiment. Tumors were initiated in KT;H11LSL-Cas9 mice with the barcoded Lenti-sgTSPool/Cre. Three weeks of treatment was initiated after 15 weeks of tumor growth. The number of mice used for each treatment arm is shown. B–D, The estimated genotype-specific treatment response (Ĝ) calculated from the inverse variance weighted average of ScoreRTN and ScoreRGM for the pharmacogenomic mapping experiment (B), negative control experiment in KT mice (C), and palbociclib repeat experiment (D). Stars represent significant effects. Close modal If |\hat{G}$| > 0, the genotype is relatively resistant to the therapy, and if |\hat{G}$| < 0, the genotype is relatively sensitive to the therapy. To be conservative, for the combined score to be called significant, we require at least one significant P value (P < 0.05), and one marginally significant P value (P < 0.1) for the two statistics ScoreRTN and ScoreRGM. ### Comparing with human cell line response database GDSC The drug-sensitivity data from human cell lines were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) database (www.cancerrxgene.org; ref. 5). Due to the limited number of LUAD cell lines, we focused on comparing the results from Pan-cancer cell lines. All 5 monotherapies used in our study were assessed by GDSC. Except for Keap1 and Rbm10, which are not reported for everolimus and paclitaxel, the GSTR of all other 51 gene–drug pairs were quantified by GDSC. The effect size and FDR-corrected P values were used for comparison. ### Analysis of clinical data for resistance to chemotherapy Despite relatively widespread genotyping, clinical treatment data and response data are extremely limited. Memorial Sloan Kettering Cancer Center (MSKCC, New York, NY) has a tremendous program to genotype patients and to collect clinical data. Most patients with oncogenic KRAS-driven lung cancer get platinum doublet therapy. Patients with metastatic or recurrent lung adenocarcinoma harboring a KRAS mutation in codons 11, 12, or 61, as detected by MSK-IMPACT (36), were reviewed. Patients who received platinum chemotherapy (carboplatin or cisplatin) with pemetrexed ± bevacizumab as first-line treatment were included (n = 216). Treatment efficacy was measured as time of first treatment with platinum doublet chemotherapy to start of next systemic therapy, or death if no subsequent therapy was received. Patients who continued on platinum doublet therapy at the last follow-up were censored. The retrospective chart review was approved by the MSK Institutional Review Board. Kaplan–Meier estimator plots of time-to-next treatment for patients with and without mutations at each of the 11 tumor suppressor genes of interest were generated. In addition, a multivariable Cox proportional hazards model analysis was performed, integrating the mutational status of the 11 genes as individual input features to assess the independent effect of co-occurring mutations. ### Data availability statement The sequencing data set generated and analyzed during the current study is available in the Gene-Expression Omnibus database (accession code: GSE146448). Other data and relevant code are available in https://github.com/lichuan199010/Tuba-seq-analysis-and-summary-statistics. ### Development of the PGx-Tuba-seq pipeline To eliminate sgRNA-sgID/barcode uncoupling due to lentiviral template switching and to minimize PCR, sequencing, and clustering errors, we made multiple improvements to our Tuba-seq experimental protocols and analysis pipeline (Fig. 1A; Tables 2 and 3; and Materials and Methods; ref. 26). We initiated lung tumors in KrasLSL-G12D/+;Rosa26LSL-Tomato;H11LSL-Cas9 (KT;H11LSL-Cas9) mice and control Cas9-negative KT mice with a pool of barcoded Lenti-sgRNA/Cre vectors targeting eleven putative tumor suppressors and four control vectors with inert sgRNAs (Lenti-sgTSPool/Cre; Fig. 1A). To eliminate template switching during lentiviral reverse transcription, we generated each vector separately and pooled each viral vector immediately prior to tumor initiation (37). Tumor suppressors were selected based on common occurrence in human lung adenocarcinomas and previously suggested roles in oncogenesis (26). Eighteen weeks after tumor initiation, the sgID-BC region from each bulk tumor-bearing lung was PCR amplified and sequenced to quantify the number of neoplastic cells in each tumor (Fig. 1A). Our new analysis pipeline essentially eliminated the impact of read errors, as assessed by two metrics, including the spurious tumors generated from spike-in barcodes with known sequences and correspondence of tumor barcodes with those from the lentiviral plasmid pool (Fig. 1B; Supplementary Fig. S1A). Quantification of the impact of tumor suppressor gene inactivation on tumor growth in KT;H11LSL-Cas9 mice using our optimized method uncovered effects that were generally consistent with our previous analyses, but with greater magnitudes of tumor suppression (Fig. 1C; Supplementary Figs. S1C and S1D and S2A–S2D; sign test for differences in magnitudes, P = 0.001; ref. 28). Consistent with the robustness of our methods, analysis of the KT mice with Lenti-sgTSPool/Cre-initiated tumor revealed no false-positive tumor-suppressive effects (Supplementary Fig. S1C and S1D). These technical improvements to the Tuba-seq method further enhance the ability of this technology to be applied to study a variety of questions in tumor progression and evolution, as well as quantification of the pharmacogenomic interactions as performed in this study. When quantifying tumor suppressor gene effects using Tuba-seq, each mouse represents an internally controlled experiment in which metrics of tumor size can be compared between tumors of each tumor suppressor genotype and tumors initiated with inert sgRNAs within the same mouse (Fig. 1C; Supplementary Fig. S1B–S1D; ref. 26). In contrast, comparing tumor size distributions between groups of mice, such as between untreated and drug-treated groups, requires methods that address the technical and biological differences among mice. To understand the statistical properties and potential biases intrinsic to this type of analysis, we rigorously modeled drug responses and genotype-specific responses. We initially performed our modeling with the assumption that cancer cells in tumors of all sizes respond equally to each treatment, whereas the treatment effects can vary by genotype. Specifically, we estimated the drug effect on inert tumors and then applied this effect to all tumors to calculate an expected distribution of tumor sizes after treatment (Fig. 2A and Materials and Methods). GSTRs were quantified by comparing the observed distribution of tumor sizes for tumors of a certain genotype after treatment with the expected distribution derived from the untreated mice. We developed two statistics to characterize GSTRs: (i) ScoreRTN—relative tumor number, which compares the relative numbers of tumors above a certain size after treatment; and (ii) ScoreRGM—relative geometric mean, which constitutes the relative change in the geometric mean of tumors from the full distribution of tumor sizes (Materials and Methods). By assessing the performance of the two statistics, we showed that both statistics are unbiased (Supplementary Fig. S3E–S3H) and exhibit substantial and similar power (Supplementary Fig. S4A–S4C), although one statistic may outperform the other if the genotype-specific response is not uniform across tumor sizes (Methods; Fig. 2B and C; Supplementary Fig. S5A and S5B). Moreover, by performing power analysis and plotting the ROC curves for both statistics across multiple sample sizes (i.e., number of mice/group), we confirmed the high sensitivity and specificity of our system (Fig. 2B and C; Supplementary Fig. S4A–S4C). We also found that relaxing the assumption that tumors of all sizes respond proportionally to treatment did not change our results substantially (Supplementary Fig. S5A and S5B). ### Complex pharmacogenomic map uncovered using the PGx-Tuba-seq pipeline We applied Tuba-seq and our statistical metrics to assess the GSTRs of 11 genotypes of lung tumors to a panel of eight single and combination therapies (Figs. 1A and 3A; Table 1). These therapies were chosen to perturb diverse signaling pathways and assess the genotype dependency of chemotherapy responses. KT;H11LSL-Cas9 mice with Lenti-sgTSPool/Cre-initiated lung tumors were treated for three weeks with one of the eight therapies followed by Tuba-seq analysis (Figs. 1A and 3A). The total cancer cell numbers estimated by Tuba-seq were highly correlated with total tumor-bearing lung weights, which varied substantially among mice even within the same groups (Supplementary Fig. S6A–S6C). Despite expected mouse-to-mouse variation, comparing the overall tumor burden and the number of tumors with inert sgRNAs in the untreated and treated mice revealed significant overall therapeutic effects for five out of the eight treatments (Supplementary Fig. S6D). We compared the tumor size profiles of treated mice with those of untreated mice and calculated the ScoreRTN and ScoreRGM (Supplementary Fig. S7A). For both statistics, we estimated the magnitudes of GSTRs and the associated P values using bootstrapping. Across all genotypes and treatments, the two statistics were well correlated in magnitude as expected under the model of proportional tumor responses (Supplementary Fig. S7B; r = 0.86, P = 10−46). Among the 88 assessed genotype–treatment pairs, 20 and 17 significant GSTRs (P < 0.05) were identified by ScoreRTN and ScoreRGM, respectively. Of these, 19 genotype–treatment interactions were significant by one statistic (P < 0.05) and at least marginally significant (P < 0.1) by the other (Supplementary Fig. S7A and S7B; Supplementary Table S1). We derived a composite measure of GSTR (⁠|\hat{G}$|⁠) with the magnitude estimated from the inverse variance weighted average of the two statistics (Materials and Methods, Fig. 3B). Analysis of genotype-specific effects across treatments highlighted similarities among tumor suppressors, including those of Lkb1 and Setd2 that we have previously suggest to have redundant tumor-suppressive effects (28). Furthermore, combination treatments clustered with their corresponding single therapies (Supplementary Fig. S7C and S7D), and an additive model shows good predictive power (Supplementary Fig. S7E and S7F). Power analysis showed that our findings were robust to the cancer cell number cutoff (Supplementary Fig. S8A), choice of inert sgRNAs (Supplementary Fig. S8B), and inaccurate estimation of drug effects (Supplementary Fig. S9A and S9B).

One of the detected GSTRs was well known in advance—the resistance of Rb1-deficient tumors to the CDK4/6 inhibitor, palbociclib. Our ability to rediscover this interaction serves as a positive control of our method and is consistent with the expectation that some pharmacogenomic interactions transcend cancer types (Supplementary Fig. S10A–S10E). This resistance is consistent with the biochemical features of this pathway (Supplementary Fig. S10F) and clinical findings in breast cancer and hepatocellular carcinoma (38–40).

To further test the performance of our experimental and statistical procedures, we performed two additional experiments. First, as a negative control for GSTR identification, we treated Cas9-negative KT mice with a combination of chemotherapy and MEK inhibition (Supplementary Fig. S11A). This treatment led to a dramatic reduction in tumor sizes compared with untreated KT mice (Supplementary Fig. S11B). Only one false-positive GSTR was identified (ScoreRTN, P = 0.03; ScoreRGM, P = 0.07) with a very weak magnitude of the effect (⁠|\hat{G}$| = 0.093, whereas the minimum magnitude of significant GSTR interactions in the main experiment was 0.108; Fig. 3C; Supplementary Fig. S11C). Furthermore, none of the individual inert sgRNAs (sgNeo1, sgNeo2, sgNeo3, and sgNT) had a significant effect by either metric for any of the eight treatments in our main pharmacogenomic mapping experiment, adding confidence in the veracity of the detected GSTRs (Fig. 3B and C). Simulations suggest that these cohort sizes have substantial albeit imperfect power (Supplementary Fig. S4A–S4C); therefore, we next attempted to rediscover the genotype–palbociclib interactions. We initiated tumors in a similar, yet somewhat smaller cohort of KT;H11LSL-Cas9 mice with Lenti-sgTSPool/Cre and repeated the palbociclib treatment. Analyses of these mice again identified Rb1 inactivation as a mediator of palbociclib resistance (Fig. 3D; Supplementary Fig. S10B). Smad4-deficient tumors, which showed modest resistance in our initial experiment, showed nominal resistance in the repeat experiment (⁠|\hat{G}$| = 0.167), although this interaction was not significant (P = 0.17 and 0.20 for ScoreRTN and ScoreRGM, respectively). Given the magnitude of this GSTR and our sample sizes, this false negative is not surprising. Assuming a true positive rate of 80%, which is considered desirable (41, 42), when identifying two genuine GSTR signals (Smad4 and Rb1, for instance) in two independent experiments, the probability of missing at least one of these findings is 1–80%4 = 59%.

### Multiple sources of evidence confirm the findings of our PGx-Tuba-seq analysis

Although most of the detected pharmacogenomic interactions we uncovered are novel, several lines of evidence derived from clinical and preclinical data are consistent with our observations. For instance, Lkb1-inactivation reduced sensitivity to mTOR inhibition in our data, which is supported by a previous in vitro study (43) and anecdotal data from the analysis of lung adenocarcinoma patient-derived primary cultures (Supplementary Fig. S12A–S12C; ref. 12). Moreover, previous studies have shown that KrasG12D;Lkb1−/− lung tumors are sensitive to phenformin (25) and resistant to MEK inhibition (23).

The ultimate goal of our study was to find genotype–treatment responses that predicted lung adenocarcinoma patient responses. Lung adenocarcinoma patients are often treated with first-line platinum-containing combination therapies. In our analysis, Keap1-inactivation led to resistance to treatments that included carboplatin, while not promoting significant resistance to the other therapies (Fig. 3B). Interestingly, KEAP1 inactivation has been previously suggested to reduce responses to chemotherapy (44–46). To further investigate the clinical impact of tumor suppressor genotype on lung adenocarcinoma responses, we queried the tumor suppressor genotype and therapeutic benefit of platinum-containing treatments (assessed as time-to-next-treatment) of 216 patients with oncogenic KRAS-driven human lung adenocarcinomas treated at Memorial Sloan Kettering Cancer Center (Materials and Methods). When each gene was assessed individually (Supplementary Fig. S13A–S13K), both KEAP1 and LKB1 mutations were associated with worse clinical outcomes (P = 6 × 10−6, Fig. 4A and P = 0.06; Supplementary Fig. S13C and S13J, respectively). However, the marginally significant effect of LKB1 mutation appears to be driven primarily by the co-occurrence of KEAP1 and LKB1 mutations (Supplementary Fig. S13L; refs. 47, 48). This finding is also well supported by our pharmacogenomic data in which Lkb1 inactivation did not confer resistance to platinum-containing treatments (Fig. 3B). We further quantified the hazard ratio of the mutational status of the 11 genes accounting for the effect of other coincident mutations. This analysis confirmed that mutations of KEAP1 correlated with a shorter time-to-next treatment, which is consistent with our Tuba-seq results as well as previous studies on the impact of KEAP1/NRF2 pathway alterations on platinum responses (Fig. 4A and B; refs. 44, 49). Our in vivo pharmacogenomic platform, in which the responses of tumors with defined genotypes can be quantified, establishes direct causal relationships between genotype and treatment responses, and enables accurate interpretation of patient data.

Figure 4.

Comparison of Tuba-seq identified GSTRs with cell line and clinical data. A, Kaplan–Meier curve (with 95% confidence interval in shading) of time-to-next treatment (months) for patients with or without KEAP1 mutations with metastatic oncogenic KRAS-driven lung adenocarcinoma to platinum-containing chemotherapy. The number of patients in each group is shown. P values were calculated from the Mantel–Haenszel test. B, Responses of patients with metastatic oncogenic KRAS-driven lung adenocarcinoma to platinum-containing chemotherapy are consistent with KEAP1 inactivation leading to resistance. KEAP1 mutations are significantly correlated with a higher hazard ratio for time-to-next treatment. C, Correlation between GSTR estimated in our study and that from the GDSC database based on cancer cell line studies. The significant GSTRs in our study are highlighted in red. D, Comparison of our identified GSTRs with data from the GDSC database. Stars represent significant cases.

Figure 4.

Comparison of Tuba-seq identified GSTRs with cell line and clinical data. A, Kaplan–Meier curve (with 95% confidence interval in shading) of time-to-next treatment (months) for patients with or without KEAP1 mutations with metastatic oncogenic KRAS-driven lung adenocarcinoma to platinum-containing chemotherapy. The number of patients in each group is shown. P values were calculated from the Mantel–Haenszel test. B, Responses of patients with metastatic oncogenic KRAS-driven lung adenocarcinoma to platinum-containing chemotherapy are consistent with KEAP1 inactivation leading to resistance. KEAP1 mutations are significantly correlated with a higher hazard ratio for time-to-next treatment. C, Correlation between GSTR estimated in our study and that from the GDSC database based on cancer cell line studies. The significant GSTRs in our study are highlighted in red. D, Comparison of our identified GSTRs with data from the GDSC database. Stars represent significant cases.

Close modal

### Comparisons with the cell line and PDX data

Although the positive and negative predictive values of cancer cell line studies are often questioned (50), the scale at which these in vitro studies can be performed has enabled the generation of drug response data across large panels of cell lines (11, 51, 52). Our study constitutes the largest in vivo survey of GSTRs; thus, we compared our findings to a study of cell line–therapeutic responses (GDSC; ref. 5) in which all five of our monotherapies were assessed (paclitaxel, palbociclib, phenformin, everolimus/rapamycin, and trametinib; ref. 5). Among the genotype–treatment pairs assessed in both studies, nine had significant effects in our analysis, but only one of these genotype–treatment pairs was significant in GDSC (RB1-palbociclib; Fig. 4C and D; Supplementary Fig. S14A and S14B). Note that in general we would not expect excellent agreement between our results and the cell line studies, given the lack of the autochthonous environment as well as the complexity of genetic backgrounds and mutation load in cell lines (50, 53).

The PRISM/DepMap compound screen has also quantified genotype-specific treatment responses (52). We tested whether mutation of each tumor suppressor gene is associated with a better or worse response for each genotype–treatment pair (the Mann–Whitney U test with FDR correction). The log viability measured by PRISM/DepMap compound screen and ScoreGSTR predicted by Tuba-seq were significantly correlated (ρ = 0.34, P = 0.01). Among the 9 significant genotype–treatment pairs predicted by Tuba-seq, 7 are in the same direction in the PRISM/DepMap compound screen data set, although only three of these effects were significant in PRISM (Supplementary Table S2). This is likely driven in part by their small sample size and the fact that in the PRISM/DepMap compound screen data set, the results are correlative and ignore all co-occurring mutations, whereas our analysis establishes a direct causal relationship.

PDX models can also be used to test for the association between genotype and drug response. Gao and colleagues conducted a very broad PDX study, generating a total of 4,759 response curves from ∼1,000 PDXs treated with 62 treatments (Supplementary Table S3; ref. 14). We used two-way ANOVA to determine whether there are any significant genotype–treatment pairs in these PDX data where the therapies overlap with our Tuba-seq results. Overall, there was no significant correlation between our ScoreGSTR and these ANOVA results (r = 0.124, P = 0.623; ρ = 0.07, P = 0.792; Supplementary Table S4). Given the large number of mutations per PDX (642 on average for the cancers used for comparison) and the small number of response curves measured per gene–drug pair (median number of treated PDXs that have the gene of interest mutated was 6; see Supplementary Table S3), the lack of correlation is not surprising. This PDX study, despite its extremely large scope, failed to identify the positive control genotype–treatment pair of RB1-mutated tumors being resistant to CDK4/6 inhibitors. These PDX results also did not uncover that KEAP1 inactivation leads to resistance to chemotherapy, which is an interaction that has been confirmed with clinical data (Fig. 4A; ref. 44).

Here we described and validated a scalable and quantitative in vivo pharmacogenomic preclinical model, which has high power to identify genotype–treatment responses using modest-size cohorts of mice. Although the number of mice required is modest, the total number of assayed tumors is large—on the order of thousands per mouse—providing the ability to assay a large number of tumor suppressors in the same experiment at a reasonable cost (Supplementary Fig. S15A–S15C). Indeed, although genetically engineered mouse models are key preclinical models to study genotype-specific treatment responses, traditional approaches are neither rigorously quantitative nor scalable, requiring impractically large numbers of mice. For instance, we estimated that with ten mice per group, the sensitivity of our approach would be > 99% to detect a genotype-specific treatment resistance that results in tumor sizes that are 50% larger than control tumors. If we had used a more traditional approach of comparing four cohorts of mice (with and without a specific tumor suppressor alteration and therapy-treated versus vehicle-treated), ∼300 mice/group would be required to achieve the same sensitivity for just one tumor suppressor genotype (Supplementary Fig. S15A–S15C). To build the pharmacogenomic map presented in this study, we would have needed to breed, initiate tumors in, and treat ∼10,000 mice instead of 58; thus, our system represents a >100-fold increase in throughput. Moreover, our power to detect effects is mostly limited by the number of mice per group and not by the number of tumors per mouse, allowing future iterations of this approach to query more genotypes per mouse.

We used one sgRNA per gene for the screening, and one may be concerned with the efficiency and off-target effects of the sgRNA. However, these sgRNAs have been extensively validated by previous studies. The ruggedness of the pharmacogenomic landscape further suggests the efficiency of the sgRNAs, with 7 of the 11 sgRNAs showing some genotype-specific treatment responses. Moreover, our pipeline is largely immune to off-target effects for sgRNAs, and such effects would not be expected to generate GSTRs (Supplementary Fig. S3 and Supplementary Methods and Discussion). Furthermore, neither differences in tumor number nor overall tumor burden across mice dramatically shifts tumor suppressive effects, suggesting that this method is not dramatically influenced by mouse-to-mouse variation (Supplementary Fig. S16A and S16B; ref. 29).

Our method is not only scalable and quantitative, but also allows the introduction of specific alterations into each tumor and the study of marginal effects of individual tumor suppressor genes in isolation, which is not possible using traditional cell lines or PDX approaches. Moreover, the use of genetically engineered mouse models allows autochthonous tumors to develop in their natural immunocompetent environment. This provides not only the ability to study immunotherapies but also the ability to recapitulate aspects of chemotherapy and targeted therapy responses that are influenced by adaptive immune responses.

The key result of this study, which had been suspected but never directly demonstrated, is that tumor suppressor genotype has a substantial impact on responses to a range of distinct therapies. The fact that this was not previously demonstrated experimentally is primarily due to the lack of appropriate systems, which underscores the need for higher-throughput quantitative preclinical models (27). Indeed, although databases like The Cancer Genome Atlas and GENIE databases provide valuable information on the mutational spectra in tumors, these databases generally lack treatment histories and cannot be used to study pharmacogenomic interactions. Prior cell line studies suggested that very few genotypes significantly affect drug responses (e.g., 0.24% of genotype–treatment pairs in GDSC), which we believe is largely due to the lack of statistical power. In contrast, we show that >20% of genotype–treatment pairs show interactions, suggesting a complex pharmacogenomic map of resistance and sensitivity of KRAS-driven lung adenocarcinoma.

There are some potential caveats for our PGx-Tuba-seq approach. We can only introduce a limited number of mutations into each tumor, reducing our ability to recapitulate the high tumor mutation burden and overall complexity of human tumors. Although we can study the genetic interaction among up to three genes, is it possible that even higher order interactions could modify the pharmacogenomic landscape. Furthermore, the extent to which our results recapitulate responses in patients remains unknown due to the lack of large-scale patient data sets. Thus, the interpretation of our results will benefit from further experimental, bioinformatic, and clinical evidence.

The complexity and rugged nature of this pharmacogenomic map has important implications for precision medicine. The complexity of human cancer genomics and the large number of potential therapies suggest that large-scale investigation of the pharmacogenomic maps in preclinical models will aid in patient selection. Our framework for in vivo functional genomic studies should easily allow a larger number of genes and additional monotherapies and combination therapies to be tested. Application to other genomic subtypes of lung cancer and potentially to other cancer types should further increase our knowledge of the pharmacogenomic determinants of therapy responses (54). We anticipate that the use of this platform to quantify the effects of additional therapies across a greater diversity of cancer genotypes will provide a cause-and-effect pharmacogenomic understanding from which novel biological hypotheses and precision treatment approaches will emerge.

H. Cai reports grants from the University of California during the conduct of the study. I.P. Winters reports grants from NIH during the conduct of the study; personal fees and other support from D2G Oncology, Inc. outside the submitted work; in addition, I.P. Winters has a patent for USPN 10,801,021 issued, licensed, and with royalties paid from D2G Oncology, Inc., a patent for USPN 10,738,300 issued, licensed, and with royalties paid from D2G Oncology, Inc., a patent for USSN 16/909,806 pending, a patent for USSN 17/017,920 pending, a patent for USSN 17/281,919 pending, a patent for USSN 63/129,161 pending, a patent for USSN 63/113,769 pending, a patent for USSN 63/132,359 pending, a patent for USSN 63/122,928 pending, a patent for USSN 63/144,866 pending, a patent for USSN 63/170,402 pending, a patent for USSN 63/173,976 pending, a patent for USSN 63/185,785 pending, and a patent for USSN 63/194,044 pending. C.M. Rudin reports personal fees from AbbVie, Amgen, AstraZeneca, Epizyme, Genentech/Roche, Ipsen, Jazz, Lilly, Syros, Bridge Medicines, Harpoon Therapeutics, and Earli outside the submitted work. D.A. Petrov reports other support from D2G Oncology outside the submitted work; in addition, D.A. Petrov is a co-founder and advisor for D2G Oncology Inc., and has a patent for compositions and methods for multiplexed quantitative analysis of cell lineages issued and licensed to D2G Oncology. M.M. Winslow reports grants from NIH during the conduct of the study; other support from D2G Oncology, Inc. outside the submitted work; in addition, M.M. Winslow has a patent for Stanford pending, issued, licensed, and with royalties paid from D2G Oncology, Inc.; and is a co-founder and advisor for D2G Oncology Inc. No disclosures were reported by the other authors.

C. Li: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. W.-Y. Lin: Conceptualization, resources, data curation, formal analysis, visualization, methodology, writing–original draft, writing–review and editing. H. Rizvi: Resources, formal analysis, validation, methodology, writing–original draft, writing–review and editing. H. Cai: Resources, data curation, formal analysis, validation, investigation, methodology, writing–review and editing. C.D. McFarland: Conceptualization, formal analysis, investigation, methodology, writing–review and editing. Z.N. Rogers: Conceptualization, methodology, writing–review and editing. M. Yousefi: Resources, data curation, investigation, methodology, writing–review and editing. I.P. Winters: Conceptualization, methodology, writing–review and editing. C.M. Rudin: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, methodology, writing–original draft, project administration, writing–review and editing. D.A. Petrov: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. M.M. Winslow: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing.

The authors thank L. Andrejka for technical support; A. Orantes and S. Mello for administrative support; members of the Stanford Animal Care staff for excellent animal care; D. Feldser, J. Sage, J. Zhang, and members of the Winslow and Petrov laboratories for helpful comments. C. Li was the Connie and Bob Lurie Fellow of the Damon Runyon Cancer Research Foundation (DRG-2331). W-Y. Lin was supported by an American Association of Cancer Research Postdoctoral fellowship (17-40-18-LIN). H. Rizvi was supported by the Druckenmiller Center for Lung Cancer Research at Memorial Sloan Kettering. H. Cai was supported by a Tobacco-Related Disease Research Program Postdoctoral Fellowship (28FT-0019). C.D. McFarland was supported by NIH K99-CA226506. M. Yousefi was supported by a Stanford Dean's fellowship, an American Lung Association senior research training grant, and NIH F32-CA236311. I.P. Winters was supported by the NSF Graduate Research Fellowship Program and NIH F31-CA210627. This work was supported by NIH R01-CA207133 (to M.M. Winslow and D.A. Petrov), NIH R01-CA231253 (to M.M. Winslow and D.A. Petrov), NIH R01-CA234349 (to M.M. Winslow and D.A. Petrov) and NIH U01-CA199215 (C.M. Rudin).

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.

1.
Miller
KD
,
Nogueira
L
,
Mariotto
AB
,
Rowland
JH
,
Yabroff
KR
,
Alfano
CM
, et al
Cancer treatment and survivorship statistics, 2019
.
CA Cancer J Clin
2019
;
69
:
363
85
.
2.
Denisenko
TV
,
Budkevich
IN
,
Zhivotovsky
B
.
Cell death-based treatment of lung adenocarcinoma
.
Cell Death Dis
2018
;
9
:
117
.
3.
Dagogo-Jack
I
,
Shaw
AT
.
Tumour heterogeneity and resistance to cancer therapies
.
Nat Rev Clin Oncol
2018
;
15
:
81
94
.
4.
Bedard
PL
,
Hansen
AR
,
Ratain
MJ
,
Siu
LL
.
Tumour heterogeneity in the clinic
.
Nature
2013
;
501
:
355
64
.
5.
Iorio
F
,
Knijnenburg
TA
,
Vis
DJ
,
Bignell
GR
,
Menden
MP
,
Schubert
M
, et al
A landscape of pharmacogenomic interactions in cancer
.
Cell
2016
;
166
:
740
54
.
6.
Wang
Y
,
Wang
Z
,
Xu
J
,
Li
J
,
Li
S
,
Zhang
M
, et al
Systematic identification of non-coding pharmacogenomic landscape in cancer
.
Nat Commun
2018
;
9
:
3192
.
7.
Qiu
Z
,
Li
H
,
Zhang
Z
,
Zhu
Z
,
He
S
,
Wang
X
, et al
A pharmacogenomic landscape in human liver cancers
.
Cancer Cell
2019
;
36
:
179
93
.
8.
Wood
DE
,
Kazerooni
EA
,
Baum
SL
,
Eapen
GA
,
Ettinger
DS
,
Hou
L
, et al
Lung cancer screening, version 3.2018, NCCN Clinical Practice Guidelines in Oncology
.
J Natl Compr Canc Netw
2018
;
16
:
412
41
.
9.
Ettinger
DS
,
Wood
DE
,
Aisner
DL
,
Akerley
W
,
Bauman
J
,
Chirieac
LR
, et al
Non-small cell lung cancer, version 5.2017, NCCN clinical practice guidelines in oncology
.
J Natl Compr Canc Netw
2017
;
15
:
504
35
.
10.
Punt
CJ
,
Koopman
M
,
Vermeulen
L
.
From tumour heterogeneity to advances in precision treatment of colorectal cancer
.
Nat Rev Clin Oncol
2017
;
14
:
235
46
.
11.
Barretina
J
,
Caponigro
G
,
Stransky
N
,
Venkatesan
K
,
Margolin
AA
,
Kim
S
, et al
The Cancer Cell Line Encyclopedia enables predictive modelling of anticancer drug sensitivity
.
Nature
2012
;
483
:
603
7
.
12.
Lee
JK
,
Liu
Z
,
Sa
JK
,
Shin
S
,
Wang
J
,
Bordyuh
M
, et al
Pharmacogenomic landscape of patient-derived tumor cells informs precision oncology therapy
.
Nat Genet
2018
;
50
:
1399
411
.
13.
Roper
N
,
Stensland
KD
,
Hendricks
R
,
Galsky
MD
.
The landscape of precision cancer medicine clinical trials in the United States
.
Cancer Treat Rev
2015
;
41
:
385
90
.
14.
Gao
H
,
Korn
JM
,
Ferretti
S
,
Monahan
JE
,
Wang
Y
,
Singh
M
, et al
High-throughput screening using patient-derived tumor xenografts to predict clinical trial drug response
.
Nat Med
2015
;
21
:
1318
25
.
15.
Ben-David
U
,
Siranosian
B
,
Ha
G
,
Tang
H
,
Oren
Y
,
Hinohara
K
, et al
Genetic and transcriptional evolution alters cancer cell line drug response
.
Nature
2018
;
560
:
325
30
.
16.
Hidalgo
M
,
Amant
F
,
Biankin
AV
,
Budinska
E
,
Byrne
AT
,
Caldas
C
, et al
Patient-derived xenograft models: an emerging platform for translational cancer research
.
Cancer Discov
2014
;
4
:
998
1013
.
17.
Ben-David
U
,
Ha
G
,
Tseng
YY
,
Greenwald
NF
,
Oh
C
,
Shih
J
, et al
Patient-derived xenografts undergo mouse-specific tumor evolution
.
Nat Genet
2017
;
49
:
1567
75
.
18.
Hoffman
RM
.
Patient-derived orthotopic xenografts: better mimic of metastasis than subcutaneous xenografts
.
Nat Rev Cancer
2015
;
15
:
451
2
.
19.
Middleton
G
,
Fletcher
P
,
Popat
S
,
Savage
J
,
Summers
Y
,
Greystoke
A
, et al
The National Lung Matrix Trial of personalized therapy in lung cancer
.
Nature
2020
;
583
:
807
12
.
20.
Kersten
K
,
de Visser
KE
,
van Miltenburg
MH
,
Jonkers
J
.
Genetically engineered mouse models in oncology research and cancer medicine
.
EMBO Mol Med
2017
;
9
:
137
53
.
21.
Sharpless
NE
,
Depinho
RA
.
The mighty mouse: genetically engineered mouse models in cancer drug development
.
Nat Rev Drug Discov
2006
;
5
:
741
54
.
22.
Singh
M
,
Murriel
CL
,
Johnson
L
.
Genetically engineered mouse models: closing the gap between preclinical data and trial outcomes
.
Cancer Res
2012
;
72
:
2695
700
.
23.
Chen
Z
,
Cheng
K
,
Walton
Z
,
Wang
Y
,
Ebi
H
,
Shimamura
T
, et al
A murine lung cancer co-clinical trial identifies genetic modifiers of therapeutic response
.
Nature
2012
;
483
:
613
7
.
24.
Schmitt
A
,
Knittel
G
,
Welcker
D
,
Yang
T-P
,
George
J
,
Nowak
M
, et al
ATM deficiency is associated with sensitivity to PARP1-and ATR inhibitors in lung adenocarcinoma
.
Cancer Res
2017
;
77
:
3040
56
.
25.
Shackelford
DB
,
Abt
E
,
Gerken
L
,
Vasquez
DS
,
Seki
A
,
Leblanc
M
, et al
LKB1 inactivation dictates therapeutic response of non-small cell lung cancer to the metabolism drug phenformin
.
Cancer Cell
2013
;
23
:
143
58
.
26.
Rogers
ZN
,
McFarland
CD
,
Winters
IP
,
Naranjo
S
,
Chuang
CH
,
Petrov
D
, et al
A quantitative and multiplexed approach to uncover the fitness landscape of tumor suppression in vivo
.
Nat Methods
2017
;
14
:
737
42
.
27.
Winters
IP
,
Murray
CW
,
Winslow
MM
.
Towards quantitative and multiplexed in vivo functional cancer genomics
.
Nat Rev Genet
2018
;
19
:
741
55
.
28.
Rogers
ZN
,
McFarland
CD
,
Winters
IP
,
Seoane
JA
,
JJ
,
Yoon
S
, et al
Mapping the in vivo fitness landscape of lung adenocarcinoma tumor suppression in mice
.
Nat Genet
2018
;
50
:
483
6
.
29.
Cai
H
,
Chew
SK
,
Li
C
,
Tsai
MK
,
Andrejka
L
,
Murray
CW
, et al
A functional taxonomy of tumor suppression in oncogenic KRAS-driven lung cancer
.
Cancer Discov
2021
. doi: 10.1158/2159-8290.
30.
Murray
CW
,
JJ
,
Tsai
MK
,
Li
C
,
Winters
IP
,
Tang
R
, et al
An LKB1-SIK axis suppresses lung tumor growth and controls differentiation
.
Cancer Discov
2019
;
9
:
1590
605
.
31.
Jackson
EL
,
Willis
N
,
Mercer
K
,
Bronson
RT
,
Crowley
D
,
Montoya
R
, et al
Analysis of lung tumor initiation and progression using conditional expression of oncogenic K-ras
.
Genes Dev
2001
;
15
:
3243
8
.
32.
Chiou
SH
,
Winters
IP
,
Wang
J
,
Naranjo
S
,
Dudgeon
C
,
Tamburini
FB
, et al
Pancreatic cancer modeling using retrograde viral vector delivery and in vivo CRISPR/Cas9-mediated somatic genome editing
.
Genes Dev
2015
;
29
:
1576
85
.
33.
L
,
Zwingman
TA
,
Sunkin
SM
,
Oh
SW
,
Zariwala
HA
,
Gu
H
, et al
A robust and high-throughput Cre reporting and characterization system for the whole mouse brain
.
Nat Neurosci
2010
;
13
:
133
40
.
34.
Pezza
JA
,
Kucera
R
,
Sun
LJNEB
.
Polymerase fidelity: what is it, and what does it mean for your PCR
.
New England Biolabs
;
2014
.
35.
MacConaill
LE
,
Burns
RT
,
Nag
A
,
Coleman
HA
,
Slevin
MK
,
Giorda
K
, et al
Unique, dual-indexed sequencing adapters with UMIs effectively eliminate index cross-talk and significantly improve sensitivity of massively parallel sequencing
.
BMC Genomics
2018
;
19
:
30
.
36.
Cheng
DT
,
Mitchell
TN
,
Zehir
A
,
Shah
RH
,
Benayed
R
,
Syed
A
, et al
Memorial Sloan Kettering-integrated mutation profiling of actionable cancer targets (MSK-IMPACT): a hybridization capture-based next-generation sequencing clinical assay for solid tumor molecular oncology
.
J Mol Diagn
2015
;
17
:
251
64
.
37.
Hill
AJ
,
McFaline-Figueroa
JL
,
Starita
LM
,
Gasperini
MJ
,
Matreyek
KA
,
Packer
J
, et al
On the design of CRISPR-based single-cell molecular screens
.
Nat Methods
2018
;
15
:
271
4
.
38.
Bollard
J
,
Miguela
V
,
Ruiz de Galarreta
M
,
Venkatesh
A
,
Bian
CB
,
Roberto
MP
, et al
Palbociclib (PD-0332991), a selective CDK4/6 inhibitor, restricts tumour growth in preclinical models of hepatocellular carcinoma
.
Gut
2017
;
66
:
1286
96
.
39.
Knudsen
ES
,
Pruitt
SC
,
Hershberger
PA
,
Witkiewicz
AK
,
Goodrich
DW
.
Cell cycle and beyond: exploiting new RB1 controlled mechanisms for cancer therapy
.
Trends Cancer
2019
;
5
:
308
24
.
40.
O'Leary
B
,
Cutts
RJ
,
Liu
Y
,
Hrebien
S
,
Huang
X
,
Fenwick
K
, et al
The genetic landscape and clonal evolution of breast cancer resistance to palbociclib plus fulvestrant in the PALOMA-3 trial
.
Cancer Discov
2018
;
8
:
1390
403
.
41.
Hong
EP
,
Park
JW
.
Sample size and statistical power calculation in genetic association studies
.
Genomics Inform
2012
;
10
:
117
22
.
42.
Wittes
J
.
Sample size calculations for randomized controlled trials
.
Epidemiol Rev
2002
;
24
:
39
53
.
43.
Xiao
P
,
Sun
LL
,
Wang
J
,
Han
RL
,
Ma
Q
,
Zhong
DS
.
LKB1 gene inactivation does not sensitize non-small cell lung cancer cells to mTOR inhibitors in vitro
.
Acta Pharmacol Sin
2015
;
36
:
1107
12
.
44.
Arbour
KC
,
Jordan
E
,
Kim
HR
,
Dienstag
J
,
Yu
HA
,
Sanchez-Vega
F
, et al
Effects of co-occurring genomic alterations on outcomes in patients with KRAS-mutant non-small cell lung cancer
.
Clin Cancer Res
2018
;
24
:
334
40
.
45.
Zhang
P
,
Singh
A
,
Yegnasubramanian
S
,
Esopi
D
,
Kombairaju
P
,
Bodas
M
, et al
Loss of Kelch-like ECH-associated protein 1 function in prostate cancer cells causes chemoresistance and radioresistance and promotes tumor growth
.
Mol Cancer Ther
2010
;
9
:
336
46
.
46.
Tian
Y
,
Wu
K
,
Liu
Q
,
Han
N
,
Zhang
L
,
Chu
Q
, et al
Modification of platinum sensitivity by KEAP1/NRF2 signals in non-small cell lung cancer
.
J Hematol Oncol
2016
;
9
:
83
.
47.
Galan-Cobo
A
,
Sitthideatphaiboon
P
,
Qu
X
,
Poteete
A
,
Pisegna
MA
,
Tong
P
, et al
LKB1 and KEAP1/NRF2 pathways cooperatively promote metabolic reprogramming with enhanced glutamine dependence in KRAS-mutant lung adenocarcinoma
.
Cancer Res
2019
;
79
:
3251
67
.
48.
Papillon-Cavanagh
S
,
Doshi
P
,
Dobrin
R
,
Szustakowski
J
,
Walsh
AM
.
STK11 and KEAP1 mutations as prognostic biomarkers in an observational real-world lung adenocarcinoma cohort
.
ESMO Open
2020
;
5
:
e000706
.
49.
Jeong
Y
,
Hellyer
JA
,
Stehr
H
,
Hoang
NT
,
Niu
X
,
Das
M
, et al
Role of KEAP1/NFE2L2 mutations in the chemotherapeutic response of patients with non-small cell lung cancer
.
Clin Cancer Res
2020
;
26
:
274
81
.
50.
Haibe-Kains
B
,
El-Hachem
N
,
Birkbak
NJ
,
Jin
AC
,
Beck
AH
,
Aerts
HJ
, et al
Inconsistency in large pharmacogenomic studies
.
Nature
2013
;
504
:
389
93
.
51.
Garnett
MJ
,
Edelman
EJ
,
Heidorn
SJ
,
Greenman
CD
,
Dastur
A
,
Lau
KW
, et al
Systematic identification of genomic markers of drug sensitivity in cancer cells
.
Nature
2012
;
483
:
570
5
.
52.
Corsello
SM
,
Nagari
RT
,
Spangler
RD
,
Rossen
J
,
Kocak
M
,
Bryan
JG
, et al
Discovering the anticancer potential of non-oncology drugs by systematic viability profiling
.
Nat Cancer
2020
;
1
:
235
48
.
53.
M
,
DiStefano
MS
,
Riesel
EA
,
Carrasco
E
,
Holzapfel
HY
,
Jones
DC
, et al
Drug response consistency in CCLE and CGP
.
Nature
2016
;
540
:
E9
E10
.
54.
Foggetti
G
,
Li
C
,
Cai
H
,
Hellyer
JA
,
Lin
WY
,
Ayeni
D
, et al
Genetic determinants of EGFR-driven lung cancer growth and therapeutic response in vivo
.
Cancer Discov
2021
;
11
:
1736
53
.