Analysis of expression differences of immune genes in non-small cell lung cancer based on TCGA and ImmPort data sets and the application of a prognostic model.

BACKGROUND: There has been little investigation carried out into the activity of immune-related genes in the prognosis of non-small cell lung cancer (NSCLC). Our study set out to analyze the correlation between the differential expression of immune genes and NSCLC prognosis by screening the differential expression of immune genes. Based on the immune genes identified, we aimed to construct a prognostic risk model and explore some novel molecules which have predictive potential for therapeutic effect and prognosis in lung cancer.
METHODS: Immune gene transcriptome data and clinical data of NSCLC samples were extracted from TCGA database, and transcription factors in the ImmPort dataset were obtained. The data were divided into two groups: normal tissues and tumor tissues. The expression levels of immune genes were compared using the edgeR algorithm, and then differential expression analysis was performed. The survival analysis was carried out by combining differential immune genes with clinical survival time, so that the immune genes influencing the prognosis of NSCLC could be determined. A risk score was calculated based on the expression levels of the immune genes related to the prognosis of NSCLC and their corresponding coefficients to construct a prognostic risk model. This model was used to calculate patient risk scores and perform clinical correlation analysis. The selected molecules were further verified by clinical samples.
RESULTS: By comparing NSCLC tissues with normal tissues, a total of 6,778 differentially expressed genes were found (P<0.05), of which 490 were differential immune-related genes. Survival analysis determined 28 differential immune genes to be associated with prognosis (P<0.05). We calculated the patient risk value based on the immune gene prognosis model. The survival curve was drawn according to the patient risk score and showed that the survival prognosis was significantly different for the high-risk and the low-risk groups (P<0.05). The area under the receiver operating characteristic (ROC) curve (AUC) was 0.723, which represented a relatively high true-positive rate. All of the results proved the reliability of our immune gene risk prognostic model. After drawing the risk curve, S100A16, IGKV4, S100P, ANGPTL4, SEMA4B, and LGR4 were found to be the high-risk immune genes in NSCLC. Clinical correlation analysis of survival-related differential immune genes revealed that in patients with lymph node metastasis, ANGPTL4 was positively correlated with T stage, S100a16 and SEMA4B were upregulated, and VIPR1 was downregulated. Further analysis revealed that VIPR1 was decreased in metastatic lung cancer compared to non-metastatic lung cancer. Furthermore, the real-time PCR detection of the clinical samples showed that S100A16 expression in lung cancer was increased, while VIPR1 expression in lung cancer was downregulated, which was consistent with the results of our bioinformatics analysis.
CONCLUSIONS: Based on big data from the TCGA and ImmPort databases, our study analyzed the relationship between differential expression of immune-related genes and clinical data, and constructed a prognostic model based on the immune genes identified. Two novel molecules, S100A16 and VIPR1, were verified to possibly have significant biological function in NSCLC. Our research may provide us with new insight into the immune genes by which the malignant biological behavior of NSCLC is mediated. 2020 Annals of Translational Medicine. All rights reserved.

Introduction

As the most common type of tumor in the world (1), lung cancer has the highest mortality rate among all cancers, and accounts for 23% of the total number of tumor-related deaths, more than any other cancer (2). Around 75% to 80% of all lung cancer cases are non-small cell lung cancer (1,3).

With the development of modern molecular biology and tumor immunology, it has gradually been recognized that tumor occurrence and progression are not only related to the internal genetic background of cancer cells, but are also dependent on the interaction between the tumor and various systems within the body, especially the immune system (4,5).

Immune-related cells and factors are involved throughout the entire process of tumor occurrence, proliferation, and development, and have a significant influence on the tumor (6-8). Particularly with the development of monoclonal antibody drugs and adoptive immunotherapy, tumor immunity has become a hot spot in tumor research. To date, there has been no in-depth study into the relationship between the screening of immune-related molecules, their expression in lung cancer, and their effect on lung cancer prognosis (7). Therefore, it is necessary to explore the characteristics of immune-related molecules in lung cancer patients, and to evaluate the relationship between the components of the immune environment of lung cancer at the cell, molecule, and gene levels, to identify new targets that can predict the therapeutic effect and prognosis of lung cancer.

In recent years, with the development of gene expression profile chips and second-generation high-throughput sequencing technology, the quantity of lung cancer expression profile data has exploded, providing a basis for a comprehensive study of differentially expressed genes and their role in lung cancer prognosis (9,10). Based on big data from the TCGA and ImmPort databases, this study analyzed the differences in the expression of immune genes in NSCLC and clinical survival time, and constructed a prognostic model based on immune genes. The feasibility of the immune gene prognosis model in clinical application was verified.

Methods

Data download

The transcriptome data and clinical data of all NSCLC in TCGA-LUAD and TCGA-LUSC were downloaded through the Genomic Data Commons (GDC) Data Portal, and 594 clinical data were obtained. The immune gene data was downloaded through the ImmPort data portal, and 2,498 immune-related genes were obtained. Transcription factor data was downloaded from the Cistorm website.

Differential expression analysis

Differential gene expression analysis

All transcriptome data were analyzed for differences using the Wilcoxon test in R software. Genes with significant differences in expression between normal and tumor tissues were screened out. The screening criteria were |logFC| >1, FDR <0.05.

Analysis of immune gene differences

The differential genes were combined with the acquired immune gene data and analyzed in R software to screen out the differential immune genes from all the differential genes.

Construction and evaluation of immune gene prognosis model

The expression of immune genes was combined with survival time and survival status to carry out a survival analysis of the differential immune genes and clinical survival time, and to identify the immune genes that affect the prognosis of NSCLC. Based on these genes, an immune gene prognosis model was constructed. The ROC and risk score curves were drawn in R to evaluate the model.

Clinical correlation analysis

The prognosis-related immune genes were combined with clinical data, and the patient’s risk value was calculated based on the immune gene prognosis model. In addition, survival-related differential immune genes were also subjected to clinical correlation analysis.

Real-time polymerase chain reaction (PCR)

The original concentrations of specific transcripts were measured by real-time PCR using an ABI PRISM 7300 sequence detector (Applied Biosystems, Foster City, CA, USA). The primer sequences are provided in .

Table S1

Sequence of primers for real-time PCR

Primer	Sequence (5' to 3')
S100A16 forward primer (homo)	ATGTCAGACTGCTACACGGAG
S100A16 reverse primer (homo)	GTTCTTGACCAGGCTGTACTTAG
VIPR1 forward primer (homo)	TCATCCGAATCCTGCTTCAGA
VIPR1 reverse primer (homo)	AGGCGAACATGATGTAGTGTACT
SEMA4B forward primer (homo)	GGCCCTCTTTGCACTTAACAG
SEMA4B reverse primer (homo)	TGTAGTTTTGACAGTCACGCTT
S100P forward primer (homo)	TCATCCGAATCCTGCTTCAGA
S100P reverse primer (homo)	AGGCGAACATGATGTAGTGTACT
ANGPTL4 forward primer (homo)	GTCCACCGACCTCCCGTTA
ANGPTL4 reverse primer (homo)	CCTCATGGTCTAGGTGCTTGT
18s forward primer	CGGCTACCACATCCAAGGAA
18s reverse primer	GCTGGAATTACCGCGGCT

Ethics statement

All of the 29 patients (including lung cancer tissues and paired adjacent tissues) involved in this study underwent treatment at the Department of Thoracic Surgery, the First Affiliated Hospital of China Medical University (Shenyang, China). All clinical investigations were approved by the hospital’s Ethical Committee and all of the patients submitted signed consent forms.

Results

Differential gene expression analysis

After screening by Wilcoxon test, 6,778 differential genes met the conditions. All of the differential genes screened were intersected with immune genes, and 490 differential immune genes were obtained. The heat map showed distinct clustering of differentially immune genes in cancerous and adjacent tissues () and the volcano plot suggested clustering of differentially immune genes of high or low expression level ().

Figure 1

A heat map (A) and a volcano map (B) of differential immune genes. The heat map abscissa represents the sample: the blue area represents normal tissue and the red area represents tumor tissue; the ordinate represents the gene. On the volcano map, the green area represents the downregulated differential genes and the red area represents the upregulated differential genes.

Prognostic related immune gene analysis

Survival analysis of differential immune genes and clinical survival time and survival status yielded 28 differential immune genes related to prognosis: S100P, S100A16, CRABP1, FURIN, FGF2, F2RL1, BIRC5, OAS1, PPARG, BTK, IGKV4-1, SEMA4B, SEMA7A, CX3CR1, CAT, GPI, IL11, INHA, INSL4, TNFSF12, TNFSF13, ADRB2, ANGPTL4, LGR4, TNFRSF11A, VIPR1, SHC3, and CBLC. Among these immune genes, 18 carry high risk, and the other 10 carry low risk (). Survival analysis to identify prognosis-related immune genes found that 9 of the 28 immune genes were significantly related to patient survival: ANGPTL4, IL11, INHA, LGR4, S100A16, S100P, SEMA4B, SHC3, and VIPR1.

Figure 2

A forest diagram of the differential immune genes. The red mark indicates that the HR value of the immune gene is greater than 1 (high risk), and the green mark indicates that the HR value of the immune gene is less than 1 (low risk).

Immune gene prognostic model

The risk scores of prognosis-related immune genes were calculated with R software () and combined with the patients’ clinical data, including survival time and survival status. A prognostic risk model based on the prognosis-related immune genes identified was constructed, and a survival curve was drawn. According to the calculated patient risk score, the patients were divided into a high-risk group and a low-risk group, which showed a significant difference in the survival prognosis of the two groups (P<0.05) (). The 5-year survival rate of the high-risk group was approximately 23.7% (95% CI: 16.25–34.7%) and the 5-year survival rate of the low-risk group was approximately 57.9% (95% CI: 48.7–68.9%). The AUC value was 0.723, which represented a higher true positive rate, thus proving the accuracy of the survival curve ().

Table 1

The immune gene prognostic model

ID	Coefficient
S100P	−0.00034
S100A16	0.001593
CRABP1	0.004666
FGF2	0.179618
IGKV4-1	−0.00037
SEMA4B	0.005702
IL11	0.163065
INHA	0.008009
INSL4	0.010981
ANGPTL4	0.004859
LGR4	0.012764
TNFRSF11A	0.240902
VIPR1	−0.08722
SHC3	−0.1968

Figure 3

The survival prognosis was significantly different between the high-risk group and low-risk group (P<0.05) (A). The AUC value is 0.723, which represents a higher true positive rate, thus proving the accuracy of the survival curve (B). AUC, area under the receiver operating characteristic curve.

The abscissa from left to right indicates that the risk value increases gradually. After drawing the risk curve, the significant molecules were identified. The risk diagram, survival state diagram, and heat map based on the risk scores all showed that the higher the risk score, the worse the patient's prognosis ().

Figure 4

Risk model diagram (A), survival state diagram (B), and heat survival diagram (C) based on risk scores.

Clinical correlation analysis of the survival-related differential immune genes revealed that in patients with lymph node metastasis, ANGPTL4 was positively correlated with T stage, S100a16 and SEMA4B were upregulated, and VIPR1 was downregulated. Further analysis revealed that VIPR1 was decreased in metastatic lung cancer compared to non-metastatic lung cancer ().

Figure 5

Immune gene ANGPTL4 had significant clinical correlation with T stage (A). S100A16, SEMA4B, VIPR1 had significant clinical correlation with N stage (B-D). VIPR1 had significant clinical correlation with M stage (E).

The survival analysis of the above genes revealed that the immune genes ANGPTL4, S100A16, SEMA4B, VIPR1, and S100P were significantly related to patient survival ().

Figure 6

The survival curves of the immune genes ANGPTL4 (A), S100A16 (B), S100P (C), SEMA4B (D), and VIPR1 (E). P<0.05.

Validation of clinical samples

To verify the expression of the selected immune-related molecules in lung cancer and adjacent tissues, 29 cases of clinical samples were analyzed. The results showed that the level of S100A16 was elevated in lung cancer compared to the adjacent tissue, thus it could possibly act as an oncogene (). In contrast, the expression of VIPR1 was found to be lower in lung cancer tissues than in the adjacent tissues, thus it could play a role as a tumor suppressor. S100A16 and VIPR1 exhibit potential to serve as clinical markers of lung cancer.

Figure 7

The expression of immune-related genes in the lung cancer samples (n=29). *, P<0.05.

According to the pathology of 29 lung cancer patients, they were divided into two groups, well-differentiated lung cancer and poorly differentiated lung cancer, and analyzed the expressions of S100A16 and VIPR1. We found that the expression of these two molecules is closely related to the pathological differentiation of lung cancer, which holds potential for lung cancer prognosis. The results showed that S100A16 was elevated in poorly differentiated lung cancer than well-differentiated lung cancer (). In contrast, the expression of VIPR1 was decreased in the poorly differentiated group compared to the well-differentiated group. Both S100A16 and VIPR1 exhibit potential to serve as clinical predictors of lung cancer.

Figure 8

The classification of S100A16 and VIPR1 expression in the lung cancer samples (n=29).

Discussion

Over the past decade, it has been recognized that tumor occurrence and progression should not only be attributed to the internal genetic background of cancer cells, but also to the interaction the tumor has with various systems within the body, especially the immune system (11,12). Immune-related cells and factors are involved throughout the whole process of tumor occurrence, proliferation, and development, and significantly influence the tumor (13). Therefore, there is an increasing need to explore the characteristics of immune-related molecules and to evaluate the function of immune genes in lung cancer (14,15). By doing this, more immune mechanisms and early warning effects of immunotherapy can be revealed. In this study, our prognostic model and survival analysis screened out immune genes that are significantly related to the prognosis and survival of lung cancer patients. Our research could provide an effective reference for clinical treatment of lung cancer. Current research has confirmed that immune genes are related to various biological behaviors, such as tumor development, metastasis, and apoptosis (16-19). At the same time, these immune genes also can be regarded as a new target of immunotherapy.

In our study, we identified the important immune-related molecules in NSCLC. Through our analysis of biological information and the prognosis analysis of lung cancer patients, we preliminarily screened five immune genes with potential early warning significance. However, in the subsequent validation of clinical samples, we finally focused on two molecules, S100A16 and VIPR1. After a further literature search, we found these two molecules to have an important pathologic effect, despite having been rarely discussed in the study of NSCLC (20,21).

S100A16 is a member of the S100 family and a partner of the Annexin family of proteins. It has been reported to be significantly upregulated in many tumors, such as HCC and breast cancer (22,23). S100A16 expression has also been associated with poor prognosis in patients with breast cancer and with the promotion of the invasive activity of breast cancer cells via interaction with cytoskeleton dynamics. In addition, the upregulation of S100A16 has been found to promote epithelial-mesenchymal transition via the Notch1 pathway in breast cancer (22). Conversely, the downregulation of S100A16 expression in oral squamous cell carcinoma specimens has been associated with poor prognosis and differentiation grade (24). Furthermore, in oral squamous cell carcinoma, S100A16 has been demonstrated to contribute to a less aggressive tumor phenotype. These findings imply that S100A16 plays distinct roles in tumor development and progression, depending on the pathologic type of cancer (21,23). However, the role of S100A16 in NSCLC remains unclear. Our findings concluded that S100A16 is a promising candidate as a prognostic marker and molecular hub in NSCLC.

Vasoactive intestinal peptide (VIP), also known as PACAP, is an important neuropeptide that controls cell physiology and mainly functions via two receptor subtypes VAPC1 and VAPC2 (25,26). VIP is a small neuropeptide involved in the relaxation of smooth muscle, exocrine, endocrine, and hydration ion flux in intestinal epithelial cells (27,28). Research also showed that the reduction of VIP may be related to the occurrence of disease, and overexpression of VIP could inhibit inflammation, therefore reducing acute lung injury in mice (29). Other studies also found that the immune response and survival rate of VIP knockout mice increased (30). VIP or VPAC1 receptor antagonist could strengthen the ability of chemotherapy by killing breast cancer cells and improving anti-viral immunity (31). These findings all indicate that, because of its potential role in tumor development and progression, VIPR can be used in tumor diagnosis and as a treatment target.

The function of VIP/PACAP and its receptor (VPAC1 and VPAC2) in normal human tumors and their role in potential new treatments were summarized by Dr. Teddy Moody (32). Some studies have reported that VIPR1 inhibits the growth of several cancers, including prostate cancer, lymph blastoma, and medulloblastoma, which suggests that VIPR1 may significantly inhibit the growth and development of cancer cells (33,34). In the present study, overexpression of VIPR1 was found to reduce the proliferative capacity of H1299 cells (35). The experimental results are consistent with those of previous studies and indicate that VIPR1 plays an important role in the inhibition of LUAD cell growth. It may provide a novel target for drug intervention in patients with CAP expressing VIPR1. The results of our research are consistent with the results of current in-vivo experiments and further explain the important role of VIPR1 in lung cancer. Our findings lay a solid foundation for further study of the clinical transformation of lung cancer.

Conclusions

In conclusion, our results present two innovative molecules, S100A16 and VIPR1, which play an important role in NSCLC and hold potential as predictors of lung cancer prognosis. Our results also provide a new insight into the immune genes that mediate the malignant biological behavior of lung cancer.

The article’s supplementary files as