Association of High Tumor Mutation Burden in Non–Small Cell Lung Cancers With Increased Immune Infiltration and Improved Clinical Outcomes of PD-L1 Blockade Across PD-L1 Expression Levels

Key Points Question Is tumor mutation burden (TMB) associated with improved outcomes of programmed cell death–1 (PD-1)/programmed death ligand–1 (PD-L1) inhibition across PD-L1 expression levels in non–small cell lung cancer (NSCLC)? Findings In this cohort study of 1552 patients with NSCLC, the group with high TMB had improved response rates and survival after receiving PD-1/PD-L1 inhibition therapy across PD-L1 expression subgroups compared with the group with low TMB. High TMB levels were associated with increased CD8-positive T-cell infiltration and distinct immune response gene expression signatures. Meaning These findings suggest that in NSCLC, a high number of nonsynonymous tumor mutations is associated with immune cell infiltration and inflammatory T-cell expression signatures, leading to increased sensitivity to PD-1/PD-L1 inhibition across PD-L1 expression subgroups.


Dana-Farber Cancer Institute cohort
Patients at the Dana-Farber Cancer Institute who consented to institutional review board-approved protocols DF/HCC 02-180, 11-104, 13-364, and/or 17-000 which allowed for conducting translational research and tumor next-generation sequencing, respectively, were included.

Memorial Sloan Kettering Cancer Center cohort
Patients at the Memorial Sloan Kettering Cancer Center were included if they had advanced NSCLC which underwent tumor next generation sequencing and if they had also consented to institutional review board-approved protocols.

Stand Up to Cancer/Mark Foundation
Patients included in the Stand Up to Cancer/Mark Foundation cohort were enrolled if they had advanced NSCLC which was treated with PD-1/PD-L1 inhibitors.

Tumor genomic profiling and somatic variant calling in the DFCI and MSKCC cohorts
Tumor genomic profiling and somatic variants were performed using clinically validated bioinformatics pipelines 1,2 . Sequence reads were aligned to reference sequence b37 edition from the Human Genome Reference Consortium using bwa (http://biobwa.sourceforge.net/bwa.shtml), and further processed using Picard (version 1.90, http://broadinstitute.github.io/picard/) to remove duplicates and Genome Analysis Toolkit (GATK) to perform localized realignment around indel sites. Single nucleotide variants were called using MuTect v1.1.4, insertions and deletions were called using GATK Indelocator, and variants were annotated using Oncotator. In the DFCI cohort, to filter out potential germline variants, the standard pipeline removed SNPs present at >0.1% in Exome Variant Server, NHLBI GO Exome Sequencing Project (ESP) (URL: http://evs.gs.washington.edu/EVS/), present in dbSNP, or present in an in-house panel of normals, but rescues those also present in the COSMIC database. For this study, variants were further filtered by removing variants present at >0.1% in the gnomAD v.2.1.1 database or were annotated as Benign or Likely Benign in the ClinVar database (PMID: 32461654, 29165669). In the MSKCC cohort, patient-matched normal DNA was used to filter out germline variants, as previously described.

DFCI OncoPanel and MSK-IMPACT
Tumor mutational burden (TMB), defined as the number of somatic, coding, base substitution, and indel mutations per megabase (Mb) of genome examined, was determined using the OncoPanel (Dana-Faber) and MSK-IMPACT (MSKCC) NGS platforms, as previously described 1 Stand Up to Cancer cohort whole exome sequencing DNA was extracted from FFPE tumor specimens and either matched normal whole blood or in cases where this was unavailable, adjacent normal FFPE specimens. Extraction was performed using the Qiagen AllPrep DNA/RNA Mini Kit (cat# 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. Mean target coverage for tumor and normal samples were 150X and 80X, respectively. Tumor mutational burden was defined as the number of non-synonymous base substitutions, indel mutations per megabase of genome examined, using an exome size of 35.8 Mb.

Tumor mutational burden normalization across different platforms
TMB distributions were harmonized between the two platforms by applying a normal transformation followed by standardization to Z-scores, as previously described 3 . Briefly, power transformations were first used to normalize cohort-specific TMB distributions; second, Tukey's ladder of powers in the rcompanion package was used to identify the optimal transformation coefficient. Third, the normalized distributions were then standardized into z scores by subtracting the transformed distribution mean and dividing by the standard deviation.

Cell subset analysis from the TCGA dataset
To perform cell type enrichment analyses, RNA sequencing data from the LUAD and LUSC TCGA cohort were deconvoluted to estimate cell subsets using the xCell package. xCell estimates the abundance scores of 64 cell types, including adaptive and innate immune cells, hematopoietic progenitors, epithelial cells, and extracellular matrix cells, based on single sample gene set enrichment analysis (ssGSEA) data 4 . Gene expression values (RSEM V2) were converted into Z-scores and used to compute cell type enrichment scores with the xCellAnalysis function. Statistical significance of differential cell type enrichment between cohorts was estimated with Wilcox Rank Sum test.

Multiplexed immunofluorescence (ImmunoProfile)
Multiplexed immunofluorescence (mIF) was performed on samples from the Dana-Farber Cancer Institute by staining 5-micron formalin-fixed, paraffin-embedded whole tissue sections with standard, primary antibodies sequentially and paired with a unique fluorochrome followed by staining with nuclear counterstain/4′,6-diamidino-2-phenylindole (DAPI), as previously described 5 . All samples were stained for PD-L1 (clone E1L3N), PD-1 [clone EPR4877(2)], CD8 (clone 4B11), FOXP3 (clone D608R), Cytokeratin (clone AE1/AE3), and DAPI (nuclear counterstain). Each sample had a single slide stained and scanned at 20x resolution by a Vectra Polaris imaging platform. Regions of Interest (ROIs) were defined for each image, and only these regions were used for quantitative image analysis currently. Within each ROI, InForm Image Analysis software (PerkinElmer/Akoya) was run to phenotype and score cells based on biomarker expression. A custom script quantified the number/percentage of cells which are positive for relevant biomarkers in specific tissue regions. Each ROI was divided into one or more of these defined regions: intra-tumoral (IT), which was defined as the region of the slide consisting of tumor beyond the tumor-stroma interface; tumor-stroma interface (TSI), which was defined as the region within 40 microns to either side of the defined border between tumor and stroma; and total (IT + TSI). Cell count was calculated per ROI and averaged (unweighted) across ROIs, reported as count per millimeter squared +/-standard error. Statistical significance of differential cell type enrichment between groups was estimated with Wilcox Rank Sum test.

Gene expression analysis
Gene expression data were downloaded from the Firehose website (TCGA Firehose Legacy version) while somatic mutation data were downloaded from cBioPortal website (cbioportal.org). The RSEM V2 values were used to represent gene expression and genes with counts less than 10 were filtered out. Gene expression profiles were analyzed according to TMB categories. Median expression within each group was used to estimate expression fold-change (FC) to minimize the possible impact of outlier samples. Gene differential expression analyses across TMB subgroups were conducted using R package DESeq2. Pvalues were corrected for multiple hypothesis testing via false discovery rate (FDR) adjustment. Fold-change threshold of an absolute value greater than 1.5 and FDR adjusted P-value threshold less than 0.1 were utilized to identify differentially expressed genes. Pathway enrichment analyses were conducted separately for up-and down-regulated genes using Molecular Signatures Database (MSigDB) collections.

Statistical analysis
Clinicopathologic data and immunotherapy response data were abstracted from the electronic medical record. Overall response rate was determined by a blinded radiologist using Response Evaluation Criteria In Solid Tumors, version (RECIST) 1.1. Progression-free survival was determined from the start date of immunotherapy until the date of disease progression or death, and overall survival was calculated from the date of diagnosis of advanced NSCLC until the date of death. All p-values are twosided and confidence intervals are at the 95% level. Overall survival among patients who never received PD-(L)1 inhibition was calculated from the date of the start of systemic therapy for advanced disease, other than immunotherapy. TMB comparisons were computed using the Mann-Whitney U test or the Kruskal-Wallis test, when appropriate. Linear correlations were evaluated using Spearman's test, and categorical variables were evaluated using Fisher's exact test. Event-time distributions were estimated using Kaplan-Meier methodology. Log-rank tests were used to test for differences in event-time distributions, and Cox proportional hazards models were fitted to obtain estimates of hazard ratios in univariate and multivariate models. The proportional hazards assumption was assessed with Schoenfeld residuals. All P-values are 2-sided and confidence intervals are at the 95% level, with significance pre-defined to be at <0.05. Multiple comparison correction was performed using the Benjamini-Hochberg procedure. Missing values were handled using inverse probability weighting (IPW) and multiple imputation approaches using R package MICE, as previously described. All statistical analyses were performed using R version 3.6.3.

TMB cut-off identification and validation
To identify and validate TMB thresholds associated with immunotherapy efficacy, an unbiased recursive partitioning algorithm was used to investigate an optimal grouping of TMB with respect to the objective response rate to immune checkpoint inhibition in a discovery cohort comprised of patients from the MSKCC cohort, using the partykit function in R, as previously described 6 . A 10-fold cross-validation method was used to train and measure the performance of the model using the caret function in R, as previously described 7 . The threshold identified was validated in two independent cohorts of patients treated with PD-(L)1 blockade in the DFCI and SU2C/Mark Foundation cohorts, following TMB harmonization across platforms, as described above and as previously described 3 . As PD-L1 tumor proportion score (TPS) is an important predictor for ICI efficacy, we applied both Inverse probability weighting (IPW) and multiple imputation approaches using R package MICE to address the potential selection bias arising from the PD-L1 TPS missingness. Variables used for multiple imputation and to calculate the weights for PD-L1 TPS missingness included sex, age, ECOG performance status, histology, smoking status, and line of therapy for ICI. IPW and multiple imputation were conducted separately in each cohort, and the multivariable analyses results were pooled based on 5 repeated complete imputed datasets.