Characterization of natural killer (NK) cells in lung adenocarcinoma and construction of an NK risk signature based on single-cell and macromolecular RNA-seg data

doi:10.21203/rs.3.rs-4840386/v1

Download PDF

Article

Characterization of natural killer (NK) cells in lung adenocarcinoma and construction of an NK risk signature based on single-cell and macromolecular RNA-seg data

https://doi.org/10.21203/rs.3.rs-4840386/v1

This work is licensed under a CC BY 4.0 License

You are reading this latest preprint version

Background/Aims: Natural killer (NK) cells play a crucial role in tumor cell apoptosis, immune milieu regulation, and angiogenesis inhibition. This study aims to analyze the NK signature in lung adenocarcinoma (LUAD) and establish an NK cell-based risk signature for predicting the prognosis of LUAD patients.

Methods: Single-cell RNA sequencing (scRNA-seq) data were obtained from the GEO database, while RNA-seq and microarray data from LUAD were simultaneously obtained from the TCGA and GEO databases. The scRNA-seq data were processed using the Seurat R package to identify NK clusters based on NK markers. Differentially expressed genes (DEGs) between normal and tumor samples were identified through differential expression analysis of LUAD-related data. Pearson correlation analysis was used to identify DEGs associated with NK clusters, followed by one-way Cox regression analysis to identify NK cell-related prognostic genes. Subsequently, Lasso regression analysis was employed to construct a risk signature based on NK cell-related prognostic genes. Finally, a column-line diagram model was constructed based on the risk signature and clinicopathological features.

Results: Based on the scRNA-seq data, we identified five Natural killer (NK)cells clusters in lung adenocarcinoma (LUAD), with four of them showing associations with prognosis in LUAD. Out of 19,495 differentially expressed genes (DEGs), a total of 725 genes significantly associated with NK clusters were pinpointed and further narrowed down to form a risk profile comprising 13 genes. These 13 genes were primarily linked to 21 signaling pathways, including vascular smooth muscle contraction, RNA polymerase, and pyrimidine metabolism. Additionally, the risk profile exhibited significant associations with stromal and immune scores, as well as various immune cells. Multifactorial analysis indicated that the risk profile served as an independent prognostic factor for LUAD, and its efficacy in predicting the outcome of immunotherapy was validated. Furthermore, a novel column-line diagram integrating staging and NK-based risk profiles was developed, demonstrating strong predictability and reliability in prognostic forecasting for LUAD.

Conclusion: The NK cell-based risk signature proves to be a valuable tool for predicting the prognosis of patients with lung adenocarcinoma (LUAD). Furthermore, a comprehensive understanding of NK cell characterization in LUAD could potentially unveil insights into the response of LUAD to immunotherapies and offer novel strategies for cancer treatment.

Biological sciences/Immunology

Health sciences/Biomarkers

natural killer cells

lung adenocarcinoma

differentially expressed genes

risk signature

columnar plots

Lung cancer is a lethal disease that causes the most frequent cancer-related deaths worldwide^[1], with adenocarcinoma accounting for more than half of all cases^[2]. Although breakthroughs in chemotherapy, radiotherapy, surgery and other novel therapeutics, the incidence and mortality of lung adenocarcinoma (LUAD) continues to increase^[3]. Therefore, it is essential to identify the relevant pathological parameters to provide foundation for clinical prognosis analysis. Omics technology is regarded as a effective tool for understanding the molecular pathogenesis of LUAD^[4–6]. Whole-genome sequencing and transcriptome sequencing has provided valuable insights into the gene amplification of RNA methyltransferase in LUAD, suggesting the activation mechanism of tumorigenesis^[4]. Single-cell RNA sequencing techniques were applied to investigate the differences in immune response between the lung squamous cell carcinoma (LUSC) and LUAD^[5].

One subset of tumor immune cells called natural killer (NK) cells is well-known for its capacity to identify and then eliminate cancer cells^[7]. Recent studies demonstrated that NK cells could identify and remove LUAD cells by engaging activating receptors with stress-induced ligands that are overexpressed on cancer cells' surfaces^[8–9]. NK cells display an potent antitumor immune cytotoxicity by regulating the tumor microenvironment via MEK/ERK and PI3K/Akt/mTOR pathways^[10–11]. NK cells have the capacity to influence the adaptive immune response to LUAD by interacting with other immune cells^[12]. Growing evidence have suggested that the status of NK cells within LUAD are correlated with patient prognosis^[13–15]. Increased NK cell infiltration and activity are frequently associated with improved survival and treatment outcomes^[14]. Therefore, NK cells interact with lung cancer in a complex and dynamic manner, and understanding these interactions is critical for predicting patients’ prognosis.

Despite some studies on NK cells in LUAD have been conducted^[16–17], the systematic characteristics of NK cells, as well as their association with LUAD prognosis and immunotherapy response, remain limited known. In this study, LUAD scRNA-seq and transcriptome data were retrieved from publicly available databases to differentiate LUAD subsets and uncover NK cell-related risk signature. The immunological environment and immunotherapy response that support the NK cell-based signature were studied, and the clinical significance of the signature was also identified. In addition, a unique nomogram that integrated the NK cell-based risk signature with clinicopathological parameters was established to make it easier to use NK cell-based factors to predict LUAD prognosis. Our study might offer researchers novel insights into the biology of LUAD, resulting in more targeted therapy and a better prognosis for those with LUAD.

2.1 Data Acquisition and Organization

ScRNA-seq data from GSE131907 were acquired from the Gene Expression Omnibus (GEO) database, encompassing two lung adenocarcinoma (LUAD) and two normal control samples. Selection criteria included cells expressing a given gene in a minimum of three cells and cells exhibiting at least 250 genes expressed. The mitochondria to rRNA ratio was assessed utilizing the PercentageFeatureSet function within the Seurat R package. Only cells expressing over 6000 genes and with a Unique Molecular Identifier (UMI) count exceeding 100 were retained, resulting in a total of 23,673 cells.

Single-cell transcriptomic and copy number variation (CNV) analyses were conducted using data from The Cancer Genome Atlas (TCGA) database, incorporating pertinent LUAD clinical data. Samples devoid of survival and outcome status were omitted, yielding 500 tumor and 59 adjacent non-tumor samples. The GSE31210 cohort, comprising 226 LUAD samples, served as a validation set; samples lacking follow-up and outcome data were excluded. Following this, ten cancer-associated pathways were identified through literature review, including WNT, PI3K, NOTCH, RAS, cell cycle, MYC, TGF-Beta, HIPPO, and NRF1 pathways^[18].

2.2 Definition of NK

In this study, we conducted a re-analysis of single-cell RNA sequencing (scRNA-seq) data from LUAD using the Seurat package^[19]. Our objective was to comprehensively characterize NK cells. We applied several data preprocessing steps to ensure data quality. First, we filtered out cells with over 6000 or under 250 expressed genes. Next, we performed log-normalization and removed batch effects from the four samples.

To reduce the dimensionality of the data, we employed principal component analysis (PCA) with 20 principal components and a resolution of 0.25 for non-linear dimensionality reduction. Subsequently, we used the FindNeighbors and FindClusters functions (dim = 40, resolution = 0.2) to cluster cells into subgroups. To visualize the clustering results, we utilized t-distributed stochastic neighbor embedding (t-SNE) with the RunTSNE function.

We then focused specifically on NK cells by identifying marker genes such as NKG7, KLRD1, KLRB1, GNLY, and GZMB. Clustering and t-SNE dimensionality reduction were performed exclusively on the NK cells. To identify marker genes specific to each NK cell cluster, we employed the FindAllMarkers function with criteria including a log fold change (logFC) threshold of 0.5, a minimum percentage threshold (minpct) of 0.35, and a corrected p-value cutoff of < 0.05.

To gain further insights into the functional characteristics of the NK cell clusters, we performed gene ontology enrichment analysis using the KEGG enrichment analysis method^[20]. This analysis was carried out using the clusterProfiler software package.

Additionally, we analyzed copy number variation (CNV) features in the NK cell clusters to distinguish between tumor and normal cells in each sample. For this analysis, we utilized the CopyKAT R software package^[21].

2.3 Identification of NK hub genes

During our analysis, we initially identified differentially expressed genes (DEGs) between tumor and normal tissues using the limma software package. We set a screening threshold of FDR < 0.05 and |log2(Fold Change)| > 1.5^[22] to determine significant DEGs.

Next, we assessed the correlation between these DEGs and the clustering of NK cells. We specifically focused on identifying key NK cell-related genes with p < 0.001 and a correlation coefficient (cor) > 0.4.

For prognostic gene screening, we conducted univariate Cox regression analysis using the survival package. Genes with p < 0.05 were selected as potential prognostic markers. To further refine the gene selection, we applied the least absolute shrinkage and selection operator (lasso) Cox regression analyses. Subsequently, multivariate Cox regression analyses using stepwise regression methods were performed.

Based on the results of the multivariate Cox model, we constructed a risk signature using the formula: risk score = Σβi * Expi. Here, i represents the genes included in the risk signature, expi denotes the gene expression, and βi signifies the coefficient obtained from the multivariate Cox model. We categorized patients into high-risk and low-risk groups after zero-mean normalization.

To evaluate the predictive performance of the risk signatures, we conducted receiver operating characteristic curve (ROC) analyses using the timeROC software package. These analyses provided an assessment of the risk signature's ability to predict patient outcomes. Similar analyses were performed in the validation cohort to validate the findings.

2.4 Immune landscape analysis

In order to delve deeper into the tumor microenvironment (TME), we utilized the CIBERSORT algorithm to evaluate the proportions of 22 distinct immune cell subtypes within the TCGA cohort. This analysis allowed us to assess the relative abundance and composition of immune cells present in the TME.

To gain additional insights into the TME, we also calculated an immunity score and stroma score using the ESTIMATE algorithm^[23]. These scores provide valuable information about the immune and stromal components within the TME. The immunity score reflects the level of immune cell infiltration, while the stroma score indicates the extent of stromal cell involvement. By utilizing these scores, we can gain a better understanding of the complex interplay between the tumor and its microenvironment.

2.5 Constructing risk characteristics and column line plots

To assess the prognostic significance of clinicopathological and risk characteristics, we conducted univariate and multivariate Cox regression analyses. In the multivariate Cox model, variables with a significance level of P < 0.05 were selected. Subsequently, we utilized the rms package^[24] to construct column-line plots, which enabled us to predict the prognosis of LUAD patients based on these selected variables.

To evaluate the predictive accuracy of the models, calibration curves were generated. These curves provide a visual representation of the agreement between the observed and predicted outcomes. Additionally, we performed decision curve analysis (DCA) to assess the reliability of the models. DCA is a useful tool for evaluating the clinical utility and net benefit of predictive models.

2.6 Responsiveness to immune checkpoint blocks

We obtained transcriptomic and clinical data from the IMvigor210 cohort of LUAD patients who were treated with anti-PD-L1 drugs (atezolizumab)^[25]. This cohort, which can be accessed at http://research-pub.gene.com/IMvigor210CoreBiologies, provided valuable insights into the response to anti-PD-L1 therapy in LUAD patients.

Additionally, we acquired transcriptomic data from the GSE78220 cohort, which consisted of pre-existing melanomas treated with anti-PD-1 checkpoint inhibition. We utilized this data to investigate the potential value of risk marker scores in predicting the effectiveness of immune checkpoint blockade therapy (ICB)^[26].

2.7 Statistical analyses

All statistical analyses were conducted using R software (v3.6.3). To assess the relationships between variables, correlation matrices were generated using either Pearson or Spearman correlation coefficients, depending on the nature of the data. The Wilcoxon test was employed to compare the two groups, providing insights into potential differences between them.

Furthermore, survival differences were evaluated using Kaplan-Meier curves and the log rank test. This allowed for the examination of potential disparities in survival outcomes between different groups or conditions. A p-value of less than 0.05 was considered statistically significant, indicating a notable result.

3.1 Screening and Identification of NK Clusters in scRNA-seq Data

The flowchart depicting the design of this study is presented in Fig. 1. After an initial screening of the scRNA-seq data, 23,673 cells were obtained. Following log-normalization and dimensionality reduction, we identified 20 subclusters. Subsequently, utilizing five marker genes (NKG7, KLRD1, KLRB1GNLY, and GZMB), we successfully identified 680 NK populations (Fig. S1A, B).

To refine our analysis, we further clustered and downscaled the cells within the 680 NK populations, ultimately identifying five distinct NK clusters (Figure S1C, D). By observing the TSNE plots of the four sample distributions (Fig. 2A), we were able to generate and focus on these five NK clusters for subsequent analysis (Fig. 2B).

Examining the expression profiles of the 737 differentially expressed genes (DEGs) within the five NK clusters, we identified the top five DEGs, which were considered as marker genes for the NKF clusters (Fig. 2C). The distribution ratio of these five clusters within each cohort is illustrated in Fig. 2D.

Moreover, through KEGG analysis (Fig. 2E), we discovered that these DEGs were significantly enriched in various pathways, including the TNF signaling pathway, T-cell receptor signaling pathway, and NK cell-mediated cytotoxicity. Additionally, based on the CNV characteristics, we observed that the five NK clusters comprised a total of 680 tumor cells and normal cells (Fig. 2F).

3.2 Exploration of Cancer-Related Pathways in NK Clusters

To gain further insights into the relationship between NK clusters and tumor progression, we conducted GSVA scoring of 10 pathways associated with tumorigenesis within the 5 NK clusters (Fig. 3A). Notably, the proportion of malignant cells in NK_0, NK_1, and NK_3 clusters was significantly higher compared to the other two clusters (Fig. 3B), and there were significant differences between each cluster.

Furthermore, we analyzed the GSVA scores of the 10 tumor-associated pathways in both malignant and non-malignant cells within each NK cluster. Interestingly, we observed significant differences in the RAS pathway within the NK_0 cluster and the NOTCH pathway within the NK_1 cluster (Fig. 3C-G).

To investigate the potential correlation between NK clusters and prognosis, we initially computed the ssGSEA scores of marker genes for each NK cluster utilizing the TCGA cohort. Strikingly, the results revealed that all five NK clusters exhibited higher scores in normal samples compared to tumor samples (Fig. 4A).

Next, we categorized the LUAD samples from the TCGA dataset into high and low NK score groups based on optimal cutoff values determined by the survminer R package. Intriguingly, within the high-NK score group, the samples belonging to NK_1, NK_2, NK_3, and NK_4 groups displayed a more favorable prognosis compared to the low-NK score group. However, it is noteworthy that the prognosis of LUAD was not associated with the NK_0 cluster (Fig. 4B-F). Although there were differences in NK_0 enrichment between LUAD and normal samples, our findings suggest that the contribution of NK_0 clusters to LUAD progression is minimal.

3.3 Identification of NK cell-related core genes

To identify core genes related to NK cells, we conducted a comprehensive analysis by screening differentially expressed genes (DEGs) between tumor tissues and normal tissues. Out of the 19,495 DEGs identified, 12,525 were up-regulated and 6,970 were down-regulated. Remarkably, 725 genes exhibited significant correlation with prognostically relevant NK clusters. Subsequently, we performed one-way Cox regression analysis to assess the prognostic value of each gene, resulting in the identification of 133 genes with prognostic significance (Fig. 5A, B).

Enrichment analysis of these genes revealed their involvement in various biological processes (BP), cellular components (CC), molecular functions (MF), and KEGG pathways (Fig. 5C). However, it is important to note that the discussion section should be completed first to provide context and specific results before filling in the details of the enrichment analysis.

To further refine the gene selection, we employed Lasso Cox regression analysis, which ultimately filtered down to 13 genes with a lambda value of 0.0311 (Fig. 5D, E). These 13 genes, namely ANOS1, ADGRE3, ADAMTS8, CCT3, CD79A, CX3CR1, DPEP2, FMO3, FBP1, OASL, RRAS, SELENOP, and SFTPC, were included in the risk signature after multifactorial Cox regression analysis using stepwise regression (Fig. 5F). The final risk signature formula consists of the following 13 genes with corresponding coefficients: RiskScore = "-0.121 * CD79A − 0.056 * DPEP2 + 0.094 * OASL + 0.363 * RRAS + 0.052 * ADAMTS8–0.282 * ADGRE3–0.039 * ANOS1 + 0.155 * CCT3 + 0.075 * CX3CR1–0.145 * FBP1–0.115 * FMO3–0.087 * SELENOP − 0.021 * SFTPC".

After Z-mean normalization, we calculated the risk scores for each sample and subsequently classified them into high and low-risk groups.

In the TCGA cohort, the area under the curve (AUC) values for 1- to 5-year survival ranged from 0.66 to 0.71. Similarly, in the GSE3141 cohort, the AUC values ranged from 0.64 to 0.72, while in the GSE31210 cohort, they ranged from 0.65 to 0.69. In the GSE37745 cohort, the AUC values ranged from 0.57 to 0.71, and in the GSE50081 cohort, they ranged from 0.70 to 0.77. Lastly, in the GSE68465 cohort, the AUC values ranged from 0.60 to 0.68 (Fig. 6A-F).

Furthermore, the Kaplan-Meier survival analyses demonstrated a similar pattern of survival outcomes in both the TCGA and GEO cohorts. High-risk patients exhibited significantly worse survival outcomes compared to low-risk patients (Fig. 6A-F).

3.4 Mutation and pathway analysis of key genes

In our investigation of the risk signature genes, we examined SNV mutations in these 13 genes. Interestingly, all genes except RRAS exhibited SNV mutations in a greater number of samples (Fig. 7A). We further analyzed the probability of co-occurrence of mutations between these key genes and the 10 most frequently mutated genes. Figure 7B demonstrates that while no significant mutation co-occurrence was observed among these 12 genes, significant co-occurrence was observed between ADGRE3 and TTN, CX3CR1 and MUC16, RYR2 and ZFHX4, ANOS1 and LRP1B, DPEP2 and CSMD3, as well as RYR2 and ZFHX4. We also used OncogenicPathways to count known oncogenic signalling pathways in the TCGA cohort (Fig. 7C).

To gain further insights into the association between these risk genes and LUAD, we investigated their correlations with various molecular features of LUAD. Remarkably, most of the genes exhibited significant negative associations with aneuploidy scores, homologous recombination defects, score alterations, fragment counts, and non-silencing mutation rates. Notably, CCTS was the only gene significantly positively associated with these molecular features (Fig. 7D).

Additionally, we explored the potential pathways associated with each risk gene. The analysis revealed significant associations between these 13 genes and a total of 21 pathways. Noteworthy pathways included the B cell receptor signaling pathway, the chemokine signaling pathway, and leukocyte migration (Fig. 7E, 7F).

3.5 Analysis of key genes in relation to immunity

Our data analysis revealed notable findings regarding the relationship between key genes and immunity. Specifically, we observed that CCT3 exhibited a significant negative correlation with matrix score, immune score, and estimated score. Conversely, the remaining 12 genes displayed significant positive correlations with these scores (Fig. 8A).

To further explore these associations, we categorized each gene based on its median expression and compared the three scores among different expression groups. Notably, for the ANOS1 gene, the high-expression group exhibited significantly lower scores compared to the low-expression group, while the high-expression groups of the remaining genes displayed significantly higher scores compared to the low-expression group (Fig. 8B).

Furthermore, our correlation analysis demonstrated significant associations between the 13 key genes and various immune cells, including B cells, T cells, and NK cells (Fig. 8C). Moreover, we observed significant differences in several immune cells between the high and low expression groups of the risk genes (Fig. 8D).

3.6 Sensitivity analysis of risk signature for PD-L1 blockade immunotherapy

The emergence of T-cell immunotherapy as a promising approach for cancer treatment has highlighted the importance of evaluating the prognostic value of risk characteristics in the context of immune checkpoint therapy. To assess this, we examined the IMvigor210 and GSE78220 cohorts.

Within the IMvigor210 cohort, consisting of 348 patients with varying responses to anti-PD-L1 receptor blockers, including complete remission (CR), partial remission (PR), stable disease (SD), and progressive disease (PD), our analysis revealed compelling results. Patients in the low-risk group experienced a significant clinical benefit and demonstrated significantly longer overall survival compared to those in the high-risk group (Fig. 9A, p = 0.0017). Furthermore, the risk scores were higher in patients with SD/PD compared to those with CR/PR (Fig. 9B). Additionally, the high-risk group exhibited a higher percentage of SD/PD compared to the low-risk group (Fig. 9C).

Intriguingly, the survival analysis demonstrated a significant difference between stage I + II patients (Fig. 9D, p = 0.041) and stage III + IV patients (Fig. 9E, p = 0.018) within different risk groups. This suggests that risk scores were particularly sensitive in distinguishing patients with different stages of LUAD.

Similarly, in the GSE78220 cohort, we observed that patients in the low-risk group had significantly longer overall survival compared to those in the high-risk group (Fig. 9F, p = 0.0067). Notably, patients with PD exhibited higher risk scores compared to those with CR/PR (Fig. 9G), and the high-risk group displayed a higher percentage of PD cases than the low-risk group.

3.7 Independent risk factor identification and histogram construction

In order to enhance the predictive performance of risk features, we conducted univariate and multivariate COX regression analyses to integrate clinicopathological factors and risk scores. The results from the multivariate COX regression analysis revealed that the risk signature (HR = 2.69, 95% CI: 1.93–3.74, P < 0.001) and metastatic status (HR = 2.06, 95% CI: 1.39–3.05, P = 0.002) were the most crucial independent prognostic factors for LUAD.

Based on these findings, we developed a histogram that incorporated both staging and risk scores (Fig. 10A-C). This histogram demonstrated effectiveness in accurately predicting actual survival outcomes, as confirmed by the calibration curve (Fig. 10D). Furthermore, the decision curve analysis (DCA) indicated that the discriminative power of the histogram in identifying high-risk patients surpassed that of the risk score and staging alone (Fig. 10E). Additionally, the time-dependent receiver operating characteristic (TimeROC) analysis demonstrated that both the risk score and the histogram exhibited higher area under the curve (AUC) values compared to other metrics in the TCGA cohort (Fig. 10F).

NK cells regulate tumor proliferation, angiogenesis and the immune system of LUAD^[27]. Studies on the pathophysiology of NK cells in LUAD have progressed to the genetic and molecular levels^[28]. Nonetheless, the characterisation of NK cells in LUAD remains unknown, and research on LUAD prognosis with respect to NK cells is comparatively limited. Our study focused on the characterization of NK cells. ScRNA-seq data were used to classify NK cells in LUAD. Eventually, we uncovered five NK cell subsets with diverse features that might regulate biology of the LUAD. Recent studies have demonstrated the prognostic significance of the NK cell-related gene signature in the assessment of LUAD outcomes^[29–30]. Consistently, our findings revealed four distinct clusters exhibiting strong associations with LUAD prognosis, as quantified by a composite score derived from DEGs encompassing all five clusters. Moreover, our results suggested that differences in the Notch signaling pathway among the identified clusters might underlie the prognostic significance of NK cells in LUAD. Notch signaling pathway plays a vital role in regulating the progression of LUAD^[31]. Suppressing the Notch signaling pathway has demonstrated efficacy in attenuating collagen synthesis and curtailing cellular invasiveness^[32]. Notably, high levels of Notch-related genes have been proven to be correlated with poor survival outcomes in LUAD patients^[33].

Considering the prognostic relevance of the four NK cell clusters, we constructed a 13-gene risk signature. This signature comprised nine protective genes(ANOS1, ADGRE3, CD79A, CX3CR1, DPEP2, FMO3, FBP1, SELENOP and SFTPC) and four risk genes(ADAMTS8, CCT3, RRAS and OAST). In our study, SNV mutations were detected in ANOS1, ADGRE3, ADAMTS8, CCT3, CD79A, CX3CR1, DPEP2, FMO3, FBP1, OASL, SELENOP, and SFTPC, while no statistically significant co-occurrence patterns were observed among these mutations. SNV mutations hold a significant position in the pathogenesis of LUAD, as they can manifest as "driver mutations" that bestow tumor cells with growth advantages, thereby fueling the onset and progression of the disease^[34–35]. While no prior studies have revealed a connection between SNV mutations in these risk genes and LUAD, our findings indicated such mutations might contribute to the progression of LUAD. We also demonstrated that the 13 genes were closely associated with 21 signaling pathways. It was proven that part of these pathways were related to onset and prognosis of LUAD. Vascular smooth muscle contraction, RNA polymerase, pyrimidine metabolism and B cell receptor signaling pathway act as effective predictors for the prognosis and therapeutic responses of LUAD patients^[36–38]. Spliceosome-related and lysosome-related genes might increase LUAD incidence rate^[39–40]. Homologous recombination pathway is up-regulated in LUAD, and increased expression of relevant genes is associated with poor overall survival^[41].

Previous research demonstrated that tumor growth could attribute to the interplay between NK cells and the tumor immune microenvironment (TIME)^[42–43]. Our findings revealed that 9 protective genes and three risk genes were positively and significantly associated with the immunological score, whereas one risk gene exhibited a negative correlation. These findings suggested potential interactions between these genes and the TIME in LUAD, hinting at their potential as therapeutic targets. The TIME is constituted by a diverse array of immune cells within the tumor islets, collectively dictating the antitumor immunological milieu of TIME^[44]. NK cells can interact in various ways with these immune cells, fostering an immunosuppressive TIME that allows tumor cells to evade immune surveillance^[45]. In the risk signature, most protective genes were negatively correlated with NK cell activation and T cell activation. Meanwhile, regulatory T cells were negatively associated with protective genes and positively associated with risk genes. It was reported that NK cells could modulate T cell apoptosis in tumors^[46]. T cells also could interact with NK cells and enhance their viability^[47]. These suggested that NK related-genes acted not only on NK cells but also on T cells, therefore regulating the whole TIME.

Despite the fact that many patients have resistance to immunotherapeutic therapy^[48], we found that relevant risk signature efficiently stratified patients based on their likelihood of benefiting from immunotherapies. Additionally, previous studies have reveal that NK cells have the potential to treat solid tumors^[49–50]. However, the therapy was less effective for certain patients, which might be a consequence to the extremely complex immunosuppressive milieu^[51]. Fortunately, our findings demonstrated that the NK cell-based signature had the capacity to predict the response to anti-PD-L1 immunotherapy. Further experimental research is warranted to elucidate the mechanisms underlying NK cell-TIME interactions in LUAD and their potential therapeutic implications in LUAD immunotherapy.

There are some limitations in our study. Firstly, the retrospective use of public data to establish NK cell clusters and the NK cell-based risk signature highlights the need for validation in prospective, multi-center LUAD cohorts in future investigations. Secondly, while a predictive nomogram was developed, additional studies are required to better investigate the specific biological correlation between the included NK cell-related genes and clinicopathological characteristics.

To sum up, we comprehensively categorized NK cell in LUAD, identifying five distinct clusters with different characteristics. The DEGs in the five clusters were abundant in vascular smooth muscle contraction, RNA polymerase, pyrimidine metabolism, the B cell receptor signaling pathway, etc. Four of the clusters showed significant correlations with LUAD prognosis and. An NK cell-based prognostic risk signature was created by 13 key genes from the four clusters. This signature appeared to be associated with the immune milieu and might be used to predict responsiveness to anti-PD-L1 therapy. Finally, we constructed a reliable nomogram incorporating the clinicopathological parameters and NK cell-related risk signature for predicting the clinical prognosis of patients with LUAD.

Conflicts of Interest:

Chuanxi Tian, Yikun Guo, Tianyi Lv, Daowen Yang, and Min Li declare that they have no competing interests.

Funding:

This work was supported by Ministry of Science and Technology of the People´s Republic of China (2022YFC3500801).

Author Contribution

Conceptualization: C.X.T., and D.W.Y.; Formal analysis and Data Curation: Y.K.G.; Writing - Original Draft: C.X.T., Y.K.G., T.Y.L.; Writing - Review & Editing: C.X.T., Y.K.G., T.Y.L. All authors read and approved the final manuscript.

Data Availability

The datasets generated and/or analysed during the current study are available in the Gene Expression Omnibus (GEO) database(GSE131907) and the Cancer Genome Atlas (TCGA) database(GSE31210).

Chen J, Yang H, Teo ASM, et al. Genomic landscape of lung adenocarcinoma in East Asians. Nat Genet. 2020;52(2):177–186. doi:10.1038/s41588-019-0569-6
Liang J, Bi G, Huang Y, et al. MAFF confers vulnerability to cisplatin-based and ionizing radiation treatments by modulating ferroptosis and cell cycle progression in lung adenocarcinoma. Drug Resist Updat. 2024;73:101057. doi:10.1016/j.drup.2024.101057
Thai AA, Solomon BJ, Sequist LV, et al. Lung cancer. Lancet. 2021; 398: 535–554. doi:10.1016/S0140-6736(21)00312-3
Haga Y, Sakamoto Y, Kajiya K, et al. Whole-genome sequencing reveals the molecular implications of the stepwise progression of lung adenocarcinoma. Nat Commun. 2023;14(1):8375. doi:10.1038/s41467-023-43732-y
Wang C, Yu Q, Song T, et al. The heterogeneous immune landscape between lung adenocarcinoma and squamous carcinoma revealed by single-cell RNA sequencing. Signal Transduct Target Ther. 2022;7(1):289. doi:10.1038/s41392-022-01130-8
Shiraishi K, Takahashi A, Momozawa Y, et al. Identification of telomere maintenance gene variations related to lung adenocarcinoma risk by genome-wide association and whole genome sequencing analyses. Cancer Commun (Lond). 2024;44(2):287–293. doi:10.1002/cac2.12498
Wu SY, Fu T, Jiang YZ, Shao ZM. Natural killer cells in cancer biology and therapy. Mol Cancer. 2020;19(1):120. doi:10.1186/s12943-020-01238-x
Schmidt L, Eskiocak B, Kohn R, Dang C, Joshi NS, DuPage M, Lee DY, Jacks T. Enhanced adaptive immune responses in lung adenocarcinoma through natural killer cell stimulation. Proc Natl Acad Sci U S A. 2019;116(35):17460–17469. doi: 10.1073/pnas.1904253116.
Park MD, Reyes-Torres I, LeBerichel J, et al. TREM2 macrophages drive NK cell paucity and dysfunction in lung cancer. Nat Immunol. 2023;24(5):792–801. doi:10.1038/s41590-023-01475-4
Zhang X, Wang K, Zhao W, et al. TRAF3IP3 at the trans-Golgi network regulates NKT2 maturation via the MEK/ERK signaling pathway. Cell Mol Immunol. 2020;17(4):395–406. doi:10.1038/s41423-019-0234-0
Chen H, Duan X, Deng X, et al. EBV-Upregulated B7-H3 Inhibits NK cell-Mediated Antitumor Function and Contributes to Nasopharyngeal Carcinoma Progression. Cancer Immunol Res. 2023;11(6):830–846. doi:10.1158/2326-6066.CIR-22-0374
Ding S, Qiao N, Zhu Q, et al. Single-cell atlas reveals a distinct immune profile fostered by T cell-B cell crosstalk in triple negative breast cancer. Cancer Commun (Lond). 2023;43(6):661–684. doi:10.1002/cac2.12429
Wu J, Zhou J, Xu Q, et al. Identification of Key Genes Driving Tumor Associated Macrophage Migration and Polarization Based on Immune Fingerprints of Lung Adenocarcinoma. Front Cell Dev Biol. 2021;9:751800. doi:10.3389/fcell.2021.751800
Liu S, Galat V, Galat Y, Lee YKA, Wainwright D, Wu J. NK cell-based cancer immunotherapy: from basic biology to clinical development. J Hematol Oncol. 2021;14(1):7. doi:10.1186/s13045-020-01014-w
Maddineni S, Silberstein JL, Sunwoo JB. Emerging NK cell therapies for cancer and the promise of next generation engineering of iPSC-derived NK cells. J Immunother Cancer. 2022;10(5):e004693. doi:10.1136/jitc-2022-004693
Zheng QW, Ni QZ, Zhu B, et al. PPDPF promotes lung adenocarcinoma progression via inhibiting apoptosis and NK cell-mediated cytotoxicity through STAT3. Oncogene. 2022;41(36):4244–4256. doi:10.1038/s41388-022-02418-3
Wan R, Bai L, Cai C, et al. Discovery of tumor immune infiltration-related snoRNAs for predicting tumor immune microenvironment status and prognosis in lung adenocarcinoma. Comput Struct Biotechnol J. 2021;19:6386–6399.doi:10.1016/j.csbj.2021.11.032
Sanchez-Vega F, Mina M, Armenia J, Chatila WK, Luna A, La KC, et al.. Oncogenic signaling pathways in the cancer genome atlas. Cell (2018) 173(2):321–337.e10. doi: 10.1016/j.cell.2018.03.035.
Butler A, Hoffman P, Smibert P, Papalexi E, Satija R. Integrating single-cell transcriptomic data across different conditions, technologies, and species. Nat Biotechnol (2018) 36(5):411–20. doi: 10.1038/nbt.4096.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an r package for comparing biological themes among gene clusters. OMICS (2012) 16(5):284–7. doi: 10.1089/omi.2011.0118.
Gao R, Bai S, Henderson YC, Lin Y, Schalck A, Yan Y, et al.. Delineating copy number and clonal substructure in human tumors from single-cell transcriptomes. Nat Biotechnol (2021) 39(5):599–608. doi: 10.1038/s41587-020-00795-2.
Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, Smyth GK, et al.. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res (2015) 43(7):e47. doi: 10.1093/nar/gkv007.
Chen B, Khodadoust MS, Liu CL, Newman AM, Alizadeh AA. Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol Biol (2018) 1711:243–59. doi: 10.1007/978-1-4939-7493-1_12.
Zhang S, Tong YX, Zhang XH, Zhang YJ, Xu XS, Xiao AT, et al.. A novel and validated nomogram to predict overall survival for gastric neuroendocrine neoplasms. J Cancer (2019) 10(24):5944–54. doi: 10.7150/jca.35785.
Mariathasan S, Turley SJ, Nickles D, Castiglioni A, Yuen K, Wang Y, et al.. TGFβ attenuates tumour response to PD-L1 blockade by contributing to exclusion of T cells. Nature (2018) 554(7693):544–8. doi: 10.1038/nature25501.
Hugo W, Zaretsky JM, Sun L, Song C, Moreno BH, Hu-Lieskovan S, et al.. Genomic and transcriptomic features of response to anti-PD-1 therapy in metastatic melanoma. Cell (2016) 165(1):35–44. doi: 10.1016/j.cell.2016.02.065.
Bald T, Krummel MF, Smyth MJ, Barry KC. The NK cell-cancer cycle: advances and new challenges in NK cell-based immunotherapies. Nat Immunol. 2020;21(8):835–847. doi:10.1038/s41590-020-0728-z
Merino A, Maakaron J, Bachanova V. Advances in NK cell therapy for hematologic malignancies: NK source, persistence and tumor targeting. Blood Rev. 2023;60:101073. doi:10.1016/j.blre.2023.101073
Ni J, Wang X, Stojanovic A, et al. Single-Cell RNA Sequencing of Tumor-Infiltrating NK Cells Reveals that Inhibition of Transcription Factor HIF-1α Unleashes NK Cell Activity. Immunity. 2020;52(6):1075–1087.e8. doi:10.1016/j.immuni.2020.05.001
Sheffer M, Lowry E, Beelen N, et al. Genome-scale screens identify factors regulating tumor cell responses to natural killer cells. Nat Genet. 2021;53(8):1196–1206. doi:10.1038/s41588-021-00889-w
Quintanal-Villalonga A, Taniguchi H, Zhan YA, et al. Multiomic Analysis of Lung Tumors Defines Pathways Activated in Neuroendocrine Transformation. Cancer Discov. 2021;11(12):3028–3047. doi:10.1158/2159-8290.CD-20-1863
Xiang H, Pan Y, Sze MA, et al. Single-cell analysis identifies NOTCH3-mediated interactions between stromal cells that promote microenvironment remodeling and invasion in lung adenocarcinoma. Cancer Res. Published online February 9, 2024. doi:10.1158/0008-5472.CAN-23-1183
Wu X, Zhao J, Yang L, et al. Next-generation Sequencing Reveals Age-dependent Genetic Underpinnings in Lung adenocarcinoma. J Cancer. 2022;13(5):1565–1572. doi:10.7150/jca.65370
Ye R, Yu Y, Zhao R, Han Y, Lu S. Comprehensive molecular characterizations of stage I-III lung adenocarcinoma with tumor spread through air spaces. Front Genet. 2023;14:1101443. doi:10.3389/fgene.2023.1101443
Szeto CW, Kurzrock R, Kato S, et al. Association of differential expression of immunoregulatory molecules and presence of targetable mutations may inform rational design of clinical trials. ESMO Open. 2022;7(1):100396. doi:10.1016/j.esmoop.2022.100396
Yu J, Lan L, Liu C, Zhu X. Improved prediction of prognosis and therapy response for lung adenocarcinoma after identification of DNA-directed RNA polymerase-associated lncRNAs. J Cancer Res Clin Oncol. 2023;149(14):12737–12754. doi:10.1007/s00432-023-05118-x
Zhang P, Pei S, Zhou G, Zhang M, Zhang L, Zhang Z. Purine metabolism in lung adenocarcinoma: A single-cell analysis revealing prognostic and immunotherapeutic insights. J Cell Mol Med. 2024;28(8):e18284. doi:10.1111/jcmm.18284
Lin A, Fang J, Cheng Q, Liu Z, Luo P, Zhang J. B Cell Receptor Signaling Pathway Mutation as Prognosis Predictor of Immune Checkpoint Inhibitors in Lung Adenocarcinoma by Bioinformatic Analysis. J Inflamm Res. 2022;15:5541–5555. doi:10.2147/JIR.S379016
Yang Y, Huang T, Fan Y, et al. Significance of Spliceosome-Related Genes in the Prediction of Prognosis and Treatment Strategies for Lung Adenocarcinoma. Biomed Res Int. 2022;2022:1753563. doi:10.1155/2022/1753563
Li H, Sha X, Wang W, et al. Identification of lysosomal genes associated with prognosis in lung adenocarcinoma. Transl Lung Cancer Res. 2023;12(7):1477–1495. doi:10.21037/tlcr-23-14
Zheng S, Yao L, Li F, et al. Homologous recombination repair rathway and RAD54L in early-stage lung adenocarcinoma. PeerJ. 2021;9:e10680. doi:10.7717/peerj.10680
Sayaman RW, Saad M, Thorsson V, et al. Germline genetic contribution to the immune landscape of cancer. Immunity. 2021;54(2):367–386.e8. doi:10.1016/j.immuni.2021.01.011
Liu Z, Ding M, Qiu P, Pan K, Guo Q. Natural killer cell-related prognostic risk model predicts prognosis and treatment outcomes in triple-negative breast cancer. Front Immunol. 2023;14:1200282. doi:10.3389/fimmu.2023.1200282
Chen Y, Liu Y, Han L. Spatial landscape of the tumor immune microenvironment. Trends Cancer. 2023;9(6):459–460. doi:10.1016/j.trecan.2023.03.006
Wang M, Zhou Z, Wang X, Zhang C, Jiang X. Natural killer cell awakening: unleash cancer-immunity cycle against glioblastoma. Cell Death Dis. 2022;13(7):588. doi:10.1038/s41419-022-05041-y
Jin WJ, Jagodinsky JC, Vera JM, Clark PA, Zuleger CL, Erbe AK, Ong IM, Le T, Tetreault K, Berg T, Rakhmilevich AL, Kim K, Newton MA, Albertini MR, Sondel PM, Morris ZS. NK cells propagate T cell immunity following in situ tumor vaccination. Cell Rep. 2023;42(12):113556. doi: 10.1016/j.celrep.2023.113556.
Besla R, Penuel E, Del Rosario G, et al. T cell-Dependent Bispecific Therapy Enhances Innate Immune Activation and Antibody-Mediated Killing. Cancer Immunol Res. 2024;12(1):60–71. doi:10.1158/2326-6066.CIR-23-0072
Desbois M, Wang Y. Cancer-associated fibroblasts: Key players in shaping the tumor immune microenvironment. Immunol Rev. 2021;302(1):241–258. doi:10.1111/imr.12982
Kim D, Jo S, Lee D, et al. NK cells encapsulated in micro/macropore-forming hydrogels via 3D bioprinting for tumor immunotherapy. Biomater Res. 2023;27(1):60. doi:10.1186/s40824-023-00403-9
Tong L, Jiménez-Cortegana C, Tay AHM, Wickström S, Galluzzi L, Lundqvist A. NK cells and solid tumors: therapeutic potential and persisting obstacles. Mol Cancer. 2022;21(1):206. doi:10.1186/s12943-022-01672-z
Liu WN, So WY, Harden SL, et al. Successful targeting of PD-1/PD-L1 with chimeric antigen receptor-natural killer cells and nivolumab in a humanized mouse cancer model. Sci Adv. 2022;8(47):eadd1187. doi:10.1126/sciadv.add1187

No competing interests reported.

SupplementaryMaterials.docx

Download PDF

Reviews received at journal
10 Sep, 2024
Reviewers agreed at journal
29 Aug, 2024
Reviewers agreed at journal
28 Aug, 2024
Reviewers invited by journal
27 Aug, 2024
Editor assigned by journal
20 Aug, 2024
Editor invited by journal
20 Aug, 2024
Submission checks completed at journal
16 Aug, 2024
First submitted to journal
01 Aug, 2024

You are reading this latest preprint version

Characterization of natural killer (NK) cells in lung adenocarcinoma and construction of an NK risk signature based on single-cell and macromolecular RNA-seg data

Status:

Version 1

Abstract

Figures

1 Introduction

2 Materials and Methods

2.1 Data Acquisition and Organization

2.2 Definition of NK

2.3 Identification of NK hub genes

2.4 Immune landscape analysis

2.5 Constructing risk characteristics and column line plots

2.6 Responsiveness to immune checkpoint blocks

2.7 Statistical analyses

3 Results

3.1 Screening and Identification of NK Clusters in scRNA-seq Data

3.2 Exploration of Cancer-Related Pathways in NK Clusters

3.3 Identification of NK cell-related core genes

3.4 Mutation and pathway analysis of key genes

3.5 Analysis of key genes in relation to immunity

3.6 Sensitivity analysis of risk signature for PD-L1 blockade immunotherapy

3.7 Independent risk factor identification and histogram construction

4 Discussion

5 Conclusion

Declarations

Conflicts of Interest:

Funding:

Author Contribution

Data Availability

References

Additional Declarations

Supplementary Files

Status:

Version 1