Totally, 1,816 samples were enrolled in this study, namely, 11 LUAD samples with scRNA-seq data, 500 LUAD samples from the TCGA, 1,007 LUAD samples from six independent GEO cohorts (https://www.ncbi.nlm.nih.gov/geo/), and 298 samples treated with immunotherapy from the IMvigor210 cohort. Single-cell RNA-sequencing data from 11 primary LUAD samples of GSE131907 were obtained from the GEO database, and were used to determine the NK cell marker genes of LUAD. The Cancer Genome Atlas (TCGA) bulk tumor transcriptomic data (FPKM normalized) and clinical information of 500 patients with LUAD were downloaded from the UCSC Xena (https://xenabrowser.net/) for identifying survival-related genes and constructing prognostic signatures. Six independent microarray datasets, namely, GSE30219 (n = 83), GSE3141 fimmu.2022.850745(n = 58), GSE50081 (n = 127), GSE26939 (n = 115), GSE72094 (n = 398), and GSE31210 (n = 226), were also obtained from the GEO database for external validation. In this study, the TCGA RNA-sequencing data of were converted into transcripts per kilobase million (TPM) values, which are more comparable between TCGA samples and microarrays (31). Transcriptomic and matched clinical data of patients who received anti-PD-L1 treatment from the IMvigor210 cohort were collected from http://research-pub.gene.com/IMvigor210CoreBiologies to explore the value of NKCMGS in speculating on the immunotherapy response (32). The study used publicly available datasets with preexisting ethics approval from original studies.

We conducted an analysis of scRNA-seq data by R packages, including “Seurat” and “SingleR” (33). To retain high-quality scRNA-seq data, three filtering measures were applied to the raw matrix for each cell: only genes that were expressed in at least 5 single cells were included, cells that expressed less than 100 genes were eliminated, and cells with more than 5% of mitochondrial genes were also removed. We first used the “Seurat” R package to normalize scRNA-seq data by the “NormalizeData” function, setting the normalization method as “LogNormalize.” Normalized scRNA-seq data were then transformed into a Seurat object, and the top 1,500 highly variable genes were identified using the “FindVariableFeatures” function. After that, we applied the “RunPCA” function of the “Seurat” R package to perform the principal component analysis (PCA) to reduce the dimension of the scRNA-seq data based on the top 1,500 genes. We used JackStraw analysis to identify significant PCs, and we selected the first 15 PCs for cell clustering analysis according to the proportion of variance explained. The “FindNeighbors” and “FindClusters” functions in the “Seurat” package were used for cell clustering analysis. The k-nearest neighbor graph was constructed based on Euclidean distance in PCA using the “FindNeighbors” function to determine the closest neighbors of each cell. Then, t-distributed stochastic neighbor embedding (t-SNE) was performed using the “RunTSNE” function. Cell clustering was demonstrated using t-SNE-1 and t-SNE-2. The “FindAllMarkers” function in the “Seurat” package was used to calculate the differentially expressed genes (DEGs) of each cluster using Wilcoxon–Mann–Whitney tests. To identify the marker genes for each cluster, the cutoff threshold values, adjusted p-value <0.01 and |log2 (fold change)| >1 were used. For cluster annotation, we performed a reference-based annotation using reference data from the Human Primary Cell Atlas (34).

A univariate Cox regression analysis was performed to evaluate the prognostic value of NK cell marker genes for OS in TCGA LUAD patients, and genes with p <0.01 were identified as prognostic genes. Next, to minimize overfitting, prognostic genes were assessed by least absolute shrinkage and selection operator (LASSO) Cox proportional hazards regression using the “glmnet” package. LASSO is a popular method for regression with high-dimensional predictors and is broadly applied to the Cox proportional hazard regression model for survival analysis (35). By using the function “cv.glmnet”, 10-fold cross-validation was conducted to select the best model. The tuning parameter λ was chosen by 1 − SE (standard error). We got a list of genes with non-zero beta coefficients. Finally, based on the genes generated by LASSO Cox regression analysis, we used a stepwise multivariate Cox regression analysis to identify the prognostic values of specific gene signatures. The risk model was constructed by a linear combination of the mRNA expression of the genes and the relevant risk coefficient. Based on the median cut-off value, the patients were classified into the low-risk or high-risk groups. To validate the prognostic power of the NKCMGS, the area under the curve (AUC) was calculated using the “survivalROC” package (36). The Kaplan–Meier method was employed for survival analysis, and the log-rank test was used to determine the statistical significance of the differences using the R package “survminer” (37). The predictive ability of the signature was validated using survival analysis and AUC in 6 independent GEO datasets.

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis by using the R package “clusterProfiler” (38). GO analysis was performed using the enrichGO function of the R package “clusterProfiler” and GO annotations were based on genome-wide annotation packages (org.Hs.eg.db) released by the Bioconductor project (39). KEGG analysis was performed using the enrichKEGG function of the R package “clusterProfiler” and “clusterProfiler” queries the latest online KEGG database through a web API to obtain the pathway data and perform functional analysis. A p-value of < 0.05 was considered significant enrichment.

The CIBERSORT algorithm, a useful method for obtaining infiltrating characteristics of 22 immune cell types with gene expression profiles (40), was applied to dissect the proportion of immune cell infiltration in high-risk and low-risk groups. A seven-metagene (HCK, IgG, Interferon, LCK, MHC-I, MHC-II, and STAT1) has been extensively used to assess the inflammatory activity in TME (41). Therefore, we conducted GSVA analysis to investigate the associations between the NKCMGS and metagenes of inflammatory activities by using the “GSVA” package (42). Heatmap plots were generated using the “ComplexHeatmap” R package from Bioconductor (43).

The ESTIMATE algorithm was employed to assess levels of stromal and immune cell infiltration using expression profiles by the “estimate” R package (44). Stromal score, immune score, ESTIMATE score, and tumor purity score were calculated using the RNA-sequencing data of the TCGA LUAD cohort and a Wilcoxon t-test was performed to compare these scores between different risk groups.

We first applied PD-L1 expression, tumor mutation burden (TMB), and TCR repertoire to predict the response to immune checkpoint blocking therapy. The PD-L1 mRNA expression of LUAD patients was collected from RNA-sequencing data of the TCGA LUAD cohort. Gene mutation data of LUAD patients were downloaded from the TCGA database and TMB was calculated using “maftools” package (45). TMB was determined as the number of somatic indels and base substitutions per million bases in the coding region of the genome detected, and was calculated as previously described (46). The richness and Shannon diversity indexes were used to characterize the diversity of the TCR repertoire. The richness measures the number of unique TCRs in the sample, while the Shannon diversity index reflects the relative abundance of the different TCRs. The richness values and Shannon diversity index valves of TCR in the TCGA LUAD patients were obtained from the Pan-Cancer Atlas study (47). Additionally, 298 urothelial carcinoma patients with both transcriptomic data and treatment response to immunotherapy from the IMvigor210 cohort were used for speculating the immunotherapy response of the signature.

Categorized variables between different risk groups were compared by the Wilcoxon t-test. Univariate and multivariate Cox regression analyses were used to investigate the prognostic value of the NKCMGS and different clinicopathological characteristics. P <0.05 was set as a significant threshold. Benjamini–Hochberg was implemented to adjust the P-value for multiple testing using the R function “p.adjust”. For data analysis and generation of figures, R software version 4.1.0 (http://www.R-project.org) was used.

Based on scRNA-seq data of GSE131907, we obtained gene expression profiles of 45,149 cells from 11 primary LUAD samples for further analysis (Figure 1A). We conducted PCA using the top 1,500 variable genes to reduce the dimensionality, and 17 cell clusters were then identified (Figure 1B). Subsequently, the cell identity of each cluster was annotated using a reference dataset from the Human Primary Cell Atlas, and cells in cluster 7 were defined as NK cells (Figure 1C). This cluster was also found to have distinct gene expression profiles, with 189 genes differentially expressed between the 17 clusters (Figure 1D), which were identified as LUAD-related NK cell marker genes (Supplementary Table 2). The functional enrichment, including GO and KEGG analysis, showed that the NK cell marker genes were mostly related to immune features, such as positive regulation of leukocyte activation, MHC protein complex, antigen binding, and hematopoietic cell lineage (Supplementary Figure 1).

To construct a prognostic signature based on the 189 NK cell marker genes, we first used the TCGA LUAD cohort as the training set to perform a univariate Cox regression analysis, and 25 NK cell marker genes were significantly related to OS (Supplementary Table 3). Next, LASSO Cox regression analysis with one standard error (SE) and 10-fold cross-validation was conducted on the 25 NK cell marker genes, and 16 genes were screened out (HPGDS, CTSG, SLC18A2, GCSAML, ADRB2, ACTG1, ACOT7, CLIC1, SELENOK, PEBP1, BEX4, BIRC3, DDIT4, TRBC1, ACAP1, and S100A10) for further analysis (Supplementary Figures 2A, B). Finally, we used stepwise multivariate Cox regression analysis to optimize the prognostic signature to only include the 7 most predictive genes: Risk score = (−0.614 × GCSAML expression) + (1.893 × ACTG1 expression) + (1.022 × ACOT7 expression) + (−1.715 × SELENOK expression) + (−1.840 × PEBP1 expression) + (1.077 × BIRC3 expression) + (−1.180 × ACAP1 expression) (Supplementary Figure 2C). The relative expression of the 7 marker genes in various clusters is shown in Supplementary Figure 3, which indicates the specificity of the expression of the signature genes. The median risk score was 0.956 by ranking the risk score from high to low, which was used to classify patients into low-risk (n = 250) and high-risk (n = 250) groups. Figure 2A exhibited the distribution of risk scores and survival status, which suggested more deaths in the high-risk group. Figure 2B shows the expression details of the 7 NK cell marker genes.

Kaplan–Meier analysis demonstrated that patients with high-risk scores had significantly inferior OS than patients with low-risk scores (Figure 2C). To assess the predictive accuracy of the risk model, time-dependent area under the ROC curves for OS was calculated, and the 1-, 3-, and 5-year AUC values were 0.710, 0.725, and 0.730, respectively (Figure 2D). The performance of the model was then evaluated using a ten-fold cross-validation procedure, and the 1-, 3-, and 5-year mean AUC values were 0.669, 0.674, and 0.652, respectively.

The predictive value of NKCMG was first assessed in TCGA LUAD patients with different genders, ages, smoking histories, and tumor stages. The results revealed that high-risk score was significantly correlated with an inferior prognosis in the male (P = 1.5 × 10⁻⁴, Supplementary Figure 4A), female (P = 3.563 × 10⁻⁵, Supplementary Figure 4B), young (P = 7.8 × 10⁻⁴, Supplementary Figure 4C) or old (P = 1.306 × 10⁻⁶, Supplementary Figure 4D), non-smoker (P = 0.016, Supplementary Figure 4E), smoker (P = 4.465 ×10⁻⁶, Supplementary Figure 4F), early stage (P =1.709 × 10⁻⁵, Supplementary Figure 4G) or advanced stage (P = 0.0094, Supplementary Figure 4H) LUAD patients. Next, we further evaluated the predictive performance of NKCMGS in the TCGA LUAD patients stratified by different molecular characteristics, including EGFR, KRAS, and TP53 mutations. Similarly, we observed that the NKCMGS showed robust predictive power in the EGFR wild-type (WT) (P = 9.923 × 10⁻⁷, Supplementary Figure 5A), EGFR mutation (MUT) (P = 0.032, Supplementary Figure 5B), KRAS WT (P = 1.23 × 10⁻⁶, Supplementary Figure 5C), KRAS MUT (P = 0.012, Supplementary Figure 5D), TP53 WT (P = 6.9 × 10⁻⁴, Supplementary Figure 5E) or TP53 MUT (P = 2.377 × 10⁻⁵, Supplementary Figure 5F) subgroup.

To validate the robustness of the NKCMGS, we included 6 independent GEO cohorts in this study, and the clinical features of these 6 GEO cohorts are shown in Table 1. We used the same formula to calculate the risk score of each patient in 6 GEO cohorts. Patients were sorted into the high-risk and low-risk groups in each cohort by the median risk score. Kaplan–Meier analysis demonstrated that the high-risk group had inferior prognosis than the low-risk group in all 6 GEO cohorts, namely, GSE30219 (Figure 3A, HR: 3.557, 95% CI: 1.856–6.818, P = 1.32 × 10⁻⁴), GSE3141 (Figure 3B, HR: 3.064, 95% CI: 1.457–6.443, P = 0.002), GSE50081 (Figure 3C, HR: 1.932, 95% CI: 1.099–3.394, P = 0.02), GSE26939 (Figure 3D, HR: 2.312, 95% CI: 1.390–3.846, P = 9.3 × 10⁻⁴), GSE72094 (Figure 3E, HR: 2.038, 95% CI: 1.387–2.995, P = 2.1 × 10⁻⁴), and GSE31210 (Figure 3F, HR: 1.555, 95% CI: 0.789–3.063, P = 0.2). The ROC curves of the risk score in the 6 validation cohorts also showed good performance (Supplementary Figure 6). Additionally, a prognostic meta-analysis was performed by R package “meta” using the random-effects model to evaluate the integrated predictive value of NKCMGS in these 6 cohorts (48). The results of the meta-analysis indicated that NKCMGS was a significant prognostic indicator for patients with LUAD (HR: 2.227, 95% CI: 1.782–2.784, P = 1.96 × 10⁻¹²) (Figure 3G).

To further investigate whether the risk score can independently affect the prognosis of LUAD, we conducted univariate and multivariate Cox regression analysis using clinical features, molecular factors, and the risk score in the TCGA LUAD patients. As expected, the results of multivariate Cox regression analysis proved that the risk score was an independent prognostic factor (HR: 1.889, 95% CI: 1.373–2.599, P = 9.37 × 10⁻⁵) (Table 2). Meanwhile, we performed a LASSO Cox regression analysis with the risk score and all these clinical features to select the most predictive variables. The results demonstrated that the risk score and tumor stage were the best predictive factors for the prognosis.

To elucidate the potential mechanism of the excellent predictive capability of NKCMGS, we further investigated biological pathways related to NKCMGS. Firstly, the correlation analysis was performed using the TCGA LUAD dataset to identify the genes that were closely correlated with the risk score (Pearson |R| >0.4, P <0.05). As shown in Supplementary Figure 7A, 100 positively correlated genes and 24 negatively correlated genes were filtered out (Supplementary Table 4). Subsequently, we performed GO and KEGG enrichment analyses using the “ClusterProfiler” package for these selected genes. GO analysis revealed that these genes were mainly implicated in the biological processes of mitotic division, namely, chromosome segregation, mitotic nuclear division, and the G2/M transition of the mitotic cell cycle (Supplementary Figure 7B). KEGG analysis also verified that these genes were closely involved in the cell cycle pathway (Supplementary Figure 7C).

As NK cells play a vital role in anti-tumor immunity, we explored the relationship of the NKCMGS with immune cell infiltration in LUAD patients. By using the ESTIMATE algorithm, we found that high-risk patients had lower immune score, stromal score, ESTIMATE score, and higher tumor purity than low-risk patients (Figures 4A–D), which suggested that the risk score was negatively correlated with the level of immune cell infiltration. Subsequently, we applied the CIBERSORT algorithm to estimate the infiltration level of different types of immune cells in the TME. The CIBERSORT analysis revealed that high-risk patients had a higher fraction of resting NK cells, M0 macrophages, M2 macrophages, activated dendritic cells, and activated mast cells, but had a lower fraction of plasma cells, CD8⁺ T cells, follicular helper T cells, regulatory T cells, resting dendritic cells, and resting mast cells (Figure 4E). In Figure 4F, the fractions of different immune cells between high- and low-risk groups are shown. We further conducted a correlation analysis between the risk score and immune cell infiltration, which showed that the risk score was positively related to macrophages and neutrophils but was negatively related to T, B, and mast cells (Supplementary Figure 8).

To figure out the relationship between NKCMGS and inflammatory activities, we explored the associations between the NKCMGS and 7 clusters of metagenes (HCK, LCK, IgG, Interferon, MHC-I, MHC-II, and STAT1), representing various inflammatory and immune activities as previously reported (41). Supplementary Figure 9A shows the expression details of these metagenes in the TCGA LUAD dataset. Next, we used Gene Sets Variation Analysis (GSVA) to calculate the expression of 7 gene clusters and the correlation between the NKCMGS and each cluster of metagenes is shown in Supplementary Figure 9B. The results showed that the risk score was negatively correlated with HCK, IgG, LCK, MHC-I, and MHC-II.

Given the important roles in anti-tumor immunity of NK cells and the promising treatment efficacy of NK cell-based immunotherapy, we explored whether the NKCMGS could predict responses of LUAD patients to immune checkpoint inhibitors. First, we analyzed the relationship between the NKCMGS and widely recognized immunotherapy biomarkers (PD-L1 expression and TMB) in the TCGA LUAD cohort. The results indicated that PD-L1 expression exhibited no significant difference between the low-risk and high-risk patients, but low-risk patients harbored a significantly lower TMB than high-risk patients (Figures 5A, B). In previous studies, the T-cell receptor (TCR) is in charge of the recognition of antigens presented by the MHC, and the repertoire analysis of TCR has been demonstrated as a useful biomarker for stratification and monitoring of patients on immunotherapy (49–51). Subsequently, we analyzed the TCR repertoire and found that the TCR richness and diversity of low-risk patients was significantly higher than that of high-risk LUAD patients (Figures 5C, D). Finally, to further explore the value of NKCMGS in predicting the immunotherapy response, 298 patients from the IMvigor210 cohort who received anti-PD-L1 treatment were enrolled in this study for analysis. Kaplan–Meier analysis showed an inferior survival rate for high-risk patients after immunotherapy (Figure 5E). Lower risk scores were associated with an objective response to anti-PD-L1 treatment (Wilcoxon test, P = 0.01; Figure 5F). The objective response rate of anti-PD-L1 treatment was significantly elevated in the low-risk group (two-sided chi-square test, P = 0.005; Figure 5G). The ROC curves showed that the combination of TMB, PD-L1, and risk score models could predict anti-PD-L1 response with 76.1% accuracy, which was superior to that of TMB (AUC = 0.728), PD-L1 (AUC = 0.569), or risk score (AUC = 0.603) alone (Figure 5H). Collectively, these findings indicate that patients with a low-risk score are more likely to benefit from immunotherapy, and the NKCMGS may be a useful biomarker to identify LUAD patients who may benefit from immunotherapy.