Single-cell analysis of human non-small cell lung cancer lesions refines tumor classification and patient stratification

Abstract

Immunotherapy is a mainstay of non-small cell lung cancer (NSCLC) management. While tumor mutational burden (TMB) correlates with response to immunotherapy, little is known about the relationship between the baseline immune response and tumor genotype. Using single-cell RNA sequencing, we profiled 361,929 cells from 35 early-stage NSCLC lesions. We identified a cellular module consisting of PDCD1⁺CXCL13⁺ activated T cells, IgG⁺ plasma cells, and SPP1⁺ macrophages, referred to as the lung cancer activation module (LCAM^hi). We confirmed LCAM^hi enrichment in multiple NSCLC cohorts, and paired CITE-seq established an antibody panel to identify LCAM^hi lesions. LCAM presence was found to be independent of overall immune cell content and correlated with TMB, cancer testis antigens, and TP53 mutations. High baseline LCAM scores correlated with enhanced NSCLC response to immunotherapy even in patients with above median TMB, suggesting that immune cell composition, while correlated with TMB, may be a nonredundant biomarker of response to immunotherapy.

Keywords: CITEseq; NSCLC; dendritic cells; high dimensional profiling; immunotherapy; macrophages; tumor cell atlas; tumor microenvironment; tumor mutational burden; tumor-associated myeloid cells.

Figure 1.. scRNA-seq and CITE-seq establish the diversity of transcriptional states in the tumor microenvironment

(A) Study overview. Resected tumor tissue and nLung were digested to single-cell suspensions, enriched for CD45⁺ cells, and subjected to single-cell assays shown. (B) Clinical data of patients indicating summary pathological stage, smoking history, histological diagnosis, and sex. (C) Expression of cell-type marker genes across immune scRNA-seq clusters (MNP, mononuclear phagocyte; pDC, plasmacytoid dendritic cell). Heatmap shows the number of unique molecular identifiers (UMI) per cell. Clusters are shown using even sampling of cells from eight patients analyzed by CITE-seq. Cells were downsampled to 2,000 UMI/cell. (D) Lineage-defining surface marker expression as measured by CITE-seq. Single cells correspond directly to (C). CITE-seq count values were first quantile-normalized across patients, then row-normalized across cells in the heatmap. (E) Cluster frequencies among immune cells in 35 tumor and 29 matched nLung samples. Clusters correspond to (C) and (D). *p < 0.05, **p < 0.01, ***p < 0.001 (Wilcoxon signed-rank test with Bonferroni correction; N = 29 matched tissue pairs). Raw and adjusted p values are tabulated in Table S6. (F) Euclidean distances between sample pairs among nLung only (nLung-nLung), tumor only (tumor-tumor), or between nLung and tumor (tumor-nLung). Distances between pairs of patient-matched samples were excluded. ***p < 0.001, Wilcoxon rank-sum test. (G) Log-ratios between cell-type frequencies in tumor and nLung in the present cohort (x axis) or that of Lambrechts et al. (2018) (y axis). Clusters were grouped by cell-type annotation. Crosses represent mean ± SEM. See also Figure S1.

Figure 2.. Intratumoral DCs comprise expanded DC3 and express an LCH-like signature

(A) Expression of genes discriminating scRNA-seq DC clusters showing number of UMI per cell, showing cells from four patients analyzed by CITE-seq with the DC panel shown in (B). (B) CITE-seq measurement of DC surface markers on cells corresponding directly to those in (A). (C) Differences between tumor and nLung of DC frequencies; *p < 0.05, **p < 0.01, ***p < 0.001 (Wilcoxon signed-rank test with Bonferroni correction; N = 25 matched tissue pairs with >50 DCs observed). (D) Average expression of LAMP3 and CD274 in DC clusters. (E) MICSSS of DC-LAMP⁺/PD-L1⁺ DC and T cells in a TLS. (F) Expression of follicular DC marker MYH11 in TLS in a section adjacent to (E). (G) Expression among putative monocyte-derived MΦ (MoMΦ), CD14⁺ monocytes, and DC of monocyte, cDC2, and MoMΦ gene signatures (see Figures S2D and S2E). Top 20 genes from each score are shown among cells evenly sampled by cell type (left) and corresponding summary scores. Cells were ordered by the ratio of monocyte:cDC2 summary scores and were downsampled to 2,000 UMI. (H) Box plots showing average expression of LCH-like signature genes in DC populations across patients. See also Figure S2.

Figure 3.. Tumors exclude AMΦ and exhibit a diversity of MoMΦ populations

(A) Average cluster expression of key monocyte, MΦ, and cDC2-like genes. (B) CITE-seq surface marker detection. Cells are evenly sampled per cluster from the four patients analyzed with the panel shown. (C) Histograms of gene module scores per cell type (see also Figures S3E–S3H). (D–F) Expression among CD14⁺ monocytes, MoMΦ, and AMΦ of cell-type-specific gene scores (see Figure S3I). Cells are plotted by AMΦ (x axis) and MoMΦ (y axis) score, and cell annotations are indicated by colored dots or contour plots (D). Cells are plotted similarly and colored by CD14⁺ monocyte score (E), or by individual gene expression (F). (G) Differences between tumor and nLung of monocyte and MΦ frequencies; *p < 0.05, **p < 0.01, ***p < 0.001 (Wilcoxon signed-rank test with Bonferroni correction; N = 29 tissue pairs). (H) Average expression of transcripts encoding secreted factors across MNP cell types. See also Figure S3.

Figure 4.. CITE-seq and TCR analysis of the lymphoid compartment

(A) Expression of genes discriminating clusters of T cells, showing cells from two patients analyzed by CITE-seq with the panel shown in (B). (B) CITE-seq measurement of T cell-surface markers on cells corresponding directly to those in (A). (C) Differences between tumor and nLung of NK and T frequencies; *p < 0.05, **p < 0.01, ***p < 0.001 (Wilcoxon signed-rank test with Bonferroni correction; N = 29 matched tissue pairs). (D and E) Phenotypic distribution of T cells among tissue-stratified clonotypes. Frequencies of unique TCRs observed by scTCR-seq in nLung (x axis) or tumor in a representative patient (D). In (E), cells were first grouped by TCR tissue tropism categories as defined in (D); for three patients, the phenotypic makeup of the cells with unique TCRs, tissue-specific TCRs, or TCRs shared across tissues is plotted for nLung (i) and tumor (ii) tissues as a percentage of cells with similarly tissue-distributed TCRs. Each patient is indicated by shape. (F) Expression of genes discriminating clusters of B cells and plasma cells, showing cells from four patients analyzed by CITE-seq with the panel shown in (G). (G) CITE-seq measurement of cell-surface markers on cells corresponding directly to those in (A). See also Figure S4.

Figure 5.. Patients stratify along the LCAM axis

(A) Correlation of lineage-normalized cell-type frequencies. Analysis includes 26 similarly processed tumors (10X Chromium V2, CD45⁺ magnetic bead enrichment). (B) Lineage-normalized LCAM^hi and LCAM^lo cell-type frequencies among samples from Mount Sinai, Lambrechts et al. (2018), and Zilionis et al. (2019) (50 tumor, 36 matched nLung). Dataset of origin is indicated by color bar below. (C) Immune lineage frequencies with columns corresponding to patient ordering in (B). (D) Scores measuring geometric averages of lineage-normalized LCAM cell-type frequencies. (E–G) MICSSS demonstrating infiltration of CD138⁺ plasma cells in an LCAM^hi lesion (right) relative to LCAM^lo (E), plasma cells arising outside of TLS (F), and colocalization of PD1⁺ T cells, CD68⁺ MΦ, and CD138 plasma cell LCAM^hi tumor stroma (G). (H–M) Log₂ ratios of ligand-receptor (LR) intensities between tumor and nLung of LCAM^hi patients, (“LR ratio;” y axis) and LCAM^lo patients (x axis). All interactions among T cells, B cells, MΦ, DC3, cDC, and monocytes, colored by indication of significance (permutation test, H). Dashed diagonal line indicates unity. (I–M) Same data as in (H), highlighting in bold LR ratios for interactions between T cell ligands and B cell receptors (I), T cell ligands and cDC receptors (J), MΦ ligands and T cell receptors (K), DC3 ligands and T cell receptors (L), and cDC ligands and T cell receptors (M). Labeled interactions are plotted in red. See also Figure S5.

Figure 6.. Tumor features related to the LCAM immune response

(A) Normalized expression of LCAM^hi and LCAM^lo bulk-RNA signature genes in TCGA lung adenocarcinoma (LUAD) dataset. Cell-type association with sets of genes for each signature is shown. Patients are sorted along y axis by ensemble LCAM score. (B) Scatterplot of the ensemble LCAM score (y axis) with the difference of CAF and normal fibroblast signature scores (x axis). (C) Scatterplot of LogTMB and ensemble LCAM score. Patients are divided into those with (black) and without (red) presence of a smoking-related mutational signature. Black and red lines indicate linear regression relationships computed over each group of patients independently (r = 0.36; p = 3.3 × 10⁻⁴ in the undetected smoking signature group; r = 0.36; p = 1.6 × 10⁻¹² in the detected signature group).

Figure 7.. LCAM is modified by driver mutations and associates with checkpoint response

(A and B) Box plots showing either the ensemble LCAM score (A) or TMB (B) among TCGA LUAD patients, divided by combinations of driver mutations. (C and D) Histograms of residuals of the regression of the ensemble LCAM score on the LogTMB, with patients stratified by TP53 (C) or KRAS (D) mutational status (two-sided t test). (E) Scatterplot of LogTMB and ensemble LCAM score, plotted for LUAD (black) or LUSC (red). (F) Correlation between TMB and the ensemble LCAM score for all TCGA cancers, showing Spearman ρ (x axis) and the −log₁₀(p_adj) after Bonferroni correction (y axis). BRCA, breast invasive carcinoma; STAD, stomach adenocarcinoma; BLCA, urothelial bladder carcinoma; UCEC, uterine corpus endometrial carcinoma; COAD, colon adenocarcinoma; THYM, thymoma; PRAD, prostate adenocarcinoma; CESC, cervical squamous cell carcinoma; ACC, adrenocortical carcinoma. (G–J) Stratification of response to immunotherapy by TMB and ensemble LCAM in the POPLAR trial. Kaplan-Meier curves showing probability of progression-free survival (PFS) among patients with high LCAM ensemble score (top 25%) versus low LCAM ensemble score (bottom 75%) in the arm receiving atezolizumab (anti-PD-L1; G) or docetaxel chemotherapy (H). Forest plots demonstrating PFS hazard ratio estimates for LCAM and TMB status in a categorical multivariate Cox proportional hazards analysis in patients treated with atezolizumab (left) or docetaxel (I). Kaplan-Meier curves showing probability of PFS among patients with high LCAM ensemble score (top 25%) versus low LCAM ensemble score (bottom 75%) in patients receiving atezolizumab who had TMB greater than the median TMB (J). See also Figure S7.