Identification of prognostic-related tumor microenvironment genes in lung adenocarcinoma and establishment of a prognostic prediction model

Introduction

Lung cancer is a malignant tumor with high incidence and mortality worldwide. According to the Global Cancer Observatory (GLOBOCAN) 2020 data, lung cancer ranks first in both incidence and mortality in China, with a 5-year survival rate of less than 20%, posing a major public health challenge (1,2). Non-small cell lung cancer (NSCLC), the most prevalent form of lung cancer, accounts for 85% of all cases. Among NSCLC subtypes, lung adenocarcinoma (LUAD) has the highest incidence, comprising 35% of cases in the Chinese population (3). In recent years, the introduction of immune checkpoint inhibitors (ICIs) has significantly improved the outcomes of NSCLC patients (4,5), highlighting the importance of studying tumor immune characteristics. Given its prevalence, LUAD has garnered considerable attention regarding its immune features and prognosis.

The tumor microenvironment (TME) is a dynamic system consisting of immune cells, stromal cells, cancer cells, and local secreted factors. Immune cells and interstitial cells are the major non-tumor components, which are considered to play an essential role in the development and prognosis of tumors (6,7). Previous study has demonstrated that stromal cells in the TME are genetically stable, making them promising therapeutic targets for minimizing the risk of drug resistance and tumor recurrence (8). Additionally, the metabolic dysregulation of tumor cells often induces metabolic stress in tumor-infiltrating immune cells, leading to impaired anti-tumor immune responses and abnormalities in immune-related genomic functions (9,10). A multi-target approach inhibiting TME-related genomes may provide a more effective cancer treatment (11). Since the TME is closely related to immune characteristics and prognosis, it is also essential to understand the correlation between the TME features and patient outcomes.

Predictive models can better assist in improving the accuracy of clinical diagnosis and treatment. Currently, various models have been developed in lung cancer to predict patient prognosis, including radiological features, clinical characteristics, tumor markers, and genomic features (12-14). While considerable research has also been conducted on immunological features and related models, most are based on cells, proteins, or imaging characteristics (15,16). These features often require advanced detection methods and analytical techniques, limiting their practicality in clinical applications. Consequently, exploring TME-related gene transcriptomic expression may provide a more applicable approach.

This study utilized bioinformatics approaches to analyze and assess the immune characteristics of LUAD patients. Key TME-related genes associated with prognosis were identified, and a prognostic prediction model was developed to aid in clinical diagnosis and treatment. Gaining deeper insights into the immune features of LUAD and the expression of associated regulatory genes is crucial for advancing the development of clinical prognostic models and guiding treatment selection strategies in the future. We present this article in accordance with the TRIPOD reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-297/rc).

Methods

Data download

The transcriptomic data (RNA-seq) and clinical information data of LUAD patients were downloaded through The Cancer Genome Atlas (TCGA) database (https://xenabrowser.net/datapages/); the information of a total of 559 samples was collected in this way, which included 58 control samples and 501 LUAD tumor samples. Among the 501 LUAD samples, males and females accounted for 46.1% and 53.9%, respectively. The median age was 66 years, and 74.8% of the patients were in stages I–II.

The Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) was used to download the GSE50081 (17) data set as a validation queue. A total of 181 samples were included in the GSE50081 cohort, among which 127 LUAD tumor samples with clinical survival and prognosis information were available for external verification. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Assessment of immune status and screening of the differentially expressed genes (DEGs) associated with TME

The transcriptomic data [Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data (ESTIMATE) algorithm] was employed to assess the immune status of LUAD patients and estimate the stromal cells and immune cells in tumors. An ESTIMATE model score was subsequently assigned, comprising an Estimate score, immune score, stromal score, and tumor purity rating. It was used to assess the degree of immune cell invasion and the immune status (18). According to the median value of the above 4 scoring dimensions, the patients were divided into two groups, and the differences in survival prognosis among sample groups with different scoring levels were explored. Further, the patients were divided into high versus low stromal score and high versus low immune score groups according to the immune score and stromal score results, respectively. Limma algorithm was used to identify the significant DEGs between the two groups. The genes with false discovery rate (FDR) <0.05 and |log₂ fold change (FC)| >0.5 were considered TME-related DEGs (TME-DEGs). Moreover, hierarchical clustering was performed on the selected DEGs, and the results were analyzed by Database for Annotation, Visualization and Integrated Discovery (DAVID; version 6.8; https://david.ncifcrf.gov/) under the 2-grade intersection of DEGs of biological annotations [including Gene Ontology (GO) function and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment] signaling pathway (19,20).

Exploration of the relationship of the TME-DEGs related to prognosis

The TME-DEGs significantly associated with survival prognosis were screened by univariate Cox regression analysis. Further, multivariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses were conducted to screen the essential genes of TME-DEGs that were independently associated with prognosis (ITME genes). Finally, the cases were divided into high and low groups according to their expression levels of ITME genes, and the difference in prognosis between the two groups was evaluated by the Kaplan-Meier curve method.

Verification and analysis of the correlation of ITME genes with immunological and clinical characteristics

The correlation between the ESTIMATE model scores and the expression level of the ITME genes in the TCGA LUAD training data set was analyzed.

Based on the gene expression levels of TCGA LUAD training data set, the online tool Tumor Immune Estimation Resource (TIMER; https://cistrome.shinyapps.io/timer/) (21) was used to calculate the association of ITME genes and infiltration level of 6 types immune cells (include B cells, CD4 + T cells, CD8 + T cells, neutrophils, macrophages, and dendritic cells).

The patients were divided into two groups according to the expression level of each ITME gene obtained by the preceding steps, and the differences in clinical information between the groups with different expression levels were compared.

Establishment and validation of a prognostic model based on ITME genes

Based on the LASSO coefficient and gene expression level of ITME genes, a risk score (RS) prediction model was established. The RS calculation formula is as follows:

RS=∑Coefgenes×Expgenes[1]

The Coef_genes represents the LASSO coefficient of the target gene, and Exp_genes represents the expression level of the target gene.

The RS values were calculated in the TCGA training data set, and then the patients were divided into high (RS score ≥ median RS value) and low (RS score < median RS value) groups. Kaplan-Meier curve analysis evaluated the association between the two groups and actual prognostic information. The correlation between the RS model, each ITME gene, and prognosis was verified by the GSE50081 data set.

Statistical analysis

All statistical analyses were performed using R 3.6.1 (https://www.r-project.org/). The “ESTIMATE package” is used to estimate immunization scores, the “limma package” is used to screening the TME-DEGs, and the “Cor package” used to verify the correlation between ITME genes and immunization scores. “Survival package” and “lars package” were used to conduct LASSO Cox regression analysis, and establish the prediction model. The ‘ROC package” used to plot receiver operating characteristic (ROC) curve for the RSs in the training set, and the area under the curve (AUC) value was calculated to evaluate the performance of the prognostic model. T-test and Fisher’s test were used for statistics, and P<0.05 was considered statistically significant. The detailed analysis process is shown in Figure 1.

Figure 1 Schematic diagram of the analysis process. DAVID, Database for Annotation, Visualization and Integrated Discovery; DFS, disease-free survival; DEGs, differentially expressed genes; ESTIMATE, an Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; KEGG, Kyoto Encyclopedia of Genes and Genomes; LASSO, least absolute shrinkage and selection operator; LUAD, lung adenocarcinoma; OS, overall survival; RFS, recurrence-free survival; TCGA, The Cancer Genome Atlas.

Results

Calculation of the ESTIMATE model scores

Five hundred and fifty-nine samples were collected from the TCGA training data set, including 58 control samples of normal individuals and 501 LUAD tumor samples of LUAD patients. The genes’ expression levels were used to calculate each patient’s Estimate score, immune score, stromal score, and tumor purity, and the differences between LUAD samples and control samples were compared (available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-1.xlsx). The results showed that the Estimate score, immune score, and stromal score in LUAD tumor samples were significantly lower than those in control samples, and tumor purity was significantly higher than in control samples (Figure 2).

Figure 2 The distribution of stromal score, immune score, ESTIMATE score, and tumor purity in tissue samples of normal control and LUAD tumor. ***, P<0.001. ESTIMATE, an Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; LUAD, lung adenocarcinoma.

To compare the differences in survival with different ESTIMATE model scores levels, the samples were divided into low- and high-level groups respectively according to the median level of the 4 scores, and the differences in overall survival (OS), relapse-free survival (RFS), and disease-free survival (DFS) were analyzed. The results showed that the OS, RFS, and DFS of the groups with high scores were all significantly better than those of the low-scoring groups (Figure 3), indicating a good prognosis for high scores patients (the corresponding clinical prognostic information and other data of the samples are shown in https://cdn.amegroups.cn/static/public/tlcr-24-297-2.xlsx).

Figure 3 Kaplan-Meier survival curves of relationship between different levels of stromal score, immune score, ESTIMATE score, and tumor purity groups and prognosis of OS, RFS, and DFS, in which the blue and red curves indicating low and high rating groups respectively. CI, confidence interval; DFS, disease-free survival; ESTIMATE, an Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; OS, overall survival; RFS, recurrence-free survival.

Screening of TME-DEGs

According to the results of the immune score and stromal score, the samples were divided into two groups (high and low stromal score group, high and low immune score group); 1,201 and 987 DEGs conformed to the threshold conditions respectively, and their expression levels were analyzed (Figure 4A, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-3.xlsx). Moreover, the DEGs sets filtered from stromal and immune score were compared and filtered further to obtain their intersections set. 66 DEGs with significantly down-regulated expression and 614 with significantly up-regulated expression were obtained (Figure 4B, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-4.xlsx).

Figure 4 DEGs related with TME. (A) Volcano map (left) and heat map (right) of DEGs in LUAD sample. In volcano map, red and blue dots indicate DEGs with significantly up-regulated and down-regulated expression respectively, and gray dots indicate DEGs expression with no significant difference; the dashed horizontal lines indicate FDR <0.05, the 2 vertical dashed lines indicate |log₂FC| >0.5. In the heat map, the horizontal axis shows the sequence of samples from low to high scores. Above is the analysis based on stromal score, and below is the analysis based on immune score. (B) Venn diagram of intersection filter by compare the DEGs sets under the 2 scores. DEGs, differentially expressed genes; FC, fold change; FDR, false discovery rate; LUAD, lung adenocarcinoma; TME, tumor microenvironment.

Based on DAVID, 680 DEGs were subjected to GO function enrichment and KEGG signaling pathway analysis. Screening was conducted of 96 biological processes (BP), 11 cellular components (CC), 24 molecular functions (MF), and 17 KEGG signaling pathways (available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-5.xlsx). In the GO function annotation, immune response of BP, external side of plasma membrane of CC, and transmembrane signaling receptor activity of MF had the most significant correlation with TME (Figure 5A). KEGG signaling pathway analysis showed that the most significant correlation was hsa04060: cytokine-cytokine receptor interaction (Figure 5B, Figure S1).

Figure 5 DEGs annotated map of relationship with biological function. (A) GO function annotated map. The horizontal axis represents the number of genes, the vertical axis represents the function item, the columns with different colors represent different classifications, and the fonts marked in red represent the item with the most significant correlation in the classification; (B) KEGG signal path enrichment analysis annotated map. The horizontal axis represents the number of DEGs participating entries, the vertical axis represents signal pathway, and the color represents degree of correlation, with more red indicating that it was the more significant one. BP, biological process; CC, cellular component; DEGs, differentially expressed genes; FDR, false discovery rate; GO, Gene Ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; MF, molecular function.

Screening of ITME genes

A total of 160 DEGs with significant survival prognosis were screened from 680 TME-DEGs by the univariate Cox regression analysis (available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-6.xlsx). Multivariate Cox regression analysis was then performed on the 160 DEGs, and 26 DEGs independently and significantly correlated with prognosis were screened (available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-7.xlsx). Finally, LASSO survival regression analysis was performed on the 26 DEGs, and 5 ITME genes were obtained (Figure 6A), all independently and significantly correlated with prognosis (Figure 6B). These 5 DEGs were ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A, important TME-DEGs associated with prognosis.

Figure 6 Correlation analysis of TME-DEGs and prognosis. (A) LASSO parameter diagram, the intersection of horizontal and vertical red dotted lines (−4.013, 0.226) in the following diagram should be the optimal parameter position. (B) Forest map of multivariate Cox regression analysis. *, P<0.05 **, P<0.01. CI, confidence interval; DEGs, differentially expressed genes; LASSO, least absolute shrinkage and selection operator; TME, tumor microenvironment.

Further, according to the expression level of each ITME gene, the patients were divided into high and low-expression groups, which were then evaluated regarding the prognostic differences in OS, RFS, and DFS. The results showed that the 5 genes were all significantly correlated with OS, RFS, and DFS, of which the samples with low ABCC2 expression level had better clinical prognosis (P<0.001, P=0.042, P=0.02). High expression levels of ECT2L (P<0.001, P=0.002, P=0.02), CD200R1 (P<0.001, P=0.02, P=0.01), ACSM5 (P<0.001, P=0.02, P=0.049), and CLEC17A (P<0.001, P=0.046, P=0.03) had better clinical prognosis (Figure 7, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-8.xlsx).

Figure 7 Kaplan-Meier survival curves of relationship between prognosis and ITME genes in TCGA training data set, in which the blue and red curves represent low and high expression sample groups respectively. The prognosis included OS, RFS, and DFS, the ITME genes included ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A. CI, confidence interval; DFS, disease-free survival; DEGs, differentially expressed genes; ITME, TME-DEGs that were independently associated with prognosis; OS, overall survival; RFS, relapse-free survival; HR, hazard ratio; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.

The correlation between ITME genes and ESTIMATE model scores

To verify the correlation between the expression levels of 5 ITME genes and ESTIMATE model scores, we analyzed the correlation between ITME genes and immune characteristics. The results showed that ECT2L, CD200R1, ACSM5, and CLEC17A have a significant positive correlation with stromal, immune, and Estimate scores, and a significant negative correlation with tumor purity scores, but the ABCC2 has an opposite significant correlation with the scores (Figure 8, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-9.xlsx), which is consistent with prognostic-related features.

Figure 8 Scatter plot of correlation between ITME genes and ESTIMATE algorithm scores. The ITME genes included ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A, the ESTIMATE algorithm scores include stromal score, immune score, estimate score, and tumor purity. DEGs, differentially expressed genes; ITME, TME-DEGs that were independently associated with prognosis; ESTIMATE, Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data; TME, tumor microenvironment.

The correlation between the ITME genes and the infiltration levels of immune cells

The TIMER resource was used to calculate the infiltration levels of immune cells of the TCGA training data set, and the correlation between the infiltration level of each type immune cell and the ITME genes was analyzed. The results showed that ECT2L, CD200R1, ACSM5, and CLEC17A were significantly positively correlated with all 6 types of immune cells, including B cells, CD8+ T cells, CD4+ T cells, macrophages, neutrophils, and dendritic cells (DC). However, ABCC2 was not significantly correlated with these immune cells, and was negatively associated with B cells, macrophages, and DC (Figure 9).

Figure 9 Scatter plot of correlation between ITME genes and immune cells level. The ITME genes included ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A, the immune cells included B cells, CD8⁺ T cells, CD4⁺ T cells, macrophages, neutrophils, and DC. DC, dendritic cell; DEGs, differentially expressed genes; ITME, TME-DEGs that were independently associated with prognosis; LUAD, lung adenocarcinoma; TME, tumor microenvironment; TPM, transcripts per million.

The correlation between the ITME genes and the clinical characteristics

In the TCGA training data set, the samples were divided into high and low groups according to the expression level of each ITME gene. The clinical information of samples in different groups was statistically analyzed. The results showed that the high expression level of ECT2L, CD200R1, ACSM5, and CLEC17A genes was significantly related to the lower T stage of LUAD patients (P<0.001). However, there was no significant correlation between ABCC2 gene expression level and clinical characteristics (Table 1, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-10.xlsx).

Establishment of a prognostic model with ITME genes

Based on the 5 ITME genes, the following RS model was established:

RS=(0.024026758)∗ExpABCC2+(−0.00052061)∗ExpECT2L+(−0.019691114)∗ExpCLEC17A+(−0.012622394)∗ExpCD200R1+(−0.004677594)∗ExpACSM5[2]

The RS values of the samples were calculated from the TCGA training dataset and divided into high and low groups based on the median value. The results of prognostic correlation analysis showed that the OS (P<0.001), RFS (P=0.009), and DFS (P=0.005) of patients in the low RS group were significantly better than those in the high RS group (Figure 10, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-11.xlsx). The ROC curve was used to analyze the accuracy of the prediction model, and the results showed that it had good prediction efficiency for 5-year OS (AUC =0.70), RFS (AUC =0.59), and DFS (AUC =0.56) (Figure 11, Figure S2).

Figure 10 Kaplan-Meier survival curves of relationship between prognosis and different RS groups based on TCGA training data set, in which the blue and red curves represent low and high RS groups respectively. The prognosis included (A) OS; (B) RFS and (C) DFS. CI, confidence interval; DFS, disease-free survival; OS, overall survival; RFS, relapse-free survival; HR, hazard ratio; RS, risk score; TCGA, The Cancer Genome Atlas.

Figure 11 ROC curve of the RS prediction model of OS based on TCGA training data set, in which the red, blue and green curves represent 1, 3, and 5 years ROC respectively. AUC, area under the curve; CI, confidence interval; OS, overall survival; ROC, receiver operating characteristic; RS, risk score; TCGA, The Cancer Genome Atlas.

Validation of a prognostic model with ITME genes

The GSE50081 cohort containing 127 LUAD samples was obtained from GEO, and its RS values were analyzed to verify the performance of the prognostic model. The results showed that the OS (P<0.001) and DFS (P=0.006) of patients in the low RS group were still significantly better than those in the high group in the GSE50081 validation data set (Figure 12, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-12.xlsx). The ROC curve invalidation data set also showed a good prediction efficiency for 5-year OS (AUC =0.72) and DFS (AUC =0.67) (Figure 13, Figure S2). ABCC2 was significantly expressed in the high RS group, whereas the other 4 genes were significantly expressed in the low RS group; the same trends were observed in the TCGA training data set (Figure 14, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-13.xlsx). Moreover, the correlation between each ITME gene and prognosis was verified in GSE50081 validation data set. The results showed that the high expression of ECT2L, CD200R1, ACSM5, and CLEC17A was significantly associated with a better prognosis, whereas high expression of ABCC2 was associated with a worse prognosis (Figure 15, available online: https://cdn.amegroups.cn/static/public/tlcr-24-297-14.xlsx). These results are also consistent with those of the TCGA training data set. Taken together, these results suggested that the prognostic model may have a good predictive ability.

Figure 12 Kaplan-Meier survival curves of relationship between prognosis and different RS groups based on GSE50081 verification data set, in which the blue and red curves represent low and high RS groups respectively. The prognosis included (A) OS and (B) DFS. CI, confidence interval; DFS, disease-free survival; OS, overall survival; HR, hazard ratio; RS, risk score.

Figure 13 ROC curve of the RS prediction model of OS based on GSE50081 verification data set, in which the red, blue and green curves represent 1, 3, and 5 years ROC respectively. AUC, area under the curve; CI, confidence interval; OS, overall survival; ROC, receiver operating characteristic; RS, risk score.

Figure 14 Expression levels of the 5 ITME genes in different RS groups of (A) TCGA training data set and (B) GSE50081 verification data set. *, P<0.05; ***, P<0.001. DEGs, differentially expressed genes; ITME, TME-DEGs that were independently associated with prognosis; RS, risk score; TCGA, The Cancer Genome Atlas; TME, tumor microenvironment.

Figure 15 Kaplan-Meier survival curves of relationship between prognosis and ITME genes in GSE50081 verification data set, in which the blue and red curves represent low and high expression sample groups respectively. The prognosis included OS and DFS, the ITME genes included ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A. CI, confidence interval; DEGs, differentially expressed genes; DFS, disease-free survival; ITME, TME-DEGs that were independently associated with prognosis; OS, overall survival; HR, hazard ratio; TME, tumor microenvironment.

Discussion

With the rapid development of medical precision, previous study has found that tumors’ occurrence and development are related to the internal genetic factors of cancer cells and the interaction of various systems in the body (22). The TME is a dynamic system composed of immune cells, stromal cells, cancer cells, and a network of complex cytokines and chemokines (23). Among them, immune-related cells and factors participate in tumorigenesis, proliferation, and development, which is particularly important (24). Therefore, exploring the characteristics of immune-related molecules and evaluating the correlation between TME-related genes and prognosis is significant. In this study, gene expression data of 501 LUAD patients in the TCGA database were analyzed, 680 DEGs significantly related to TME were obtained, and subjected to LASSO and Cox multivariate analysis. Finally, 5 TME genes independently associated with prognosis were screened out, among which ECT2L, CD200R1, ACSM5, and CLEC17A were positively related to pre-treatment and post-treatment, and ABCC2 was negatively associated with prognosis. The prognosis risk prediction model constructed by the above 5 genes showed that the prognosis of high- and low-risk patients was significantly different in the training set and verification set, which may be an effective indicator to predict the prognosis of LUAD patients. In addition, the correlation between the expression of the 5 TME genes and the level of immune cell infiltration and clinical features was analyzed. Since these genes are significantly related to immunity and prognosis, they may be of great significance in diagnosing and treating LUAD patients, which is worth further exploration.

Based on the mechanism correlation between TME and immunity, significant efforts have been made in the past 5 to 10 years to explore the characteristics of TME associated with anti-tumor immune efficacy and prognosis prediction of lung cancer (24). At present, the “immune score” (IS) based on the number of T lymphocytes in tumors has been included in the National Comprehensive Cancer Network (NCCN) guidelines of Colon Cancer (25-27), and the TME analysis characterized by the presence of tumor-infiltrating lymphocytes (TILs) has also been used to predict the treatment and prognosis of patients (28). Although these scores may be robust biomarkers in clinical practice, it is difficult to detect the infiltration of immune cells in patients with limited tissue accessibility, especially in advanced-stage patients. Gene expression analysis can replace cell and protein characteristics to estimate TME at the transcriptomic level, which has caused a tremendous upsurge in the scientific community (29). With the emergence of high-throughput sequencing technology [next-generation sequencing (NGS)], there have been multiple studies exploring TME and prognosis through RNA-seq construction algorithms (28,30,31). The ESTIMATE model is an algorithm used to estimate the immune score and tumor purity by inferring the level of invasion matrix and immune cells in tumor tissue from gene expression data. The estimation ability of this method has been verified in large independent data sets, and many solid tumor-related studies have shown a correlation between the ESTIMATE model score and prognosis (18,32-34). In a retrospective exploration of the ORIENT-11 data, it was also emphasized that the ESTIMATE model is a robust predictive factor for immune efficacy (35). However, the ESTIMATE model score will still be interfered by partial gene expression due to its dependence on infiltrating stromal cells and immune cells (18). Therefore, we used the data of LUAD patients in the TCGA training data set to optimize the screening of TME-related genes. Finally, 5 TME genes (ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A) independently and significantly related to the prognosis were screened as ITME genes, were used to construct a prognosis risk prediction model. Moreover, the 5 ITME genes were significantly correlated with the ESTIMATE model score. The RS prediction model built in this study has a more significant difference in prognosis between high and low groups of LUAD patients to the ESTIMATE scoring model, especially in terms of OS (P value: 0.002 vs. <0.001), which may have a better prediction efficiency.

Aside from the prediction model, a single abnormal gene expression is also related to the prognosis of tumor patients. It has been shown that the expression of 15 characteristic genes, such as ATP1B1, TRIM14, and FAM64A, is significantly related to the prognosis of NSCLC patients, and may be used as a prognostic marker of NSCLC patients (36). Similarly, RRAGB, RSPH9, RPS6KL1, RXFP1, RRM2, and RTL1 are associated with 10-year survival in LUAD patients (37). Our research results show that ABCC2, ECT2L, CD200R1, ACSM5, and CLEC17A are significantly correlated with OS, DFS, and RFS of LUAD patients, among which ECT2L, CD200R1, ACSM5, and CLEC17A are positively correlated. ACSM5 is a gene encoding acyl coenzyme A synthetase protein. Previous studies have shown that the high expression of ACSM5 is significantly related to a better prognosis of thyroid and breast cancer (38,39), consistent with our findings in LUAD. CD200R1 is a gene encoding the receptor of OX-2 membrane glycoprotein, expressed in various immune cells (40). The CD200/CD200R1 signaling pathway inhibits anti-tumor response by regulating the function of macrophages and T cells, and its high expression may be related to better immune level (41,42). However, in a study focusing on patients with early NSCLC, the expression of CD200R1 protein was associated with poor immune microenvironment and poor prognosis (43), was contrary to our findings. However, whether the expression of CD200R1 protein has the same prognostic impact on patients with advanced NSCLC has not been reported, and further exploration is needed. CLEC17A encodes a unique glycan- binding receptor, often expressed on the surface of B cells (44). Previous study has shown that it is associated with the migration of epithelial tumor cells (45), but the correlation between its high gene expression and prognosis has not been reported.

Similarly, the protein encoded by ECT2L, as a guanine exchange factor, has only been reported to be expressed in patients with early acute lymphoblastic leukemia (46). The correlation between this gene’s high expression and solid tumors’ prognosis has not been reported. ABCC2 is a gene encoding a protein that is a member of the adenosine triphosphate (ATP) binding cassette (ABC) transporter superfamily, and is the only gene negatively related to prognosis. The protein encoded by this gene is mainly involved in transporting body substrates (drugs, poisons, and endogenous compounds) (47). Previous studies have shown that the high expression of ABCC2 is associated with drug resistance to chemotherapy (47,48). In a hepatobiliary function study, high macrophage infiltration will lead to NF-KB signal transduction, interfere with the binding of the ABCC2 gene and promoter, and cause bile transport disorder (49). Its correlation with chemotherapy resistance and the antagonism of macrophages may be the reasons for its negative correlation with immunity and prognosis.

In addition, this study also analyzed the correlation between the expression of 5 genes and the level of immune cell infiltration. The results showed that ECT2L, CD200R1, ACSM5, and CLEC17A are significantly positively correlated with 6 kinds of immune cells, whereas ABCC2 was not significantly correlated with them, even exhibiting a negative trend. It has been shown that the expression of ABCC2 gene is associated with poor immune status in LUAD patients (50), and several previous studies have shown that the infiltration level of T cells, B cells, and macrophages is positively correlated with the prognosis of LUAD patients (51-53). These results further verified the correlation between the above 5 genes and the immunity and prognosis of LUAD patients.

This paper also has the following shortcomings: (I) It only conducted analysis and model establishment verification based on the public cohort, with no further clinical verification; (II) the mechanism of the 5 genes screened related to immunity and tumor prognosis is still limited, and our study did not further explore the mechanism.

Conclusions

This study screened 5 independent prediction genes related to LUAD immunity and prognosis through a TCGA data set, and constructed a prognostic prediction model. The model has good prediction efficiency and is expected to provide help for clinical diagnosis and therapy.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-297/rc

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-297/prf

Funding: This research was funded by the Guangdong Medical Science and Research Foundation (No. B2021023) and Guangzhou Planned Project of Science and Technology (No. 202102080609).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-24-297/coif). L.B. received speakers’ fees from Astra-Zeneca, MSD, Roche, Takeda and travel fees from Takeda and Sanofi, outside the current work. The other authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin 2021;71:209-49. [Crossref] [PubMed]
2. Li W, Wu M, Wang Q, et al. A comparative genomics analysis of lung adenocarcinoma for Chinese population by using panel of recurrent mutations. J Biomed Res 2020;35:11-20. [Crossref] [PubMed]
3. Li D, Shi J, Dong X, et al. Epidemiological characteristics and risk factors of lung adenocarcinoma: A retrospective observational study from North China. Front Oncol 2022;12:892571. [Crossref] [PubMed]
4. Pakkala S, Owonikoko TK. Immune checkpoint inhibitors in small cell lung cancer. J Thorac Dis 2018;10:S460-7. [Crossref] [PubMed]
5. Xiong A, Wang J, Zhou C. Immunotherapy in the First-Line Treatment of NSCLC: Current Status and Future Directions in China. Front Oncol 2021;11:757993. [Crossref] [PubMed]
6. Jia D, Li S, Li D, et al. Mining TCGA database for genes of prognostic value in glioblastoma microenvironment. Aging (Albany NY) 2018;10:592-605. [Crossref] [PubMed]
7. Belluomini L, Dodi A, Caldart A, et al. A narrative review on tumor microenvironment in oligometastatic and oligoprogressive non-small cell lung cancer: a lot remains to be done. Transl Lung Cancer Res 2021;10:3369-84. [Crossref] [PubMed]
8. Quail DF, Joyce JA. Microenvironmental regulation of tumor progression and metastasis. Nat Med 2013;19:1423-37. [Crossref] [PubMed]
9. Zhong Y, She Y, Deng J, et al. Deep Learning for Prediction of N2 Metastasis and Survival for Clinical Stage I Non-Small Cell Lung Cancer. Radiology 2022;302:200-11. [Crossref] [PubMed]
10. Li X, Wenes M, Romero P, et al. Navigating metabolic pathways to enhance antitumour immunity and immunotherapy. Nat Rev Clin Oncol 2019;16:425-41. [Crossref] [PubMed]
11. Roma-Rodrigues C, Mendes R, Baptista PV, et al. Targeting Tumor Microenvironment for Cancer Therapy. Int J Mol Sci 2019;20:840. [Crossref] [PubMed]
12. Ren H, Xiao Z, Ling C, et al. Development of a novel nomogram-based model incorporating 3D radiomic signatures and lung CT radiological features for differentiating invasive adenocarcinoma from adenocarcinoma in situ and minimally invasive adenocarcinoma. Quant Imaging Med Surg 2023;13:237-48. [Crossref] [PubMed]
13. Zhang G, Wang X, Jia J, et al. Development and validation of a nomogram for predicting survival in patients with surgically resected lung invasive mucinous adenocarcinoma. Transl Lung Cancer Res 2021;10:4445-58. [Crossref] [PubMed]
14. Koyama J, Morise M, Furukawa T, et al. Artificial intelligence-based personalized survival prediction using clinical and radiomics features in patients with advanced non-small cell lung cancer. BMC Cancer 2024;24:1417. [Crossref] [PubMed]
15. Tong H, Sun J, Fang J, et al. A Machine Learning Model Based on PET/CT Radiomics and Clinical Characteristics Predicts Tumor Immune Profiles in Non-Small Cell Lung Cancer: A Retrospective Multicohort Study. Front Immunol 2022;13:859323. [Crossref] [PubMed]
16. Zhang J, Liu X, Huang Z, et al. T cell-related prognostic risk model and tumor immune environment modulation in lung adenocarcinoma based on single-cell and bulk RNA sequencing. Comput Biol Med 2023;152:106460. [Crossref] [PubMed]
17. Der SD, Sykes J, Pintilie M, et al. Validation of a histology-independent prognostic gene signature for early-stage, non-small-cell lung cancer including stage IA patients. J Thorac Oncol 2014;9:59-64. [Crossref] [PubMed]
18. Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun 2013;4:2612. [Crossref] [PubMed]
19. Huang da W. Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc 2009;4:44-57. [Crossref] [PubMed]
20. Huang da W. Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res 2009;37:1-13. [Crossref] [PubMed]
21. Li T, Fan J, Wang B, et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res 2017;77:e108-10. [Crossref] [PubMed]
22. Ladbury CJ, Rusthoven CG, Camidge DR, et al. Impact of Radiation Dose to the Host Immune System on Tumor Control and Survival for Stage III Non-Small Cell Lung Cancer Treated with Definitive Radiation Therapy. Int J Radiat Oncol Biol Phys 2019;105:346-55. [Crossref] [PubMed]
23. Graves EE, Maity A, Le QT. The tumor microenvironment in non-small-cell lung cancer. Semin Radiat Oncol 2010;20:156-63. [Crossref] [PubMed]
24. Giatromanolaki A, Koukourakis IM, Balaska K, et al. Programmed death-1 receptor (PD-1) and PD-ligand-1 (PD-L1) expression in non-small cell lung cancer and the immune-suppressive effect of anaerobic glycolysis. Med Oncol 2019;36:76. [Crossref] [PubMed]
25. Fridman WH, Pagès F, Sautès-Fridman C, et al. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer 2012;12:298-306. [Crossref] [PubMed]
26. Pagès F, Mlecnik B, Marliot F, et al. International validation of the consensus Immunoscore for the classification of colon cancer: a prognostic and accuracy study. Lancet 2018;391:2128-39. [Crossref] [PubMed]
27. Galon J, Costes A, Sanchez-Cabo F, et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science 2006;313:1960-4. [Crossref] [PubMed]
28. Brambilla E, Le Teuff G, Marguet S, et al. Prognostic Effect of Tumor Lymphocytic Infiltration in Resectable Non-Small-Cell Lung Cancer. J Clin Oncol 2016;34:1223-30. [Crossref] [PubMed]
29. Genova C, Dellepiane C, Carrega P, et al. Therapeutic Implications of Tumor Microenvironment in Lung Cancer: Focus on Immune Checkpoint Blockade. Front Immunol 2021;12:799455. [Crossref] [PubMed]
30. Li B, Severson E, Pignon JC, et al. Comprehensive analyses of tumor immunity: implications for cancer immunotherapy. Genome Biol 2016;17:174. [Crossref] [PubMed]
31. Sturm G, Finotello F, Petitprez F, et al. Comprehensive evaluation of transcriptome-based cell-type quantification methods for immuno-oncology. Bioinformatics 2019;35:i436-45. [Crossref] [PubMed]
32. Ji Y, Zhang T, Yang L, et al. The effectiveness of three-dimensional reconstruction in the localization of multiple nodules in lung specimens: a prospective cohort study. Transl Lung Cancer Res 2021;10:1474-83. [Crossref] [PubMed]
33. Yu X, Chen YA, Conejo-Garcia JR, et al. Estimation of immune cell content in tumor using single-cell RNA-seq reference data. BMC Cancer 2019;19:715. [Crossref] [PubMed]
34. Huang H, Cai X, Lin J, et al. A novel five-gene metabolism-related risk signature for predicting prognosis and immune infiltration in endometrial cancer: A TCGA data mining. Comput Biol Med 2023;155:106632. [Crossref] [PubMed]
35. Sun D, Liu J, Zhou H, et al. Classification of Tumor Immune Microenvironment According to Programmed Death-Ligand 1 Expression and Immune Infiltration Predicts Response to Immunotherapy Plus Chemotherapy in Advanced Patients With NSCLC. J Thorac Oncol 2023;18:869-81. [Crossref] [PubMed]
36. Yang D, Ma X, Song P. A prognostic model of non small cell lung cancer based on TCGA and ImmPort databases. Sci Rep 2022;12:437. [Crossref] [PubMed]
37. Xie H, Xie C. A Six-Gene Signature Predicts Survival of Adenocarcinoma Type of Non-Small-Cell Lung Cancer Patients: A Comprehensive Study Based on Integrated Analysis and Weighted Gene Coexpression Network. Biomed Res Int 2019;2019:4250613. [Crossref] [PubMed]
38. Mo XB, Zhang YH, Lei SF. Genome-wide identification of m(6)A-associated SNPs as potential functional variants for bone mineral density. Osteoporos Int 2018;29:2029-39. [Crossref] [PubMed]
39. Yazdani B, Jazini M, Jabbari N, et al. Altered expression level of ACSM5 in breast cancer: An integrative analysis of tissue biomarkers with diagnostic potential. Gene Reports 2021;22:100992.
40. Broderick C, Hoek RM, Forrester JV, et al. Constitutive retinal CD200 expression regulates resident microglia and activation state of inflammatory cells during experimental autoimmune uveoretinitis. Am J Pathol 2002;161:1669-77. [Crossref] [PubMed]
41. Coles SJ, Wang EC, Man S, et al. CD200 expression suppresses natural killer cell function and directly inhibits patient anti-tumor response in acute myeloid leukemia. Leukemia 2011;25:792-9. [Crossref] [PubMed]
42. Bisgin A, Meng WJ, Adell G, et al. Interaction of CD200 Overexpression on Tumor Cells with CD200R1 Overexpression on Stromal Cells: An Escape from the Host Immune Response in Rectal Cancer Patients. J Oncol 2019;2019:5689464. [Crossref] [PubMed]
43. Yoshimura K, Suzuki Y, Inoue Y, et al. CD200 and CD200R1 are differentially expressed and have differential prognostic roles in non-small cell lung cancer. Oncoimmunology 2020;9:1746554. [Crossref] [PubMed]
44. Graham SA, Jégouzo SA, Yan S, et al. Prolectin, a glycan-binding receptor on dividing B cells in germinal centers. J Biol Chem 2009;284:18537-44. [Crossref] [PubMed]
45. Breiman A, López Robles MD, de Carné Trécesson S, et al. Carcinoma-associated fucosylated antigens are markers of the epithelial state and can contribute to cell adhesion through CLEC17A (Prolectin). Oncotarget 2016;7:14064-82. [Crossref] [PubMed]
46. Zhang J, Ding L, Holmfeldt L, et al. The genetic basis of early T-cell precursor acute lymphoblastic leukaemia. Nature 2012;481:157-63. [Crossref] [PubMed]
47. Ito K, Oleschuk CJ, Westlake C, et al. Mutation of Trp1254 in the multispecific organic anion transporter, multidrug resistance protein 2 (MRP2) (ABCC2), alters substrate specificity and results in loss of methotrexate transport activity. J Biol Chem 2001;276:38108-14. [Crossref] [PubMed]
48. Huisman MT, Smit JW, Crommentuyn KM, et al. Multidrug resistance protein 2 (MRP2) transports HIV protease inhibitors, and transport can be enhanced by other drugs. AIDS 2002;16:2295-301. [Crossref] [PubMed]
49. El Kasmi KC, Vue PM, Anderson AL, et al. Macrophage-derived IL-1β/NF-κB signaling mediates parenteral nutrition-associated cholestasis. Nat Commun 2018;9:1393. [Crossref] [PubMed]
50. Li Q, Liu XL, Jiang N, et al. A new prognostic model for RHOV, ABCC2, and CYP4B1 to predict the prognosis and association with immune infiltration of lung adenocarcinoma. J Thorac Dis 2023;15:1919-34. [Crossref] [PubMed]
51. Ohtaki Y, Ishii G, Nagai K, et al. Stromal macrophage expressing CD204 is associated with tumor aggressiveness in lung adenocarcinoma. J Thorac Oncol 2010;5:1507-15. [Crossref] [PubMed]
52. Liu X, Wu S, Yang Y, et al. The prognostic landscape of tumor-infiltrating immune cell and immunomodulators in lung cancer. Biomed Pharmacother 2017;95:55-61. [Crossref] [PubMed]
53. Schulze AB, Evers G, Görlich D, et al. Tumor infiltrating T cells influence prognosis in stage I-III non-small cell lung cancer. J Thorac Dis 2020;12:1824-42. [Crossref] [PubMed]