Pan-cancer profiling of tumor-infiltrating natural killer cells through transcriptional reference mapping

Abstract

The functional diversity of natural killer (NK) cell repertoires stems from differentiation, homeostatic, receptor-ligand interactions and adaptive-like responses to viral infections. In the present study, we generated a single-cell transcriptional reference map of healthy human blood- and tissue-derived NK cells, with temporal resolution and fate-specific expression of gene-regulatory networks defining NK cell differentiation. Transfer learning facilitated incorporation of tumor-infiltrating NK cell transcriptomes (39 datasets, 7 solid tumors, 427 patients) into the reference map to analyze tumor microenvironment (TME)-induced perturbations. Of the six functionally distinct NK cell states identified, a dysfunctional stressed CD56^bright state susceptible to TME-induced immunosuppression and a cytotoxic TME-resistant effector CD56^dim state were commonly enriched across tumor types, the ratio of which was predictive of patient outcome in malignant melanoma and osteosarcoma. This resource may inform the design of new NK cell therapies and can be extended through transfer learning to interrogate new datasets from experimental perturbations or disease conditions.

Fig. 1. NK cell differentiation at the transcriptional level.

a, Integration process of scRNA-seq data of NK cells from 12 donors and 4 different laboratories using scVI showing a UMAP based on the scVI latent representation, followed by a UMAP based on the diffusion map components. b, AUCell scores of gene signatures for CD56^bright and CD56^dim NK cell subsets. c, UMAP representation of five sorted subsets from a donor with an adaptive expansion (left) and a donor without an adaptive expansion (right). d, Heatmap depicting accuracy of our prediction model for subset annotation tested on the held-out 15% of cells from the subset-specific dataset (two donors). e, UMAP of the scANVI representation of both bulk and sorted NK cells, showing original annotation of NK cells (12 donors, left) and subset labels predicted (right) using the scANVI model trained with sorted subset data (2 donors). f, Dot plots showing the top three up- and downregulated genes between all pairs of subsets (x and y axes) as identified by the differential expression module in scANVI. These top genes were then visualized across all NK cell subsets within the differentiation spectrum (x axis), to highlight the continuous nature of NK cell differentiation. g, Diffusion map representation showing the predicted subset labels for the bulk data (top) and depicting Leiden clustering of the 12 donor NK cell dataset (bottom). h, Heatmap showing distribution of our annotated 12 donor NK cell subsets over the 5 Leiden clusters. i, Frequency (freq.) of annotated late CD56^dim and adaptive NK cell subsets in donors with and without an adaptive NK cell expansion. Int., intermediate. Source data

Fig. 2. GRNs defining conventional and adaptive NK cell fates.

a, UMAP representation highlighting the starting cell (blue) with the lowest value CD56^dim signature score and the two terminal cells (orange) as predicted by Palantir. b, UMAP representation of the data from the sorted subsets (two donors) showing the RNA velocity vector field as a stream plot and the inferred pseudotime. c,d, PAGA graph with directionality and transitions from RNA velocity analysis for the sorted subsets (2 donors) (c) and subset-inferred bulk donors (12 donors) stratified based on the presence or absence of adaptive expansion (d). e, Gene trends clustered into five overall trends of expression along pseudotime, showing expression of KLRC2, CD52 and IL32 in both terminal fates (pink, conventional fate; orange, adaptive fate). f, Inferred GRNs where dominant TFs for each trend are highlighted. g, Selection of regulons showing differential expression over pseudotime within the conventional and adaptive fate.

Fig. 3. Pan-cancer atlas of healthy tissue-resident and solid TiNK cells.

a, Graphic overview of healthy tissue datasets included in the analysis, with the number of donors denoted in brackets. b, UMAP representation showing integration of all healthy tissue datasets. c, Graphic overview of solid tumor datasets included in the analysis, with the number of donors denoted in brackets. d, UMAP representation showing integration of all solid tumor datasets. e, Scoring of tissue-residency signatures (sig.) in PB-NK cell subsets, as well as CD56^bright- and CD56^dim-annotated TrNK and TiNK subsets. f, UMAP representation showing integration of subset-annotated PB-NK, TrNK and TiNK cells. g,h, PAGA graphs (g) and connectivity heatmap (h) showing connectivity of PB-NK, TrNK and TiNK subsets across all tissues/tumor types, with individual tissues/tumor types highlighted (g). The scale represents gene set activity computed by AUCell (e). Panels a and c created with BioRender.com.

Fig. 4. Cellular composition of pan-cancer cell atlas and subset distribution of TiNK cells.

a, Heatmap depicting changes in immune subset proportion in tumor samples compared with healthy tissue samples at the pan-cancer level and within individual tumor types. b,c, Proportions of major immune subsets within healthy tissue (b) and tumor samples (c). d,e, Predicted subset annotations of CellTypist-identified NK cells in healthy tissue (d) and tumor samples (e) compared with annotated PB-NK cells. f,g, Frequency of CD56^bright NK cells (f) and relative frequency of subsets (g) identified by flow cytometry in a cohort of healthy blood donors (n = 19) and central tumor samples from patients with NSCLC (n = 25), from 23 independent experiments. Data were analyzed using two-sample Student’s t-test with Bonferroni’s correction (a) and a two-tailed Mann–Whitney U-test (f): *P < 0.05, ^**P < 0.01, ^***P < 0.001, ^****P < 0.0001. The bar graph in f represents the mean ± s.d, with the actual P value indicated (P < 0.0001). Source data

Fig. 5. Distinct cellular states of NK cells identified in pan-cancer atlas.

a, UMAP depicting neighborhood (Nhood) groups identified by Milo and computed using the scVI representation. b, Beeswarm plot depicting differential abundance of neighborhoods (TiNK versus Ref-NK enriched). Colored neighborhoods are differentially abundant at a false recovery rate (FDR) of 0.1. c, Pie charts showing distribution of NK subsets across neighborhood groups annotated using our annotation model (Fig. 1). d, Expression of dominant TF regulons of NK cell differentiation across NK cell states (neighborhood groups). e, Expression of TF regulons uniquely expressed across cellular states. f, Graphic representation of cellular states. g,h, Volcano plots depicting DEGs between group 1 versus group 2 (g) and group 3 versus group 4/5/6 (h) cellular states. Differential expression analysis was performed using the findNhoodGroupMarkers method within the miloR package. Counts were aggregated per sample; groups were compared using edgeR and the adjusted P values were used for the plots. i, Scoring of pathway gene signatures in NK cell states. Func., function; homeo., homeostasis. j–n, Dot plots depicting selected genes belonging to stress response (j), immune suppression (k), metabolism (l), cytotoxicity (m) and chemokine/cytokine secretion (n). o, Pie charts depicting distribution of NK cell states in blood, tissues and tumors. Volcano plots: log(fold-change) cutoff at 0.5, P < 0.05. The scale represents regulon activity (d and e) or gene set activity (i) computed by AUCell.

Fig. 6. Intercellular communication of distinct cellular states in the TME.

a, Selected predicted incoming signaling pathways involving TiNK cells common across tumor type, identified by CellChat. b, Violin plots showing expression of receptors for the MIF, COLLAGEN and LAMININ communication pathway in NSCLC. c, Circle plot depicting predicted incoming signaling via CD74, CXCR4 and CD44 expression (NSCLC). d, Violin plots showing expression of receptor and ligand for the MHC-I communication pathway in NSCLC. e, Selected predicted outgoing signaling pathways involving TiNK cells across tumor type. f, Heatmap depicting interaction role of individual cell populations in CCL, PARs and IFN-II signaling pathways in NSCLC based on network centrality analysis. g, Dot plot depicting GZMA expression (ligand for PARs) and IFNG expression (ligand for IFN-II) in NK cell states across tumor type. h, Frequency of granzyme A⁺CD56^dim NK cells in healthy blood donors (n = 6) and patients with NSCLC (n = 11) from 9 independent experiments. i, Geometric mean fluorescence intensity (gMFI) of granzyme A⁺CD56^dim NK cells in healthy blood donors and patients with NSCLC (n = 7) from six independent experiments. j,k, Degranulation (CD107a) (j) and granzyme B release (k) of CD56^bright and CD56^dim NK cells against A549 target cells pre-treated with (black box) and without (white box) IFN-γ (24 h) in the presence (black box) or absence (white box) of α-NKG2A antibody (E:T 1:1, 4 h, n = 5 biologically independent replicates from one experiment). Data were analyzed using two-way analysis of variance followed by Šidák’s multiple-comparison test (j and k), two-tailed Mann–Whitney U-test (h) or two-tailed Wilcoxon’s test (i). All bar graphs represent the mean ± s.d, with the actual P values indicated. Source data

Fig. 7. Distinct cellular states in spatial RNA-seq and association with patient outcome.

a, Deconvoluted spatial RNA-seq images from SKCM, NSCLC and GBM at the level of immune populations. b, Pie charts depicting compositional analysis of major immune populations from scRNA-seq datasets and spatial-seq datasets for SKCM, NSCLC and GBM samples. c,d, Annotation of CD56^bright and CD56^dim NK cell subsets (c) and the six cellular states of NK cells (d) in SKCM. e–g, Dot plots depicting selected genes belonging to NK cytotoxicity (e), stress response (f) and immunosuppression (g) scored across NK cell states in spatial-seq data from SKCM. h, Kaplan–Meier survival curves showing association of high/low group 1/3 gene signatures with patient outcome across tumor types. Survival analysis was performed using Cox’s proportional hazards model; P values were computed using the log(rank) test.

Extended Data Fig. 1. Peripheral blood NK cell subsets and sorting strategy.

(a, b) AUCell scores of gene signatures for CD56^bright and CD56^dim NK regulones (a) and proteomes (b). (c, d) Sorting strategy for phenotypically defined functional PB-NK cell subsets sequenced in one donor with (c) and one without (d) an adaptive expansion. (e) Heatmap depicting similarity between our five annotated transcriptional NK cell subsets (y-axis) and the Meta-NK defined NK subsets (x-axis). The scale represent gene set activity calculated by AUCell (a-b, e).

Extended Data Fig. 2. RNA velocity.

(a) Graphical depiction of inferring RNA velocity based on spliced vs unspliced transcripts. (b) RNA velocity plots for ZEB2 and BCL11B transcripts stratified by subset annotation in donors with an adaptive expansion.

Extended Data Fig. 3. Healthy tissue dataset annotation using CellTypist.

(a) Heatmap depicting expression of signature genes of the main immune populations annotated by CellTypist across all tissue samples. (b-f) UMAP representation showing integration of all healthy tissue datasets, prostate (b), lung (c), pancreas (d), skin (e), breast (f), with individual cell subtypes annotated using CellTypist.

Extended Data Fig. 4. Solid tumor dataset annotation using CellTypist.

(a-g) UMAP representation showing integration of all solid tumor datasets, PRAD (a), NSCLC (b), SKCM (c), PAAD (d), BRAC (e), GBM (f), SARC (g), with individual cell subtypes annotated using CellTypist.

Extended Data Fig. 5. Tissue-residency scoring of NK cells.

(a, b) UMAP representation showing integration of all healthy tissue (a) and solid tumor (b) datasets, with lymphocytes populations visualized. (c, d) IL7R expression in annotated NK cells (CD56^bright, CD56^dim) and ILCs (ILC2, ILC3) in tissues (c) and tumors (d). (e, f) Dotplots depicting expression of genes defining the literature-TR and atlas-TR signatures in CD56^bright and CD56^dim subsets in healthy blood and across all tissue types (e) and stratified by individual tissues (f). (g) Tissue-residency scoring (atlas-TR) of CD56^bright and CD56^dim annotated NK cells in individual tissue and tumor types. The scale represents gene set activity calculated by AUCell (g).

Extended Data Fig. 6. Phenotyping of NK cells in NSCLC patient samples.

(a) Gating strategy for CD56^bright and CD56^dim NK cells in healthy donor (PBMC) and NSCLC samples (Tumor). (b-e) Representative plots of CD56^dim NK cells (b) and quantification of NGK2A (c), KIR (d) and CD57 (e) expression in PBMC (n = 19) and NSCLC samples (n = 25), from 23 independent experiments. Data were analyzed using two-tailed Mann-Whitney test (c-e). All bar graphs represent the mean ± s.d. Actual p values are indicated. Source data

Extended Data Fig. 7. Characterization of cellular states of NK cells identified in pan-cancer cell atlas.

(a) Beaswarm plot depicting differential abundance of TiNK or Ref-NK (PB-NK, TrNK) enriched neighborhoods, clustered based on subset annotation of individual neighborhoods. (b) TiNK fraction of cells in neighborhoods within each neighborhood group. The boxplot indicates the median with the interquartile range (IQR), whiskers extend to the farthest point within 1.5 times the IQR from the box. n is the number of neighborhoods in each group: group 1, n = 382; group 2, n = 1261; group 3, n = 1239; group 4, n = 871; group 5, n = 1427; group 6, n = 1752. (c-e) Volcano plots depicting differentially expressed genes (DEGs) between Group 4 vs. Group 3/5/6 (c), Group 5 vs. Group 3/4/6 (d), Group 6 vs. Group 3/4/5 (e). Differential expression analysis was performed using the findNhoodGroupMarkers method within the miloR package. Counts were aggregated per sample; groups were compared using edgeR and the adjusted p-values were used for the plots. (f-i) Gene set enrichment analysis (GSEA) for DEGs identified between Group 1 vs Group 2 (f), Group 3 vs Group 4/5/6 (g), Group 5 vs Group 3/4/6 (h) and Group 6 vs Group 3/4/5 (i). Volcano plots: log fold change cutoff at 0.5, p < 0.05. GSEA plots: p value cutoff 0.5 (red line).

Extended Data Fig. 8. Intercellular communication in TME.

(a-g) Scatterplot depicting incoming and outgoing interaction strength of individual cell types in BRAC (a), PAAD (b), PRAD (c), NSCLC (d), SARC (e), SKCM (f), GBM (g) as identified by CellChat. (h) Violin plots showing expression of ligands for the CCL (CCL3, CCL5) communication pathway in NSCLC.

Extended Data Fig. 9. In vitro validation of IFNG-HLA-E-KLRC1 axis in NSCLC.

(a) Representative histogram of HLA-E expression of A549 cells pre-treated (24 h) with and without IFNγ. (b-d) Viability (b), frequency (c) and geometric MFI (d) of HLA-E + A549 cells (n = 4, biological replicates from two independent experiments). (e) Gating strategy and representative contour plots for functional readout of CD56^bright and CD56^dim NK cells against A549 target cells pre-treated with and without IFNγ (24 h) in presence and absence of α-NKG2A antibody (E:T 1:1, 4 h). Data were analyzed using two-tailed Mann-Whitney test (b-d). All bar graphs represent the mean ± s.d. Actual p values are indicated. Source data

Extended Data Fig. 10. Deconvolution of TCGA datasets.

Distribution of CD56^bright and CD56^dim NK cells in deconvoluted TCGA datasets. The boxplots indicates the median with the interquartile range (IQR), whiskers extend to the farthest point within 1.5 times the IQR from the box. For each plot and each subset n is the number of patients for each tumor type: SARC, n = 88; PAAD, n = 183; NSCLC, n = 600; BRAC, n = 1231; SKCM, n = 473; PRAD, n = 554; GBM, n = 175.