Identification prognosis-associated immune genes in colon adenocarcinoma

Colorectal cancer is the fourth of the new cases and deaths among 36 cancers in 2018 [1]. It is also the second most frequent cancer diagnosed in women and the third most in men [2]. Several clinical trials have shown that the enhancement of tumor immune response can lead to long-term clinical response and benefit. Many expression profiling studies have proved that there is a relationship between the immune genes’ (IGs) characteristics and the better prognosis or efficacy to therapy, either chemotherapy or immunotherapy in tumors [3,4]. It is also essential to understand the immune environment’s contribution and the molecular subtypes of different tumor indications, but it lacks adequate wisdom to promote a reasonable combination of immunotherapy [5,6].

Whole genome-wide gene expression clusters such as The Cancer Genome Atlas (TCGA) have been constructed to classify and identify genomic abnormalities in large populations worldwide to detect the effects of the genetic composition of tumors on clinical prognosis [7]. With the human genome map finished, gene expression profiling has been cumulatively accepted by clinical diagnostic criteria. Besides, tumor microenvironment (TME) has significantly influenced therapeutic response and clinical outcomes [8,9]. Transcription and genomic data stemmed from many tumor samples used to research the TME, and immune infiltration assessment illustrated molecular subtypes of acute myeloid leukemia, bladder cancer, and cervical cancer [10–12]. Several separate reports about immune infiltration and lncRNA were associated with colon adenocarcinoma (COAD) [13,14]. Immunoscore, the primary significant prognostic ability, and robustness in COAD were validated [15]. However, few studies on the combination of IG and immune infiltration in the prognosis model of COAD patients.

The present paper mainly attempted to construct the prognosis model based on prognosis-related IGs (PRIGs), explore the prognostic value, and inquiry immune infiltration in COAD, which provide a novel orientation to the potential targets. Besides, we validated the PRIGs signature in GSE39582, GSE12945, and GSE103479 cohorts. In summary, we identified a new essential marker (risk score) to predict COAD patients’ prognosis. Comparison with the previous research, our study provided a new risk model as an independent prognostic biomarker in risk stratification for COAD patients.

In the present study, both the transcriptome and clinical data of COAD were obtained from the TCGA database (Data Release 20.0, Release Date: 11 November 2019, https://tcga-data.nci.nih.gov/tcga/). The main filter criteria for our data were as follows: (1) The keywords for cases were colon [Primary Site], ‘TCGA [Program],’ ‘TCGA-COAD [Project],’ ‘Adenomas and Adenocarcinomas [Disease Type].’ (2) The keywords for files for RNA‐sequencing data were ‘Transcriptome Profiling [Data Category],’ ‘Gene Expression Quantification [Data Category],’ ‘RNA-Seq [Experimental Strategy],’ ‘HTSeq - FPKM [Workflow Type]’). (3) The keywords for files for clinical data were ‘Clinical [Data Category],’ ‘BCR XML [Data Format]’). The expression of mRNA was obtained from the RNA‐sequencing data. The data access (GSE39582, GSE12945, and GSE103479), as the validation data, was obtained from Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database. Gene chips’ raw data were normalized using the Robust Multi-Array Average (RMA) algorithm provided by ‘limma’ [16]. Perl and R-package ‘sva’ were used to merge microarray data and reduce heterogeneity among the GSE39582, GSE12945, and GSE103479 datasets [17]. We use R software (version 3.6.2) to extract and annotate the expression of immune genes in TCGA and GEO databases.

Data extraction and integration were performed by Perl software. Screening differentially expressed genes, IGs, transcription factors (TFs) was undertaken by the Wilcoxon Test [16], using R-package ‘limma’ ‘edgeR.’False-discovery rate (FDR) < 0.05 and Log2 | (fold change, FC) | > 1.5 were set as the cutoffs. Bidirectional hierarchical clustering was analyzed and heat map was drawn by R-package ‘pheatmap’ (https://cran.r-project.org/web/packages/pheatmap/). R-package ‘ggplot2’ was used to draw a volcano map. Correlation analysis between differentially expressed IGs (DEIGs) and TFs (DETFs) was performed using correlation tests based on R software (cor > 0.4, P<0.001). The protein–protein interaction (PPI) network was constructed by Cytoscape software (version 3.7.2) [18]. R-package ‘Survival’ was used to conduct univariate cox regression analysis (P<0.01). Forest map was drawn by the R language ‘ggforest function.’

We constructed an optimal prognosis model of PRIGs via the Akaike Information Criterion. Genes expression and survival analysis were evaluated by the Kaplan–Meier method and the log-rank test (P<0.05). The risk score of each patient was calculated by the following formula: $Riskscore=Σj=1nCoefj×Xj$⁠, where, Coefj denotes the coefficient and Xj indicates the expression levels of each IG [19]. The survival data were obtained from the clinical data and combined with the previously acquired expression profiling data. The median risk score was selected as a cutoff to divide COAD cohorts into high-risk and low-risk groups. The evaluation indicator of survival analysis was overall survival (OS), which refers to the period from the date of diagnosis until death from any cause. The survival curve was drawn according to the high- and low-risk value using R-package ‘survival’ and ‘survminer.’ R-package ‘survival ROC’ was used to draw the Receiver Operating Characteristic (ROC) curve for the difference and calculate the value of Area Under Curve (AUC). The risk curve, survival state diagram, and heat map were drawn based on the patient risk score. PRIGs were recognized by univariate and multivariate Cox regression analyses. The results were verified in the GSE39582, GSE12945, and GSE103479 cohorts.

We used the expression of PRIGs to draw a nomogram via the R-package ‘Hmisc,’ ‘lattice,’ ‘Formula,’ ‘ggplot2,’ ‘foreign,’ and ‘rms.’ Moreover, calibration curves were used to assess the consistency between actual and predicted survival rates. Furthermore, the consistency index (C-index) was calculated to evaluate the model prognosis capability. The values of 0.5 and 1.0 represent a random probability and an excellent performance for predicting survival with the model, respectively. The GSE39582, GSE12945, and GSE103479 cohorts were used to verify the risk score’s prognostic value and the nomogram.

Based on the R software, the functional enrichment analysis of PRIGs was performed to identify Gene Ontology (GO) categories of biological processes (BPs) and molecular functions (MFs). The Database for Annotation, Visualization and Integrated Discovery (DAVID) database was used to conduct enrichment analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [20]. FDR < 0.05 was invoked as the cutoff. R-package ‘GOplot’ was used to integrate expression data and functional enrichment analysis. Moreover, Gene Set Enrichment Analysis (GSEA) was also used to analyze the differences between the high-risk and low-risk groups. The number of permutations was set to 1000 and FDR < 0.25 was acknowledged as statistically significant [21].

TCGA-COAD cohort included 385 patients. Patients without survival time and less than 90 days were excluded, and 334 patients were obtained as the training group. GSE39582 dataset included 566 stage I–IV COAD patients, GSE12945 contained 62 stages I–IV colorectal cancer patients, and GSE103479 included 156 stage I–IV colorectal cancer patients. Using the same exclusion criteria of the training group, 587 colorectal cancer patients were selected from the GEO cohorts as the test group. The detailed clinical characteristics of patients are listed in Supplementary Table S1.

A total of 437 transcriptome data were obtained, 39 (8.9%) were derived from healthy samples and 398 cases (91.1%) came from tumor samples. We screened out 3833 differentially expressed genes (DEGs) included 2610 up-regulated and 1223 down-regulated (Figure 1A,D); 332 differentially expressed immune genes (DEIGs) contained 120 up-regulated and 212 down-regulated (Figure 1B,E). We also obtained 37 differentially expressed transcription factors (DETFs) including 25 up-regulated and 12 down-regulated (Figure 1C,F).

We filtered out 24 PRIGs by univariate Cox regression analysis. CXCL3 was low-risk PRIG (Hazard Ratio (HR) < 1, P<0.01), others were high-risk PRIGs (HR > 1, P<0.01); the detailed information is shown in Supplementary Figure S1 and Table S2. Twelve PRIGs were selected by multivariate regression analysis, including SLC10A2, CXCL3, NOX4, FABP4, ADIPOQ, IGKV1-33, IGLV6-57, INHBA, UCN, VIP, NGFR, and TRDC (Table 1). The prognosis model was constructed based on these 12 PRIGs. The coefficient of PRIGs is shown in Supplementary Figure S2. Each patient’s risk score was computed based on each gene’s expression level and risk coefficient. A risk score was used to predict prognosis.

The calculation formula of the risk score was as follows: Risk score = 0.1966 × expression of TRDC + (−0.3930) × expression of NGFR + 0.0473 × expression of VIP + 0.4528 × expression of UCN + 0.1227 × expression of INHBA + 0.0027 × expression of IGLV6-57 + 0.0625 × expression of IGKV1-33 + (−0.3076) × expression of ADIPOQ + 0.0692 × expression of FABP4 + (−1.1528) × expression of NOX4 + (−0.0249) × expression of CXCL3 + 0.6105 × expression of SLC10A2.

According to the median risk score, the patients were divided into high-risk and low-risk groups. A heat map visualized the gene expression profiles of the high- and low-risk groups. As the risk score escalated, the patients’ survival time decreased, and the death numbers increased (Figure 2A). In the training group, the Venn diagram and R software were used to select the validation gene by intersecting with the 24 prognostic-related IGs screened from the TCGA database, and finally 5 genes were verified (Figure 2B). OS was significantly lower in the high-risk subgroup than in the low-risk subgroup (Figure 2C, P=1.294e−05), the 1-, 3-, and 5-year AUC value of ROC was 0.792, 0.744, and 0.797, respectively (Figure 2D). Similarly, in the test group, OS was also apparently lower in the high-risk subgroup than in the low-risk subgroup (Figure 2E, P=2.178e−05), and the 1-, 3-, and 5-year AUC value of ROC was 0.625, 0.646, and 0.713, respectively (Figure 2F). Univariate and multivariate independent prognostic analyses manifested that the risk score was an independent prognostic predictor both in the training and test groups (Tables 2 and 3, Supplementary Figure S3).

Nomogram is a useful tool quantifying an individual’s risk in a clinical setting by integrating multiple risk factors [22,23]. We exploited the nomogram to predict 3- and 5-years OS probabilities by integrating the 12 PRIGs signature for TCGA-COAD and verified in the GEO datasets (Figure 3A). Each factor is assigned a score in proportion to its contribution to the survival risk. We calculated each patient’s total points by summarizing the number of points for all PRIGs and estimating survival rates via drafting a vertical line between each prognosis axis and the total point axis, which may help appropriate practitioners make clinical decisions for COAD patients. The calibration curve showed that the actual survival rate matched the predicted 3- and 5-years survival rate well; the C-index was 0.77 in the TCGA data (Figure 3B–D). The nomogram was verified in the GEO cohorts, the C-index was 0.72, and 3- and 5- year calibration curves were shown in Figure 3E–G.

The PPI networks between DIGs and DTFs were made up of ten DEIGs and four DETFs. All of them had a positive relationship (correlation coefficient = 0.4, P<0.001, Figure 4A). GO enrichment analysis of the 12 PRIGs manifested in Figure 4B,C, a total of 36 GO terms of BP, 26 GO terms of MF were identified to be significant (Supplementary Material S1, FDR < 0.05). These 12 genes were involved in several essential BPs, including hormone secretion, hormone transport, lipid transport, lipid localization, signal release, macrophage differentiation regulation, negative regulation of blood pressure, and negative regulation of secretion, glucose homeostasis, brown fat cell differentiation. The TOP-10 terms of BP and MF are shown in Figure 4C.

The KEGG pathway analysis showed that PRIGs are mainly enriched in the cytokine–cytokine receptor interaction, Peroxisome Proliferators-activated Receptor (PPAR) signaling pathway (P<0.05). GSEA analysis showed that changed genes were enriched in several common pathways. A total of 111/178 genesets were up-regulated in phenotype high-risk group, 66 gene sets were significant at FDR < 25%. The top-eight of the high-risk group are shown in Figure 4D, including hedgehog signaling pathway (Enrichment score, NES = 2.09, P=0.000), basal cell carcinoma (NES = 1.98, P=0.000), axon guidance (NES = 1.87, P=0.004), notch signaling pathway (NES = 1.86, P=0.010), Gonadotropin-releasing hormone (GnRH) signaling pathway (NES = 1.84, P=0.002), glycosaminoglycan biosynthesis chondroitin sulfate (NES = 1.83, P=0.006), extracellular matrix (ECM) receptor interaction (NES = 1.77, P=0.025), and dorsoventral axis formation (NES = 1.77, P=0.011).

The diagnostic strategy for evaluating the tumor immune infiltration has a comprehensive clinical application value. Although the overlapping biological features of tumor immunophenotypes have been described [3], gene composition lacks consensus, the cell specificity of gene expression signal is not apparent, and there is no systematic analysis of the gene’s prognostic value in various tumor types. Until recently, a comprehensive meta-analysis confirmed the prognostic role of the IG in cancer as a whole [24]. In the current work, we singled out 24 PRIGs and constructed a prognosis model based on 12 PRIGs, which was entirely accurate and had better predictive value. A short bioinformatics protocol flow chart listed in Figure 5.

First, univariate regression analysis demonstrated that 22 genes had a relationship with poor prognosis, especially SLC10A2, PLCG2, and PTH1R. Previous reports have mentioned that secondary bile acids increase the risk of COAD, and SLC10A2 plays a crucial role in intestinal bile acid reabsorption. This transporter is the primary mechanism for the uptake of intestinal bile acids by apical cells in the distal ileum. Bile acid malabsorption is related to epidermal growth factor receptor expression and M3 muscarinic [25,26]. PLCG2 is associated with immune responses and PLCG2 mutations have been found in the chronic lymphocytic leukemia patient [27,28]. However, its role in COAD’s occurrence and development is not very clear, and is worth further study. Researchers found that PTH1R was associated with both miRNA and DNA methylation in colorectal liver metastasis patients [29].

Second, whether the IG expression can be used as a prognostic marker is a vital research topic. In the present study, our prognosis model is considered to be valuable. The 1-, 3-, and 5-year AUC values of ROC were 0.792, 0.744, and 0.797, respectively, which was more accurate than previous reports [30,31]. The 12 IGs may be ideal prognostic markers. Five-year OS was ∼75% in the low-risk group and 50% in the high-risk group, consistent with the previous report [32]. High-risk genes in the model, including SLC10A2, NOX4, FABP4, ADIPOQ, IGKV1-33, IGLV6-57, INHBA, UCN, VIP, NGFR, and TRDC. NADPH oxidase complex can positively predict relapse in stage II left-side COAD, NOX4 encodes a member of the NOX-family of enzymes that functions as the catalytic subunit [33]. FABP4 plays a crucial role in regulating lipid metabolism and adipocyte differentiation associated with the tumor suppressor, PTEN [34]. FABP4 and ADIPOQ have been reported to increase the risk of colorectal neoplasia related to the COAD survival probability [35,36]. Dudgeon et al. found that UCN can inhibit COAD cell growth via inducing apoptosis by PUMA and a p53 target [37]. VIP promotes the initiation and progression of COAD by the generation of a local proinflammatory environment [38]. NGFR hypermethylation can be observed through COAD development [39]. The risk score is closely correlated with the COAD malignant clinicopathological characteristics and is an independent prognostic factor. The prediction ability of 12 immune genes is also cross-verified in three GEO data sets, which shows that our cohort analysis based on TCGA and GEO data has predictive value. Moreover, the roles of IGKV1-33, IGLV6-57, and TRDC in the development of COAD are not very clear, which deserve further study.

Subsequently, we constructed a nomogram to predict individual clinical outcomes. Nomogram is a stable and reliable tool to quantify individual risk by combining and describing risk factors, which has been used for tumor prognosis, including colorectal cancer [40,41]. By summarizing all points of view, the model provides individuals with numerical possibilities for clinical outcomes, such as OS, relapse, and drug noncompliance. Except for standard clinicopathological features, such as genetic markers, can also be included in predictive nomogram models to predict clinical outcomes [41,42]. The combination of autophagy gene features and prognostic factors has a better prognostic value than a single application [43]. Besides, the calibration curve showed that the nomogram accurately predicted 3- and 5-years survival rates. Consistently, in the current study, we illustrated that a nomogram, including the 12 PRIGs signature, could better predict the 3- and 5-years survival possibility of COAD patients. These results suggest that the 12-PRIGs prognostic model maybe have a convinced value in adjusting COAD patients’ treatment plans.

Furthermore, GO enrichment and KEGG pathway analysis show that PRIG participates in regulating several crucial processes and signaling pathways. ADIPOQ, INHBA, UCN, and VIP participate in hormone secretion and hormone transport BPs. SLC10A2, FABP4, ADIPOQ, and INHBA involved in lipid transport and lipid localization BPs [44–46]; INHBA, CXCL3, and NGFR participate in the regulation of the cytokine–cytokine receptor interaction pathway; FABP4 and ADIPOQ involved in the regulation of PPAR signaling pathway, which coincided with the previous report, PPAR signaling pathway might affect the recurrence of COAD patients [47,48]. The GSEA analysis further suggested that the pro-tumor inflammatory reaction might be more prominent in stage IIB/IIC of COAD [49]. The blockade of PD-L1 increased FABP4 expression in tissue-resident memory T (Trm) cells, promoting lipid uptake by Trm cells [50]. Previous reports indicate that the hedgehog signaling pathway might affect the recurrence of COAD patients [51]. Dandawate et al. reported that the Notch signaling pathway involved regulating colon cancer cells’ proliferation [51]. Our results are consistent with these reports.

From the above discussion, the conclusion can be reached that 12 IGs were more significantly associated with the COAD prognosis and the prognosis model was accurate. The prognostic model based on these genes is validated in three independent GEO cohorts. Further research might be necessary about these genes in the clinical application and may provide a new perspective into the COAD’s pathogenesis. Some of the antecedently overlooked genes maybe have a chance to become additional biomarkers for COAD. Finally, the present study enhanced our cognition of the sophistication of the occurrence and development of COAD and may detect novel therapeutic targets.

The authors certify that all the original data in this research could be obtained from public database. Both the transcriptome and clinical data of COAD were obtained from the TCGA database (https://tcga-data.nci.nih.gov/tcga/). The data of training group are available at NCBI GEO (GSE39582 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE39582, GSE12945 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12945, and GSE103479 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE103479). Other data that were used to support the findings of the present study are included within the supplementary information files. All the raw data in the present study are available from the first author or corresponding author upon request.

The authors declare that there are no competing interests associated with the manuscript.

The authors declare that there are no sources of funding to be acknowledged.

D.h.M. and Y.d.M. designed the study. X.p.M. downloaded the data. Y.d.M. and Y.Y. analyzed the data. Y.d.M. and J.t.W. drafted the article. All authors approved the paper.

We thank the TCGA and GEO databases for the availability of the data.

AUC

Area Under Curve

BP

Biological Process

COAD

Colon Adenocarcinoma

C-index

Consistency Index

DAVID

The Database for Annotation, Visualization and Integrated Discovery

DEG

Differentially Expressed Gene

DEIG

Differentially Expressed Immune Gene

DETF

Differentially Expressed Transcription Factor

ECM

Extracellular Matrix

GEO

Gene Expression Omnibus

GnRH

Gonadotropin-releasing Hormone

GO

Gene Ontology

GSEA

Gene Set Enrichment Analysis

HR

Hazard Ratio

IG

Immune Gene

KEGG

Kyoto Encyclopedia of Genes and Genomes

MF

Molecular Function

NES

Enrichment Score

OS

Overall Survival

PPAR

Peroxisome Proliferators-activated Receptor

PPI

Protein–Protein Interaction

PRIG

Prognosis-related IGs

RMA

Robust Multi-Array Average

ROC

Receiver Operating Characteristic

TCGA

The Cancer Genome Atlas

TF

Transcription Factor

TME

Tumor Microenvironment

TrmR

Resident memory T cell