Dexin Zhang1*
Yujie Zhao2- 1Respiratory Department of the Second Affiliated Hospital of Xi’an Jiaotong University Medical College, Xi’an, China
- 2Regional Marketing Department, Yuce Biotechnology Co, Ltd., Dabaihui Center, Shenzhen, China
Introduction: This research explored the immune characteristics of natural killer (NK) cells in lung adenocarcinoma (LUAD) and their predictive role on patient survival and immunotherapy response.
Material and methods: Molecular subtyping of LUAD samples was performed by evaluating NK cell-associated pathways and genes in The Cancer Genome Atlas (TCGA) dataset using consistent clustering. 12 programmed cell death (PCD) patterns were acquired from previous study. Riskscore prognostic models were constructed using Least absolute shrinkage and selection operator (Lasso) and Cox regression. The model stability was validated in Gene Expression Omnibus database (GEO).
Results: We classified LUAD into three different molecular subgroups based on NK cell-related genes, with the worst prognosis in C1 patients and the optimal in C3. Homologous Recombination Defects, purity and ploidy, TMB, LOH, Aneuploidy Score, were the most high-expressed in C1 and the least expressed in C3. ImmuneScore was the highest in C3 type, suggesting greater immune infiltration in C3 subtype. C1 subtypes had higher TIDE scores, indicating that C1 subtypes may benefit less from immunotherapy. Generally, C3 subtype presented highest PCD patterns scores. With four genes, ANLN, FAM83A, RHOV and PARP15, we constructed a LUAD risk prediction model with significant differences in immune cell composition, cell cycle related pathways between the two risk groups. Samples in C1 and high group were more sensitive to chemotherapy drug. The score of PCD were differences in high- and low-groups. Finally, we combined Riskscore and clinical features to improve the performance of the prediction model, and the calibration curve and decision curve verified that the great robustness of the model.
Conclusion: We identified three stable molecular subtypes of LUAD and constructed a prognostic model based on NK cell-related genes, maybe have a greater potential for application in predicting immunotherapy response and patient prognosis.
1 Background
Lung cancer is a leading cause of cancer mortality in the world (Hirsch et al., 2017). Statistics reported that in 2022 in the United States will die from cancer, and approximately 350 of them die from lung cancer every day (Siegel et al., 2022). Adenocarcinoma (lung adenocarcinoma, LUAD) is currently the predominant histologic type, which accounts for approximately 50% of all lung cancer cases, and is notable for its high incidence, high mortality, and poor prognosis (Succony et al., 2021). Currently, surgery is recommended for early-stage lung cancer and is considered the most effective treatment option, while those with advanced disease are often further supplemented with radiotherapy, chemotherapy, targeted therapy, and immunotherapy (Hoy et al., 2019). Regardless of the interventions used, the overall 5-year survival of LUAD patients remains below 20% (Duma et al., 2019). Therefore, it is necessary to develop current understanding on the pathogenesis of LUAD to provide a theoretical basis for reducing the occurrence of LUAD, improving the treatment of LUAD and its prognosis.
The development of LUAD involves external environment, gene mutation, tumor immunity, and family genetics, and is a multistep, cascade process (Suster and Mino-Kenudson, 2020). As a component of the tumor microenvironment, tumor immune cells are present in all stages of LUAD and play an important role in shaping tumor development (Saab et al., 2020). For example, tumor-associated macrophages can accelerate tumor progression by promoting tumor angiogenesis, metastasis and immune escape. Regulatory T cells inhibit anti-tumor immune responses, thereby promoting the development of immunosuppressive tumor microenvironments and promoting cancer progression (Hsieh et al., 2012). Cytotoxic CD8+ memory T cells kill tumor cells by recognizing specific antigens on them and stimulating an immune response (Arneth, 2019). Dendritic cells are antigen-presenting cells, which are an important bridge between innate and adaptive immunity. Dendritic cells can not only induce cellular immunity and humoral immunity, but also activate natural killer (NK) cells and NK T cells (Sadeghzadeh et al., 2020). NK cells are anti-tumor immune cells that kill cancer cells in the body, but in the tumor microenvironment NK cells are generally reduced in number and impaired in function (Russell et al., 2022). Basic experiments and clinical studies together have shown that NK cells are in the first line of defense against tumors and do not require pre-stimulation to cause NK cells to migrate to the lesion and play an immunomodulatory role (Guillerey, 2020). Phenotypically, NK cell subpopulations display potent antitumor immune cytotoxicity via MEK/ERK and PI3K/Akt/mTOR pathways upon stimulation through cytokines such as interleukin (IL) (Valipour et al., 2019). Although patient’s immune system can recognize neoantigens produced by tumors with high mutational load (immunogenic “hot” tumors), in terms of its mutational load, lung cancer is immunogenic, only moderately, to some extent. Therefore, the highly complex interaction between LUAD and NK cells is a major challenge to improve immunotherapy.
Studies on the pathogenesis of NK cells in LUAD have delved into the genetic-molecular field, and it is mostly believed that the development of LUAD is the result of a multigene, multistage involvement (Crinier et al., 2020). However, the genetic landscape and immune profile of NK cells in LUAD are unclear, also the prognosis and immune efficacy prediction of LUAD based on NK cells have not been reported. This study first identified stable molecular subtypes of LUAD by consistent clustering of NK cell-associated genes, and further compared clinicopathological, mutational, immunological, and pathway characteristics among the subtypes. Then, we constructed a risk model and a clinical prognostic model, which can be used to evaluate personalized treatment for LUAD patients.
2 Materials and methods
2.1 Source of clinical information and gene expression profile data of NK cells
The clinical information and mRNA transcriptome data of LUAD patients were downloaded from the TCGA GDC API (Colaprico et al., 2016). To verify the accuracy of the results, we also downloaded the clinical and mRNA gene expression data of LUAD patients from the Gene Expression Omnibus database (GEO) database (Toro-Domínguez et al., 2019), including GSE72094, GSE31210 datasets. The TCGA dataset contained 500 LUAD samples as the training set, while the GSE72094 and GSE31210 datasets contained 398 and 226 LUAD samples, respectively, as the validation set.
To ensure the quality and reliability of the downloaded data, quality control was performed, and the inclusion and exclusion criteria were (Hirsch et al., 2017) to remove samples with incomplete clinical information; (Siegel et al., 2022); to remove samples with unknown survival time and survival status; (Succony et al., 2021); to remove probes with one probe matching to multiple genes, and the mean value was taken as the expression value of that gene when multiple probes matched to one gene.
NK cell-associated genes were obtained from three parts, including the ImmPort official website (https://www.immport.org/resource), the MSigDB database (Molecular Signatures Database, https://ngdc.cncb.ac.cn/databasecommons/database/id/1077) and the LM22 database (Newman et al., 2015), containing 134 cell-associated genes, 18 NK cell-associated pathways, and 79 NK cell-associated genes, respectively.
2.2 Subtyping of LUAD patients based on NK cell-associated genes
A total of 213 NK cell-associated genes and 18 NK cell-associated pathways were obtained from the three databases, and we used the single sample gene set enrichment analysis (ssGSEA) method to evaluate these 213 NK cell-associated genes and 18 NK cell-associated pathways in the TCGA and GEO datasets, respectively. The samples were then clustered by ConsensusClusterPlus using these pathway scores in the TCGA and GEO cohorts, and the “K-M” algorithm and “1-Pearson correlation” as the metric distance (Azman et al., 2006). We conducted 500 bootstraps, with each one including 80% patients of in the training set and 20% those of the validation set. Finally, based on the cumulative distribution function (CDF), the optimal number clusters were decided, and the optimal classification and the sample molecular subtyping was obtained by calculating the consistency matrix and the consistency cumulative distribution function (Zhang et al., 2021a).
2.3 Immunological features and pathway analysis among different molecular subtypes
We obtained the molecular characteristics of LUAD genomic alterations from published literature, including LOH, Aneuploidy Score, tumor mutation burden (TMB), purity, and ploidy, Homologous Recombination Defects, Intratumor heterogeneity. The relative abundance of 22 immune cells were calculated using CIBERSORT R package. At the same time, we used the ESTIMATE algorithm R package to calculate the proportion of immune cells and finally compared the inflammatory and immune activity scores (Chakraborty and Hossain, 2018; Chen et al., 2018).
We performed gene set enrichment analysis (GSEA) on all NK cell-associated genes in the Hallmark database, and then used the ssGSEA method to calculate the pathway scores for both TCGA and GEO datasets in the GSVA package (Barbie et al., 2009). A false discovery rate (FDR) of <0.05 in this study was considered statistically significant.
2.4 Drug sensitivity analysis between molecular subtypes
Immune checkpoint inhibitor (ICI)-based therapy has become one of the standard treatments for advanced lung cancer (Zhang et al., 2021b). We first assessed the expression of genes associated with immunotherapy, such as CTLA4, PD-L1, and PD-1, among various molecular subtypes to determine whether there were differences in immunotherapy responsiveness among them. Next, we used the TIDE software (http://tide.dfci.harvard.edu/) to assess the potential clinical effects of immunotherapy in our defined molecular subtypes. Greater likelihood of immune escape was correlated with a higher TIDE prediction score, suggesting that patients may benefit less from immunotherapy (Jiang et al., 2018). Finally, we performed drug sensitivity prediction for LUAD in the “pRRophetic” package (Geeleher et al., 2014).
2.5 Identification of key NK cell-related genes among molecular subtypes
The differentially expressed genes among different molecular typing were calculated by the “limma” package, using FDR <0.05 and |log2FC| > 1 as the statistical difference criteria, and visualized the differentially expressed genes by “pheatmap” and “ggplot2” R packages in a heatmap and volcano map. Then, all genes with statistically significant differences were enriched using the “clusterProfiler” package.
Next, we performed univariate Cox regression analysis for differentially expressed genes between molecular subtypes, and then reduced the prognosis-related genes by Lasso regression (Sun et al., 2021), which can better solve the problem of multicollinearity in regression analysis by compressing some coefficients and setting some coefficients to zero at the same time. With the gradual increase of lambda, we selected the number of factors when the coefficients of independent variables tended to zero. Then, we used the AIC deficit pool information criterion through stepwise regression, which has the advantage of considering the statistical fit of the model and the number of parameters used to fit it, and at the same time indicates a sufficient fit of the model obtains with fewer parameters (Zhang, 2016).
2.6 Construction and validation of the prognostic model
We calculated the NK cell-related prognostic RiskScore for each sample according to the formula defined by the sample RiskScore (below), and normalized it (Nie et al., 2021).
After that, LUAD patients were divided into high- and low-risk groups based on the relationship between RiskScore and 0, where those with RiskScore >0 were considered as having a high risk and those with RiskScore <0 were considered as having a low risk. Finally, the survival differences between the two groups were compared by log-rank test. In order to verify the robustness of the model, we performed immune signature analysis, survival curve, and drug treatment difference analysis for the patients in the two groups.
2.7 Improvement of prognostic models and survival prediction in LUAD patients
To more accurately quantify the risk assessment and survival probability of LUAD patients, we combined the RiskScore with other clinicopathological characteristics of LUAD patients and constructed a nomogram using the “nomogramEx” R package. To validate the accuracy of the model, a calibration curve was plotted by the “PredictABEL” function to visualize the goodness-of-fit. This was followed by decision curve analysis (DCA) to describe the change in net benefit as the threshold probability changed under the intervention of the predicted value by the model (Van Calster et al., 2016; Van Calster et al., 2018).
2.8 Programmed cell death (PCD) analysis
12 PCD patterns (apoptosis, necroptosis, pyroptosis, ferroptosis, cuproptosis, entotic cell death, netotic cell death, parthanatos, lysosome-dependent cell death, autophagy-dependent cell death, alkaliptosis, and oxeiptosis) have been taken from the previous research (Zou et al., 2022). ssGSEA analysis based on the expression data of PCD related genes using the R package GSVA. Spearman analysis was conducted to know the relationship among PCD patterns, clinical features, RiskScore in LUAD samples.
2.9 Statistical analysis
Unless otherwise specified, all statistical tests were bilateral and conducted using R software (version 4.1.3, https://www.r-project.org/), and p < 0.05 was considered statistically significant.
3 Results
3.1 Molecular subtyping of LUAD based on NK cell-associated genes
We first calculated the NK cell-related genes showing close relationship with LUAD survival chance by univariate Cox regression analysis, and screened 63 prognostically significant genes (p < 0.05, Figure 1A), including the prognostic (Protective) genes SHC1, TICAM1, PVR, RAET1E, RAC1 (HR > 1), and KLRB1, CD160, KIR3DL2, CLEC12B, and KIR2DL1 (HR < 1). Then, we used these 63 differential genes for consistent clustering, and determined the best cluster number according to CDF. And we could see from Figures 1B, C that Cluster = 3 had more stable clustering results, hence, k = 3 was selected to obtain three molecular subtypes (C1, C2, and C3) (Figure 1D). Then, we performed survival analysis of patients with these three molecular subtypes using the K-M survival method, and the results identified a significant difference in prognosis among the three molecular subtypes, with C1 patients having the worst prognosis and C3 patients having the optimal prognosis (Figure 1E). The results were also validated in the GSE72094 dataset (Figure 1F). Meanwhile, we found that the “Risk” genes were high-expressed in the C1 subtype and the “Protective” genes were high-expressed in the C3 subtype in the heat map (Figure 1G). These results suggested that the molecular subtyping based on NK cell-related genes was reasonable, and there were significant differences in gene expression and prognosis among patients with different subtypes.
FIGURE 1. Molecular subtyping based on natural killer cell-associated genes. (A) Forest plot of prognostically significant natural killer cell-associated genes in the TCGA-LUAD cohort. (B) CDF curves of the TCGA-LUAD cohort. (C) CDF Delta area curves of the TCGA-LUAD cohort. (D) heat map of sample clustering at consensus k = 3 in the TCGA-LUAD cohort. (E) KM curves of the relationship between overall survival (OS) prognosis of the three subtypes in the TCGA-LUAD cohort. (F) Prognostic differences between the three molecular subtypes in the GSE72094 cohort. (G) Heatmap of prognosis significant natural killer cell genes expression in different subtypes of TCGA-LUAD.
3.2 Genetic landscape between molecular subtypes of LUAD
To explore the differences in specific gene expression profiles among different molecular subtypes, we compared the molecular profiles among C1, C2, and C3 subtypes of LUAD samples, respectively, and it is obvious from Figure 2A that purity, and ploidy, TMB, Aneuploidy Score, LOH, Homologous Recombination Defects expression were the highest in C1 but the lowest in C3, which was consistent with previous studies (Thorsson et al., 2018). In addition, we compared the differences between the molecular subtyping of published studies and that in this study. Here it was found that the C3 subclass occupied the most of the C3 subtypes we defined, suggesting that the C3 subtype was the major subtype of LUAD (Figure 2B). In addition, a significant correlation between molecular subtypes and gene mutations was detected after analyzing the correlation between gene mutations and molecular subtypes, and observed. TTN, MUC16, CSMD3, and RYR2 were the most widely mutated genes in LUAD (Figure 2C), and this finding indicated that the development of LUAD was closely related to the above-mentioned gene mutations.
FIGURE 2. Genomic alterations in the molecular subtypes of the TCGA cohort. (A) Comparison of differences in Homologous Recombination Defects, Aneuploidy Score, Fraction Altered, Number of Segments, and Tumor mutation burden in the TCGA cohort molecular subtypes. (B) Comparison of the three molecular subtypes with immune molecular subtypes. (C) Somatic mutations in the three molecular subtypes (chi-square test). *p < 0.05; **p < 0.01; ***p < 0.001; ****p < 0.0001.
3.3 Pathways enrichment analysis among the molecular subtyping of LUAD
To investigate pathway differences in LUAD among different molecular subtypes, we performed GSEA enrichment analysis among molecular subtypes. As shown in Figure 3A, we enriched a total of 33 significant pathways in the TCGA-LUAD dataset, including MYC_TARGETS_V2, E2F_TARGETS, NFLAMMATORY_ RESPONSE, MYOGENESIS, INTERFERON_GAMMA_RESPONSE, MYC_TARGETS_V1, GLYCOLYSIS, G2M_CHECKPOINT, EPITHELIAL_MESENCHYMAL_TRANSITION, ALLOGRAFT_REJECTION, suggesting that these NK cell genes were mainly associated with cell cycle and immunity in C1 and C3. Additionally, pathways different between C1 and C3 subtypes, between C2 and C3 subtypes, between C1 and C2, were analyzed (Figure 3B). Overall, the cell cycle pathway was activated in C1 patients, while the immune-related pathway was suppressed, therefore we hypothesized that these NK cell genes might play an important role in the cell cycle pathway as well as in the tumor microenvironment. To validate these results, we presented the pathway differences between C1 and C2, and C2 and C3 as radar plots, and the results showed that they both had significant consistency in cell cycle (MYC_TARGETS_V2, MTORC1_SIGNALING, MYC_TARGETS_V1) and immune-related pathways (G2M_CHECKPOINT, E2F_TARGETS, UNFOLDED_PROTEIN_ RESPONSE) (Figure 3C).
FIGURE 3. Pathway analysis between molecular subtypes. (A) Bubble plots of GSEA results for C1 vs. C3 subtypes in two lung adenocarcinoma cohorts. (B) Bubble plots of GSEA results for different molecular subtypes compared to each other in the TCGA-LUAD cohort. (C) Radar plots of C1 vs. C2 and C2 vs. C3 activation pathways in the TCGA-LUAD cohort.
3.4 Immune characteristics among different molecular typologies of LUAD
The immune system plays a dual role in the development of LUAD, as it can recognize and destroy tumor cells, while tumor cells can also evade host immune attack by forming a complex immunosuppressive network under the pressure of immune selection using the immune system’s own negative regulatory mechanisms, thus the TME is in a constant state of change (Anichini et al., 2020; Spella and Stathopoulos, 2021). To explore the immune landscape among different molecular subtypes of LUAD, we first assessed the differences in the components of immune cells in the TCGA-LUAD cohort using the CIBERSORT algorithm and observed that most immune cells (B cells, T cells, NK cells, etc.) were significantly different (p < 0.05) (Figure 4A). We then used the ESTIMATE algorithm to assess immune cell infiltration, and the results showed that StromalScore, ImmuneScore and EstimateScore were significantly different between C1, C2, and C3 (p < 0.05)), with ImmuneScore accounting for the largest proportion of C3 types, suggesting a higher degree of immune infiltration in C3 subtypes (Figure 4B). Similarly, we obtained results in the GSE72094-LUAD cohort that were consistent with the TCGA-LUAD cohort (Figures 4C, D). In addition, we assessed the inflammatory activity of C3, C2, C1, except for IgG, the remaining six out of 7 metagenes clusters (HCK, Interferon, LCK, MCH I, MCH II, and STAT1) showed significantly different enrichment scores, with the C4 subtype having higher inflammatory activity (Figure 4E). The findings were consistent in the GSE72094-LUAD cohort (Figure 4F).
FIGURE 4. Proportions of immune cell components in the two lung adenocarcinoma cohorts. (A) Differences in 22 immune cell scores between different molecular subtypes in the TCGA-LUAD cohort. (B) Differences in ESTIMATE immune infiltration between different molecular subtypes in the TCGA-LUAD cohort. (C) Differences in the GSE72094 cohort 22 immune cell scores between different molecular subtypes. (D) Differences in ESTIMATE immune infiltration between different molecular subtypes in the GSE72094 cohort. (E) Differences in seven inflammation-related gene cluster scores across molecular subtypes in TCGA-LUAD cohort. (F) Differences in gene cluster scores between different molecular subtypes in seven inflammatory-related genes in GSE72094 cohort.
3.5 Differences in immunotherapy between molecular subtypes
In recent years, immunotherapy has led to new opportunities in the treatment of small cell lung cancer. Clinical trials of some immune checkpoint inhibitors have demonstrated their efficacy and safety in LUAD (Hua et al., 2021). Based on this, we first evaluated the expression of the representative molecules of immunotherapy (PD-1, PD-L1, CTLA4) among the three molecular subtypes, and observed that PD-1, PD-L1, and CTLA4 were significantly more expressed in C3 subtype (p < 0.05) (Figure 5A). We also applied the “T-cell-inflamed GEP score” to assess the predictive potential of different molecular subtypes to cancer immunotherapy, and the results also showed that the score was highest in C3 (Figure 5B). Considering that IFN-γ is a cytokine that plays a key role in immunomodulation and immunotherapy, we downloaded the GOBP_RESPONSE_TO_INTERFERON_GAMMA gene set from the GO database for ssGSEA analysis, and found that the IFN-γ response was significantly enhanced in the C1 subtype (Figure 5C). We also compared the differences in INFG gene expression in the three subtypes and found that INFG was noticeably high-expressed in the C3 subtype (Figure 5D). Moreover, CYT score, which reflects the cytotoxic effect, was significantly higher in the C3 subtype than in the other subtypes (Figure 5E). In addition, the TIDE prediction data indicated that the C1 subtype had a higher TIDE score, suggesting that the C1 subtype was less likely to benefit from immunotherapy (Figure 5F). The estimated IC50 of docetaxel, vincristine, paclitaxel and cisplatin among 3 subtypes showed that C1 was more sensitive to the four chemotherapy drugs (Figure 5G). The above results indicated that predicting immunotherapy for LUAD based on NK cell-related genes was a practical approach.
FIGURE 5. Differences in immunotherapy/chemotherapy treatment between molecular subtypes. (A) Differences in expression of immune checkpoint genes between molecular subtypes. (B) difference in “T cell inflamed GEP score” between molecular subtypes. (C) Differences in “response to IFN-gamma” between different molecular subtypes. (D) Differences in INFG gene expression in different isoforms. (E) Differences in “Cytolytic activity” between different molecular subtypes. (F) Differences in TIDE scores between subtypes. (G) Box plots of estimated IC50 of docetaxel, vincristine, paclitaxel and cisplatin in TCGA-LUAD.
3.6 The analysis of PCD patterns among molecular subtypes
The ssGSEA analysis calculated the score of 12 PCD patterns in each sample in TCGA dataset and GSE72094 dataset. We found that 9 PCD scores had differences among 3 subtypes in both two datasets (Figures 6A, B). In TCGA dataset, Stage, Gender, especially, Age had closely associated to PCD patterns (Figure 6C), but in GSE72094 dataset, clinical features had litter associated to PCD patterns (Figure 6D). Autophagy score were increased in early Stage, the scores of Pyroptosis, Autophagy, Necroptosis and Oxeiptosis were enhanced in Male samples, and samples with age > 60 had higher Pyroptosis, Entotic. cell.death scores in TCGA dataset (Figure 6E). In GSE72094 dataset, Oxeiptosis score was highest in StageⅢ, and Ferroptosis and Necroptosis scores were greater in patients with age>60 (Figure 6F).
FIGURE 6. The PCD characteristic among 3 molecular subtypes. (A) the ssGSEA analysis of 12 PCD patterns among 3 molecular subtypes in TCGA-LUAD dataset. (B) the ssGSEA analysis of 12 PCD patterns among three molecular subtypes in GSE72094 dataset. (C) The spearman analysis between clinical feature and PCD in TCGA-LUAD dataset. (D) The spearman analysis between clinical feature and PCD in GSE72094 dataset. (E) The ssGSEA analysis of PCD in TCGA-LUAD samples with Stage, Gender and Age. (F) The ssGSEA analysis of PCD in GSE72094 samples with Stage, and Age.
3.7 Establishment of LUAD risk model
We first calculated the NK cell-related genes significantly differentially expressed among the three molecular subtypes by the limma package, significant expression differences of NK cell-related genes among C1, C2, and C3 were detected, including 11 upregulated genes and 180 downregulated genes (Supplementary Figures S1A, B). Differentially expressed downregulated genes were related to immune-related pathways, as shown by the results of enrichment analysis (Supplementary Figure S1C). Genes with upregulated level were related to inflammatory and immune pathways (Supplementary Figure S1D). 173 genes with high prognostic impact (p < 0.05), including 159“Protective” and 14“Risk” genes, were identified from those genes by conducting one-way Cox regression analysis (Supplementary Figure S2A). Further, we observed the trajectory of each gene with lambda using Lasso analysis, and the model was optimal when lambda = 0.0382, which corresponded to 9 differential genes (Supplementary Figures S2B, C). After that, we reduced the genes to four, namely, ANLN, FAM83A, RHOV, and PARP15, by the stepAIC method in the MASS package (Supplementary Figure S2D).
Then, we calculated the Riskscore score for each TCGA-LUAD patient using these four genes and the above formula (Figure 7A). We classified those RiskScore with 0 ≤ as low-risk group and with RiskScore >0 as high-risk group. Then, we performed a prognostic classification ROC analysis in the “timeROC” package for analyzing 1-year, 2-year, 3-year, and 5-year prognostic prediction classification efficiency, and we found that the model had a high AUC (0.71, 0.69, 0.7, and 0.67) (Figure 7B). The results of survival analysis showed that patients in the low-risk group developed a significantly better prognosis (p < 0.001) (Figure 7C). To confirm the robustness of this clinical prognostic model, we validated it in the GSE72094 and GSE31210 cohorts and used the same approach to calculate the RiskScore of patients (Figures 7D–G).
FIGURE 7. Risk modeling and validation. (A) RiskScore, survival time vs. survival status and expression of necroptosis-related genes in TCGA-LUAD dataset. (B) ROC curves with AUC for RiskScore classification in the TCGA-LUAD dataset. (C) Distribution of KM survival curves for RiskScore in the TCGA-LUAD dataset. (D,E): ROC curves and KM survival curves for RiskScore in the GSE72094 cohort. (F,G): ROC curves and KM survival curves of RiskScore in the GSE31210 cohort.
3.8 Pathological characteristics of high- and low-risk groups
To investigate the reliability of this risk model classification method, we first compared the clinical characteristics of patients in both high- and low-risk groups. The results showed that the RiskScore scores of patients with Stage III-IV, M Stage, N Stage, T Stage were significantly higher than Stage I-II ones. In addition, we also found that male patients had a higher RiskScore (Figure 8A). Also, we compared the differences in RiskScore by molecular subtype and found that the RiskScore for the C1 subtype with poorer prognosis was significantly higher than C3 with a better prognostic outcome (Figure 8B), and that the majority of the samples with high RiskScore were “C1” patients (Figure 8C). In addition, we also compared whether there was a prognostic difference in the—high- and low-risk groups between the different clinicopathological characteristics subgroups in the TCGA-LUAD cohort. Across different clinical subgroups, the risk grouping performed equally well, pointing to the reliability of the grouping (Figure 8D). This finding also applied to the GSE72094-LUAD cohort (Supplementary Figure S3).
FIGURE 8. Performance of RiskScore in TCGA-LUAD cohort with different clinicopathological characteristics. (A) Differences between RiskScore between different clinicopathological subgroups in the TCGA-LUAD cohort. (B) Differences in RiskScore between different molecular subtypes in the TCGA-LUAD cohort. (C) Differences between molecular subtypes and RiskScore subgroups in the TCGA-LUAD cohort. (D) KM curves between high- and low-risk groups in the TCGA-LUAD cohort between different clinicopathological subgroups.
3.9 Immune infiltration and pathway characteristics of low-risk and high-risk patients
We compared the relative abundance of 22 immune cell types in the two subgroups of the TCGA-LUAD cohort in high- and low-risk groups. We discovered that the majority of immune cells (B cells, macrophages, T cells, and mast cells) were significantly different in high- and low-risk groups (p 0.05, Figure 9A). It is worth noting that activated NK cells had no significance between high- and low-group. We also examined the connection between the RiskScore and 22 immune cell components (Figure 9B). Also, we assessed the immune cell infiltration using the ESTIMATE method. The three scores were significantly different between two risk groups (p < 0.05), and the low-Riskscore group had higher immune infiltration (Figure 9C). The relationship between biological function in different samples with RiskScore was analyzed by “ssGSEA” analysis and found that the high risk group was significantly enriched to some cell cycle-related pathways, such as HALLMARK_SPERMATOGENESIS, and HALLMARK_REPAIR, SPERMATOGENESIS, HALLMARK_DNA_REPAIR, ALLMARK_MYC_TARGETS_V2, HALLMARK_UNFOLDED_PROTEIN_RESPONSE, etc. (Figure 9D). Further, we selected functional pathways with correlations greater than 0.4, from which we could see that RiskScore showed positive correlation with cell cycle-related pathways, such as HALLMARK_MYC_TARGETS_V1, HALLMARK_MTORC1_SIGNALING, HALLMARK_E2F_TARGETS (Figure 9E).
FIGURE 9. Immune infiltration/pathway characteristics between RiskScore subgroups. (A) Proportion of immune cell components in the TCGA cohort. (B) Correlation analysis of 22 immune cell components in the TCGA cohort with the RiskScore. (C) Proportion of immune cell components in the TCGA cohort calculated by ESTIMATE software. (D) The top 10 pathways with the most significant differences between high- and low-groups. (E) Results of correlation analysis between the RiskScore and KEGG pathways with correlation greater than 0.4.
3.10 Differences in immunotherapy/chemotherapy for patients in high- and low-risk groups
First, we used the “T-cell-inflamed GEP score” to assess the predictive potential of the different RiskScore subgroups in cancer immunotherapy. The results showed that the “T-cell-inflamed GEP score” was elevated in the low-risk group, but the difference was not statistically significant (Figure 10A), however, in the low-risk group the IFN-γ response was noticeably elevated (Figure 10B). The CYT score, which reflects cytotoxic effects, was elevated in the low-risk group, showed no statistically significant differences (Figure 10C). The expression of representative molecules of immunotherapy (CTLA4, PD-L1, and PD-1) was calculated in the risk groups and showed that CTLA4 was significantly more expressed in the low-risk group (p < 0.05), while the difference in PD-1 and PD-L1 expression was not significant (Figure 10D). We looked at the connection between RiskScore and medication response in cancer cell lines to better understand the impact of RiskScore on drug response. We found 49 substantially linked relationships between RiskScore and drug sensitivity in the Genomics of Drug Sensitivity in Cancer (GDSC, http://cancer.sanger.ac.uk/cell_lines#) database using Spearman correlation analysis. Of these 49 pairs, 15 pairs were significantly associated with Riskscore correlations, such as Vinorelbine, Sabutoclax, Vinblastine, Entinostat, Vincristine, and Sorafenib (Figure 10E). We found that these drugs mainly target the EGFR signaling and TNKS2 pathways through the study on the signaling pathways of the genes targeted by these drugs (Figure 10F). In addition, we also explored the response of different molecular subtypes in the TCGA-LUAD cohort to the traditional chemotherapeutic agents Docetaxel, Vinorelbine, Paclitaxel and Cisplatin, and found that overall patients in the high-risk group were more sensitive to all the four chemotherapeutic agents (Figure 10G), suggesting that patients in the high-risk group may benefit from these four drugs.
FIGURE 10. Differences in immunotherapy/chemotherapy between RiskScore subgroups. (A) Difference in “T cell inflamed GEP score” between molecular subtypes. (B) Difference in “response to IFN-γ" between molecular subtypes. (C) Differences in “Cytolytic activity” between molecular subtypes. (D) Differences in expression of immune checkpoint genes between molecular subtypes. (E) 15 pairs drugs were significantly associated to RiskScore. (F) 15 pairs drugs mainly target EGFR signaling and TNKS2 pathways. (G) IC50 box plots of docetaxel, vincristine, paclitaxel and cisplatin in TCGA-LUAD dataset.
3.11 PCD characteristics in high- and low-risk groups
We also determine the PCD characteristics in high- and low-risk groups using ssGSEA analysis. 6 of 12 PCD styles had differences between high- and low-risk groups in TCGA dataset (Figure 11A). In GSE72094 dataset, 10 PCD patterns scores presented differentiation in high- and low-risk groups (Figure 11B). Moreover, the differences of 9 PCD scores between high- and low-groups was observed in GSE31210 dataset (Figure 11C). RiskScore as well as four model genes were obviously related to PCD patterns (Figure 11D).
FIGURE 11. The ssGSEA analysis of 12 PCD patterns in high- and low-groups. (A) ssGSEA analysis of 12 PCD in high- and low-group in TCGA-LUAD dataset. (B) ssGSEA analysis of 12 PCD in high- and low-group in GSE72094 dataset. (C) ssGSEA analysis of 12 PCD in high- and low-group in GSE31210 dataset. (D) the association among RiskScore, model genes and 12 PCD patterns.
3.12 RiskScore combined with clinicopathological features to further improve prognostic models and survival prediction
Univariate and multifactorial Cox regression analyses revealed RiskScore as the most significant prognostic factor (Figures 12A, B). We created a nomogram (Figure 12C) combining RiskScore and other clinicopathological traits for the risk assessment and prediction of survival probability for LUAD patients. The model results revealed that RiskScore had the biggest influence on survival prediction. The prediction calibration curves at the three calibration points of 1, 3, and 5 year(s) nearly overlapped with the standard curve, which indicated that the nomogram plot had excellent prediction performance. We further assessed the prediction accuracy of the model using the calibration curve (Figure 12D). We also used DCA (Decision curve) to test the model’s dependability, and it was shown that RiskScore and Nomogram performed much better than the extreme curve and had the strongest ability to predict survival among other clinicopathological factors (Figure 12E).
FIGURE 12. Establishment of nomogram. (A,B): Univariate and Multivariate Cox analysis of RiskScore and clinicopathological characteristics; (C): The nomogram model; (D): Calibration curves for 1, 3, and 5 years for the nomogram; (E): Decision curves for the nomogram.
4 Discussion
Lung cancer is currently the most aggressive malignancy in the world, of which LUAD is the most common histological subtype of primary lung cancer, accounting for 64% of peripheral lung cancers, and has been reclassified from invasive precancerous lesions to invasive adenocarcinoma (Denisenko et al., 2018; Hutchinson et al., 2019). Despite the current advances in the treatment of LUAD, the median survival is only 8.6 months and immune escape is considered one of the main factors leading to treatment failure in LUAD (Yotsukura et al., 2021). In contrast to the remarkable efficacy of immune checkpoint inhibitor (ICI) in metastatic melanoma, Hodgkin’s lymphoma, and bladder cancer, not all patients with LUAD are sensitive to ICI (Zhang et al., 2020). Mechanisms of immune escape that lack adaptive immune response include hypoxia-driven immunosuppressive factors, anti-apoptotic pathways, chronic inflammation, metabolic damage, and immune cells such as regulatory T (Treg) cells, tumor-associated M2 macrophages (TAM), myeloid-derived suppressor cells (MDSC) (Yu et al., 2021). Recent studies have shown that T and NK cell dysfunction and depletion or deficiency of antitumor-specific effector cells are involved in LUAD immune escape (Hong et al., 2019), and although the exact mechanism is unclear, it points to new ideas for the study of immune escape in LUAD and provides new targets for immunotherapy in LUAD.
LUAD is usually resistant to chemotherapy and/or radiotherapy and leads to the development of distant metastases (Jiang et al., 2021). NK cell dysfunction and failure in patients with LUAD could be caused by immune escape mechanisms mediated by lung cancer cells or tumor microenvironment, leading to failure of immunotherapy. The reason for this is related to tumor upregulation of inhibitory ligands (e.g., HLA-C molecules) and recognition by autoinhibitory KIR receptors carrying ITIM motifs (Daëron et al., 2008). Cellular experiments showed that other inhibitory receptors, for instance, KLRG-1, LAG-3, CD94/NKG2A, TIM3, TIGIT, and their ligands were also frequently upregulated on NK cells from LUAD patients (Lee et al., 1998; Nayyar et al., 2019), which was consistent with our study, where we found significantly different NK cell-related gene expression in different subtypes. CTLA-4 (ipilimumab) improved clinical prognosis of patients with LUAD (Paulsen et al., 2017) in addition to the common PDL-1 inhibitors (avelumab, atezolizumab, durvalumab) and PD-1 (camrelizumab, spartalizumab, nivolumab, pembrolizumab). Our study identified the expression patterns of PD-1/PD-L1 and CTLA-4 in different subtypes, confirming a possible immune escape mechanism of NK cells in LUAD and providing a new perspective for blocking immune dysregulation.
The tumor microenvironment (TME) consists of associated fibroblasts (CAF), tumor cells, other immune cells, and endothelial cell constituents (ECs) (Vitale et al., 2019). Ghiringhelli F et al. showed that suppressive immune cells such as Treg cells, CTLA-4+ regulatory, and that N2 neutrophils and M2 macrophages can disrupt the anti-lung cancer activity of NK cells (Domagala-Kulawik et al., 2014). Similarly, our data showed significant differences in the proportion of NK cells, B cells, and T cell content between different molecular subtypes, suggesting that other immune cells may impair the cytotoxic and migratory activity of NK cells with numerical and functional advantages, and thus causing NK cell depletion (Bi and Tian, 2017). But we found that activated NK cells had no differences between high- and low-group, maybe caused by insufficient samples.
Changes in NK cell counts, including peripheral blood, circulation and TME in healthy individuals, can be used as prognostic markers in patients with head and neck and lung tumors (Lin et al., 2017; Lin et al., 2020; Zhong et al., 2021). We constructed the prognosis model by NK cell-related genes (ANLN, FAM83A, RHOV, and PARP15), which is a powerful tool to assist clinical decision-making with effective prediction of patient survival and drug sensitivity. ANLN is an actin-binding protein, and previous studies have demonstrated that ANLN is associated with actin cytoskeleton dynamics (Xu et al., 2019). Xu J et al. showed that ANLN overexpression promotes distant metastasis of lung cancer cells and is associated with epithelial mesenchymal transformation (EMT) of LUAD cells transformation (EMT) in LUAD cells. Similar to previous bioinformatic analyses, our study found that upregulated FAM83A in LUAD tissues, which was relate to LUAD prognosis (Suzuki et al., 2005; Deng et al., 2021). Knockdown of FAM83A inhibited proliferation, migration and invasion of LUAD cells. In addition, the lncRNA FAM83A-AS1 regulates FAM83A expression by acting as a competing endogenous RNA for miR-495-3p (Wang et al., 2021). These results suggested that FAM83A plays an oncogenic role in LUAD and that FAM831-AS1 can regulate FAM83 expression by taking up miR-495-3p. Similar to FAM83A, invasion, migration and proliferation of LUAD cells could be stimulated by RHOV overexpression, while knockdown of RHOV inhibits the functionalistic behavior of the cells. In addition, RHOV knockdown inhibits metastasis and LUAD tumor growth of nude mice, which may be related to RHOV activation of the JNK/c-Jun signaling pathway (Zhang et al., 2021c). There are fewer basic studies on PARP15 in LUAD, and genomic data with large sample sizes suggested that RHOV is a useful marker for immunotherapy and survival in LUAD (Han et al., 2020). The above studies revealed a novel regulatory mechanism of NK cells in LUAD tumor development, which may be a new biomarker and therapeutic target for LUAD.
Docetaxel, Vinorelbine, Paclitaxel and Cisplatin are currently widely used chemotherapy drugs for lung cancer, which cause cell cycle arrest (Clegg et al., 2001; Dasari and Tchounwou, 2014). However, resistance can develop, leading to further tumor development and side effects such as myelosuppression, drug nephritis, nausea, vomiting, hearing loss and polyneuropathy, which will significantly reduce the patient’s quality of life (Dasari and Tchounwou, 2014). Acquired chemotherapy resistance is a major problem faced by clinicians and a major cause of treatment failure. Regardless of the type of resistance, loss of tumor sensitivity to the drug leaves very little time for therapy to correct, with the goal of improving patient survival. Patients’ clinical outcomes can be significantly improved by personalizing treatment regimens and predicting the effects of drug therapy. The results of this study showed that patients in C1 subtype and high-risk group were more sensitive to and benefited from four chemotherapy drugs. We speculated that may be the number of NK cells affects drug sensitivity.
Although this study reveals the immune signature of NK cell-related genes in LUAD and confirms the role in prognosis and immunotherapy of LUAD, the following limitations remain: (Hirsch et al., 2017): The wide variety, rapid development of bioinformatics tools can help predict potential key molecules and pathways, narrow the scope and improve the efficiency of the study, but the final findings should be validated based on real genetic data in basic and clinical settings; (Siegel et al., 2022); The database used to conduct functional and signaling pathway enrichment analysis has comprehensive and complete data, but its slow updates may have some unpredictable effects on the results; (Succony et al., 2021); The results were based on extrapolation of the raw signal algorithm and should be supported by further laboratory and clinical evidence.
5 Conclusion
Based on NK cell-related genes, we identified three stable molecular subtypes of LUAD, which differed significantly in terms of immunity, pathways, prognosis and drug sensitivity among different molecular subtypes. Based on NK cell-related genes, this study developed a prognostic model, which was highly robust and had a greater potential for application in predicting immunotherapeutic response and patient prognosis.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Materials, further inquiries can be directed to the corresponding author.
Author contributions
All authors contributed to this present work: DZ, designed the study. YZ, acquired the data. YZ, drafted the manuscript. DZ, revised the manuscript. All authors read and approved the manuscript.
Funding
This study was supported by Clinical Study on Diagnosis and Treatment of Peripheral Pulmonary Nodules by Bronchoscopic Navigation and Thoracic Wall Navigation (No. S2023-YF-YBSF-0407).
Conflict of interest
Author YZ was employed by Yuce Biotechnology Co, Ltd.
The remaining author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1156230/full#supplementary-material
References
Anichini, A., Perotti, V. E., Sgambelluri, F., and Mortarini, R. (2020). Immune escape mechanisms in non small cell lung cancer. Cancers 12 (12), 3605. doi:10.3390/cancers12123605
Arneth, B. (2019). Tumor microenvironment. Med. Kaunas. Lith. 56 (1), 15. doi:10.3390/medicina56010015
Azman, J., Frković, V., Bilić-Zulle, L., and Petrovecki, M. (2006). [Correlation and regression]. Acta medica Croat. cas. Hravatske akad. Med. znan. 60 (1), 81–91.
Barbie, D. A., Tamayo, P., Boehm, J. S., Kim, S. Y., Moody, S. E., Dunn, I. F., et al. (2009). Systematic RNA interference reveals that oncogenic KRAS-driven cancers require TBK1. Nature 462 (7269), 108–112. doi:10.1038/nature08460
Bi, J., and Tian, Z. (2017). NK cell exhaustion. Front. Immunol. 8, 760. doi:10.3389/fimmu.2017.00760
Chakraborty, H., and Hossain, A. (2018). R package to estimate intracluster correlation coefficient with confidence interval for binary data. Comput. methods programs Biomed. 155, 85–92. doi:10.1016/j.cmpb.2017.10.023
Chen, B., Khodadoust, M. S., Liu, C. L., Newman, A. M., and Alizadeh, A. A. (2018). Profiling tumor infiltrating immune cells with CIBERSORT. Methods Mol. Biol. Clift. NJ) 1711, 243–259. doi:10.1007/978-1-4939-7493-1_12
Clegg, A., Scott, D. A., Sidhu, M., Hewitson, P., and Waugh, N. (2001). A rapid and systematic review of the clinical effectiveness and cost-effectiveness of paclitaxel, docetaxel, gemcitabine and vinorelbine in non-small-cell lung cancer. Health Technol. Assess. Winch. Engl. 5 (32), 1–195. doi:10.3310/hta5320
Colaprico, A., Silva, T. C., Olsen, C., Garofano, L., Cava, C., Garolini, D., et al. (2016). TCGAbiolinks: An R/bioconductor package for integrative analysis of TCGA data. Nucleic acids Res. 44 (8), e71. doi:10.1093/nar/gkv1507
Crinier, A., Narni-Mancinelli, E., Ugolini, S., and Vivier, E. (2020). SnapShot: Natural killer cells. Cell. 180 (6), 1280–1280.e1. doi:10.1016/j.cell.2020.02.029
Daëron, M., Jaeger, S., Du Pasquier, L., and Vivier, E. (2008). Immunoreceptor tyrosine-based inhibition motifs: A quest in the past and future. Immunol. Rev. 224, 11–43. doi:10.1111/j.1600-065X.2008.00666.x
Dasari, S., and Tchounwou, P. B. (2014). Cisplatin in cancer therapy: Molecular mechanisms of action. Eur. J. Pharmacol. 740, 364–378. doi:10.1016/j.ejphar.2014.07.025
Deng, F., Xu, Z., Zhou, J., Zhang, R., and Gong, X. (2021). ANLN regulated by miR-30a-5p mediates malignant progression of lung adenocarcinoma. Comput. Math. methods Med. 2021, 9549287. doi:10.1155/2021/9549287
Denisenko, T. V., Budkevich, I. N., and Zhivotovsky, B. (2018). Cell death-based treatment of lung adenocarcinoma. Cell. death Dis. 9 (2), 117. doi:10.1038/s41419-017-0063-y
Domagala-Kulawik, J., Osinska, I., and Hoser, G. (2014). Mechanisms of immune response regulation in lung cancer. Transl. lung cancer Res. 3 (1), 15–22. doi:10.3978/j.issn.2218-6751.2013.11.03
Duma, N., Santana-Davila, R., and Molina, J. R. (2019). Non-small cell lung cancer: Epidemiology, screening, diagnosis, and treatment. Mayo Clin. Proc. 94 (8), 1623–1640. doi:10.1016/j.mayocp.2019.01.013
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PloS one 9 (9), e107468. doi:10.1371/journal.pone.0107468
Guillerey, C. (2020). NK cells in the tumor microenvironment. Adv. Exp. Med. Biol. 1273, 69–90. doi:10.1007/978-3-030-49270-0_4
Han, L., Shi, H., Luo, Y., Sun, W., Li, S., Zhang, N., et al. (2020). Gene signature based on B cell predicts clinical outcome of radiotherapy and immunotherapy for patients with lung adenocarcinoma. Cancer Med. 9 (24), 9581–9594. doi:10.1002/cam4.3561
Hirsch, F. R., Scagliotti, G. V., Mulshine, J. L., Kwon, R., Curran, W. J., Wu, Y. L., et al. (2017). Lung cancer: Current therapies and new targeted treatments. Lancet (London, Engl. 389 (10066), 299–311. doi:10.1016/S0140-6736(16)30958-8
Hong, G., Chen, X., Sun, X., Zhou, M., Liu, B., Li, Z., et al. (2019). Effect of autologous NK cell immunotherapy on advanced lung adenocarcinoma with EGFR mutations. Precis. Clin. Med. 2 (4), 235–245. doi:10.1093/pcmedi/pbz023
Hoy, H., Lynch, T., and Beck, M. (2019). Surgical treatment of lung cancer. Crit. care Nurs. Clin. N. Am. 31 (3), 303–313. doi:10.1016/j.cnc.2019.05.002
Hsieh, C. S., Lee, H. M., and Lio, C. W. (2012). Selection of regulatory T cells in the thymus. Nat. Rev. Immunol. 12 (3), 157–167. doi:10.1038/nri3155
Hua, Z. D., Liu, X. B., Sheng, J. H., Li, C., Li, P., Cai, X. Q., et al. (2021). UBE2V2 positively correlates with PD-L1 expression and confers poor patient survival in lung adenocarcinoma. Appl. Immunohistochem. Mol. Morphol. AIMM 29 (8), 585–591. doi:10.1097/PAI.0000000000000928
Hutchinson, B. D., Shroff, G. S., Truong, M. T., and Ko, J. P. (2019). Spectrum of lung adenocarcinoma. Seminars ultrasound, CT, MR 40 (3), 255–264. doi:10.1053/j.sult.2018.11.009
Jiang, P., Gu, S., Pan, D., Fu, J., Sahu, A., Hu, X., et al. (2018). Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558. doi:10.1038/s41591-018-0136-1
Jiang, T., Fang, Z., Tang, S., Cheng, R., Li, Y., Ren, S., et al. (2021). Mutational landscape and evolutionary pattern of liver and brain metastasis in lung adenocarcinoma. J. Thorac. Oncol. official Publ. Int. Assoc. Study Lung Cancer 16 (2), 237–249. doi:10.1016/j.jtho.2020.10.128
Lee, N., Llano, M., Carretero, M., Ishitani, A., Navarro, F., López-Botet, M., et al. (1998). HLA-E is a major ligand for the natural killer inhibitory receptor CD94/NKG2A. Proc. Natl. Acad. Sci. U. S. A. 95 (9), 5199–5204. doi:10.1073/pnas.95.9.5199
Lin, M., Liang, S. Z., Shi, J., Niu, L. Z., Chen, J. B., Zhang, M. J., et al. (2017). Circulating tumor cell as a biomarker for evaluating allogenic NK cell immunotherapy on stage IV non-small cell lung cancer. Immunol. Lett. 191, 10–15. doi:10.1016/j.imlet.2017.09.004
Lin, M., Luo, H., Liang, S., Chen, J., Liu, A., Niu, L., et al. (2020). Pembrolizumab plus allogeneic NK cells in advanced non-small cell lung cancer patients. J. Clin. investigation 130 (5), 2560–2569. doi:10.1172/JCI132712
Nayyar, G., Chu, Y., and Cairo, M. S. (2019). Overcoming resistance to natural killer cell based immunotherapies for solid tumors. Front. Oncol. 9, 51. doi:10.3389/fonc.2019.00051
Newman, A. M., Liu, C. L., Green, M. R., Gentles, A. J., Feng, W., Xu, Y., et al. (2015). Robust enumeration of cell subsets from tissue expression profiles. Nat. methods 12 (5), 453–457. doi:10.1038/nmeth.3337
Nie, J., Shan, D., Li, S., Zhang, S., Zi, X., Xing, F., et al. (2021). A novel ferroptosis related gene signature for prognosis prediction in patients with colon cancer. Front. Oncol. 11, 654076. doi:10.3389/fonc.2021.654076
Paulsen, E. E., Kilvaer, T. K., Rakaee, M., Richardsen, E., Hald, S. M., Andersen, S., et al. (2017). CTLA-4 expression in the non-small cell lung cancer patient tumor microenvironment: Diverging prognostic impact in primary tumors and lymph node metastases. CII 66 (11), 1449–1461. doi:10.1007/s00262-017-2039-2
Russell, É., Conroy, M. J., and Barr, M. P. (2022). Harnessing natural killer cells in non-small cell lung cancer. Cells 11 (4), 605. doi:10.3390/cells11040605
Saab, S., Zalzale, H., Rahal, Z., Khalifeh, Y., Sinjab, A., and Kadara, H. (2020). Insights into lung cancer immune-based biology, prevention, and treatment. Front. Immunol. 11, 159. doi:10.3389/fimmu.2020.00159
Sadeghzadeh, M., Bornehdeli, S., Mohahammadrezakhani, H., Abolghasemi, M., Poursaei, E., Asadi, M., et al. (2020). Dendritic cell therapy in cancer treatment; the state-of-the-art. Life Sci. 254, 117580. doi:10.1016/j.lfs.2020.117580
Siegel, R. L., Miller, K. D., Fuchs, H. E., and Jemal, A. (2022). Cancer statistics, 2022. CA a cancer J. Clin. 72 (1), 7–33. doi:10.3322/caac.21708
Spella, M., and Stathopoulos, G. T. (2021). Immune resistance in lung adenocarcinoma. Cancers 13 (3), 384. doi:10.3390/cancers13030384
Succony, L., Rassl, D. M., Barker, A. P., McCaughan, F. M., and Rintoul, R. C. (2021). Adenocarcinoma spectrum lesions of the lung: Detection, pathology and treatment strategies. Cancer Treat. Rev. 99, 102237. doi:10.1016/j.ctrv.2021.102237
Sun, J., Yue, W., You, J., Wei, X., Huang, Y., Ling, Z., et al. (2021). Identification of a novel ferroptosis-related gene prognostic signature in bladder cancer. Front. Oncol. 11, 730716. doi:10.3389/fonc.2021.730716
Suster, D. I., and Mino-Kenudson, M. (2020). Molecular pathology of primary non-small cell lung cancer. Archives Med. Res. 51 (8), 784–798. doi:10.1016/j.arcmed.2020.08.004
Suzuki, C., Daigo, Y., Ishikawa, N., Kato, T., Hayama, S., Ito, T., et al. (2005). ANLN plays a critical role in human lung carcinogenesis through the activation of RHOA and by involvement in the phosphoinositide 3-kinase/AKT pathway. Cancer Res. 65 (24), 11314–11325. doi:10.1158/0008-5472.CAN-05-1507
Thorsson, V., Gibbs, D. L., Brown, S. D., Wolf, D., Bortone, D. S., Ou Yang, T. H., et al. (2018). The immune landscape of cancer. Immunity 48 (4), 812–830.e14. doi:10.1016/j.immuni.2018.03.023
Toro-Domínguez, D., Martorell-Marugán, J., López-Domínguez, R., García-Moreno, A., González-Rumayor, V., Alarcón-Riquelme, M. E., et al. (2019). ImaGEO: Integrative gene expression meta-analysis from GEO database. Bioinforma. Oxf. Engl. 35 (5), 880–882. doi:10.1093/bioinformatics/bty721
Valipour, B., Velaei, K., Abedelahi, A., Karimipour, M., Darabi, M., and Charoudeh, H. N. (2019). NK cells: An attractive candidate for cancer therapy. J. Cell. physiology 234 (11), 19352–19365. doi:10.1002/jcp.28657
Van Calster, B., Nieboer, D., Vergouwe, Y., De Cock, B., Pencina, M. J., and Steyerberg, E. W. (2016). A calibration hierarchy for risk models was defined: From utopia to empirical data. J. Clin. Epidemiol. 74, 167–176. doi:10.1016/j.jclinepi.2015.12.005
Van Calster, B., Wynants, L., Verbeek, J. F. M., Verbakel, J. Y., Christodoulou, E., Vickers, A. J., et al. (2018). Reporting and interpreting decision curve analysis: A guide for investigators. Eur. Urol. 74 (6), 796–804. doi:10.1016/j.eururo.2018.08.038
Vitale, I., Manic, G., Coussens, L. M., Kroemer, G., and Galluzzi, L. (2019). Macrophages and metabolism in the tumor microenvironment. Cell. metab. 30 (1), 36–50. doi:10.1016/j.cmet.2019.06.001
Wang, G., Li, X., Yao, Y., Jiang, Z., Zhou, H., Xie, K., et al. (2021). FAM83A and FAM83A-AS1 both play oncogenic roles in lung adenocarcinoma. Oncol. Lett. 21 (4), 297. doi:10.3892/ol.2021.12558
Xu, J., Zheng, H., Yuan, S., Zhou, B., Zhao, W., Pan, Y., et al. (2019). Overexpression of ANLN in lung adenocarcinoma is associated with metastasis. Thorac. cancer 10 (8), 1702–1709. doi:10.1111/1759-7714.13135
Yotsukura, M., Asamura, H., Motoi, N., Kashima, J., Yoshida, Y., Nakagawa, K., et al. (2021). Long-term prognosis of patients with resected adenocarcinoma in situ and minimally invasive adenocarcinoma of the lung. J. Thorac. Oncol. official Publ. Int. Assoc. Study Lung Cancer 16 (8), 1312–1320. doi:10.1016/j.jtho.2021.04.007
Yu, Y., Wang, Z., Zheng, Q., and Li, J. (2021). GREB1L overexpression correlates with prognosis and immune cell infiltration in lung adenocarcinoma. Sci. Rep. 11 (1), 13281. doi:10.1038/s41598-021-92695-x
Zhang, C., Zhang, G., Sun, N., Zhang, Z., Zhang, Z., Luo, Y., et al. (2020). Comprehensive molecular analyses of a TNF family-based signature with regard to prognosis, immune features, and biomarkers for immunotherapy in lung adenocarcinoma. EBioMedicine 59, 102959. doi:10.1016/j.ebiom.2020.102959
Zhang, H., Luo, Y. B., Wu, W., Zhang, L., Wang, Z., Dai, Z., et al. (2021). The molecular feature of macrophages in tumor immune microenvironment of glioma patients. Comput. Struct. Biotechnol. J. 19, 4603–4618. doi:10.1016/j.csbj.2021.08.019
Zhang, Q., Tang, L., Zhou, Y., He, W., and Li, W. (2021). Immune checkpoint inhibitor-associated pneumonitis in non-small cell lung cancer: Current understanding in characteristics, diagnosis, and management. Front. Immunol. 12, 663986. doi:10.3389/fimmu.2021.663986
Zhang, D., Jiang, Q., Ge, X., Shi, Y., Ye, T., Mi, Y., et al. (2021). RHOV promotes lung adenocarcinoma cell growth and metastasis through JNK/c-Jun pathway. Int. J. Biol. Sci. 17 (10), 2622–2632. doi:10.7150/ijbs.59939
Zhang, Z. (2016). Variable selection with stepwise and best subset approaches. Ann. Transl. Med. 4 (7), 136. doi:10.21037/atm.2016.03.35
Zhong, R., Chen, D., Cao, S., Li, J., Han, B., and Zhong, H. (2021). Immune cell infiltration features and related marker genes in lung cancer based on single-cell RNA-seq. Clin. Transl. Oncol. official Publ. Fed. Span. Oncol. Soc. Natl. Cancer Inst. Mexico 23 (2), 405–417. doi:10.1007/s12094-020-02435-2
Keywords: lung adenocarcinoma, NK cells, programmed cell death, immunity, survival, prediction model, chemotherapy drug
Citation: Zhang D and Zhao Y (2023) Identification of natural killer cell associated subtyping and gene signature to predict prognosis and drug sensitivity of lung adenocarcinoma. Front. Genet. 14:1156230. doi: 10.3389/fgene.2023.1156230
Received: 01 February 2023; Accepted: 20 March 2023;
Published: 07 April 2023.
Edited by:
Min Sun, Hubei University of Medicine, ChinaCopyright © 2023 Zhang and Zhao. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Dexin Zhang, dexin1994@mail.xjtu.edu.cn