Genomic and transcriptomic analysis of checkpoint blockade response in advanced non-small cell lung cancer

1The introduction of PD-1/PD-L1 inhibitors in the management of advanced non-small cell lung cancer (NSCLC) has led to a major paradigm shift in the treatment of the disease. Following multiple studies demonstrating improved overall survival, these agents have garnered approval either alone– or in combination with chemotherapy, or CTLA4 blockade. However, with responses observed in only one in five unselected patients–, improved predictors of response are needed to identify patients most likely to benefit.

2Given that the potential for long-term disease control is only realized in a minority of patients, extensive effort has been dedicated to identifying biomarkers of response and resistance. The dominant biomarkers to date are PD-L1 protein expression on tumor cell membranes and tumor mutational burden (TMB)–, which may underlie the generation of neoantigens that can serve as targets for immune recognition and targeting.

3While additional features have begun to emerge including potential roles for mutation clonality, an inflamed microenvironment, and alterations in individual genes such as EGFR, and STK11 (ref. ), further identification and integration of relevant predictors have been hindered by the absence of large, multi-omic, NSCLC-specific patient cohorts.

4Here we describe findings from the first integrative analysis of the Stand Up To Cancer-Mark Foundation (SU2C-MARK) NSCLC cohort, a dataset of 393 patients treated with checkpoint inhibitors in the advanced-stage setting. We performed whole exome sequencing (WES) and RNA sequencing (RNA-seq) along with detailed clinical response assessments, enabling the composite assessment of genomic and transcriptomic biomarkers of response and resistance. Collectively, these richly annotated data will be a resource to the field in furthering both the basic and applied investigation into the role of PD-1/PD-L1 agents in advanced NSCLC.

5We analyzed formalin-fixed paraffin-embedded (FFPE) tumor samples collected before receipt of checkpoint blockade (defined as the first line of therapy in which a patient received a PD-1/PD-L1 agent) from a total of 393 patients with advanced NSCLC across nine cancer centers (Table 1 and Fig. 1a). The majority of these patients were treated with single-agent therapy (81%), with additional subsets receiving combination therapy including either CTLA4 blockade (17%) or chemotherapy (1%). Both tumor and matched normal specimens (from blood, or in rare cases, adjacent normal tissue) underwent WES; for a subset of patients, tumor tissue was additionally profiled by whole transcriptome RNA-seq. After stringent quality control (Methods), a total of 309 WES and 152 RNA-seq specimens were included for analysis. The primary outcome was best overall response (BOR) determined by a dedicated review of clinical imaging and quantified using RECIST v1.1 criteria.

6As is typical for patients with NSCLC, the SU2C-MARK cohort consisted predominantly of adenocarcinoma (73%) and squamous cell carcinoma (20%), with smaller contributions from large cell neuroendocrine (LCNE) carcinoma (2%) and other histologies (4%; Extended Data Fig. 1a). Among patients with annotated PD-L1 staining assessments (224/393 available, 43% missing), 25% had a tumor proportion score (TPS) of less than 1%, 33% had PD-L1 TPS 1–49%, and 42% had PD-L1 TPS ≥ 50%. As expected, higher PD-L1 TPS was associated with an increased response rate to checkpoint blockade (Extended Data Fig. 1b). Thus, our dataset reflected the histologic and biomarker compositions typically observed in unselected, real-world NSCLC cohorts,.

7To better understand the relationship between mutational drivers and response, we assessed the prevalence of known drivers in lung cancer across our three response categories: partial or complete response (PR/CR), stable disease (SD) and progressive disease (PD; Fig. 1b and Extended Data Fig. 2a). Consistent with prior reports–, nonsynonymous TMB associated with response category (P = 6 × 10−9), with median TMB 14.0 mut/MB among those with PR/CR, compared to 9.0 mut/MB for SD, and 7.4 mut/MB for PD (Fig. 1c). Initial examination of the cohort was also consistent with previously observed driver associations,,, such as EGFR alteration or KRAS/STK11 comutation being a negative predictor of checkpoint blockade response (Extended Data Fig. 1c,d).

8To facilitate a more comprehensive analysis, we performed logistic regression, testing the relationship between 49 known lung cancer drivers, and response (that is, CR/PR versus SD/PD; Methods). In all, six genes achieved significance or near significance, defined as a false discovery rate (FDR) threshold of 10% or 25%, respectively (Fig. 1d). In this analysis, mutations in ATM appeared to be most favorable with respect to checkpoint blockade response (logistic regression FDR q = 0.04, OR = 3.5, CI95% (1.5, 8.0)), while EGFR alterations were least favorable (q = 0.12, OR = 0.29, CI95% (0.11, 0.79)). Given the strong association between ATM and response in our cohort, we tested this association in an independent cohort of patients with NSCLC treated with PD-(L)1 blockade and profiled by MSK-IMPACT. In this external cohort, ATM alterations were associated with improved overall survival following checkpoint blockade (P = 0.03; Extended Data Fig. 2b). As this association was not seen at the cohort-wide level (P = 0.45), these results suggest a predictive rather than simply prognostic role for ATM alteration.

9We next explored relationships between copy number alterations and response in the cohort (Extended Data Fig. 3a). Among focal events, only focal amplification of 5p15.33, the cytoband containing TERT, achieved significance, and was associated with reduced response to immunotherapy (q = 0.07, OR = 0.59, CI95% (0.40, 0.87); Fig. 1e and Extended Data Fig. 3b). Of note, this association was not reproduced in the MSK-IMPACT cohort, which may be a function of the more limited sensitivity of amplifications in panel data (Extended Data Fig. 3c). Taken together, these results suggest that in addition to the aggregate metric of TMB, individual driver events may also define favorable and unfavorable NSCLC subsets for checkpoint blockade.

10To better understand how the determinants of immune recognition in our cohort related to response, we calculated the neoantigen burden for each exome in the SU2C-MARK cohort (Methods). Total neoantigen burden was significantly associated with response (q = 4 × 10−5, OR = 8.8, CI95% (4.2,19); Fig. 1f). As clonal neoantigens have been suggested to be more effective targets of immune recognition, we additionally examined the role of clonal and subclonal neoantigen burden, along with total subclone count (Methods). Indeed, clonal neoantigen burden was also significantly associated with response (q = 2 × 10−4, OR = 5.4, CI95% (2.7,11)), whereas neither subclonal neoantigen burden nor total subclone count was significant (q = 0.7 and q = 0.6, respectively; Fig. 1f).

11As different mutational processes may have different propensities for neoantigen generation, we also evaluated the mutation burden attributable to distinct mutational signatures (Extended Data Fig. 4a,b; Methods). Of the three dominant signatures, smoking was most strongly associated with response (q = 5 × 10−5), consistent with its association with clonal neoantigens, while aging (q = 0.05) and APOBEC (q = 0.01) were more weakly associated with response (Fig. 1f). We additionally observed a significant response association for indels (q = 2 × 10−5), which are suspected to be particularly immunogenic given their potential to generate new reading frames,.

12Previous studies have suggested that compromised antigen presentation, due to loss of heterozygosity (LOH) at HLA loci, decreased total unique HLA alleles, or loss of B2M may enable immune evasion in certain cancer types. We did not observe an association of any of these factors measured before therapy and response in this cohort (Fig. 1f), potentially suggesting disease-specific variation in mechanisms of resistance.

13To further assess for variation in immune infiltrate, we used MiXCR to identify B- and T-cell clonotypes from rearranged VDJ reads in our WES data (Methods). Of these subsets, T-cell receptor (TCR) burden was associated with response but did not reach statistical significance (q = 0.3). Thus, among our expanded set of exome-derived features, tumor-intrinsic markers reflective of TMB as well as clonal mutation burden emerged as top predictors of response.

14We next focused on the identification of transcriptional predictors of response. Using limma voom, we performed a genome-wide analysis of differentially expressed genes between responders (PR/CR) and nonresponders (SD/PD; Fig. 2a and Methods). Initial assessment of these results identified three related genes that achieved cohort-wide significance (padj < 0.05; Methods): PSME1, PSME2 and PSMB9. These genes are notable for their prominent role in the function of the immunoproteasome (further described below), a noncanonical peptide processing complex thought to promote differential and enhanced antigen presentation in the setting of proinflammatory cytokines. Examination of the broader collection of genes achieving nominal significance (nominal P < 0.05) revealed additional interferon-gamma (IFN-γ)-induced transcripts including TAP1 (a cytosolic peptide transporter in the antigen presentation pathway) and CD274 (which encodes PD-L1), inflammatory chemokines such as CXCL9, CXCL10 and CXCL11, and lymphocyte receptor genes (for example, CD3D and CD7), potentially surrogates for immune infiltration (Fig. 2a and Extended Data Fig. 5a). Top genes associated with nonresponse appear to span both developmental and immune-related pathways. AUTS2 and TCF7L1, interacting transcription factors within the Wnt/B-catenin signaling axis, are postulated to have roles in both stem cell and immune signaling. Another nonresponse-associated gene, PDLIM3, is a member of a protein family thought to negatively regulate NF-κB-mediated inflammatory responses. KALRN, a guanine nucleotide exchange factor expressed in stromal and myeloid cells, has been associated with inflammation in the context of atherogenesis (Extended Data Fig. 5a).

15To systematically identify differentially expressed pathways, we performed gene set enrichment analysis (GSEA) using the Hallmark Gene Sets (Fig. 2b). Top response-associated pathways included ALLOGRAFT_REJECTION, INTERFERON_GAMMA_RESPONSE and DNA_REPAIR, which has previously been observed as a predictor of checkpoint blockade response in urothelial carcinoma,. Pathways associated with resistance were diverse, with EPITHELIAL_MESENCHYMAL_TRANSITION, WNT_BETA_CATENIN_SIGNALING and TGF_BETA_SIGNALING gene sets all significantly associated with nonresponse (Fig. 2b). Taken together, these top genes and gene sets from bulk RNA-seq suggest the relevance of both immune and nonimmune components to the biology of checkpoint blockade.

16Given the remarkable convergence of all three genes (PSME1, PSME2 and PSMB9) on components of the proteasome/immunoproteasome system responsible for peptide generation, we expanded our exploration of genes specific to this antigen presentation pathway. Notably, PSME1 and PSME2 encode for the IFN-γ inducible PA28ɑβ complex that binds and enhances peptide processivity of both the constitutive and immunoproteasome. PSMB9 (LMP2) encodes the β1i IFN-γ inducible subunit that together with β2i (PSMB10) and β5i (PSMB8) represent the three inducible subunits whose incorporation transforms the constitutive proteasome into a specialized immunoproteasome with distinct peptide cleavage patterns. Hence, as all the inducible components of this complex (PSMB8, PSMB9, PSMB10, PSME1 and PSME2) are known to be downstream of IFN-γ, which itself was nominally associated with response in our analysis (IFNG P = 0.001; log2 fold change 1.1), we evaluated the response association of these components alongside canonical IFN-γ targets (HALLMARK_INTERFERON_GAMMA_RESPONSE) as well as a comprehensive list of proteasome components (GOCC_PROTEASOME_COMPLEX; Fig. 2c). Notably, immunoproteasome components were enriched in terms of the significance of association with response relative to both IFN-γ targets more broadly, as well as all proteasome components (P = 9 × 10−9 and P = 2 × 10−5, respectively; Fig. 2c and Extended Data Fig. 5b).

17Although the inducible subunits of the immunoproteasome were highly correlated with one another, increases in their expression could only partly be explained by elevated levels of IFN-γ (Extended Data Fig. 5c). Given that experimental evidence suggests they may also be induced by tumor necrosis factor-α (TNF-α), we evaluated whether higher levels of TNF may also contribute to upregulation of these components. Indeed, a linear combination of IFNG and TNF demonstrated improved model fit for immunoproteasome subunit expression (R2 = 0.31 for the combined model compared to 0.19 for the univariate model; Fig. 2d). Thus, immunoproteasome subunits appear to be singularly important predictors of response—even among the broader class of IFN-γ-induced transcripts—perhaps owing to their role as integrators of multiple cytokine cascades, enabling more efficient generation of peptide epitopes for HLA-I presentation.

18Given that both individual gene and pathway level analysis highlighted key roles for immune signaling, we aimed to better delineate discrete immune cell subsets in our bulk transcriptome data using previously identified signatures derived from single-cell RNA data (Methods). Of the 11 signatures we evaluated, exhausted CD8+ T-cells showed the strongest positive association with response, while the monocyte/macrophage and dendritic cell signatures were most strongly associated with resistance (Fig. 2e).

19As a growing body of work suggests that distinct myeloid subsets may have differing roles in antitumor immunity,, we investigated more specific subsignatures related to these cell types. Using a marker list derived from a comprehensive single-cell RNA-seq study of infiltrating myeloid cells in human and mouse lung cancers, we identified the hMono3 and hN3 subtypes as being particularly associated with resistance to checkpoint blockade (Extended Data Fig. 6). Notably, the hMono3 subtype is characterized by high expression of S100A8, a cytokine-like protein that can drive the accumulation of myeloid-derived suppressor cells. The neutrophil hN3 subtype is defined by high expression of CXCR2, which has been shown to inhibit CD8 T-cell activation within the lung cancer microenvironment. Thus, our focused analysis of immune subsets identified plausible mechanistic connections between myeloid infiltration and decreased response to checkpoint blockade.

20To identify M signatures relevant to immunotherapy response beyond individual cell types, we applied Bayesian non-negative matrix factorization (B-NMF) to our top 770 differentially expressed genes, yielding three distinct M signatures as follows: M-1, M-2 and M-3 (Fig. 3a,b; Methods). Because these signatures were derived from bulk sequencing, they are expected to reflect the complete microenvironmental signature, inclusive of both tumor and nontumor (that is, immune and stromal) sources. GSEA of these signatures revealed M-1 to be associated with epithelial–mesenchymal transition (a gene set that includes wound healing and fibrosis) and M-2 to be associated with allograft rejection/IFN-γ response, consistent with an inflamed immune environment (Fig. 3c). M-3 had a weak association with cell cycle-related E2F targets, potentially reflecting a proliferative tumor signature, which in conjunction with the relative depletion of infiltrating myeloid and lymphoid cells, most resembles the previously reported immune desert phenotype (Fig. 3d and Extended Data Fig. 7). Notably, the response rate to checkpoint blockade varied across these subtypes, with increased response rates observed in M-2 relative to M-1 and M-3 (P = 0.06; Fig. 3e). Overall, these results suggest that there may be at least two distinct transcriptional states associated with checkpoint blockade resistance in NSCLC.

21Having explored aggregate microenvironmental states, we next turned our attention to tumor intrinsic expression factors that may have a relationship with response. To define relevant tumor intrinsic (TI) lung cancer subtypes, we assembled a large reference collection of over 1,000 transcriptomes (TCGA-LCNE) representing the three predominant NSCLC histologies, namely adenocarcinoma, squamous cell carcinoma and large cell neuroendocrine carcinoma (Fig. 4a and Methods). To define signatures of individual subtypes in this collection, we first performed B-NMF across this cohort, converging on a robust four-cluster solution (Fig. 4b and Extended Data Fig. 8a). Of these TI clusters, TI-1 and TI-2 contained predominantly adenocarcinomas, TI-3 was composed largely of squamous cell carcinomas, and TI-4 was primarily large cell neuroendocrine carcinomas (Extended Data Fig. 8b). Notably, unlike our M signatures above—which were derived solely from the subset of genes with significant response associations and were enriched for immune and stromal components—our TI signatures emerged from the unsupervised factorization of primary lung cohorts spanning three distinct histologies, explaining the high concordance between our TI subtypes and existing histologic categories.

22To understand these signatures in more detail, we explored the expression of canonical markers of adenocarcinoma and squamous differentiation, namely NAPSA (which encodes Napsin A) and TP63 (which encodes both p63 and p40), respectively (Extended Data Fig. 8c). While TI-2 and TI-3 showed the expected lineage marker preferences, TI-1 samples showed weak expression of both markers. Decreased expression of lung lineage markers has previously been described in a subtype of poorly differentiated adenocarcinomas in which markers for adjacent gut lineages (neighboring endodermal territories during development) can become activated. Indeed, a comparison of these subtypes to immunohistochemical markers of various endodermal lineages revealed enrichment in these gut-specific marker genes in TI-1 samples, such as TFF1, FGA and CPS1 (Fig. 4c). TI-1 samples were also notable for an elevated TMB relative to the well-differentiated TI-2 adenocarcinoma subtype and the TI-3 squamous subtype (Extended Data Fig. 8d).

23Having established a reference collection of TI expression signatures, we applied these signatures to RNA-seq data from the SU2C-MARK cohort and assessed their association with response to checkpoint inhibitors. Notably, the dedifferentiated TI-1 cluster was most closely associated with response (Fig. 5a), consistent with the elevated mutational burden in this subtype as well as its stronger association with the M-2 ‘immune-activated’ subtype (Fig. 5b and Extended Data Fig. 8e). Indeed, patients with both immune-activated (M-2) and dedifferentiation (TI-1) signatures had the highest response rates to checkpoint blockade (67% ORR; Fig. 5c). Thus, TI states and immune M signaling may independently and additively govern responses in NSCLC.

24Having evaluated a broad set of clinical, genomic and transcriptomic features relevant to checkpoint blockade response in NSCLC, we set out to better understand the relationships between these predictors. Combining the top predictive features from each analysis, we generated a cross-correlation matrix to better understand how they relate to each other as well as to previously published signatures relevant to tumor biology and immune response (Fig. 6 and Methods),–. Notably, three strong correlation blocks could be observed, with consistent response associations within each subset. The first correlation block (C1) appeared to reflect a canonical ‘wound healing’ microenvironment, including immunosuppressive myeloid and stromal signatures. The second correlation block (C2) reflected the more classic cytokine and immune milieu associated with ‘immune activation/exhaustion,’ including both infiltrating immune signatures and proteasome subunits. The third correlation block (C3) consisted of features related to mutational burden, presumably all proxies for neoantigen abundance and consequent enhanced immune recognition.

25The remaining nine features were somewhat loosely correlated as a fourth cluster (C4) enriched for single-gene alterations with potentially distinct immunobiologies. Notably, this cluster included EGFR mutations, which interestingly showed minimal association with the immune signatures but a moderate anticorrelation with mutational burden features, suggesting the intrinsic resistance of this subtype may predominantly be driven by insufficient neoantigens (Fig. 6 and Extended Data Fig. 9a).

26To evaluate whether the additional genomic predictors identified in this study could augment existing biomarker-defined subsets of NSCLC, we selected the top two significant predictors from each cluster and evaluated their potential to further stratify progression-free survival (PFS) in three clinically relevant subgroups: TMB > 10 mut/MB (favorable; n = 27), PD-L1 TPS ≥ 50% (favorable; n = 34) and PD-L1 TPS ≤ 1% (unfavorable; n = 18). Following FDR correction, we identified multiple near-significant and significant associations (q < 0.25 and 0.1, respectively; Extended Data Fig. 9b,c and Methods), particularly when evaluating features from the immune activation/exhaustion and wound healing clusters (dedifferentiated TI-1 in PD-L1 TPS ≤ 1% q = 0.23; immune-activated M-2 in PD-L1 TPS ≤ 1% q = 0.16; macrophage/monocytes in PD-L1 TPS ≥ 50% q = 0.06; hMono3 in PD-L1 TPS ≥ 50% q = 0.11). Therefore, the presence of these factors may augment prediction based on standard clinical variables.

27Given that the predictors identified in this study were derived from bulk specimens, they likely reflect contributions from multiple distinct cell types within the tumor microenvironment. To gain additional insight into the specific cellular components that may be driving response and resistance, we explored these predictors in the context of published single-cell sequencing data from NSCLC within mixed tumor environments that may be contributing to these signals in bulk data. Evaluation of the marker expression from the 13 cancer-related clusters revealed a straightforward mapping to several TI subtypes described earlier, including one cluster (cluster 12) which mapped to our dedifferentiated TI-1 subtype (Fig. 7a and Extended Data Fig. 10a; Methods).

28Deconvolution of the unfavorable wound healing (C1) predictors suggested that the EMT and TGF-β signatures predominantly reflected fibroblasts and endothelial cells as opposed to a mesenchymal epigenetic state per se within the tumor cells; conversely, some of the dominant single-gene transcriptional predictors such as AUTS2 and TCF7L1 demonstrated substantial tumor intrinsic expression (Fig. 7b and Extended Data Fig. 10b). Similarly, analysis of the favorable predictors in the immune activation/exhaustion cluster (C2) revealed that while immunoproteasome subunits are expressed in most cell types, CXCL9 may be predominantly expressed by myeloid sources, and CXCL11 may be primarily derived from endothelial cells (Fig. 7b and Extended Data Fig. 10b). Finally, our favorable dedifferentiated (TI-1) and immune-activated (M-2) predictors, while correlated at the bulk level, did appear to identify distinct subpopulations (cancer cells and T-cells, respectively) at the single-cell level, consistent with our labeling of these signatures as tumor intrinsic versus microenvironmental (Fig. 7b). Taken together, these findings suggest the presence of rich, interacting ecosystems that may broadly underlie response and resistance to checkpoint blockade and provide a collection of specific signaling pathways and cell types that may be promising targets for future intervention.

29Comprehensive identification of predictors of checkpoint blockade response in patients with NSCLC has been limited by the availability of large, well-annotated patient cohorts with matched genomic data, particularly within individual cancer types. Here we present a joint analysis of the SU2C-MARK cohort, a collection of nearly 400 patients with NSCLC, enabling the identification of diverse molecular predictors of immunotherapy response.

30Among the top genomic features identified were ATM mutation and TERT amplification. Given emerging literature associating ATM loss with the release of cytosolic DNA and activation of the cGAS/STING pathway in other cancer types–, it is conceivable that a similar mechanism underlies the association observed in our cohort between ATM loss and response. Although less well characterized in the context of immunotherapy, TERT amplification may serve a protective function against telomere crisis, thereby forestalling a parallel mechanism, which has been linked to cGAS/STING activation and subsequent sensitization to checkpoint blockade in mouse models.

31Transcriptomic analysis in the SU2C-MARK cohort was notable for the identification of immunoproteasome subunit genes as key predictors of response, with greater enrichment than general IFN-γ targets or proteasome subunits. These findings are consistent with those described in melanoma, where a supervised signature consisting specifically of PSMB8 and PSMB9 was found to be predictive of immune checkpoint blockade response. We speculate that enhanced peptide supply to MHC-1 via increased expression of the PA28ɑβ complex and immunoproteasome may result in superior CD8+ T-cell responses. In addition, the altered cleavage specificity of the immunoproteasome relative to the constitutive proteasome—particularly in terms of preferences for branched-chain amino acids and chymotrypsin-like target sites— may confer increased antigen quality in addition to quantity in immunoresponsive tumors.

32Higher level organization of the strongest genes associated with response and resistance identified microenvironmental signatures previously associated with relevant immune states such as the immune-activated (M-2) signature and immune desert (M-3) signature. The wound healing (M-1) signature, although less well described in the context of lung cancer, does match the TGF-β transcriptional signature thought to drive T-cell exclusion in bladder cancer. While the immune desert (M-3) signature was somewhat more enigmatic, the top-weighted genes appear to be largely tumor intrinsic, suggesting they may directly reflect a tumor state unfavorable to immune invasion. Consistent with this notion, one of the top-weighted genes in the signature, DSC3, is a component of intercellular desmosome junctions that can act as barriers to immune infiltration.

33In addition to features such as these global immune states that may have pan-cancer relevance, we also describe a dedifferentiated (TI-1) NSCLC-specific subtype identified independently in both bulk and single-cell data using unsupervised approaches. A similar subtype has been described in mouse lung cancer models featuring a decreased expression of classic lung lineage markers as well as enhanced expression of developmentally adjacent endodermal lineages. The correlation between this tumor intrinsic state and our immune-activated (M-2) signature could represent an underlying differentiation state more susceptible to immune recognition (for example, via the presentation of oncofetal antigens), or conversely, a cell state change in response to an inflammatory cytokine milieu. Establishing the direction of causality between these signatures may have important implications for further therapeutic intervention.

34Finally, integrative analysis of our genomic features along with previously reported signatures relevant to immune and tumor biology supported the notion of a complex interplay between distinct signaling pathways (for example, CXCL9 versus TGF-β signaling) and distinct cell types (for example, myeloid cells versus fibroblasts), shedding light on some of the multifaceted interactions underlying checkpoint blockade responsiveness. Particularly noteworthy in this respect is the recognition that a number of features identified here may be truly tumor intrinsic predictors, which aside from a handful of specific driver events, or defects in antigen presentation have been somewhat elusive in NSCLC. It is our hope that the SU2C-MARK cohort continues to serve as a rich resource for further unraveling the complex architecture of relevant genomic predictors, and for generating deeper insights into the biology of antitumor immunity.

35All patients in the SU2C-MARK cohort consented through umbrella sequencing protocols approved under local institutional review board protocols at their respective cancer centers (Dana-Farber Cancer Institute 02-180, Massachusetts General Hospital 13-416, MD Anderson PA13-0589, Memorial Sloan Kettering 12-245, Columbia University IRB-AAA05706, University of California Davis LCRP-001, Yale 1411014879, Johns Hopkins IRB00100653). All samples in this study were from patients treated with anti-PD(L)1 therapy either as a single agent or in combination with other agents between 2009 and 2019. Although this cohort predominantly corresponds to standard-of-care therapy, a subset of patients from MSKCC treated with dual checkpoint blockade was derived from sequencing of specimens collected during the course of Checkmate 012 (NCT01454102; ref. ).

36Samples collected typically correspond to the first standard-of-care confirmation of metastatic disease, and therefore reflect a timepoint before receipt of any advanced therapy. Response data were assessed using RECIST v1.1 criteria through a dedicated radiologist review of standard-of-care clinical restaging studies (or in a subset of cases, imaging obtained while on a trial protocol). Confirmed BOR was determined using radiographic data following the first line of therapy involving a PD(L)1-based agent. PFS and overall survival were defined from the date of treatment start with a PD(L)1 agent until the first evidence of radiographic/clinical progression or date of death, respectively, and censoring was based on the date of last follow-up. To facilitate further analyses, WES and RNA-seq specimens were divided into two cohorts with cohort 1 corresponding to roughly the first 80% of available samples. Of note, a subset of these samples has been described previously in institution-specific collections,.

37Informed consent was obtained under the institutional protocols listed above. Patients were not compensated for their participation. In all, the cohort consisted of 393 patients undergoing checkpoint blockade therapy. Patients in the cohort ranged in age from 29 to 90 years. In total 182 patients were male and 207 patients were female. Additional details on the cohort distribution are described in Extended Data Fig. 1a.

38WES of DNA was performed at the Genomics Platform of the Broad Institute of Harvard and MIT as described previously,, with the exception of samples previously sequenced at Johns Hopkins and Yale University. In brief, DNA was extracted from FFPE tumor specimens and either matched normal whole blood, or in cases where this was unavailable, from adjacent normal FFPE specimens. Extraction was performed using the Qiagen AllPrep DNA/RNA Mini Kit (80204). A single aliquot of 150–500 ng input DNA in 100 μl TE buffer was used for library generation. Library preparation was performed using the Kapa HyperPrep kit, and quantification was performed using PicoGreen. Adapter ligation was performed using the TruSeq DNA exome kit from Illumina per manufacturer’s instructions. Sequencing of pooled libraries was performed using a HiSeq2500 with 76 bp paired-end reads. The mean target coverages for tumor and normal samples were 150× and 80×, respectively.

39Initial alignment of all samples to the hg19 genome was performed using the Broad Picard pipeline (v2.4.1), specifically with bwa 0.5.9 (ref. ). The Broad Cancer Genome Analysis group somatic mutation pipeline was run in the cloud platform Firecloud/Terra. Specifically, the first-pass quality control was performed by assessing sample contamination using ContEst and identifying potential sample swaps using the Picard CrossCheckFingerprints tool (using software versions from the GATK 4.0.5.1 release). Somatic single nucleotide variants (SNVs) and indels were called using a combination of MuTect, MuTect2 (ref. ) and Strelka. Recovery of somatic variants filtered due to tumor contamination in the matched normal was performed using DeTiN v1.7 (ref. ) followed by annotation with Oncotator v1.9 (refs.,). Adjacent SNV events were merged to di-nucleotide variants (DNVs), and filtering was performed using OxoG and FFPE Orientation Bias filters as well as removal of events observed in a panel of normals composed of TCGA and Illumina Capture Exome normals. Finally, a BLAT realignment filter was implemented to eliminate potentially spurious variants resulting from mismapped reads. To meet quality control criteria for inclusion in the exome cohort, samples were required to have mean and median target coverage >50×, contamination <5% and tumor purity >10% as assessed by ABSOLUTE (v1.5). Comparison of MutSig2CV driver analysis from the SU2C-MARK cohort agreed well with previously published results for TCGA lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) cohorts (Extended Data Fig. 2a).

40TMB was calculated as the natural log of nonsynonymous SNVs, DNVs and indels in a sample divided by the size of the Illumina exome capture territory in megabases. Signatures for the SU2C-MARK cohort were determined using the SignatureAnalyzer Bayesian NMF (v1.2) method–. In brief, we pooled TCGA LUAD, TCGA LUSC and SU2C-MARK cohort samples to improve our power for detection of rare signatures and performed unsupervised signature extraction using 20 random initializations. Thirteen runs converged to a seven-signature solution, so the k = 7 solution with maximum posterior probability was selected for downstream analysis. Assessment of cosine similarity between the seven signatures identified and the previously described COSMIC signatures was used to assign labels to each, with the three dominant signatures representing aging, APOBEC and smoking. Signature attributable mutation burden was calculated as the relative projection strength for each signature in a given sample. Dominant signatures identified across the cohort are shown in Extended Data Fig. 4. Log values of the mutation, signature and clonal/subclonal burdens were calculated using a pseudocount of one event per MB.

41Potential neoantigens were identified by first running POLYSOLVER (v1.0) to identify MHC Class I alleles from matched normal WES data. Predicted binding affinity for all possible 9mer and 10mer peptide sequences overlapping single and di-nucleotide somatic variants was assessed using NetMHCPan-4.0 (refs.,,). Neoantigens with percentile ranks of two or less for any Class I allele in the same patient were counted as predicted binders.

42Somatic copy number alterations were assessed from WES using the GATK4 CNV pipeline on Firecloud/Terra (corresponding to GATK v4.0.8.0). A copy number panel of normals (n = 820 samples) was generated from a collection of FFPE as well as fresh frozen samples filtered to have less than 1% of tumor in normal contamination. GATK CNV bin length was set to zero, and read counts were processed using the hg19 Illumina Capture Exome (ICE) targets with padding of 250 bases. Intervals were filtered for having less than 1% of samples with zero coverage by setting—maximum-zeros-in-interval-percentage to 1. The minimum total allele count for informative heterozygous SNPs was set to 10. GISTIC2.0 was used to process the allelic somatic copy number data to identify recurrent copy number altered regions across the cohort. The continuous copy number output values (rather than binned value) for focal and gene-specific events from GISTIC were used as inputs for downstream analysis. Comparison of significant recurrent alterations showed good consistency between the SU2C-MARK cohort and prior TCGA publications (Extended Data Fig. 3).

43Tumor purity and ploidy were estimated using ABSOLUTE (v1.5),. Specifically, somatic mutation and copy number data were used as inputs, and purity/ploidy solutions were evaluated manually. In general, solutions were selected with a preference for describing the observed data well at modeled integer copy numbers, being parsimonious (for example, diploid as opposed to genome doubled), appropriately fitting full deletions and having an alpha/2 line centered within the highest somatic SNV allelic fraction peak (where alpha represents the model purity). A gene-specific integer copy number and LOH were inferred from integer copy number segmented output from total or allele-specific copy number analysis, respectively. Samples with less than 10% purity were excluded as a filtering step during WES quality control as above.

44The subclonal architecture was inferred from ABSOLUTE input using PhylogicNDT (v1.0). Mutation clonality across single samples was modeled using a Dirichlet process, enabling the assignment of mutations to discrete subclones with imputed cancer cell fractions (CCFs). Variants assigned to clusters with CCF over 0.85 were classified as clonal, while the remainder were deemed subclonal. Subclone count was based on the total number of unique subclones identified by 1D Phylogic analysis.

45Rearranged reads corresponding to T- and B-cell receptors were identified from WES data using MiXCR v3.0 (ref. ). Primary BAM files were processed with the ‘analyze shotgun’ pipeline, and reads corresponding to TCR or Ig clonotypes with productive rearrangements (that is, those leading to in-frame rearrangements without stop codons) were summed to give a total TCR or Ig read count per sample. To infer relative T- or B-cell abundance, these read counts were normalized by calculating T-cell and B-cell burden, defined as (rearranged receptor count reads + 1)/(aligned reads/106). Natural log of this burden metric was used during the significance assessment.

46In total, 106 features derived from whole exome and transcriptome analysis were evaluated (Supplementary Table 30). Features reflecting mutation burden (for example, TMB, Neoantigens, etc.) were log-transformed before evaluation. Mutation and copy number features were filtered to include only those present in at least 5% of the cohort. Each feature was assessed in a univariate logistic regression model of BOR, binned as responders (PR/CR) versus nonresponders (SD/PD). FDR calculation was performed using the Benjamini–Hochberg method, with features categorized as significant (FDR < 0.1) or near-significant (FDR < 0.25).

47RNA-seq data were processed using the GTEx RNA-seq pipeline with the use of the GENCODE v19 reference transcriptome, followed by quality control evaluation using the RNA-SeQC2 (v1.0) pipeline,, generating both expression data as transcripts per million (TPM) as well as quality metrics. Specifically, this pipeline uses STAR (v1.0) alignment with the following settings: alignIntronMax = 1,000,000, alignIntronMin = 20, alignMatesGapMax = 1,000,000, alignSJDBoverhangMin = 1, alignSoftClipAtReferenceEnds = True, chimJunctionOverhangMin = 15, chimMainSegmentMultNmax = 1, chimSegmentMin = 15. Alignment is then followed by: (1) omission of reads that are unmapped, have secondary alignments (0 × 100 flag) or have the quality control fail flag (0 × 200), and (2) filtering for high-quality exonic reads that uniquely map as pairs (0 × 2 flag) and have fewer than six mismatches to ultimately generate gene-level expression data as well as associated quality metrics. Using the median exon TPM (CV), the number of genes detected, and other measures, we selected the highest quality samples (n = 152) for subsequent analysis.

48To analyze differentially expressed genes, we restricted our search to protein-coding transcripts, and those minimally expressed at a log2TPM of 0.5 or higher in at least 30% of our samples. Using the BOR groupings of responders (PR/CR) versus nonresponders (SD/PD), we then used the R package limma voom to identify genes differentially expressed with respect to response.

49Using the signed, log-transformed P values from the differential expression results, we performed enrichment analyses using the ‘fgsea’ package (v3.16) and the Hallmark Gene Sets from the Molecular Signatures Database (MSigDB).

50Using existing literature, we derived metagenes for clinically important features. Starting with groups of genes associated with a certain feature (for example, genes expressed according to B-cell abundance), we took the mean of the log2-transformed TPMs in our cohort, then compared samples to each other by z-scoring those averages. These analyses include metagenes for different groups of leukocytes, which we use as a proxy for the level of immune infiltration indicated by RNA-seq. Additionally, we used previously published markers of developmental lineage, and NSCLC subtypes to better understand the developmental identity of each sample. We also defined an additional gene set for neuroendocrine identity using markers from a published characterization of large-cell neuroendocrine lung cancer. For cell type-specific characterization, we used metagenes from single-cell studies of lung cancer developmental subtypes and immune infiltrate,.

51We applied the B-NMF algorithm,,, to organize the significantly differentially expressed gene set from cohort 1 of our RNA-seq data (n = 123) into three distinct clusters, that is, our M subtypes. We first filtered our log2(TPM + 1) gene expression matrix to keep only genes with differential expression P value < 0.05 and absolute log-fold change > 0.5, thus limiting our analysis to genes potentially involved in response. We further filtered out genes with sparse or low expression, that is more than 10% NA or zero values, or in the bottom 10% of mean expression. We transformed the values to fold changes by subtracting the median for each gene, then obtained the Spearman correlation matrix of these fold changes, and performed hierarchical clustering while varying the number of clusters (K) from 2 to 10 and repeating 500 iterations for each K value. We then obtained consensus matrices for each K (calculating the number of times samples clustered together in the 500 iterations), summed these matrices across all K values, and normalized the resulting matrix by the number of iterations. Using B-NMF with a half-normal prior, this matrix was used to decide on the optimal number of clusters. Using this empirically determined value of K, we then applied the B-NMF algorithm to the original log2TPM gene expression matrix. In this case, the gene expression matrix is approximated by W*H, where H is the cluster membership matrix and W is the gene weight matrix. We used the W matrix to narrow down the genes most closely associated with each cluster, keeping only genes in the top 50% of normalized weights for each cluster, as well as those with the largest difference between within-cluster versus outside-cluster expression. Using this reduced marker gene list, we classified the remaining samples in cohort 2 into our three-cluster scheme. Of note, given re-annotation of the RNA sample from patient SU2CLC-DFC-DF0732 as having been post-treatment, this specimen was removed from our analysis (Supplementary Note and Supplementary Fig. 1). We used this same procedure to define the TI subtypes, with the exception of initially filtering to keep high-variance genes (instead of keeping genes of interest from the differential expression analysis, as in the M subtypes). We similarly used TI marker genes to classify additional samples.

52A collection of the top clinical, genomic and transcriptomic predictors identified in the SU2C-MARK cohort or published previously as relevant to antitumor immunity were first compared across samples in the SU2C-MARK cohort. Unsupervised hierarchical clustering was performed on this set, identifying four broad clusters that were ultimately designated wound healing, immune activation/exhaustion, neoantigens and others. As validation of these predictor classes, recalculation of these features was performed in TCGA data by combining publicly available mutation calling and RNA-seq data for the TCGA LUAD and TCGA LUSC cohorts (combined n = 1018). Unsupervised hierarchical clustering was again performed, and feature membership was compared to assignments made earlier from analysis within the SU2C-MARK dataset. As with other sections described here, integrative analysis was performed with Python (v3.7) and R (v3.4).

53Using single-cell data from a previously published NSCLC cohort, we performed preprocessing, integration and Leiden clustering in Scanpy (v1.9.1) to identify distinct cell types. For preprocessing, we filtered counts to cells with at least 200 genes, and then filtered out genes that were observed in fewer than 50 cells. Further filtering was performed on cells with between 1,000 and 8,000 genes, total counts between 3,000 and 100,000, percent of mitochondrial counts less than 15%, and percent of ribosome counts less than 20%. Cell cycle effects were regressed out using Scanpy, and samples were then integrated using Harmony (harmony2019). The cell type of the Leiden clusters was annotated based on gene markers described in Wu et al. as well as canonical IHC cancer subtype markers. Clusters were assigned one of the cell types alveolar, B cell, cancer, endothelial, epithelial, fibroblasts, myeloid or T cell based on these expression markers. Metagene expression level was calculated as the mean expression of the gene markers that comprised the metagene. For signatures M-2 and TI-1, the top ten genes by weight were selected. Of note, in some cases, single genes or individual genes in a signature did not pass filtering or were not detected, and therefore were not plotted/included in a given metagene.

54Single-feature survival analysis was performed using progression-free and overall survival data with censoring as described above. For the MSK impact cohort, patients with alterations in ATM found on panel sequencing who also had received checkpoint blockade therapy were included in the cohort. For integrative analysis across the feature list, the top two genomic features from each correlation cluster were selected for PFS analysis as follows: the monocyte/macrophage score, the hMono3 score, dedifferentiated signature TI-1, immune-activated signature M-2, TMB, TMB indel, ATM Mutation and TERT amplification. Participants were binned into high and low categories for each feature (using 0 as a cut point for z-score features, cluster identity for signatures, median for mutation burden features and the presence or absence of alteration/copy gain for single-gene features). FDR values were subsequently computed from the nominal P values obtained via the log-rank test using the Benjamini–Hochberg method. A complete list of the log-rank test results including median PFS for each subgroup is provided in Supplementary Table 31.

55This study was designed as a retrospective immunogenomic analysis of biospecimens from NSCLC patients receiving checkpoint blockade in the advanced setting. As such, no statistical method was used to predetermine the sample size. Patients who did not have at least one pretreatment whole exome or RNA-seq sample that passed QC following library construction or alignment were excluded from the analysis (as described above). There was no randomization or stratification performed during described analyses, and investigators were not blinded to participant outcomes during primary data analysis.

56Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.

57Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at 10.1038/s41588-023-01355-5.