High-throughput Phenotyping of Lung Cancer Somatic Mutations

Abstract

Recent genome sequencing efforts have identified millions of somatic mutations in cancer. However, the functional impact of most variants is poorly understood. Here we characterize 194 somatic mutations identified in primary lung adenocarcinomas. We present an expression-based variant-impact phenotyping (eVIP) method that uses gene expression changes to distinguish impactful from neutral somatic mutations. eVIP identified 69% of mutations analyzed as impactful and 31% as functionally neutral. A subset of the impactful mutations induces xenograft tumor formation in mice and/or confers resistance to cellular EGFR inhibition. Among these impactful variants are rare somatic, clinically actionable variants including EGFR S645C, ARAF S214C and S214F, ERBB2 S418T, and multiple BRAF variants, demonstrating that rare mutations can be functionally important in cancer.

Figure 1. Expression-Based Variant Impact Phenotyping

(A) Overview of the experimental pipeline from reagent generation to eVIP analysis. (B) Dot-plot and heat map representation of replicate consistency (WT vs. WT or variant vs. variant) comparisons and signature identity (WT vs. variant) comparisons. Correlation is measured by a weighted connectivity score (wtcs). *, adjusted p < 0.05. n.s., adjusted p > 0.05. (C) Schematic of the decision tree-based eVIP algorithm. The first test (“impact test”) outputs a single Bonferroni-adjusted p value indicating likelihood of mutation impact. For impactful mutations, the next tests are used to determine the directional impact of the mutations. For mutations found to be non-impactful, a “power test” assesses whether the two signatures are similar to one another due to real signal, or due to noise. See also Figure S1, Table S1, and File S1.

Figure 2. eVIP Classifications for 110 Lung Cancer Somatic Mutations

(A) A “sparkler” plot representation of eVIP predictions. Each variant allele with an eVIP prediction is represented as a point. The x-axis indicates the –log₁₀ (adjusted p) of a Kruskal-Wallis test comparing wild-type and mutant ORF replicate consistency and signature identity. The y-axis is the “impact direction score,” the absolute value of which is equal to the –log₁₀ (adjusted p) of a Wilcoxon test directly comparing wild-type and mutant ORF replicate consistency. The sign of the impact direction score is positive if the mutant ORF replicate consistency is greater than the wild-type replicate consistency and negative if the mutant ORF replicate consistency is less than the wild-type ORF replicate consistency. The line connecting each point and the graph origin is drawn to emphasize that longer distance from the origin implies more confidence in the prediction. (B) Enrichment of impactful variants in genes found to be significantly mutated cancer genome studies. GOF, COF and LOF predictions are all considered impactful. ***, p < 0.0001 by Fisher's Exact test. (C) Distribution of eVIP calls in known oncogenes, known tumor suppressor genes (TSGs), or genes of unknown function. (D-F) Gene-specific sparkler plots for known oncogenes (D), known tumor suppressor genes (E), and genes of unknown function (F). The tested variant allele is labeled and colored based on the eVIP prediction. Coloring is as in (A). See also Table S2, Table S3, Table S4, and File S2.

Figure 3. Hypomorphic and Dominant-negative KEAP1 Variants Identified by Gene Expression Profiling

(A) Hierarchical clustering of KEAP1 expression signatures in A549 cells. The similarity matrix was computed using the weighted connectivity score (wtcs) as the similarity metric. (B) Two-class comparison of NRF2 ORF signatures versus control (EGFP, HcRed, Luciferase) signatures across three KEAP1 wild-type cell lines (AALE, SALE, H1299). The top transcripts up- or down-regulated by NFE2L2 expression were determined by a signal-to-noise statistic. (C) Expression of the direct NRF2 transcriptional target TXNRD1 is shown as a biosensor of NRF2 and KEAP1 activity. Upper panels show NRF2 variants and KEAP1 wild-type and neutral variants. The lower panels show three classes of loss-of-function KEAP1 variants identified by hierarchical clustering of expression profiles generated in four cellular contexts: A549 (KEAP1 mutant) and three KEAP1 wild-type cell lines, H1299, SALE, and AALE. See also Figure S2 and File S3.

Figure 4. Orthogonal Assays in Different Cellular Contexts Identify Oncogenic EGFR Pathway Mutations

(A) Cell viability of PC9 cells after mutant allele library infection and 72 hours of treatment with 3 μM erlotinib or DMSO. Data shown are the average robust Z scores of two replicates per ORF condition. A dashed line indicates the threshold used to select ORFs for validation (Z > 2). (B) Relationship between mutation frequency in COSMIC (x-axis) and ability to rescue cell viability in erlotinib (y-axis). Two dashed lines indicate the Z score thresholds used to select ORFs for validation (Z > 2) and the threshold at which all ORFs retested in validation (Z > 3). See also Figure S3 and Table S5.

Figure 5. A Multiplex In Vivo Tumor Formation Screen for Identification of Activated Oncogenes

(A) Left, experimental schematic showing screening workflow. Right, pie charts showing the median corrected reads per million (RPM) per ORF in pre-injection and tumor samples for all pools in the tumor experiment. Because ORFs were assayed in different pools, the proportion of each ORF on the pie chart may not reflect the actual relative levels of oncogenic activity. (B) One-sided volcano plot showing distribution of ORFs from all pools. Each datapoint represents data generated from all pre-injection and tumor replicates for a given ORF. The log2 fold-change (x-axis) was calculated by comparing the median corrected RPM of each ORF from the pre-injection samples to the median corrected RPM in tumor samples. (C) Plot showing relationship between ORF variants across functional and expression-based screens. All alleles with both tumor formation and EGFR epitasis data were plotted; alleles not analyzed by eVIP are indicated by “ND” (not determined). The colored shapes indicate the predicted mutation impact as assessed by eVIP in A549 cells. See also Figure S4, Table S6, and Table S7.

Figure 6. Integration of Expression Signatures and Functional Data Identifies Rare Activating Mutations in the RAS/MAPK pathway

(A) Heat map showing hierarchical clustering of expression signatures in A549 cells. Similarity of signatures was compared using the weighted connectivity score (wtcs) and transformed to a rankpoint distribution with the best correlation set to 1 and best anti-correlation set to −1. Asterisks indicate engineered kinase-dead ARAF variants (not found in human cancer). (B) Western blot of A549 cells expressing wild-type BRAF or BRAF variants, or control vector (HcRed). The primary antibodies used are indicated on the left. p-ERK, phosphorylated Thr202/Thr204 ERK1/2. (C) Hierarchical clustering of a similarity matrix of all A549 expression signatures in the study. Similarity of signatures was computed as in (A). A bar chart above the heatmap indicates the average robust Z score achieved by each respective ORF in the PC9 EGFR epistasis screen. (D) Expression signature analysis in A549 cells showing all ORFs ranked from left to right in order of correlation to the KRAS G12C signature. Canonical alleles of EGFR/RAS pathway oncogenes and rare alleles induce signatures highly similar to KRAS G12C (left, gray shading), with the exception of ARAF V145L, determined to be neutral by eVIP, and the three apparently inactive BRAF mutants described in the text (right side, gray text). An additional track shows allele performance in the EGFR epistasis assay (black bars).

Figure 7. Sensitivity of Rare Oncogenic Mutations to MEK Inhibition

(A) Viability of PC9 stable isogenic cell lines after treatment with the small molecules shown for 96 hours. 36-point dose response curves were performed and curves plotted in GraphPad Prism. Data shown are the ratio of the area-under-the-curve (AUC) of the mutant PC9 line vs. the PC9-luciferase negative control. (B) Viability of PC9 stable isogenic cell lines after treatment with the small molecules shown in combination with erlotinib. Inhibitors were delivered in a 1:1 molar ratio across a 36-point dose-response curve. Data shown are the ratio of the AUC of the mutant PC9 line in the combination treatment vs. the AUC of the mutant PC9 in erlotinib only. Note that EGFR K754E appears to be the least sensitive to trametinib inhibition but this is actually due to EGFR K754E conferring the least resistance to trametinib, as shown in panel (C). (C) Dose-response curves of PC9 stable isogenic cell lines expressing EGFR, ERBB2, or KRAS variants. Data shown are the same curves used to generate the heat maps in panel (B). See also Figure S5.

Figure 8. Summary of Validated Rare Functional Alleles of Oncogenes

Stick plot representation of predicted protein sequence and domain structure of oncogenes harboring rare, functional alleles. Canonical hotspot mutations are shown for reference (gray). Rare variants shown in red scored in eVIP, the EGFR epistasis screen, and the tumor formation screen, and additionally are sensitive to MEK inhibition with trametinib.