Bioinformatic analysis of RNA-seq data from TCGA database reveals prognostic significance of immune-related genes in colon cancer.

The tumor immune microenvironment is of crucial importance in cancer progression and anticancer immune responses. Thus, systematic exploration of the expression landscape and prognostic significance of immune-related genes (IRGs) to assist in the prognosis of colon cancer is valuable and significant. The transcriptomic data of 470 colon cancer patients were obtained from The Cancer Genome Atlas database and the differentially expressed genes were analyzed. After an intersection analysis, the hub IRGs were identified and a prognostic index was further developed using multivariable Cox analysis. In addition, the discriminatory ability and prognostic significance of the constructed model were validated and the characteristics of IRGs associated overall survival were analyzed to elucidate the underlying molecular mechanisms. A total of 465 differentially expressed IRGs and 130 survival-associated IRGs were screened. Then, 46 hub IRGs were identified by an intersection analysis. A regulatory network displayed that most of these genes were unfavorable for the prognosis of colon cancer and were regulated by transcription factors. After a least absolute shrinkage and selection operator regression analysis, 14 hub IRGs were ultimately chose to construct a prognostic index. The validation results illustrated that this model could act as an independent indicator to moderately separate colon cancer patients into low- and high-risk groups. This study ascertained the prognostic significance of IRGs in colon cancer and successfully constructed an IRG-based prognostic signature for clinical prediction. Our results provide promising insight for the exploration of diagnostic markers and immunotherapeutic targets in colon cancer.

1. Introduction

Colon cancer, a tumor of the large intestine, is the most common human malignancy in the digestive system.[ Many factors, such as age, gender, dietary habits, geography, and genetic background, are involved in the occurrence and development of colon cancer.[ Although most patients suffering from colon cancer can get better if the cancer is detected early, the young-onset colon cancer incidence rate is increasing and numerous patients are diagnosed at an advanced stage in developing countries.[ Thus, it is important and essential to investigate early diagnostic and prognostic biomarkers as well as underlying molecular mechanisms of differentially expressed genes (DEGs) or proteins in colon cancer.

Due to the rapid development of large-scale sequencing technology, many studies have focused on exploring valuable molecules in human colon cancer, including long noncoding RNAs,[ alternative splicing events,[ tumor-infiltrating immune cells (TIICs),[ and immunoscores.[ The tumor immune microenvironment is a battleground for tumor cells and the immune system during the neoplastic process and plays an important role in the proliferation, metastasis, and immune escape of tumor cells.[ The composition, content, properties, and function of immune cells in the tumor microenvironment are closely associated with the clinical outcomes of multiple tumors.[ Recently, researchers have found that immune-related genes (IRGs) differentially expressed in cancers can reflect the immune status and display considerable promise in the prognosis of cancer patients.[ The current studies have demonstrated that IRGs display high prognostic performance in predicting the outcomes of colorectal cancer.[ Although a large number of deaths from rectal cancer are misclassified as colon cancer, these cancers are not similar to each other. For example, the incidence rate of colon cancer is approximately 2.5 times higher than that of rectal cancer, whereas rectal tumors are more common in people aged younger than 50 years and have a better prognosis for patients.[ Therefore, it is necessary to investigate the expression profiles of IRGs and develop an independent prognostic signature for colon cancer prediction.

This study aimed to estimate the prognostic value of IRGs in colon cancer and develop an independent prognostic signature for outcome prediction using a series of bioinformatic methods. The transcriptomic RNA-seq data of 470 colon cancer patients from The Cancer Genome Atlas (TCGA) database,[ and their corresponding clinicopathological information were obtained. Then, the differentially expressed IRGs (DE-IRGs) that were also associated with the overall survival (OS) of patients were screened for the development of an independent indicator. Thus, it is of great significance for further discovery of diagnostic and prognostic markers of colon cancer and for understanding the clinical significance of the tumor immune microenvironment.

2. Materials and Methods

2.1. Data collection

The workflow of this study is presented in Figure S1 (Supplemental Digital Content, http://links.lww.com/MD/G959). The transcriptomic RNA-seq data and clinicopathological information for colon cancer patients were downloaded from the TCGA database. The RNA-seq data including 470 primary colon cancer tissues and 41 normal tissues are shown in Supplemental Digital Content (Table S1, Supplemental Digital Content, http://links.lww.com/MD/G959). The clinical information included age, sex, TNM stage, and OS. The primary tumor characteristics and clinical information are shown in Supplemental Digital Content (Table S2, Supplemental Digital Content, http://links.lww.com/MD/G959). In addition, a list of IRGs were downloaded from the Immunology Database and Analysis Portal (ImmPort) database.[

2.2. Analysis of DE-IRGs

The genes expressed in colon cancer and normal tissues were analyzed using package language R (v3.3.2) and Bioconductor. A false discovery rate of <0.05 and a log2|fold change|>1 as the cutoff values were set for the identification of DEGs. Next, the DE-IRGs were screened out from these genes for further analyses. Simultaneously, a univariate Cox analysis was deployed to select survival-associated IRGs by assessing the relationships between IRGs and the clinical outcomes of colon cancer patients. The hazard ratio (HR) and P value were calculated, and the difference was considered significant at P < .05. Subsequently, the hub IRGs were determined by intersection analysis of DE-IRGs and survival-associated IRGs.

2.3. Analysis of IRG characteristics

The potential biological function of the survival-associated IRGs was analyzed by Gene Ontology function and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The protein–protein interaction network was performed based on the STRING online database (https://stringdb.org/) and visualized using Cytoscape software version 3.7.1.[ The molecular characteristics of the hub IRGs including gene mutations and copy number variations were derived from cBioPortal (http://www.cbioportal.org/).[ In addition, the Cistrome Cancer web resource (http://cistrome.org/CistromeCancer/) was used to analyze the regulatory network between the hub IRGs and transcription factors.[

2.4. Construction of prognostic signature

A least absolute shrinkage and selection operator (LASSO) Cox regression analysis was conducted to screen candidate IRGs from the identified hub IRGs for the development of a risk model.[ A Kaplan–Meier test was performed to illustrate the survival probability of the constructed risk model and the prognostic validity was assessed by creating a receiver operating characteristic curve. According to the risk score, high- and low-risk groups for patients with colon cancer were separated. Then, the discriminatory capability of the model was evaluated in colon cancer patients according to the risk scores. Univariate and multivariate Cox regression analyses referring to age, sex, TNM stage, and risk score were performed to assess the constructed prognostic model. In the end, the abundance of TIICs, including B cells, CD4+ T cells, CD8+ T cells, dendritic cells, macrophages, and neutrophils, was analyzed and their relationships with the risk score were visualized.[

3. Results

3.1. Expression of IRGs in colon cancer

The RNA-seq data of the colon cancer cohort were obtained from the TCGA database. After a contrastive analysis, a total of 6420 DEGs, including 4503 upregulated and 1917 downregulated genes were identified (Fig. 1A, C). Further characterization revealed that a total of 465 genes, containing 179 upregulated and 286 downregulated genes, were assigned to IRGs (Fig. 1B, D).

Figure 1.

DEGs in colon cancer. Heatmap of 6420 DEGs (A) and 465 DE-IRGs (B) between normal and tumor tissues. Volcano plot of 6420 DEGs (C) and 465 DE-IRGs (D). Blue dots indicate upregulated genes, red dots indicate downregulated genes, and black dots mean genes without significant differences. DEG = differentially expressed gene, DE-IRG = differentially expressed immune-related gene, FDR = false discovery rate.

3.2. Identification of OS-associated IRGs in colon cancer

To identify possible prognostic IRGs, a univariate Cox analysis was conducted and 130 survival-associated IRGs were identified. Then, functional enrichment analyses were conducted and are shown in Table 1 and Figure 2. The results showed that the primary molecular function terms were “growth factor activity,” “cytokine activity,” and “receptor binding”; the primary biological process terms were “positive regulation of cell proliferation,” “innate immune response,” and “semaphorin-plexin signaling pathway”; and the primary cellular component terms were “extracellular space,” “extracellular region,” and “semaphorin receptor complex” (Table 1). In addition, KEGG pathway demonstrated that these OS-associated IRGs were involved in several processes of the tumor immune response, such as “cytokine-cytokine receptor interaction,” the “B-cell receptor signaling pathway,” and the “Ras signaling pathway” (Fig. 2).

Table 1

GO terms of OS-associated IRGs in colon cancer.

Ontology	ID	Description	P adjust	Count
Molecular function	GO:0008083	Growth factor activity	2.25E–13	16
	GO:0005125	Cytokine activity	2.89E–07	11
	GO:0005102	Receptor binding	6.85E–07	14
	GO:0008009	Chemokine activity	8.71E–07	7
	GO:0017154	Semaphorin receptor activity	9.04E–07	5
	GO:0046934	Phosphatidylinositol-4,5-bisphosphate 3-kinase activity	3.57E–06	7
	GO:0030215	Semaphorin receptor binding	1.53E–05	5
	GO:0045499	Chemorepellent activity	2.97E–05	5
	GO:0016814	Hydrolase activity, acting on carbon–nitrogen bonds, in cyclic amidines	4.63E–05	4
	GO:0008201	Heparin binding	1.00E–04	8
Biological process	GO:0008284	Positive regulation of cell proliferation	7.88E–14	24
	GO:0045087	Innate immune response	1.18E–11	21
	GO:0071526	Semaphorin–plexin signaling pathway	4.55E–11	9
	GO:0050853	B-cell receptor signaling pathway	9.91E–11	10
	GO:0030335	Positive regulation of cell migration	3.96E–10	14
	GO:0070374	Positive regulation of ERK1 and ERK2 cascade	2.73E–09	13
	GO:0001525	angiogenesis	4.15E–08	13
	GO:0006954	Inflammatory response	4.70E–08	16
	GO:0050919	Negative chemotaxis	1.04E–07	7
	GO:0043406	Positive regulation of MAP kinase activity	1.49E–07	8
Cellular component	GO:0005615	Extracellular space	3.54E–18	42
	GO:0005576	Extracellular region	7.83E–18	45
	GO:0002116	Semaphorin receptor complex	4.97E–07	5
	GO:0005886	Plasma membrane	6.99E–07	51
	GO:0009897	External side of plasma membrane	8.89E–06	10
	GO:0009986	Cell surface	4.14E–05	14
	GO:0072562	Blood microparticle	5.33E–05	8
	GO:0005887	Integral component of plasma membrane	2.15E–04	22
	GO:0042571	Immunoglobulin complex, circulating	2.26E–04	4
	GO:0031093	Platelet alpha granule lumen	4.14E–04	5

Figure 2.

KEGG pathways of OS-associated IRGs in colon cancer. HTLV = human T-cell leukemia virus type 1, IRG = immune-related gene, KEGG = Kyoto Encyclopedia of Genes and Genomes, MAPK = mitogen-activated protein kinase, NF-kappa B = nuclear factor-kappa B, OS = overall survival.

3.3. Characterization of hub IRGs

To develop a prognostic model based on OS-associated IRGs, hub IRGs that actively participated in the progression of colon cancer were further characterized (Fig. 3A). The results displayed that a total of 46 DE-IRGs were identified as hub IRGs, and their prognostic values are shown in Fig. 3B. Furthermore, the protein–protein interaction analyses of these hub genes showed that LEP, CXCL1, and CD19 are at the core of the interaction network (Fig. 3C).

Figure 3.

Characterization and analyses of hub IRGs in colon cancer. (A) The intersection of DE-IRGs and OS-associated IRGs. (B) Prognostic value of hub IRGs. (C) Protein–protein interaction of hub IRGs. DE-IRG = differentially expressed immune-related gene, IRG = immune-related gene, OS = overall survival.

Owing to their potential prognostic significance, the gene mutations and copy number variations of these hub IRGs were analyzed. As shown in Figure 4, gene mutations occurred at an approximately 30.58% rate. At the same time, missense mutations were found to be the most ordinarily occurring type in 29 hub IRGs and the PLCG2 gene had the highest mutation frequency. For gene copy number variation, the FABP4, ADIPOQ, and CCL28 genes were the most frequent amplifications, whereas the OXTR, JAG2, and UCN genes were the most frequent deletions.

Figure 4.

Mutation frequency of hub IRGs. IRG = immune-related gene.

3.4. Regulation of transcription factors on hub IRGs

Transcription factors are of crucial importance in the regulation of gene expression, and a regulatory network can be applied to elucidate the potential regulatory mechanisms of hub IRGs. The expression landscape of transcription factors was analyzed, and a total of 68 factors were differentially expressed in colon cancer (Fig. 5A). Then, a regulatory network between these transcription factors and hub IRGs was constructed, and their relationships are illustrated in Figure 5B. As a result, 17 transcription factors participated in the positive regulation of 14 IRGs, and most of these IRGs were unfavorable for the prognosis of colon cancer, except for BIRC5.

Figure 5.

Regulation of transcription factors on hub IRGs. (A) Transcription factors differentially expressed in colon cancer. (B) Regulatory network between transcription factors and hub IRGs. Triangles represent transcription factors. Green and red dots represent IRGs with favorable and poor prognosis for colon cancer, respectively. Red and green lines represent positive and negative regulation, respectively. IRG = immune-related gene.

3.5. Development of a IRG-based prognostic index

LASSO regression analysis confirmed that 14 hub IRGs could be selected for the construction of a prognostic model (Fig. 6). As shown in Figure 7, this model displayed a strong potential to predict the survival outcome of colon cancer patients. The prognostic accuracy was verified using time-dependent receiver operating characteristic curve and the area under the curve was 0.827, suggesting moderate capacity for survival prediction. According to the risk scores calculated by the prognostic model, 391 colon cancer patients could be well separated into low- and high-risk groups. The distribution of survival status in different groups, risk score curves, and the heatmap of the IRGs used for the construction of the prognostic index are illustrated in Figure 8. The formula based on the expression level of hub IRGs was as follows: (0.6784 × SLC10A2expr) + (0.0111 × FABP4expr) + (0.2968 × FGF2expr) + (–0.0806 × CCL28expr) + (0.0090 × IGKV1.6expr) + (0.0274 × IGKV1.8expr) + (0.2138 × EMS1expr) + (0.0406 × STC2expr) + (0.4381 × UCNexpr) + (0.1991 × UST2expr) + (0.0876 × VIPexpr) + (–4.7863 × GLP2Rexpr) + (0.1909 × IL1RL2expr) + (0.1034 × TRDCexpr).

Figure 6.

LASSO regression analysis of hub IRGs. IRG = immune-related gene, LASSO = least absolute shrinkage and selection operator.

Figure 7.

The prediction capability of the prognostic index. (A) The survival probability over time for the constructed model. (B) ROC curves of the constructed model. AUC = area under curve, ROC = receiver operating characteristic.

Figure 8.

The separating capacity of the prognostic index. (A) Survival status in low- and high-risk groups. (B) Rank of the prognostic index and distribution of different groups. (C) Heatmap of IRGs for the construction of the prognostic signature. IRG = immune-related gene.

3.6. Validation of the prognostic index

It is valuable to develop an independent predictor with clinical utility. Univariate and multivariate Cox regression analyses were performed to compare the prognostic value between risk score and other clinical indices, such as age, sex, and stage. The results indicated that the constructed prognostic index was superior to other clinical parameters and could act as an independent predictor for outcome prediction of colon cancer patients (Fig. 9).

Figure 9.

Comparison of the prognostic value between risk score and some clinical indices. (A) Univariate Cox regression analysis. (B) Multivariate Cox regression analysis.

Furthermore, the relationships between the prognostic model and 6 types of TIICs were analyzed. Accompanied by the increasing risk score, the abundance of macrophages, neutrophils, CD8+ T cells, dendritic cells, and CD4+ T cells was also increased (P < .05; Fig. 10), suggesting that the expression of these IRGs can reflect the immune status of the tumor microenvironment in colon cancer patients.

Figure 10.

The correlation between the risk score and TIICs. (A) B cells, (B) CD4+ T cells, (C) CD8+ T cells, (D) DCs, (E) macrophages, and (F) neutrophils. DC = dendritic cell, TIIC = tumor-infiltrating immune cell.

4. Discussion

The immune system plays a dual role in malignancies, which can launch an effective antitumor response or promote tumor progression and metastasis.[ The immune equilibrium is the middle phase between immunosurveillance and immune escape, during which tumor cells may produce variants and acquire the capacity to avoid immune elimination.[ Thus, it is necessary for cancer patients to rebuild immune equilibrium and maintain immune homeostasis through immunotherapy.[ During tumorigenesis and progression, an immunosuppressive tumor microenvironment is established, and many suppressive proteins and cytokines, including indoleamine-2,3-dioxygenase, programmed death-1, vascular endothelial growth factor (VEGF), interleukin-10 and transforming growth factor-β1, are produced by tumor and regulatory immune cells.[ These factors can lead to immune nonresponse or immune tolerance in cancer patients as deficiencies in antigen presentation and T-cell activation.[ Notably, TIICs are one of the important components in the tumor microenvironment, and their composition is highly associated with cancer prognosis, including colon cancer.[ With the significance of the immune system, many studies have focused on reinvigorating preexisting anticancer immune responses by rebuilding an immunostimulatory tumor microenvironment.[

At present, the TCGA and GEO databases provide numerous RNA-sequencing datasets and many computational methods have been developed for bioinformatics modeling and biomedical discovery.[ For more efficient usage of large datasets, deep learning has become the method of choice for many genomics modeling tasks, such as the prediction of the effects of genetic variation on gene regulatory mechanisms.[ IRGs have been indicated to be differentially expressed in the tumor microenvironment and to be associated with the immune status of cancer patients.[ Currently, some prognostic signatures based on single or multiple IRGs have been developed for the prediction of multiple cancers, including papillary thyroid cancer,[ gastric cancer,[ breast cancer,[ hepatocellular carcinoma,[ and laryngeal squamous cell carcinoma.[ At the same time, some researchers suggested that an immune-related prognostic model could be deployed for estimating the prognosis of colorectal cancer patients.[ To investigate the prognostic significance of IRGs in colon cancer, the DE-IRGs and survival-associated IRGs from 391 colon cancer datasets were screened in this study. After intersection analysis, a total of 46 hub IRGs were ascertained to be markedly correlated with the OS of colon cancer patients. Among these genes, a LASSO analysis further screened out 14 hub IRGs, including SLC10A2, CCL28, ESM1, STC2, IGKV1.6, IGKV1.8, UTS2, GLP2R, VIP, UCN, IL1RL2, FABP4, FGF2, and TRDC, which were suitable for the construction of a prognostic index. The area under the curve was 0.827 and the patients with high- and low-risk scores could be well distinguished. Thus, the constructed prognostic signature displayed a moderate prognostic capacity for colon cancer patients.

Functional enrichment analysis revealed that the cytokine–cytokine receptor interaction is the most significant KEGG pathway, which is similar to other studies in papillary thyroid cancer,[ gastric cancer,[ and colorectal cancer.[ Once again, it has been confirmed that cytokines are of crucial significance in tumor development and the immune response. Many cytokines produced in the tumor microenvironment can impair the function of immune cells and promote tumor progression.[ Correspondingly, a robust anticancer immune response requires the coordination of numerous stimulatory and inhibitory cytokines;[ thus, many immunotherapy efforts are focused on enhancing the efficacy in combination with agents that target cytokines and their receptors.[ Except for the “cytokine-cytokine receptor interaction” pathway, there are great differences in KEGG pathways in colon cancer when compared with the analysis of colorectal cancer.[ Compared with the 10- or 18-gene signature-based risk score for colorectal cancer,[ 14 hub IRGs were deployed for construction of prognostic index of colon cancer. Moreover, only the FABP4, UCN, and VIP genes were simultaneously involved in these prognostic models, further validating the differences between colon cancer and rectal cancer.

To explore the underlying molecular mechanisms, we analyzed the molecular characteristics of gene mutations and copy number variations and constructed a regulatory network between the differentially expressed transcription factors and hub IRGs. The results showed that approximately 30.58% of gene mutations occurred in the hub IRGs, and amplification and deletion events were induced. Simultaneously, 17 transcription factors participated in regulating the expression of IRGs. Moreover, most of these regulated IRGs are unfavorable for the prognosis of colon cancer. Therefore, these factors play a vital role in IRG expression and anticancer immunity.

As mentioned above, the compositions and fractions of TIICs were shown to be associated with the prognosis of multiple tumors.[ Thus, the relationships between the prognostic index and TIICs were analyzed. In colon cancer, the risk score had a positive correlation with some immune cells (Fig. 10), and this is similar to those in hepatocellular carcinoma,[ but contrary to some immune cells in papillary thyroid cancer.[ TIICs and IRGs in the tumor microenvironment differ in various tumors, and their actions are still being investigated. Our preliminary observations could provide a perspective for further research in the future. Certainly, there are some limitations to the prognostic index in providing guidance in clinical practice. The transcriptomics analyses can only determine some aspects of the immune status in the tumor microenvironment rather than global alterations.[ Moreover, further bioinformatic analyses, such as homology modeling and protein–protein docking, and experimental verifications in cell lines and clinical samples are required to confirm these findings.[

5. Conclusion

This study comprehensively analyzed the immunogenomic landscape of survival-associated IRGs and accomplished the goal of constructing an independent prognostic index for the survival prediction in colon cancer patients. In this study, the survival-associated IRGs were first screened and a functional enrichment analysis was conducted to elucidate their function in tumor immunity. After the intersection with DE-IRGs, the hub genes were screened for further analysis for their prognostic value, protein–protein interactions and genomic alterations. Next, 14 hub IRGs were chosen to build a prognostic model and further validation showed its feasibility as an independent predictor. Therefore, an independent prognostic signature is successfully constructed based on IRGs for the prognosis of colon cancer patients. Our results provide an alternative that could yield an effective prognosis and personalized immunotherapies for colon cancer.

Author contributions

Conceptualization: Zuquan Hu, Zhu Zeng, Shichao ZhangData curation: Yan Ouyang, Yun WangFormal analysis: Yan Ouyang, Jiangtao Huang, Yun WangMethodology: Jiangtao Huang, Fuzhou Tang, Shichao ZhangFunding acquisition: Yan Ouyang, Zuquan Hu, Zhu Zeng, Shichao ZhangProject administration: Zuquan Hu,Zhu Zeng, Shichao ZhangSupervision: Yun Wang, Fuzhou Tang

Validation: Jiangtao Huang, Yun Wang, Fuzhou Tang, Zuquan HuVisualization: Jiangtao Huang, Shichao ZhangWriting – original draft: Yan Ouyang, Jiangtao Huang

Writing – review & editing: Zuquan Hu, Zhu Zeng, Shichao Zhang