Analysis of the tumor immune environment identifies an immune gene set-based prognostic signature in non-small cell lung cancer.

BACKGROUND: The tumor immune environment plays a critical role in lung cancer initiation and prognosis. Therefore, understanding how the tumor immune environment impacts the overall survival (OS) of patients with advanced lung cancer post immunotherapy is of great importance. In this article, we aimed to identify the immune components of lung cancer and develop an immune prognostic signature to predict OS.
METHODS: Differentially expressed immune-related genes were calculated between tumor and normal tissues using expression data from The Cancer Genome Atlas (TCGA) database. Then univariate Cox regression analysis was conducted to select prognosis-related genes and the prognostic risk model was constructed by multivariate Cox regression analysis. Patient risk scores were calculated, and a clinical correlation analysis was performed within the risk model. In addition, immune cell infiltration patterns were identified to find the immune cell subtypes related to prognosis.
RESULTS: A gene model consisting of 12 immune-related genes was used as our signature. The model showed that the high-risk group experienced a shorter survival time, with an area under the receiver operating characteristic (ROC) curve (AUC) of 0.733. High-risk immune genes, such as S100 calcium binding protein A16 (S100A16) and angiopoietin-like 4 (ANGPTL4), were associated with more malignant clinical manifestations. Further, we discovered that extensive infiltration of B cells, dendritic cells, and mast cells indicated a favorable prognosis.
CONCLUSIONS: The signature developed in this paper could be an effective model for estimating OS in lung cancer patients, and the immune cell infiltration analysis of the tumor immune microenvironment could shed light on more effective treatment in clinical practice. 2022 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer is the leading cause of cancer-related morbidity and mortality, with non-small cell lung cancer (NSCLC) as the most prevalent form (1). Despite recent advances in lung cancer diagnosis and treatment, the survival time for patients with lung cancer is still short, and the majority of patients present with advanced disease (2). Predicting which patients are at high risk and investigating potential therapeutic targets are therefore necessary to promote the overall survival.

The immune system plays a key role in the monitoring and destruction of cancer cells; however, tumor cells can bypass this natural defense, and tolerance can form during disease progression (3). It has long been proposed that tumor immune environments are associated with tumor initiation and can influence therapeutic response (4). Tumor-reprogrammed immune cells are involved in both tumor-promoting and tumor-suppressing activities, which may predict clinical outcome (5). In this regard, immunotherapy focus on restoring the anti-tumor host immunity could promote occult cancerous cells destruction. Furthermore, patients who response to immunotherapeutic intervention would have better survival rates (6). Recently, immunotherapy has been recommended as the first-line therapy for patients with advanced lung cancer (7). However, the majority of patients with lung cancer do not benefit from this treatment, and so a deeper understanding of tumor immune environments is needed.

Gene expression profiling can reflect the interactions between tumor and immune cells and provide clues for finding predictive biomarkers and new targets (8). To this end, we calculated differentially expressed genes between tumor and normal tissues to uncover potential important players in tumor immune environments. An immune-related gene signature was then established to predict the prognosis of patients with lung cancer. Additionally, we described immune cell infiltration patterns in NSCLC. In this study, we aimed to provide the comprehensive description of immune characteristics of lung cancer, offering reference information for improving NSCLC treatments. We present the following article in accordance with the TRIPOD reporting checklist (available at https://dx.doi.org/10.21037/atm-21-6043).

Methods

Data acquisition and processing

The gene expression data of 497 lung adenocarcinoma tumor tissue samples and 54 normal lung tissue samples and the corresponding clinical information were acquired from The Cancer Genome Atlas (TCGA) database. A Wilcoxon test was used to analyze differentially expressed genes in the TCGA sample using the “limma” package in R (The R Foundation for Statistical Computing, Vienna, Austria). To identify differentially expressed genes in lung cancer, the cutoff threshold in TCGA was set as |log2 fold change (FC)| ≥1.0, and the false discovery rate (FDR) was set at <0.01. Immune-related genes and transcription factors (TFs) were extracted from the above-mentioned differentially expressed genes and analyzed for differential expression using the “limma” package. The corresponding clinical information of the patients with lung cancer was collected and used for the subsequent analyses. The external validation cohort consisted of the expression data and the comparative clinical data of 45 patients with squamous cell carcinoma, acquired from the Gene Expression Omnibus (GEO) database (GSE50081, n=181). The study was conducted in accordance with the Declaration of Helsinki (as revised in 2013).

Co‑expression network construction and functional enrichment analyses

A network of immune-related genes and TFs was created using Cytoscape v3.7.1 (Cytoscape Consortium, San Diego, CA, USA). Upregulated genes were represented by red circles, downregulated genes were represented by green ellipses, and TFs were represented by triangles. To identify the biological function of differentially expressed immune-related genes, we used the “clusterProfiler” package in R to perform Gene Ontology (GO) analyses. The functional enrichment of upregulated immune-related genes was represented by red bar graphs, and the functional enrichment of downregulated immune-related genes was represented by blue bar graphs.

Prognostic model construction

We used the “survival” package in R to perform a univariate Cox regression analysis for all differentially expressed immune-related genes and screened for significant candidate genes. Subsequently, based on the screened candidate genes, a prognostic model containing 12 immune-related genes was constructed, combining the effects of each gene. The coefficients of each gene were determined by multivariable Cox regression model. The risk score for each patient was determined using the following formula: risk score = Exp1 × Coe1 + Exp2 × Coe2 +……+ Expi × Coei.

According to the median value of the risk score, the patients with lung cancer were classified into high-risk and low-risk groups. Kaplan-Meier analysis was used to construct a survival curve. Then, a logrank test was applied to compare the overall survival (OS) of the 2 subgroups. Subsequently, a receiver operating characteristic (ROC) curve was drawn to evaluate the performance of the prognostic model using the “survivalROC” package in R. The correlation between the immune-related genes that constituted the prognostic model and various clinical factors were analyzed using the “beeswarm” package in R.

Quantitative real-time polymerase chain reaction (qRT-PCR)

Total RNA was extracted using TRIZOL (Invitrogen Inc., USA), and subjected to qRT-PCR using the following primers: VIPR1 (TCATCCGAATCCTGCTTCAGA; AGGCGAACATGATGTAGTGTACT), S100A16 (ATGTCAGACTGCTACACGGAG; GTTCTTGACCAGGCTGTACTTAG), GAPDH (GGCTCATGACCACAGTCCAT; GACGGACACATTGGGGGTAG). Gene expression was normalized to the level of GAPDH within each sample using the relative ∆∆Ct method.

Unsupervised clustering and analysis of immune cell-related genes

The abundance of tumor infiltrating immune cells was calculated using the CIBERSORT website (https://cibersort.stanford.edu/). Based on the abundance of immune cell infiltration provided by CIBERSORT, we used the “ConsensusClusterPlus” package in R to perform unsupervised clustering of the patients with lung cancer. We then visualized programmed cell death protein 1 (PD-1) infiltration between the 2 cluster subgroups with the “ggplot2” and “ggpubr” packages in R.

Statistical analysis

We used Kaplan-Meier analysis to construct survival curves using the “survival” and “survminer” packages in R and a log rank test to evaluate the significance of differences between the 2 subgroups. Wilcoxon signed-rank tests were performed to explore quantitative variables. Significance was defined as P<0.05. All statistical analyses were performed using R version 4.0.3.

Results

Differentially expressed immune genes and TFs in lung cancer

We used data from the TCGA dataset to calculate differentially expressed immune genes between tumor and normal tissues. A total of 489 differentially expressed genes were identified, in which 321 genes were upregulated in tumor tissues (). Differentially expressed TFs were then selected (), and a network of TFs and genes was constructed (). Among TFs, regulatory molecules enhancer of zeste homolog 2 (EZH2) (9) and transcription factor 21 (TCF21) (10) are involved in the immunosuppressive microenvironment, and chromobox homolog (CBX7) (11) is associated with lung carcinogenesis. GO analysis was carried out to explore the potential functions of these differentially expressed genes (). For upregulated genes in tumors, immunoglobulin complex, humoral immune response, complement activation, and lymphocyte-mediated immunity were the most populated subcategories, indicating that the immune system plays a key role in tumor progress. For downregulated genes, signaling receptor activator activity, receptor ligand activity, and cell chemotaxis were the most abundant activities, indicating that homeostasis-related terms are disrupted in the tumor microenvironment and potentially robust immunity may be suppressed. Our enrichment analysis results were consistent with pathological changes in lung cancer.

Figure 1

A heat map (A) and a volcano map (B) of differentially expressed immune genes; a heat map (C) and a volcano map (D) of differentially expressed transcription factors; network of transcription factors and genes (E). FDR, false discovery rate.

Figure 2

GO and pathway analysis of differentially expressed immune genes. The blue mark indicates downregulated genes, while the red mark indicates upregulated genes. GO, Gene Ontology.

Immune signature for lung cancer prognosis

Nineteen significant prognostic immune genes were identified using a univariate Cox regression analysis with P<0.05 (). Related genes then underwent a Cox proportional hazards regression to generate the best gene model. Finally, a gene model with 12 immune genes (Figure S1) was created. We calculated the risk score for each patient according to the coefficient value of the 12 genes. Subsequently, the patients were divided into high-risk and low-risk groups based on the median risk score. A longer survival time was found in the low-risk patients when compared with that of the high-risk patients (), and the risk score distribution between the 2 groups was shown in . To further evaluate the prognostic accuracy of the risk scores, time-dependent ROC curves were plotted (), yielding an area under the curve (AUC) of 0.733. To investigate the stability of our model, we assessed the prognostic value of the immune signature in the GEO database. The patients in the GEO cohort were similarly divided into a high-risk group (n=90) and a low-risk group (n=91). The survival of the patients in the low-risk group was significantly better than that of the patients in the high-risk group (P<0.001; Figure S2). Furthermore, we validated the representative genes expression in 8 matched pairs of primary lung tumors and the adjacent normal samples. Consistent with the previous findings, the protective gene vasoactive intestinal peptide receptor 1 (VIPR1) showed decreased expression level in tumors compared with adjacent normal samples, while the risky gene S100A16 expression was up-regulated in tumors (Figure S3).

Figure 3

A forest diagram of differentially expressed immune genes. The red mark indicates that the HR value of the immune gene is greater than 1, and the green mark indicates that the HR value is less than 1. HR, hazard ratio.

Figure 4

Survival analyses of two group patients. The survival was significantly different between the high-risk group and low-risk group (A) (P<0.05). The risk model and survival state diagram based on risk scores (B). The area under the ROC curve (C). ROC, receiver operating characteristic.

Association with clinicopathologic factors

A univariate Cox regression analysis was conducted using the risk score and various clinical characteristics. The analysis showed that clinical stage, degree of tissue involvement (T), lymphatic involvement (N), tumor metastasis (M), and risk score were correlated with the OS of patients (). A subsequent multivariate Cox regression analysis showed that only the risk score and the stage were associated with the OS of patients with lung cancer (P<0.05; ), proving the accuracy of the risk score in predicting survival outcomes. In addition, we assessed the relationship between the genes in our model and clinicopathologic factors. We found that genes related to poor survival, such as S100A16 and angiopoietin-like 4 (ANGPTL4), were significantly associated with T, N, and stage, and that the expression of higher-risk genes was correlated with more malignant lung cancer. In contrast, increased expression of the protective gene VIPR1 was correlated with a decreased degree of malignancy (N; ).

Figure 5

Univariate Cox regression analysis showed that clinical stage, T, M, N, and risk score were correlated with the overall survival of patients with lung cancer (A) (P<0.05). Multivariate Cox regression analysis indicated that stage and risk score were associated with overall survival of patients with lung cancer (B). T, degree of tissue involvement; N, lymphatic involvement; M, tumor metastasis.

Figure 6

The association of survival-relevant immune gene expression and clinicopathologic factors. T, degree of tissue involvement; N, lymphatic involvement.

Immune cell infiltration patterns

To illustrate the immune cell infiltration patterns in tumors, we used a gene expression-based approach to estimate the abundances of specific cell types. As shown in , most immune cell populations exhibited a positive correlation without any anticorrelation, indicating the reliability of the gene sets. Unsupervised clustering showed that the tumor samples were predominantly separated into 2 clusters; more memory B cells, mast cells, and activated dendritic cells were found in cluster A, while more CD8 T cells, Treg T cells, resting natural killer (NK) cells, and M0 macrophages (an intermediate macrophages with differentiative ability) were found in cluster B (). The patients in cluster A had better OS compared with that of the patients in cluster B () and lower PD-1 expression ().

Figure 7

Immune characteristics of lung cancers. Correlation of the immune cells in lung cancer (A); unsupervised clustering separates the patients into 2 clusters (B); the survival prognosis (C) and PD-1 expression (D) between patients in 2 clusters. PD-1, programmed cell death protein 1. ***, P<0.01.

Discussion

The immune environment plays a critical role in tumor initiation and therapeutic response (4). Lung cancer typically harbors extensive genomic alterations and is responsive to immunotherapy (12). Immune cells influence patient survival and the effectiveness of treatments. Previous studies have shown that levels of tumor-infiltrating lymphocytes are significantly correlated with prognosis in NSCLC (13). The neutrophil-to-lymphocyte ratio has also been recognized as an indicator of poor prognosis (14). Furthermore, research has found that lung cancer exhibits increased inflammatory gene transcripts and that circulating T cell subpopulations can predict immune responses (15). Traditional methods to assess immune cell infiltration mainly rely on immunohistochemistry, but the number of cell types and the characteristics of each immune cell are limited from section staining, confined by the quantity and quality of antibodies. The gene expression profile of bulk tumors provides more information in this regard. The immune environment, including immune cells and proteins expressed by immune-related genes, has a significant impact on predicting the clinical outcomes of cancers (16). Furthermore, gene expression profiling of lung cancer offers useful information for cancer diagnosis, prognosis as well as new therapy development (17). Research has shown that more than 3,000 protein-coding genes change in NSCLC after tumor resection, especially genes in pathways associated with the innate immune response suppression (18). In this study, we collected differentially expressed immune genes between tumor and normal tissues and constructed a prediction model based on prognostic genes. In our model, the patients in the high-risk group had a poor prognosis with shorter survival time.

According to previous studies, the genes used in our risk prediction model are tightly associated with lung cancer prognosis. Some of the risk genes in our model participate in lung cancer progression. Higher S100A16 messenger RNA (mRNA) expression in lung cancer is significantly correlated with poor survival (19). Retinol binding protein 2 (RBP2) performs an oncogenic function in lung tumorigenesis and progression (20) and induces epithelial-mesenchymal transition (21). Fibroblast growth factor 2 (FGF2) promotes lung cancer cell proliferation (22), and some clinical trials have evaluated the efficacy of its inhibitors on lung cancer treatment (23). In addition, interleukin 11 (IL11) acts on leukocytes to induce pro-tumorigenic and pro-metastatic neutrophils (24). In contrast, some of the tumor suppressor genes of our model, such as VIPR1, have been shown to have protective effects (25). Research has shown that VIPR1 significantly inhibits the growth, migration, and invasion of lung cancer cells (26). Moreover, VIPR1 expression is decreased in metastatic lung cancer compared with that in non-metastatic lung cancer (27).

In our study, we described the immune cell infiltration patterns in lung cancer and classified patients into 2 groups based on immune cell population. The patients in cluster B had worse OS with more M0 macrophages and Treg cells. It has been reported that M0 macrophages are associated with decreased survival (28). Treg cells, as the most notable immunosuppressive cells, can help cancer cells escape from immunological surveillance, and therapies based on Treg blockade have benefited a group of patients with lung cancer (29). The patients in cluster A exhibited more infiltration of B cells, dendritic cells, and mast cells. Consistent with previous studies, the high infiltration of B cells, dendritic cells, and mast cells indicated a favorable prognosis (30). Furthermore, the expression of PD-1 was higher in cluster B. PD-1 is an inhibitory checkpoint molecule which can inhibit T cell activation. Immune checkpoint blockade therapy targeting PD-1 has yielded promising clinical responses in patients with lung cancer (31). Therefore, we propose that the patients in cluster B would be more likely to benefit from immunotherapy.

In recent years, increasing attention has been paid to tumor immune microenvironments, and immune prognostic signatures have been constructed for various cancer types. For NSCLC, a gene model consisting of 30 immune-related genes has been shown to have a strong association with recurrence and clinical stage (32). Immune-associated risk score prognostic models have also been established (33), with an AUC of 0.702 and 0.723, respectively. Our prognostic model showed a similar true positive rate compared with that of previous studies. In addition to an immune gene signature, we established a connection between genes and TFs and the related pathway. We also selected representative genes of different immune cells to identify immune cell infiltration patterns and pave a path to better targeting of prognosis-related immune cell subsets. Our retrospective analysis provides comprehensive information about lung cancer immune environments, although further studies are required to confirm our findings.

Conclusions

In summary, we constructed a risk score prognostic model based on 12 immune-related genes. The survival analysis revealed that a high-risk score was significantly associated with poor survival outcomes. In addition, we described the immune cell infiltration patterns and found that high infiltration of B cells, dendritic cells, and mast cells indicated a favorable prognosis.

The article’s supplementary files as