PDF

Int J Med Sci 2020; 17(13):1879-1896. doi:10.7150/ijms.45813 This issue Cite

Research Paper

A Prognostic Model Based on the Immune-related Genes in Colon Adenocarcinoma

Yuan-Lin Sun, Yang Zhang, Yu-Chen Guo, Zi-Hao Yang, Yue-Chao Xu

Department of Gastrointestinal Surgery, The First Hospital, Jilin University, Changchun 130021, Jilin Province, China.

✉ Corresponding author: Yue-Chao Xu, MD, PhD, Department of Gastrointestinal Surgery, The First Hospital, Jilin University, Changchun 130021, Jilin Province, China. Tel.: +86 13804332828; E-mail: ycxuedu.cn.

Abstract

Background: Immune-related genes (IRGs) are critically involved in the tumor microenvironment (TME) of colon adenocarcinoma (COAD). Here, the study was mainly designed to establish a prognostic model of IRGs to predict the survival of COAD patients.

Methods: The Cancer Genome Atlas (TCGA), Immunology Database and Analysis Portal (ImmPort) database, and Cistrome database were utilized for extracting data regarding the expression of immune gene- and tumor-related transcription factors (TFs), aimed at the identification of differentially expressed genes (DEGs), differentially expressed IRGs (DEIRGs), and differentially expressed TFs (DETFs). Univariate Cox regression analysis was subsequently performed for the acquisition of prognosis-related IRGs, followed by establishment of TF regulatory network for uncovering the possible molecular regulatory association in COAD. Subsequently, multivariate Cox regression analysis was conducted to further determine the role of prognosis-related IRGs for prognostic prediction in COAD. Finally, the feasibility of a prognostic model with immunocytes was explored by immunocyte infiltration analysis.

Results: A total of 2450 DEGs, 8 DETFs, and 79 DEIRGs were extracted from the corresponding databases. Univariate Cox regression analysis revealed 11 prognosis-related IRGs, followed by establishment of a regulatory network on prognosis-related IRGs at transcriptional levels. Functionally, IRG GLP2R was negatively modulated by TF MYH11, whereas IRG TDGF1 was positively modulated by TF TFAP2A. Multivariate Cox regression analysis was subsequently performed to establish a prognostic model on the basis of seven prognosis-related IRGs (GLP2R, ESM1, TDGF1, SLC10A2, INHBA, STC2, and CXCL1). Moreover, correlation analysis of immunocyte infiltration also revealed that the seven-IRG prognostic model was positively associated with five types of immunocytes (dendritic cell, macrophage, CD4 T cell, CD8 T cell, and neutrophil), which may directly reflect tumor immune state in COAD.

Conclusions: Our present findings indicate that the prognostic model based on prognosis-related IRGs plays a crucial role in the clinical supervision and prognostic prediction of COAD patients at both molecular and cellular levels.

Keywords: immune-related genes (IRGs), differential expressed analysis, prognostic model, Cox regression analysis, colon adenocarcinoma (COAD)

Introduction

Colorectal cancer (CRC) is one of the primary lethal malignancies which can be divided into colon cancer and rectal cancer based on the corresponding primary tumor sites [1]. The incidence cases of CRC are estimated to exceed 2.2 million and the death cases are predicted to exceed 1.1 million deaths by 2030 [2]. Colon adenocarcinoma (COAD) is considered as the most common pathological type of CRC [3]. The incidence and mortality of COAD accounts for 6.1% and 5.8%, respectively, ranking the fourth and fifth among all types of new cancer cases [4]. Although the diverse types of advanced therapies, including surgery, adjuvant radiation therapy or chemotherapy, and targeted molecular therapy, are currently used to treat CRC, the poor prognosis and survival of patients demand prompt solutions due to delayed diagnosis and adverse drug effects [5]. In consideration of the dilemma, it is urgent to explore new biomarkers, which might exert an impact on the prognosis of COAD.

At present, immunotherapy is wisely applied in the individualized treatment in a variety of tumors. Of note, interferons (IFNs) have long been demonstrated to play diverse roles, including antimicrobial and antiviral response, cell cycle progression, apoptosis and mediators of other cytokines [6-8]. More recently, James et al. have revolutionarily concluded that blocking CTLA-4 would prime T-cells to attack cancer cells [9], and Salem M et al. have revealed the removal of cancer cells by knocking out GARP of Treg cells [10]. Piero et al. have estimated that CDX2 could be recognized as a biomarker for malignant tumors in clinical surveillance and prognosis due to the fact that patients with stage II and stage III colon cancer and a lack of CDX2 suffered neoplasm recurrence and subsequent death [11]. These findings have illustrated the therapeutic importance of immune systems in COAD. Recently, Li et al. and Peng et al have expounded the prognostic value of immune-related genes (IRGs) to non-small cell lung cancer and papillary thyroid cancer, respectively [12, 13]. However, it remains unknown of the clinical correlation and prognostic evaluation of IRGs in COAD.

The present study was designed to explore the potential correlation between the clinical prognosis and immune-related genes (IRGs), which were molecular biomarkers that can be further applied to individualized and targeted therapies. In particular, Cox regression analysis was performed to construct an IRGs-based prognostic model. The visual regulatory network formed by DETFs and prognosis-related IRGs demonstrated the potential mechanisms at a molecular level. Additionally, immune infiltration analysis concerning seven-IRG prognostic model as well as immunocyte accumulation might shed new light on the function of immunocytes for prognostic prediction in COAD.

Methods

Patient samples and data collection

The transcriptome RNA-sequencing data were retrieved and downloaded from The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/); a total of 473 tumor tissues and 41 healthy tissues were included. The data on IRGs were downloaded from the Immunology Database and Analysis Portal (ImmPort) (https://www.immport.org/) [14]. Moreover, the cancer transcriptional genes data were available in Cistrome Project (http://www.cistrome.org/).

Functional enrichment analysis for DEIRGs

To integrate and analyze the IRGs among transcriptome RNA-sequencing data and IRGs, the cancer transcriptional genes data, and the “limma package” in R software [15] was used for data management. Subsequently, DEIRGs were screened, and the cut-off value assigned for false discovery rate (FDR) <0.05 and |log₂ fold change| (|logFC|) >1. DEIRGs were used to elucidate the underlying molecular mechanisms based on the findings from Database for Annotation, Visualization and Integrated Discovery (DAVID), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. For GO analysis conducted using “GO plot packages” (https://cran.rproject.org/web/packages/GOplot/citation.html) [16], GOCircle and GOChord graphs were obtained on the basis of a precondition and by setting a P <0.05 to indicate statistical significance.

Prognosis-related IRGs and DETFs

After incorporating overall survival (OS) and DEIRGs of patients using Perl language, we performed univariate Cox regression analysis for the preliminary screening of prognosis-related IRGs based on a P <0.05. Afterwards, prognosis-related IRGs were categorized into high-risk prognosis-related IRGs (HR>1) and low-risk IRGs (HR<1) based on hazard ratios (HRs) derived from univariate analysis. TFs bound with specific DNA sequences and are thereby critically involved in the direct regulation of gene expression. In this study, relevant TFs were acquired from Cistrome database (http://www.cistrome.org/), which is an online tool, and from TCGA by integrating chromatin profiling data and cancer genomics data, respectively, which in turn facilitated the exploration of regulatory correlations of TFs with cancer transcriptomes [17]. For further investigating the association of DETFs, which were extracted from the Cistrome database, with prognosis-related IRGs obtained on the basis of univariate Cox regression analysis, the regulatory network was constructed [18, 19] on the basis of a co-expression analysis according to Cor filter >0.04 and a P value filter <0.001 and the TF regulatory network was visualized using Cytoscape 3.7.2.

Prognostic model construction

Multivariate Cox regression analysis was further conducted on prognosis-related IRGs that were obtained based on univariate analysis, aiming at establishing a prognostic model and subsequently calculating the riskScore. Moreover, COAD patients were divided into two groups based on their riskScore, namely the low-risk and high-risk group. Kaplan-Meier (K-M) survival analysis with survival R packages could be used for obtaining the survivorship curve, which revealed the difference in terms of survival between the low-risk and high-risk group by demonstrating their respective survival periods and their corresponding survival rates. Furthermore, the risk plot R package was used to plot the risk curve to reveal the difference in survival periods between the low-risk and high-risk group. And the receiver operating characteristic (ROC) curve along with the area under the curve (AUC) could be obtained by analyzing the prognostic value of prognostic model.

 Figure 1 

Flow diagram of the study.

The correlation between the prognostic model and immunocytes

Tumor Immune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/), a publicly available resource for analyzing and visualizing the abundance of tumor-infiltrating immune cells, includes 10,897 samples encompassing 32 types of cancer identified from TCGA and can be used for estimating the abundance of six TIIC subsets (CD4 T cells, CD8 T cells, B cells, neutrophils, macrophages, and dendritic cells) [20]. In this study, we carried out immunocyte infiltration analysis with the Immune Estimation file retrieved from TIMER to explore the underlying correlation between the prognostic model and immunocytes.

Further verification on IRGs in the prognostic model

Oncomine (http://www.oncomine.org) [21] and TIMER (https://cistrome.shinyapps.io/timer/) [22] were both employed for further verifying differences in the expression of IRGs between tumor tissues and healthy tissues in the prognostic model. Moreover, Human Protein Atlas (HPA) (https://www.proteinatlas.org/) [23] was used to extract immunochemical images, which showed the distribution and protein expression of IRGs in the prognostic model.

Results

Identification of DEGs, DEIRGs, and DETFs

In this study, procedures are shown in the Figure 1 wherein 514 cases were collected, including 473 cases of COAD tissues and 41 cases of non-COAD tissues. The clinical features of the subjects are displayed in Table 1. A total of 2450 DEGs (Table S1), 79 DEIRGs (Table 2), and 8 DETFs (Table 3) were identified from the genes extracted from COAD tissue and non-COAD tissue samples based on the cut-off values (|LogFC|>1, FDR <0.05). Among these genes, there were 1765 upregulated DEGs, 685 downregulated DEGs (Figure 2A, B), 48 upregulated DEIRGs, 31 downregulated DEIRGs (Figure 2C, D), 6 upregulated DETFs, and 2 downregulated DETFs (Figure 2E, F).

 Figure 2 

DEGs, DEIRGs, and DETFs. (A) Volcano plot revealing clusters of upregulated and downregulated DEGs. (B) The distinction between DEG expression in tumor and normal tissues revealed by a heatmap. (C) Volcano plot demonstrated clusters of upregulated and downregulated DEIRGs. (D) Heatmap showing the distinction between the expression of DEIRGs in tumor and normal tissues. (E) Volcano plot showing clusters of upregulated and downregulated DETFs. (F) Discrimination between DETFs expression in tumor and normal tissues revealed by a heatmap.

GO enrichment and KEGG pathway analysis

To validate the biological characteristics of DEIRGs, functional relationship analyses were performed, including GO and KEGG pathway analysis. DEIRGs were significantly enriched in five GOs (P <0.05), with the most enriched term, “GO: 0005576 extracellular region” in CC category (Figure 3A, B). The visual presentation of the association of 79 DEIRGs and relevant GO terms is shown in Figure 3C. Figure 3D showed the nine KEGG pathways enriched DEIRGs with statistical significance (P <0.05). Moreover, it was apparent from Figure 3E that the nine KEGG pathways sorted by the number of enriched DEIRGs were in order as follows: Cytokine-cytokine receptor interaction, Neuroactive ligand-receptor interaction, Rheumatoid arthritis, IL-17 signaling pathway, cAMP signaling pathway, Amoebiasis, Natural killer (NK) cell mediated cytotoxicity, Salmonella infection and Renin-angiotensin system (RAS). The nine statistically significant signaling pathways in the KEGG database (P <0.05) are shown in Table 4, and a network (Figure 3F) used for visualizing the interaction between signaling pathways and DEIRGs (P <0.05) was constructed. Based on the network visualization, it was evident that hsa04060 (cytokine-cytokine receptor interaction) was frequently used to validate the KEGG pathway.

TFs-IRGs regulatory network

Univariate Cox regression analysis was performed on DEIRGs for screening prognosis-related IRGs using “survival package” of R software (P<0.05). Consequently, 11 types of genes (Figure 4A) were identified, including six high-risk IRGs as well as five low-risk IRGs. Co-expression analysis (Cor filter>0.4 and P value filter<0.001) was performed by incorporating seven prognosis-related IRGs and DETFs, which revealed two types of DETFs (TFAP2A and MYH11) and two types of low-risk prognosis-related IRGs (GLP2R and TDGF1) for regulatory network construction (Figure 4B). IRG GLP2R was negatively modulated by TF MYH11, whereas IRG TDGF1 was positively modulated by TF TFAP2A.

Seven-IRG prognostic model

Multivariate Cox regression analysis further screened seven prognosis-related IRGs, including four high-risk IRGs and three low-risk IRGs, out of the 11 prognosis-related IRGs identified from univariate Cox regression analysis. The expression of risk coefficients was combined to calculate the riskScore, which was supportive of the prognostic model construction. This riskScore was calculated as the sum of the expression quantities of selected IRGs when multiplied with their corresponding coefficients, as represented by the following formula: riskScore = the expression quantity of SLC10A2*(0.9319) + the expression quantity of CXCL1*(-0.2474) + the expression quantity of ESM1*(0.4490) + the expression quantity of INHBA*(0.2186) + the expression quantity of GLP2R*(-2.1451) + the expression quantity of STC2* (0.1822) + the expression quantity of TDGF1* (-0.2097). The above results of multivariate Cox regression analysis are shown in Table 5. Patients were categorized into two groups, namely the high-risk group and the low-risk group, using median riskScore (0.984) as a cut-off value. The Kaplan-Meier (KM) survival curves of both high-risk and low-risk group are displayed in Figure 5A, in which the red curve represents the high-risk group (N=195) and the blue curve represents the low-risk group (N=196). And it was clear that the survival rate of patients in the high-risk group was significantly lower than that in the low-risk group (P =2.567e-04). The AUC of ROC was approximately 0.715 (Figure 5B), and the IRG-based prognostic model was relatively accurate. RiskScore curve, reflecting patient distribution in both high-risk and low-risk groups, are displayed in Figure 5C. Similarly, survival status plot, illustrating the survival status of patients in both high-risk and low-risk group, is shown in Figure 5D. A heatmap revealing the expression of seven prognosis-related IRGs is shown in Figure 5E.

Independent prognostic analysis

Multivariate Cox regression analysis revealed that seven types of prognosis-related IRGs correlated with the prognosis of COAD patients. Univariate independent prognostic analysis (Figure 6A) and multivariate independent prognostic analysis (Figure 6B) revealed that age, T staging, and riskScore were significantly independent prognostic factors (Table 6) (P <0.05).

 Figure 3 

The functional enrichment analysis of DEIRGs. (A) and (B) The outer area of GOchord indicated the expression of DEIRGs, where the maximal aggregation of red dots (upregulated DEIRGs) and blue dots (downregulated DEIRGs) was located on the GO term: extracellular region (GO: 0005576). (C) The GOcircle revealed the relationship of 79 DEIRGs with the corresponding GO terms. (D, E) The bar plot and bubble plot showed the nine KEGG pathway-enriched DEIRGs. (F) KEGG pathways and the corresponding DEIRGs. The green rectangles represented the KEGG pathways, the red ellipses indicated the upregulated DEIRGs, and the blue ellipses indicate the downregulated DEIRGs.

 Figure 4 

Univariate Cox regression analysis and TF regulatory network. (A) Forest plot demonstrated the risk classification of 11 prognosis-related IRGs, including six high-risk prognosis-related IRGs (HR>1) and five low-risk prognosis-related IRGs (HR<1). The red squares and green squares indicate the high-risk and low-risk prognosis-related IRGs, respectively. (B) The network revealed the relationship of DETFs with prognosis-related IRGs. Purple rhombuses represent DETFs, blue circles represent prognosis-related IRGs, and the red and green lines represent positive and negative correlation, respectively.

Seven-IRG prognostic model and clinical features

Some clinical features of COAD, including age, sex, tumor stage, T staging, N staging, and M staging, were assessed for their probable correlation with prognosis-related IRGs via “beeswarm R” packages of R software, as shown in Table 7 (P <0.05). The median values of the selected gene expression were used as cut-off values, and the amounts of medians were found to be directly proportional to the specific clinical features. The median values in T1-2 staging were lower than those in T3-4 staging among ESM1 expression, INHBA expression, STC2 expression and riskScore (Figure 7E,F,O and I). In terms of the N staging, the median values of CXCL1 expression were relatively higher in N0 staging than in N1-2 staging (Figure 7A). Meanwhile, the median values of INHBA expression and riskScore were relatively higher in N1-2 staging than in N0 staging (Figure 7G, J). Moreover, in terms of M staging, the median values of CXCL1 expression were relatively higher in M0 staging than in M1 staging (Figure 7B), and for the tumor stage, the median values of INHBA expression and riskScore were significantly elevated in stage Ⅲ-Ⅳ than in stage I-Ⅱ (Figure 7H, K). Meanwhile, the median values of CXCL1 expression were higher in stage I-Ⅱ than in stage Ⅲ-Ⅳ (Figure 7C). The patient's age also played an important role, and the median values of CXCL1 expression were significantly higher in patients aged >65 years than in patients aged <65 years (Figure 7D); the trend was exactly the opposite for riskScore (Figure 7L). Additionally, the median values of SLC10A2 expression were also statistically different between T staging and M staging (P <0.05) (Figure 7M, N).

 Figure 5 

The prognostic model based on IRGs in COAD. (A) OS curves for patients in high-risk (red curve) and low-risk group (blue curve). (B) ROC curve suggesting the feasibility of our prognostic model. (C) Patients of high-risk (red dots) and low-risk group (green dots), and the distribution of their corresponding riskScore. (D) Patients in high-risk (red dots) and low-risk group (green dots), and their corresponding survival status. (E) Discrimination of the expression of seven prognosis-related IRGs between high-risk and low-risk group, as revealed by a heatmap.

 Figure 6 

Univariate and multivariate independent prognostic analysis. (A, B) Forest plots of univariate and multivariate independent prognostic analysis.

Immunocyte infiltration analysis

Figure 8 revealed a positive correlation of the riskScore of the seven-IRG prognostic model with immunocytes, including dendritic cells, neutrophils, CD8 T cells, macrophages, and CD4 T cells; macrophages were found to have the most significant relationship with riskScore (Cor=0.306, P=6.264e-10).

External verification of seven IRGs using online databases

In the seven-IRG prognostic model, five IRGs were upregulated and the remaining two IRGs were downregulated in COAD. Furthermore, Oncomine was used to externally validate the discrimination regarding the expression of seven IRGs between tumor and normal tissues. As displayed in Figure 9A-E and 10B-F, the expression of STC2, ESM1, INHBA, CXCL1 and TDGF1 was significantly increased in tumor tissues than in normal tissues, whereas GLP2R expression was higher in normal tissues than in tumor tissues (Figure 9F and 10G). Similarly, SLC10A2 expression was also increased in normal tissues than in tumor tissues (Figure 10A). However, we failed to extract corresponding information on SLC10A2 from Oncomine. The distribution and expression of SLC10A2, STC2, and GLP2R at protein level are shown in Figure 10G-L, whereas the distribution and expression of CXCL1, INHBA, ESM1, or TDGF1 remained inaccessible in HPA.

 Figure 7 

Comparison of the median expression of seven IRGs based on clinical features. (A-D) A comparison of the median expression of CXCL1 based on N staging, M staging, tumor stage, and age (≤65/>65 years) with statistical significance, respectively. (E) A comparison of the median expression of ESM1 in terms of T staging with statistical significance. (F-H) A comparison of the median expression of INHBA in terms of T staging, N staging, and tumor stage with statistical significance, respectively. (I-L) A comparison of the median riskScore in terms of T staging, N staging, tumor stage, and age (≤65/>65 years) with statistical significance, respectively. (M, N) A comparison of the median expression of SLC10A2 in terms of T staging and N staging with statistical significance, respectively. (O) A comparison of the median expression of STC2 in terms of T staging with statistical significance.

 Figure 8 

Correlation of the seven-IRG prognostic model with immunocyte infiltration. The relationship of the seven-IRG prognostic model with immunocytes, including (A) B cell, (B) CD4 T cell, (C) CD8 T cell, (D) Dendritic cell, (E) Macrophage, and; (F) Neutrophil, was revealed by scatter diagrams.

 Figure 9 

External verification of IRGs expression in tumor and normal tissues from Oncomine and HPA. (A-F) The expression of six IRGs in tumor and normal tissues with Oncomine, where columns subscript 1 and 2 indicated “tumor tissues” and “normal tissues”, respectively. (G-L) The comparison of protein expression of SLC10A2, STC2 and GLP2R between tumor and normal tissues.

Discussion

The researches on the predictive effects of biomarkers and their expression on cancer prognosis are a hotspot [24-27]. CTHRC1 has been estimated as a peritoneal metastasis-related gene for prognostic prediction in CRC [28, 29]. However, it is more dominant to study how the immune-related molecular mechanisms underlying the prognosis of COAD. The formation of the inflammatory environment caused by the deficiency of p53 might become the accelerant of colorectal tumors [30], which may be a direct reference to the investigation, progression and prognosis of COAD related to the immune microenvironment. In addition, some studies have shown that better performance of immunoscore in evaluating high-risk recurrence and metastasis of CRC [31, 32]. In this study, we mainly aimed to construct the IRGs-related prognostic model, which were screened out based on immune microenvironment.

By analyzing signaling pathways as well as functional enrichment on DEIRGs, the nine KEGG pathways enriched DEIRGs were shown as follows: Cytokine-cytokine receptor interaction, Neuroactive ligand-receptor interaction, Rheumatoid arthritis, IL-17 signaling pathway, cAMP signaling pathway, Amoebiasis, NK cell mediated cytotoxicity, Salmonella infection and RAS. The above nine KEGG pathways were related to inflammatory response, most of which are validated to have relationship with progression and therapies of colon cancer, including cAMP signaling pathway, Salmonella infection, IL-17 signaling pathway, cytokine-cytokine interaction receptor, NK cell mediated cytotoxicity and RAS. Cytokines along with correlated pathways may be involved in COAD progression in early stage [33]. Tseng et al. have shown that the level of IL-17 signaling was enhanced in CRC tissues than normal tissues [34]. Additionally, Chin et al. revealed that IL-17 could enhance the DNA binding capacity of NF-κB to stimulate CCR6 expression, potentially involved in the mechanisms of CRC migration [35]. cAMP signaling stimulation has been demonstrated to suppress tumor cell migration, including melanoma [36, 37], breast cancer [37] and colon cancer [37, 38] cells. Prospective cohort study has previously demonstrated that high activity of NK cell-mediated cytotoxicity is related to attenuated tumor risk [39], and NK cells have been revealed to be significantly decreased in primary colorectal lesion and liver metastases [40]. In a population-based study on patients with Salmonella infection, Mughini-Gras et al. have shown Salmonella infection as a risk factor for low-grade colon tumors [41]. To further investigate the mechanisms, Lu et al. have reported that β-catenin signals could be stimulated following Salmonella infection, thereby enhancing colonic tumorigenesis [42]. RAS is vitally involved in tumor angiogenesis and tumor cell growth [43, 44]. Nakamura et al. have revealed that a combination of RAS inhibitor ARB and anti-PD-L1 antibody showed synergistic anti-tumor effects [44]. Herein, in this study, the biological functions of DEIRGs were comprehensively investigated in COAD populations, which could possibly offer promising foundation to illustrate possible molecular regulatory mechanisms.

Multivariate Cox regression analysis revealed seven types of IRGs, including four types of high-risk IRGs (SLC10A2, STC2, ESM1 and INHBA) and three low-risk IRGs (CXCL1, TDGF1 and GLP2R), can be used to construct prognostic model of COAD, which also showed favorable feasibility for the result that AUC being 0.715. Wnt signaling pathway activation, a signal for CRC carcinogenesis, has recently been reported to cause decreased level of SLC10A2 [45]. STC2 has recently been reported as an independent prognostic biomarker, whose expression is related to colon cancer progression [46]. The overexpression of ESM1 in serum of CRC patients has been revealed Ji et al. [47]. INHBA expression has been previously reported to be increased in CRC tissues than normal tissues, and high INHBA expression might be used as an independent prognostic factor for lymph node involvement in CRC [48]. In terms of CXCL1, increased CXCL1 expression in colon cancer cells has recently been reported, and the carcinogenesis-promoting effect of CXCL2-CXCR2 axis is mediated by Gαi-2 and Gαq/11 [49]. Miyoshi et al. have demonstrated that CRC patients with high TDGF1 expression following surgery are significantly more likely to suffer palindromia and have poorer DFS in comparison with those with low expression [50]. GLP2R is regarded to be a pivotal gene type in CRC progression via the modulation of colonic epithelial integration [51]. Both univariate and multivariate independent prognostic analysis showed that age, T staging and riskScore were independent prognostic factors. In this study, seven types of specific IRGs were incorporated to explore the prognostic significance, compared with a previous analogous study. At present, the seven-IRG prognostic model could reflect their correlation with immunocyte infiltration according to the riskScore, which could be used as a type of evaluation indicator of immunologic status. Accumulative studies have demonstrated that immune system components, such as CD8+ and CD45RO+ memory T cells with specific cytokine signatures and possibly B cells, might be prognostic biomarkers, which is associated with tumor evolution with different presence [52-56]. Tumor-associated macrophages (TAMs) infiltration has been reported to be associated with OS of patients with prostate cancer and breast cancer [57, 58]. Here, notably, riskScore of prognostic model is significantly associated with macrophages, expanding the research vision concerning the relationship of immunocyte with progression and prognosis of COAD. In recent years, the application of nanomedicine based on immune microenvironment, has been identified as a type of novel and individualized anti-cancer therapy. Amphiphilic nanoparticles have been developed as an adjuvant agent, which can be combined with EphA2-derived peptide to assist the immune system to fight against liver metastasis in CRC [59]. Kuai et al. have demonstrated that nanodiscs-based phospholipids and apolipoprotein A1 (ApoA1)-mimetic peptides could be used for antigen presentation on dendritic cells to further elevate antitumor immunity [60], and the prognostic model at the center of IRGs might be used as a favorable evaluation of therapeutic efficiency.

 Figure 10 

External verification of IRGs expression based on TIMER. (A-G) The expression of seven IRGs in tumor and normal tissue, where red box plots and blue box plots suggested tumor tissues and normal tissues, respectively. “***” indicated P <0.001.

A previous research has shown a TFs-miRNA-targeted genes network for exploration of COAD progression [61]. To further integrally explore the possible molecular regulatory mechanisms, an IRGs-TFs network was established to elucidate the correlation between MYH11 and TFAP2A, TFs screened out, with IRGs (GLP2R and TDGF1). Nevertheless, studies concerning the regulatory mechanisms underlying TFs and IRGs in COAD are still in progress. Relevant studies have shown that MYH11 can foster the pathogenesis of leukemia and NSCLC [62, 63]. Zhang et al. have demonstrated that TFAP2A can enhance anlotinib resistance via the promotion of tumor-triggered angiogenesis in lung tumor cells [64]. In our research, IRG GLP2R was positively modulated by TF MYH11, while IRG TDGF1 was negatively modulated by TF TFAP2A. Of note, the TFs-IRGs network could be used to novelly present the prospective molecular mechanisms underlying the tumorigenesis and progression of COAD.

Strengths and limitations

In this research, IRGs were synthesized for constructing a prognostic model to predict prognosis in COAD patients using bioinformatics. Survival analysis as well as independent prognostic analysis confirmed that the prognostic model was clinically feasible. Additionally, in-depth molecular regulatory correlation analysis, including immunocyte infiltration analysis as well as TFs-IRGs regulatory network, broadened the horizon on researches concerning COAD progression. Multi-database (Oncomine, TIMER and HPA) analysis was utilized to verify the seven prognosis-related IRGs on the basis of RNA expression and protein level. Nevertheless, this study had certain limitations. First, the applicability of the prognostic model needs to be validated for a larger sample population in future studies. Second, we will continue to advance molecular mechanism experiments on the role of IRGs in COAD progression.

Conclusions

Collectively, DEGs, DEIRGs, and DETFs were retrieved from TCGA, Immport database, and Cistrome database. Univariate and multivariate Cox regression analysis of DEIRGs facilitated the construction of a seven-IRG prognostic model (GLP2R, INHBA CXCL1, STC2, SLC10A2, TDGF1, and ESM1), which could be reliably used for prognostic prediction in COAD patients. Furthermore, both protein and RNA expression of seven IRGs were verified in tumor and normal tissues from Oncomine, TIMER and HPA. In addition, the transcriptional regulatory network and underlying exploration of the association of seven-IRG prognostic model with immunocyte infiltration could potentially uncover new biomolecular interactions in COAD progression, which could be utilized as a potential biomarker for personalized treatment in COAD patients.

Abbreviations

COAD: Colon adenocarcinoma; DEGs: Differentially expressed genes; DEIRGs: Differentially expressed immune-related genes; DETFs: Differentially expressed transcription factors; TCGA: The Cancer Genome Atlas Project; IMMPORT: Immunology Database and Analysis Portal; TIMER: Tumor Immune Estimation Resource; GO: Gene ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DAVID: The Database for Annotation, Visualization and Integrated Discovery; AUC: Area under the curve; ROC: Receiver operating characteristic; FDR: False discovery rate.

Supplementary Material

Supplementary table S1.

Acknowledgements

Authors' contributions

Xu YC and Sun YL designed the study; Zhang Y and Yang ZH collected the data; Zhang Y and Guo YC analyzed the data; Sun YL drafted the manuscript; Xu YC revised the manuscript. All authors read and approved the final manuscript.

Ethics Committee Approval and Patient Consent

Consent for participation from all patients was acquired from public database.

Competing Interests

The authors have declared that no competing interest exists.

References

1. Fidler MM, Soerjomataram I, Bray F. A global view on cancer incidence and national levels of the human development index. Int J Cancer. 2016;139:2436-46

2. Arnold M, Sierra MS, Laversanne M, Soerjomataram I, Jemal A, Bray F. Global patterns and trends in colorectal cancer incidence and mortality. Gut. 2017;66:683-91

3. Mutch MG. Molecular profiling and risk stratification of adenocarcinoma of the colon. J Surg Oncol. 2007;96:693-703

4. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin. 2018;68:394-424

5. Muro K. Systemic chemotherapy for metastatic colorectal cancer -Japanese Society for Cancer of the Colon and Rectum (JSCCR) Guidelines 2016 for treatment of colorectal cancer. Nihon Shokakibyo Gakkai Zasshi. 2017;114:1217-23

6. Harada H, Taniguchi T, Tanaka N. The role of interferon regulatory factors in the interferon system and cell growth control. Biochimie. 1998;80:641-50

7. Gough DJ, Levy DE, Johnstone RW, Clarke CJ. IFNgamma signaling-does it mean JAK-STAT? Cytokine Growth Factor Rev. 2008;19:383-94

8. Slattery ML, Lundgreen A, Bondurant KL, Wolff RK. Interferon-signaling pathway: associations with colon and rectal cancer risk and subsequent survival. Carcinogenesis. 2011;32:1660-7

9. Barber DL, Wherry EJ, Masopust D, Zhu B, Allison JP, Sharpe AH. et al. Restoring function in exhausted CD8 T cells during chronic viral infection. Nature. 2006;439:682-7

10. Salem M, Wallace C, Velegraki M, Li A, Ansa-Addo E, Metelli A. et al. GARP Dampens Cancer Immunity by Sustaining Function and Accumulation of Regulatory T Cells in the Colon. Cancer Res. 2019;79:1178-90

11. Dalerba P, Sahoo D, Paik S, Guo X, Yothers G, Song N. et al. CDX2 as a Prognostic Biomarker in Stage II and Stage III Colon Cancer. N Engl J Med. 2016;374:211-22

12. Li B, Cui Y, Diehn M, Li R. Development and Validation of an Individualized Immune Prognostic Signature in Early-Stage Nonsquamous Non-Small Cell Lung Cancer. JAMA Oncol. 2017;3:1529-37

13. Lin P, Guo YN, Shi L, Li XJ, Yang H, He Y. et al. Development of a prognostic index based on an immunogenomic landscape analysis of papillary thyroid cancer. Aging (Albany NY). 2019;11:480-500

14. Bhattacharya S, Dunn P, Thomas CG, Smith B, Schaefer H, Chen J. et al. ImmPort, toward repurposing of open access immunological assay data for translational and clinical research. Sci Data. 2018;5:180015

15. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W. et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43:e47

16. Walter W, Sanchez-Cabo F, Ricote M. GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics. 2015;31:2912-4

17. Mei S, Meyer CA, Zheng R, Qin Q, Wu Q, Jiang P. et al. Cistrome Cancer: A Web Resource for Integrative Gene Regulation Modeling in Cancer. Cancer Res. 2017;77:e19-e22

18. Liang J, Cui Y, Meng Y, Li X, Wang X, Liu W. et al. Integrated analysis of transcription factors and targets co-expression profiles reveals reduced correlation between transcription factors and target genes in cancer. Funct Integr Genomics. 2019;19:191-204

19. Lin P, Guo Y-N, Shi L, Li X-J, Yang H, He Y. et al. Development of a prognostic index based on an immunogenomic landscape analysis of papillary thyroid cancer. Aging. 2019;11:480-500

20. Li T, Fan J, Wang B, Traugh N, Chen Q, Liu JS. et al. TIMER: A Web Server for Comprehensive Analysis of Tumor-Infiltrating Immune Cells. Cancer Res. 2017;77:e108-e10

21. Rhodes DR, Kalyana-Sundaram S, Mahavisno V, Varambally R, Yu J, Briggs BB. et al. Oncomine 3.0: genes, pathways, and networks in a collection of 18,000 cancer gene expression profiles. Neoplasia. 2007;9:166-80

22. Liu G-M, Xie W-X, Zhang C-Y, Xu J-W. Identification of a four-gene metabolic signature predicting overall survival for hepatocellular carcinoma. Journal of cellular physiology. 2020;235:1624-36

23. Pontén F, Schwenk JM, Asplund A, Edqvist PHD. The Human Protein Atlas as a proteomic resource for biomarker discovery. J Intern Med. 2011;270:428-46

24. Rosenwald A, Wright G, Wiestner A, Chan WC, Connors JM, Campo E. et al. The proliferation gene expression signature is a quantitative integrator of oncogenic events that predicts survival in mantle cell lymphoma. Cancer cell. 2003;3:185-97

25. Sotiriou C, Wirapati P, Loi S, Harris A, Fox S, Smeds J. et al. Gene expression profiling in breast cancer: understanding the molecular basis of histologic grade to improve prognosis. J Natl Cancer Inst. 2006;98:262-72

26. Qu L, Wang ZL, Chen Q, Li YM, He HW, Hsieh JJ. et al. Prognostic Value of a Long Non-coding RNA Signature in Localized Clear Cell Renal Cell Carcinoma. Eur Urol. 2018;74:756-63

27. Whitfield ML, Sherlock G, Saldanha AJ, Murray JI, Ball CA, Alexander KE. et al. Identification of genes periodically expressed in the human cell cycle and their expression in tumors. Mol Biol Cell. 2002;13:1977-2000

28. Tan F, Liu F, Liu H, Hu Y, Liu D, Li G. CTHRC1 is associated with peritoneal carcinomatosis in colorectal cancer: a new predictor for prognosis. Med Oncol. 2013;30:473

29. Wan Q, Tang J, Han Y, Wang D. Co-expression modules construction by WGCNA and identify potential prognostic markers of uveal melanoma. Exp Eye Res. 2018;166:13-20

30. Schwitalla S, Ziegler PK, Horst D, Becker V, Kerle I, Begus-Nahrmann Y. et al. Loss of p53 in enterocytes generates an inflammatory microenvironment enabling invasion and lymph node metastasis of carcinogen-induced colorectal tumors. Cancer Cell. 2013;23:93-106

31. Mlecnik B, Bindea G, Angell HK, Maby P, Angelova M, Tougeron D. et al. Integrative Analyses of Colorectal Cancer Show Immunoscore Is a Stronger Predictor of Patient Survival Than Microsatellite Instability. Immunity. 2016;44:698-711

32. Mlecnik B, Bindea G, Kirilovsky A, Angell HK, Obenauf AC, Tosolini M. et al. The tumor microenvironment and Immunoscore are critical determinants of dissemination to distant metastasis. Sci Transl Med. 2016;8:327ra26

33. Saadi A, Shannon NB, Lao-Sirieix P, O'Donovan M, Walker E, Clemons NJ. et al. Stromal genes discriminate preinvasive from invasive disease, predict outcome, and highlight inflammatory pathways in digestive cancers. Proc Natl Acad Sci U S A. 2010;107:2177-82

34. Tseng J-Y, Yang C-Y, Liang S-C, Liu R-S, Yang S-H, Lin J-K. et al. Interleukin-17A modulates circulating tumor cells in tumor draining vein of colorectal cancers and affects metastases. Clin Cancer Res. 2014;20:2885-97

35. Chin C-C, Chen C-N, Kuo H-C, Shi C-S, Hsieh MC, Kuo Y-H. et al. Interleukin-17 induces CC chemokine receptor 6 expression and cell migration in colorectal cancer cells. Journal of cellular physiology. 2015;230:1430-7

36. Dua P, Gude RP. Pentoxifylline impedes migration in B16F10 melanoma by modulating Rho GTPase activity and actin organisation. European journal of cancer (Oxford, England: 1990). 2008;44:1587-95

37. Dong H, Claffey KP, Brocke S, Epstein PM. Inhibition of breast cancer cell migration by activation of cAMP signaling. Breast Cancer Res Treat. 2015;152:17-28

38. Murata K, Sudo T, Kameyama M, Fukuoka H, Muka M, Doki Y. et al. Cyclic AMP specific phosphodiesterase activity and colon cancer cell motility. Clin Exp Metastasis. 2000;18:599-604

39. Imai K, Matsuyama S, Miyake S, Suga K, Nakachi K. Natural cytotoxic activity of peripheral-blood lymphocytes and cancer incidence: an 11-year follow-up study of a general population. Lancet. 2000;356:1795-9

40. Halama N, Braun M, Kahlert C, Spille A, Quack C, Rahbari N. et al. Natural killer cells are scarce in colorectal carcinoma tissue despite high levels of chemokines and cytokines. Clin Cancer Res. 2011;17:678-89

41. Mughini-Gras L, Schaapveld M, Kramers J, Mooij S, Neefjes-Borst EA, Pelt Wv. et al. Increased colon cancer risk after severe Salmonella infection. PLoS ONE. 2018;13:e0189721

42. Lu R, Wu S, Zhang YG, Xia Y, Liu X, Zheng Y. et al. Enteric bacterial protein AvrA promotes colonic tumorigenesis and activates colonic beta-catenin signaling pathway. Oncogenesis. 2014;3:e105

43. Egami K, Murohara T, Shimada T, Sasaki K-I, Shintani S, Sugaya T. et al. Role of host angiotensin II type 1 receptor in tumor angiogenesis and growth. J Clin Invest. 2003;112:67-75

44. Nakamura K, Yaguchi T, Ohmura G, Kobayashi A, Kawamura N, Iwata T. et al. Involvement of local renin-angiotensin system in immunosuppression of tumor microenvironment. Cancer Sci. 2018;109:54-64

45. Aymeric L, Donnadieu F, Mulet C, du Merle L, Nigro G, Saffarian A. et al. Colorectal cancer specific conditions promote Streptococcus gallolyticus gut colonization. Proc Natl Acad Sci U S A. 2018;115:E283-E91

46. Zhang C, Chen S, Ma X, Yang Q, Su F, Shu X. et al. Upregulation of STC2 in colorectal cancer and its clinicopathological significance. Onco Targets Ther. 2019;12:1249-58

47. Ji NY, Kim Y-H, Jang YJ, Kang YH, Lee CI, Kim JW. et al. Identification of endothelial cell-specific molecule-1 as a potential serum marker for colorectal cancer. Cancer Sci. 2010;101:2248-53

48. Okano M, Yamamoto H, Ohkuma H, Kano Y, Kim H, Nishikawa S. et al. Significance of INHBA expression in human colorectal cancer. Oncol Rep. 2013;30:2903-8

49. Chen M-C, Baskaran R, Lee N-H, Hsu H-H, Ho T-J, Tu C-C. et al. CXCL2/CXCR2 axis induces cancer stem cell characteristics in CPT-11-resistant LoVo colon cancer cells via Gαi-2 and Gαq/11. Journal of cellular physiology. 2019;234:11822-34

50. Miyoshi N, Ishii H, Mimori K, Sekimoto M, Doki Y, Mori M. TDGF1 is a novel predictive marker for metachronous metastasis of colorectal cancer. International journal of oncology. 2010;36:563-8

51. Lu Y, Kweon SS, Tanikawa C, Jia WH, Xiang YB, Cai Q. et al. Large-Scale Genome-Wide Association Study of East Asians Identifies Loci Associated With Risk for Colorectal Cancer. Gastroenterology. 2019;156:1455-66

52. Jochems C, Schlom J. Tumor-infiltrating immune cells and prognosis: the potential link between conventional cancer therapy and immunity. Exp Biol Med (Maywood). 2011;236:567-79

53. Galon J, Costes A, Sanchez-Cabo F, Kirilovsky A, Mlecnik B, Lagorce-Pagès C. et al. Type, density, and location of immune cells within human colorectal tumors predict clinical outcome. Science. 2006;313:1960-4

54. Fridman WH, Pages F, Sautes-Fridman C, Galon J. The immune contexture in human tumours: impact on clinical outcome. Nat Rev Cancer. 2012;12:298-306

55. Schmidt M, Hellwig B, Hammad S, Othman A, Lohr M, Chen Z. et al. A comprehensive analysis of human gene expression profiles identifies stromal immunoglobulin κ C as a compatible prognostic marker in human solid tumors. Clin Cancer Res. 2012;18:2695-703

56. Nielsen JS, Sahota RA, Milne K, Kost SE, Nesslinger NJ, Watson PH. et al. CD20+ tumor-infiltrating lymphocytes have an atypical CD27- memory phenotype and together with CD8+ T cells promote favorable prognosis in ovarian cancer. Clin Cancer Res. 2012;18:3281-92

57. Tsutsui S, Yasuda K, Suzuki K, Tahara K, Higashi H, Era S. Macrophage infiltration and its prognostic implications in breast cancer: the relationship with VEGF expression and microvessel density. Oncology reports. 2005;14:425-31

58. Lissbrant IF, Stattin P, Wikstrom P, Damber JE, Egevad L, Bergh A. Tumor associated macrophages in human prostate cancer: relation to clinicopathological variables and survival. International journal of oncology. 2000;17:445-51

59. Yamaguchi S, Tatsumi T, Takehara T, Sasakawa A, Yamamoto M, Kohga K. et al. EphA2-derived peptide vaccine with amphiphilic poly(gamma-glutamic acid) nanoparticles elicits an anti-tumor effect against mouse liver tumor. Cancer Immunol Immunother. 2010;59:759-67

60. Kuai R, Ochyl LJ, Bahjat KS, Schwendeman A, Moon JJ. Designer vaccine nanodiscs for personalized cancer immunotherapy. Nat Mater. 2017;16:489-96

61. Wei S, Chen J, Huang Y, Sun Q, Wang H, Liang X. et al. Identification of hub genes and construction of transcriptional regulatory network for the progression of colon adenocarcinoma hub genes and TF regulatory network of colon adenocarcinoma. J Cell Physiol. 2020;235:2037-48

62. Clayton EA, Rishishwar L, Huang TC, Gulati S, Ban D, McDonald JF. et al. An atlas of transposable element-derived alternative splicing in cancer. Philos Trans R Soc Lond B Biol Sci. 2020;375:20190342

63. Ma Q, Xu Y, Liao H, Cai Y, Xu L, Xiao D. et al. Identification and validation of key genes associated with non-small-cell lung cancer. J Cell Physiol. 2019;234:22742-52

64. Zhang L-l, Lu J, Liu R-Q, Hu M-J, Zhao Y-M, Tan S. et al. Chromatin accessibility analysis reveals that TFAP2A promotes angiogenesis in acquired resistance to anlotinib in lung cancer cells. Acta Pharmacol Sin. 2020

Author contact

Corresponding author: Yue-Chao Xu, MD, PhD, Department of Gastrointestinal Surgery, The First Hospital, Jilin University, Changchun 130021, Jilin Province, China. Tel.: +86 13804332828; E-mail: ycxuedu.cn.

Received 2020-3-9
Accepted 2020-7-7
Published 2020-7-19

Citation styles

©2024 Ivyspring International Publisher. Terms of use