Distinct patterns of somatic genome alterations in lung adenocarcinomas and squamous cell carcinomas

Comparison of somatically altered genes

In order to compare the somatic profiles of lung ADC and lung SqCC and to identify novel genetic alterations, we studied 660 lung adenocarcinoma/normal paired exome sequences (including 274 previously unpublished cases, 227 previously described from The Cancer Genome Atlas (TCGA)⁶, and 159 cases from the Imielinski et al⁸ cohort) and 484 lung SqCC/normal paired exome sequences (including 308 previously unpublished cases and 176 previously described from TCGA⁷; Supplementary Tables 1–4). Similarly to previous studies^6,7, we observed a median somatic mutation rate of 8.7/Mb and 9.7/Mb for lung ADCs and SqCCs, respectively. After excluding genes with lower median expression (6.16 log₂ FPKM for ADCs and 6.27 log₂ FPKM for SqCCs; Supplementary Fig. 1; Online Methods), we identified 38 genes as significantly mutated in ADC and 20 genes significantly mutated in SqCC using MutSig2CV¹⁰ (q-value < 0.1; Supplementary Table 5,6). Only 6 genes, TP53, RB1, ARID1A, CDKN2A, PIK3CA, and NF1, were significantly mutated in both tumor types and of these, the frequency of TP53, CDKN2A, and PIK3CA mutation was significantly higher in SqCC tumors (p < 0.01; Fisher’s exact test; Figure 1a). Likewise, only 11 out of 42 focal amplification peaks were identified as altered in both tumor types (Figure 1b), while 13 out of 50 focal deletion peaks were altered in both tumor types (Figure 1c). Interestingly, when compared to 19 other tumor types from TCGA¹⁰, the lists of significantly mutated genes for lung ADC and SqCC had a greater overlap with significantly mutated gene lists from other tumor types (>13% overlap; FDR q-value < 0.1) than to each other (12% overlap; p = 0.105; Supplementary Fig. 2), consistent with previous pan-cancer analyses¹¹. Recurrently mutated and amplified genes in lung SqCC most closely resembled the alterations in head and neck squamous cell carcinoma (HNSC) and bladder cancer (BLCA), two other epithelial cancers with epidemiological associations with smoking (>25% overlap; Supplementary Fig. 2). Within these overlapping genes, TP53, CDKN2A, and FAT1 are specifically enriched in HPV- HNSC¹². In contrast, the significantly mutated genes in lung ADC were most similar to glioblastoma (GBM) and colorectal cancer (CRC) (FDR q-value < 0.1). While lung ADC and lung SqCC did share several focal deletion peaks, five of these peaks were putative fragile sites (shown in green in Figure 1c). Taken together, these results suggest that the somatic drivers of carcinogenesis may be largely distinct between lung adenocarcinoma and lung squamous cell carcinoma.

Mutational signatures in lung cancer

Various carcinogenic and cancer-related processes contribute to the mutational patterns observed in tumors^13,14. Previous large-scale studies of lung cancer genomes have identified signatures associated with non-smoking and smoking cases^6,8,15; here we extend these findings based on the improved statistical power of our larger sample set. Using non-negative matrix factorization^13,16 (NMF, Online Methods), we identified six mutational signatures in this cohort, many of which are strongly correlated with previously defined signatures in the COSMIC database^13,17 (Supplementary Figs. 3–5; Supplementary Table 7). These included a UV-related signature of C>T at TpCpC or CpCpC (COSMIC Signature 7, abbreviated SI7), a smoking-related signature of C>A transversions (SI4), a mismatch repair (MMR) signature of C>T at GpCpG (SI15/SI6), two APOBEC-related signatures of C>G or C>T at TpCpT or TpCpA (SI13 and SI2), and a final signature with a moderate correlation to COSMIC signature 5 (SI5) with putative “molecular clock” properties¹⁸ (Supplementary Fig. 5). In addition to identifying mutational signatures, NMF also estimates the number of mutations contributed by each signature within each tumor. The estimated number of SI4 (i.e. smoking-related) mutations per Mb displayed a bimodal pattern in lung ADC but not in lung SqCC (Figure 2a). Furthermore, the rate of SI4 mutations per Mb was able to classify tumors into those from never vs. ever smokers substantially better in lung ADCs (AUC=0.87; Supplementary Fig. 6) than in lung SqCCs (AUC=0.62) suggesting that the smoking statuses for the 18 never smokers with lung SqCC may be inaccurate. 87% of lung ADCs from never smokers were categorized as transversion-low (TV-L; ≤0.696 of SI4 per Mb; p = 8.5 × 10⁻³⁷; Fisher’s exact test; Figure 2b; Supplementary Fig. 6). However, only 45% of transversion-low lung ADCs were from patients who were never smokers (Figure 2b). Within each tumor, we also derived the fraction of estimated mutations for a signature by dividing the number of estimated mutations for that signature by the sum of estimated mutations from all signatures. Lung SqCCs displayed significantly higher overall rates of SI5 mutations per Mb compared to all lung ADCs (p < 0.001; Wilcoxon rank-sum test). However, lung ADCs from never smokers displayed the highest fraction of estimated mutations from this signature on average (Figure 2c; Supplementary Fig. 7). In lung SqCC, we also observed moderate associations of tumor stage with SI5 activity and total mutation rate (p < 0.01; Supplementary Fig. 8).

The mutational profiles of three lung SqCCs (~1% of lung SqCCs) exhibited a pattern of UV-related mutations (SI7) commonly observed in melanoma and displayed a significantly higher mutation rate of somatic single nucleotide variants (SSNVs) and somatic dinucleotide polymorphisms (DNPs) compared to the other lung tumors (p < 0.01) but not higher rates of indels (p > 0.05; Figure 2d). One of these patients (TCGA-18–3409) had a previous history of basal cell carcinoma in the forehead raising the possibility that metastasis from the skin to the lung had occurred. The other two lung SqCCs with this signature may also represent squamous cell skin carcinomas metastatic to the lung. Mutational profiles for another 7 tumors (4 lung ADCs and 3 lung SqCCs) exhibited an MMR-like signature (SI15/SI6) commonly observed in microsatellite instable (MSI) colorectal carcinomas (Figure 2e)¹³. These tumors had significantly higher rates of both SSNVs and short indels (p < 0.00). They also displayed lower expression levels of the mismatch repair gene MLH1 (p = 0.011) suggesting a potential etiology for this signature in lung.

Novel significantly mutated genes

Comparing the significantly mutated genes to other tumor types from the TCGA Pan-Cancer study¹⁰ revealed that there were several genes significantly mutated exclusively in lung ADC including STK11, RBM10, KEAP1, RAF1, RIT1, and MET (MutSig2CV q-value < 0.1; Figure 3a; Supplementary Table 5). NFE2L2, KDM6A, RASA1, NOTCH1, and HRAS were significantly mutated in lung SqCC but not in other cancer types (excluding HNSC and BLCA) (Figure 3b; Supplementary Table 6). Genes that reached modest statistical significance in lung ADC that have been observed previously in lung cancer or in other tumor types include AKT1 with a recurrent mutation at p.E17K, CDK4 with a recurrent mutation at p.R24L, and DNMT3A (p < 0.005; Supplementary Table 5). Novel significantly mutated genes exclusive to lung ADC and which are absent in other tumor types include PPP3CA, which is the catalytic subunit for the calcium-dependent phosphatase, calcineurin. The mutations in PPP3CA clustered in the autoinhibitory domain near the C-terminus suggesting they may be gain-of-fuction alterations (Figure 4a). In addition, mutations in the autoinhibitory domain also tended to co-occur with activating KRAS mutations (p = 0.033) suggesting a potential relationship between K-Ras and calcineurin signaling pathways. Significantly mutated methyltransferase genes included MLL3 (KMT2C) and SETD2. A novel gene in this class was the H3K79 methyltransferase DOT1L, which was mutated in 3% of lung adenocarcinomas with enrichment for truncating mutations (Figure 4a). Recurrent mutations in lung adenocarcinoma have been previously reported in splicing factors such as U2AF1 and loss of function mutations in the RNA binding protein RBM10⁸. In the current dataset, a cap methyltransferase, FTSJD1 (also known as CMTR2), was significantly mutated and enriched for frame shift mutations (Figure 4a). We also examined genes for other known proteins in this class and found recurrent mutations in SF3B1¹⁹ and SNRPD3 (Supplementary Fig. 9). EGFR mutations were enriched in females, and SMARCA4 mutations were enriched in males (FDR q-value < 0.1; Supplementary Table 8). RBM10 mutations were modestly enriched in males as previously reported (q-value = 0.219)⁶. Novel significantly mutated genes in lung SqCC that were enriched for frame shift mutations (p < 0.001) included RASA1, whose protein product is p120GAP²⁰ (Figure 4b). CUL3, whose protein product is a known interaction partner of KEAP1, also reached statistical significance in the lung SqCC cohort²¹ (Figure 4b). RB1 mutations were enriched in females, whereas PASK mutations were exclusive to males (FDR q-value < 0.1; Supplementary Table 9). We did not observe significant associations between mutation status and patient survival or tumor stage after correction for multiple hypothesis testing (Supplementary Tables 10–13). Controlling for tumor stage did not reveal additional significant associations between mutation status and survival.

Previous studies have shown that joint analysis of different tumor types can yield additional statistical power to detect low-frequency events even if the tumor types are from vastly different tissues of origin and/or etiologies¹⁰. Additionally, although the individual drivers may be distinct between two tumor types, pathways such as MAP kinase are often altered similarly in both. We therefore hypothesized that combining the lung ADC and SqCC tumor cohorts (i.e. Pan-Lung) would reveal additional recurrent somatic pathway alterations common to both. We found 14 genes to be significantly mutated in the Pan-Lung cohort but not significantly mutated in either individual tumor type (q-value < 0.1; Supplementary Fig. 10; Supplementary Table 14). Many of these genes are involved in epigenetic regulation or immune-related pathways. KLF5, a transcription factor critical for lung development²² contained a novel recurrent mutation in the zinc finger domain, which was observed in both ADCs and SqCCs (Figure 4c). A regulator of KLF5, the E3 ubiquitin ligase FBXW7²³, was also significantly mutated in the lung SqCC and Pan-Lung cohorts but did not co-occur with KLF5 mutations. A super-enhancer duplication associated with increased KLF5 expression has also been recently reported in HNSC by our group²⁴. The paralogs EP300 and CREBBP had a mutational hotspot region within the histone acetyltransferase (HAT) domain. All missense mutations in the HAT domains and other loss-of-function mutations outside this domain were non-overlapping within these two proteins. For sites with sufficient sequencing depth in the RNA-seq (power > 95%), we observed a somatic SNV validation rate of 88%.

Novel somatic copy number alterations

With a larger sample size, we had better resolution to detect novel copy number changes and ascertain the putative target genes of focal amplifications and deletions. For some peaks that still contained many genes, we inferred the most likely target gene by examining the same peak in a Pan-Cancer copy number analysis across 11 tumor types²⁵ that included a subset of lung cancers from this set. The most significantly focally amplified genes in lung ADC were NKX2-1, MYC, TERT, MCL1, and MDM2 (Figure 5a; Supplementary Table 15), while peaks at SOX2, CCND1, WHSC1L1/FGFR1, MYC, and EGFR were among the top for lung SqCC (Figure 5b; Supplementary Table 16). Amplification peaks previously described in other tumor types but less characterized before in the lung tumors included KAT6A, ZNF217, and MYCL1 for lung ADC (Figure 5a) and IGF1R, KDM5A, PTP4A1/PHF3, and MYCL1 for lung SqCC (Figure 5b). CCND3 was specifically amplified in lung ADC while an amplification peak near MIR21/TUBD1 (Figure 5c) was also observed in breast cancer²⁵. MIR21 expression has been shown to be a prognostic factor for early stage ADC^26,27. Likewise, novel amplification peaks for lung SqCC included YES1, a Src family non-receptor protein kinase, and MIR205 (Figure 5d). Expression of MIR205 has been used to distinguish lung squamous cell carcinomas from other NSCLC types²⁸ suggesting that amplification of this microRNA may represent a lineage specific alteration similar to that of SOX2 amplification. Finally, combined Pan-Lung copy number analysis revealed additional amplification peaks around MAPK1 (Figure 5a–d; Supplementary Table 17).

Focal deletion peaks in lung ADC included the chromatin modifier genes SMARCA4 and ARID2 (Supplementary Fig. 11; Supplementary Table 18), which were also significantly mutated and enriched for loss-of-function mutations. Novel lung SqCC focal deletions observed in other tumor types included ZMYND11, CREBBP, ROBO1, USP22, and KDM6A (Supplementary Fig. 11; Supplementary Table 19). B2M (Beta2-microglobulin), a component of the MHC complex, was focally deleted in both tumor types, enriched for loss-of-function mutations in both tumor types (p < 0.01), and was significantly mutated in the Pan-Lung analysis (FDR q-value = 0.006). Combined Pan-Lung copy number analysis revealed another focal deletion peak around TRAF3 (Supplementary Table 20), which was also reported in HNSC¹². In general, mRNA expression was significantly associated with copy number levels for target genes (Supplementary Figs. 12 and 13). We did not observe substantial batch effects within or across tumor types in both the mRNA expression and CNV data (Supplementary Figs. 14 and 15).

Identifying Ras/Raf/RTK drivers in lung ADC

In lung ADC, mutually exclusive alterations have been characterized in members of the receptor tyrosine kinase (RTK)/Ras/Raf signaling pathways. These alterations are of particular interest because of the dramatic responses that have been observed in response to RTK inhibitors in clinical trials such as those for lung ADC patients harboring EGFR mutation or ALK or ROS1 translocations²⁹. However, many lung ADCs do not exhibit a known activating mutation in the pathway raising the possibility that additional genes with low frequency somatic events are yet to be identified. To further understand the somatic landscape of this pathway, we first characterized alterations among the known pathway members and then identified novel genes with mutually exclusive alterations. Novel alterations in known pathway genes included a recurrent in-frame-insertion in MAP2K1 and a fusion of MET with its neighboring gene, CAPZA2 (Figure 6; Supplementary Table 21)³⁰. Previously reported TRIM24-NTRK2 and KIF5B-MET fusions³⁰ were observed in tumors without other known activating alterations. Interestingly, another NTRK2 fusion with TP63 was also found in a lung SqCC (Figure 6; Supplementary Table 21). As observed previously, high MET and ERBB2 amplifications were enriched in tumors without other known activating alterations in this pathway (p < 0.01; Supplementary Fig. 16)⁶. A single lung adenocarcinoma (TCGA-49–4512) contained an activating EGFR kinase domain duplication³¹. By manual review, we found additional canonical mutations in KRAS, EGFR, or ERBB2 in 17 tumors and complex indels in EGFR or MET in 11 tumors, some of which have been previously reported^6,8,32 (Supplementary Table 22).

Lung ADCs that had an activating SSNV, indel, amplification or gene fusion in a known RTK/Ras/Raf driver^6,33,34 were designated “oncogene positive” (n=418) while the remaining lung ADCs were considered “oncogene negative” (n=242). For the purposes of this analysis, we did not include NF1-altered tumors in the oncogene positive group as mutations in this gene are not entirely mutually exclusive with alterations in other Ras/Raf/RTK related genes. To identify additional potential drivers in this pathway, we determined if genes that are significantly mutated in any of the MutSig2CV analyses (Supplementary Tables 5, 6, and 14) or that are important in regulation of the Ras pathway³⁵ were enriched in oncogene negative samples using a Fisher’s exact test. In total, 15 genes were significantly enriched among oncogene negative samples including known Ras pathway members SOS1 and RASA1, and Rho kinase pathway members VAV1 and ARHGAP35 (q-value < 0.1; Figure 7a,c; Supplementary Table 23). SOS1 is a guanine nucleotide exchange factor (GEF) bound to the RTK complex and assists in the activation of Ras proteins³⁶. Recurrent p.N233Y mutations were observed in the autoinhibitory domain (DH) of SOS1 in 4 lung ADCs and the mutation p.D309Y in the same region has been reported in Noonan syndrome^37,38 (Supplementary Fig. 17). Similarly, VAV1 is a GEF for the Rho family GTPases. Interactions between the calponin homology (CH), Acidic (Ac), and pleckstrin homology (PH) domains are important for autoinhibition of the catalytic Dbl homology domain³⁹. The p.S67Y mutation is located near the interface of the CH, Ac and PH domains and mutagenesis at this site has been shown to increase overall GEF activity³⁹ (Supplementary Fig. 17). RASA1 and ARHGAP35 (p190RhoGAP) are GTPase activating proteins (GAPs) for the Ras and Rho kinases, respectively, and are each enriched for loss-of-function mutations (p < 0.01). We also identified amplifications peaks near FGFR1/WHSC1L1 (8p11.21), PDGFRA/KIT/KDR (4q12), and MAPK1 (chr22q11) that were only significant in the oncogene negative tumor set (q-value < 0.25; Figure 7b,c). In total, 499 (76%) lung ADCs displayed an alteration in known or putative Ras/Raf/RTK driver genes (Figure 7c). Moreover, 193 out of 227 (85%) lung ADCs that previously underwent secondary expert pathological review and had RNA-seq data available for fusion analysis⁶ contained a predicted activating alteration in the RTK/Ras/Raf pathway.

Novel co-occurrences included MET amplifications and NF1 mutations (p = 0.019; Supplementary Figure 16). Additionally, high EGFR amplification significantly overlapped with activating EGFR mutations (p = 1.9 × 10⁻⁸)^40,41 and STK11 mutations significantly overlapped with activating KRAS mutations (p = 1.1 × 10⁻⁶; Figure 7c)^42,43. Furthermore, 28 lung ADCs that remain oncogene-negative for the RTK/Ras/Raf pathway harbor STK11 mutations (Figure 7c), suggesting the possibility of an additional hitherto-unrecognized KRAS-related genome alteration complementary to STK11 mutation in these cancer samples.

Assessment of neoantigen load and recurrence

With increasing interest in the use of immune checkpoint inhibitors in lung cancer^44,45, we comprehensively analyzed the potential immunogenic properties of the mutational landscape. Within each patient, we evaluated the ability of each somatic missense mutation to be processed and presented to immune cells by any one of the patient-specific HLA alleles^46,47. We then assessed the association between the number of immunogenic mutations (i.e. neoepitopes or neoantigens) and clinical characteristics and identified the most common neoepitopes observed in lung cancer. Both nonsynonymous mutation and neoepitope counts were not significantly different between lung ADCs and lung SqCCs from ever smokers (Figure 8a,b). However, these counts were significantly lower in lung ADCs from never smokers compared to lung ADCs from ever smokers (p < 0.001; Wilcoxon rank-sum test; Figure 8a,b) and associated with overall smoking history in lung ADCs but not lung SqCCs (p < 0.001; Kruskal-Wallis test; Supplementary Fig. 18). Mutations predicted to be neoepitopes in at least 4 tumors included PIK3CA p.E542K, NFE2L2 p.E79Q, BRAF p.G466V, and EGFR p.G719A and several mutations in TP53, including p.V157F, p.G154V, p.R175G, and p.P278A (Figure 8c). A gene not previously implicated in lung cancer, C3orf59 (also known as MB21D2), contained a recurrent mutation at p.Q311E with predicted neoepitope properties (Figure 8c). Overall, 47% of lung ADC and 53% of lung SqCC samples had at least 5 predicted neoepitopes suggesting a great potential for immunotherapy.

Sample collection and pathology review

Sample collection and DNA sequencing were performed for the Imielinski et al and TCGA cohorts as previously described^6–8. All specimens were obtained from patients with appropriate consent and with approval from the relevant institutional review boards. All patients were treatment-naïve with the exception of four patients with lung SqCC and three with lung ADC who received neoadjuvant treatment prior to resection (Supplementary Table 2). Initial pathological review was performed at the contributing tissue source sites (TSS) where each tumor was given an initial histological classification. After shipment of the frozen tissue to the Biospecimen Core Resource (BCR), one or two additional frozen sections were cut and stained with H&E to confirm the histological classification of the original TSS. 159 of the lung ADCs from Imielinski et al, 289 of the lung ADCs from TCGA, and 213 of the lung SqCCs from TCGA had also undergone additional histological review by an expert pathology committee led by Dr. William Travis (MSKCC) in previous studies^6–8. Nucleic acid extraction and molecular quality control were performed at the BCR.

DNA-sequencing, alignment, and mutation calling

Exome capture was performed using the Agilent SureSelect Human All Exon 50MB kit followed by Illumina paired-end sequencing. Reads were processed using the Picard pipeline⁶. This pipeline utilizes BWA for read alignment, Picard tools for marking duplicates, and the Genome Analysis Tool Kit (GATK) for realignment around small insertions and deletions (indels) as well as base quality recalibration⁵². Contamination in tumor exomes was estimated using ContEst⁵³. Only tumors with <5% contamination, an available SNP6.0 array for copy number analysis, and a valid ABSOLUTE⁵⁴ solution were considered in the final analysis. The final sample set included 227 previously described lung ADCs from the TCGA⁶, 274 newly reported lung ADCs from the TCGA, and 159 lung ADCs from the Imielinski et al⁸ cohort, together with 176 previously described lung SqCCs from TCGA⁷, and 308 newly reported lung SqCCs from TCGA. Somatic single nucleotide variants (SSNVs) and indels were called using MuTect⁵⁵ and Indelocator (www.broadinstitute.org/cancer/cga/indelocator), respectively. These algorithms compare the tumor to the matched normal in order to exclude germline variants. Somatic calls were excluded if found in a panel of over 2,900 normal exomes as previously described¹⁰. Coding mutation patterns can be viewed for individual genes at pubs.broadinstitute.org/panlung.

Identification of significantly mutated genes

Significantly mutated genes were identified using MutSig2CV which combines p-values from tests for high mutational frequency relative to the background mutation rate (pCV), clustering of mutations within the gene (pCL), and enrichment of mutations within evolutionarily conserved sites (pFN)¹⁰. For 660 lung adenocarcinomas, we had 100% power to detect genes mutated in 10% of patients and 73% power for genes mutated in 5% of patients assuming a mutation rate of 8.7/Mb¹⁰. For 484 lung squamous cell carcinomas, we had 100% power to detect genes mutated in 10% of patients and 41% power for genes mutated in 5% of patients assuming a mutation rate of 9.7/Mb¹⁰. In order to reduce the number of hypotheses tested in the MutSig2CV analysis, we excluded genes that exhibited low expression across tumors with relatively high purity. The median log₂ FPKM value for each gene was obtained for 185 ADCs or 238 SqCCs which had a purity estimate from ABSOLUTE of >50% and available RNA-seq data (Supplementary Fig. 1). For each tumor type, a mixture model of two normal distributions was fit in R using the mclust package v4.2. Genes with 95% probability of belonging to the cluster with higher expression were considered in the multiple hypothesis correction of the MutSig2CV combined p-values. One gene, TRERF1, was excluded from the final results as closer inspection of its mutations revealed a recurrent frameshift deletion that was likely a false positive as all of these mutations had low allelic fractions (<1.5%) and had no supporting reads in matching RNA-seq data. A one-sided Fisher’s exact test was used to determine if the proportion of loss-of-function mutations (including nonsense, frameshift, and de novo start out-of-frame mutations) to other mutations for a given gene was significantly higher compared to the proportion of loss-of-function mutations to other mutations across all other genes.

Identification of recurrent copy number changes

DNA was hybridized onto Affymetrix SNP 6.0 arrays and normalized as previously described⁶. Segmentation was performed using Circular Binary Segmentation algorithm⁵⁶ followed by Ziggurat Deconstruction to infer the length and amplitude each segment. Recurrent focal SCNA peaks were identified using GISTIC2.0⁵⁷. A peak was considered focally amplified or deleted within a tumor if the GISTIC2.0-estimated focal copy number ratio was greater than 0.1 or less than −0.1, respectively. Purity and ploidy were estimated using ABSOLUTE⁵⁴. Two peaks were considered the same across tumor types if 1) the known target gene of each peak was the same or 2) the genomic location of the peaks overlapped +/− 1 Mb and each of the overlapping peaks had less than 25 genes and was smaller than 10 Mb.

RNA sequencing for expression and fusion analyses

Of the 1,144 tumors examined in this study, 495 lung ADCs and 476 lung SqCCs also had corresponding RNA-seq data from TCGA. RNA reads were generated, aligned to the hg19 genome assembly with Mapsplice⁵⁸, and normalized with RSEM⁵⁹ to Fragment per Kilobase per Million (FPKM) expression estimates as previously described⁶. Expression values less than 1 FPKM were set to 1 and all data were log₂ transformed. Exon skipping of MET exon 14 was identified with juncBASE⁶⁰ as previously described⁶. Lists of fusions were obtained from previous studies^6,3061. Fusions for additional tumors were identified with the PRADA pipeline⁶². For plotting of exonic expression of fusion transcripts, exon expression levels were counted and normalized to reads per kilobase per million (RPKM) as previously described⁶. Expression for an individual exon was first Z-score transformed across all tumors within each tumor type. Subsequently, all exons for a gene were Z-score transformed again within each tumor. Transcript annotations used for this analysis included ENST00000397752 for MET, ENST00000361183 for CAPZA2, ENST00000302418 for KIF5B, ENST00000323115 for NTRK2, ENST00000343526 for TRIM24, and ENST00000354600 for TP63.

Identification of mutational signatures

Non-negative matrix factorization (NMF) was used to deconvolute a K × G matrix of mutation catalogues into a K × N matrix of mutational processes and an N × G matrix of mutational exposures (where G is the number of lung exomes, K is the number of mutational states, and N is the number of estimated mutational processes)¹⁶. Code to run NMF was obtained from MATLAB Central (see URLs) and run using the nnmf function from the MATLAB Statistics Toolbox. We used 6 mutation types with 16 different trinucleotide contexts and 2 transcriptional strands for a total of 192 mutational states. The number of possible signatures was varied from 1 to 10 and signature stability was assessed via bootstrapping as previously described¹⁶. Within each tumor, the fraction of estimated mutations for a signature was derived by dividing the number of estimated mutations for that signature by the sum of estimated mutations from all signatures.

Predicting Immunogenicity

HLA alleles were called with POLYSOLVER⁴⁶ for all lung cancer exomes. Within each tumor, epitope predictions were made between confidently called alleles and single amino acid missense mutations. Separate lists were made consisting of wildtype and mutant peptides of length 8, 9, 10 and 11 amino acids since these are the possible peptide lengths known to be presented by human MHC class I molecules⁶³. We then predicted MHC binding affinity for each of the peptide as described previously⁴⁷. First, the proteasome processing score was calculated using the NetChop program⁶⁴. Then, we used the NetMHC⁶⁵, NetMHCpan⁶⁶, SMM⁶⁷, and SMMPMBEC⁶⁸ methods to predict the MHC binding affinity values for each peptide and used the median affinity value across all algorithms as a composite measure of binding strength. We also defined the neoepitope ratio for each mutant-wildtype peptide pair as the mutant median affinity value divided by wildtype median affinity value. This value was found to be a reliable comparator of the relative immunogenicities of the mutant versus wildtype peptide sequences⁴⁷. Peptide pairs were further considered if the mutant peptide displayed a processing score ≥0.7, a median affinity value ≥0.01, a neoepitope ratio ≥1, and the mRNA transcript of the gene was expressed in the RNA-seq data for that tumor (top 15,000 expressed genes within each tumor). Since epitope binding is HLA-dependent, the previous steps were performed for each of the called MHC-I alleles. After that, only those peptides predicted to be the best epitopes for each mutation were considered.

Statistical comparisons

Nonparametric tests such as the Wilcoxon rank-sum test (comparison between 2 groups) or the Kruskal-Wallis test (comparison between more than two groups) were used for continuous variables unless otherwise noted. The Fisher’s exact test was used when comparing 2 categorical variables. In total, longitudinal data on survival were available for 481 patients with lung ADC and 473 patients with lung SqCC from TCGA. The Cox proportional hazards model was used to examine associations between patient survival and mutation status, with and without controlling for stage. Correction for multiple hypothesis testing was performed with the Benjamini-Hochberg procedure.