Identification and validation of an immune-related gene signature predictive of overall survival in colon cancer

Construction of TME, screening of DEGs and functional enrichment analysis

Detailed clinical characteristics of the training and validation cohorts are listed in Supplementary Table 1. The ESTIMATE algorithm was applied to estimate immune and stromal scores. The association of immune and stromal scores with clinical characteristics showed that the immune scores were different in pathological stage groups, tumor M staging groups, and tumor location groups. The immune score of stage I was significantly higher than stage IV (P = 0.008). M0 was significantly higher than M1 (P = 0.002), and the right side was significantly higher than the left side (P < 0.001) (Figure 1A). White people and non-adenocarcinomas had higher stromal scores than other races and adenocarcinomas (Race: P = 0.001; Histology classification: P = 0.012) (Figure 1B). Other variables had no significant differences in the distribution of immune scores and stromal scores (Supplementary Figure 1A, 1B). Immune scores range from -899.82 to 2959.54, and stromal scores range from -2171.21 to 1943.63. Colon cancer patients were divided into high score group and low score group by the optimal cut-off value (Supplementary Figure 1C, 1D). According to the log-rank test results, the OS of the high immune score group was significantly higher than that of the low immune score group (P = 0.018) (Figure 1C). There was no significant difference in OS between the high and low stromal score groups (P = 0.390) (Figure 1D).

Figure 1. Association of stromal and immune scores with colon cancer clinical characteristics and prognosis in TCGA. (A) Significant differences in the distribution of immune scores among different tumor stage, metastasis, and tumor location groups. (B) Significant differences in the distribution of stromal scores among different race, tumor location, and histology classification groups. (C) Kaplan-Meier survival curves of high and low immune score groups. (D) Kaplan-Meier survival curves of high and low stromal score groups.

Immune- or stromal-related DEGs were identified by comparing the RNA-expression of colon cancer patients with high and low immune (or stromal) scores. A cluster analysis screened out immune-related DEGs with high score and low score groups, as displayed in the heat map (Figure 2A). A volcano plot revealed significantly differentially expressed genes (Figure 2C). A total of 1,076 genes were identified as immune-related DEGs, which contain 120 up-regulated genes and 956 down-regulated genes. The heat map and volcano plot of stromal related DEGs are shown in Figure 2B, 2D. 1,199 stromal-related genes were screened, including 153 up-regulated genes and 1,046 down-regulated genes. A Venn diagram displayed 844 intersecting immune- and stromal-related DEGs, including 55 up-regulated and 789 down-regulated genes (Figure 2E, 2F).

Figure 2. Comparison between gene expression profiles of immune and stromal scores in TCGA. (A, B) Heat maps showing expression profiles for immune score and stromal score-related DEGs. (C, D) Volcano plots showing up-regulated and down-regulated DEGs related to immune score and stromal score. (E, F) Venn diagrams showing the intersection of immune score and stromal score related up-regulated /down-regulated DEGs. (G) Histogram showing the top ten Gene Ontology terms in BP, CC, and MF. (H) Bubble chart exhibiting top fifteen KEGG analysis terms.

Top terms of GO analysis included immune response, inflammatory response, cell adhesion, and innate immune response in BP; plasma membrane, extracellular exosome, and extracellular region in CC; and calcium ion binding, receptor activity, and serine-type endopeptidase activity in MF (Figure 2G). The results of KEGG enrichment were also related to immune responses, including cytokine-cytokine receptor interaction, phagosome, chemokine signaling pathway (Figure 2H). Collectively, the results indicated that the enriched GO terms and KEGG pathways were mainly related to immune response.

Construction, verification, and subgroup analysis of IRG risk score model

Thirty genes related to the prognosis of colon cancer from 1096 DEGs were screened by univariate Cox analysis. The LASSO regression analysis model further identified 18 genes associated with OS (Supplementary Figure 2A, 2B). Six significant independent prognostic genes were selected by multivariate Cox regression analysis. Among them, AXIN2, CCL22, and CLEC10A were protective genes whose high expression was associated with higher survival probability. CRIP2, RUNX3, and TRPM5 were dangerous genes whose high expression was associated with a lower probability of survival (Supplementary Figure 3A, 3B). The prognostic gene risk score model = (-0.273 × expression value of AXIN2) + (-0.372 × expression value of CCL22) + (-0.299 × expression value of CLEC10A) + (0.344 × expression value of CRIP2) + (0.324 × expression value of RUNX3) + (0.341 × expression value of TRPM5) (Figure 3A). In accordance with the median risk score, all colon cancer patients from the TCGA and GEO cohorts were divided into high-risk groups and low-risk groups, respectively. The difference in OS between the two risk score groups was determined to be significant by a log-rank test (TCGA: P < 0.001, GEO: P = 0.019) (Figure 3B, 3E). The distribution of risk score, survival status, and gene expression among patients in the training set and validation set are given in Figure 3C, Figure 3F. The clustering heat maps showed that the expression of prognostic genes AXIN2, CCL22, and CLEC10A was up-regulated in the high-risk group, while the expression of CRIP2, RUNX3, and TRPM5 was down-regulated in the high-risk group (Figure 3C, 3F). The AUC for 1-, 3-, and 5-year OS were 0.795, 0.757, and 0.728 in the TCGA cohort, respectively (Figure 3D). The AUC values for 1-, 3-, and 5-year OS were 0.715, 0.685, and 0.666 in GEO cohort, respectively (Figure 3G).

Figure 3. Prognostic analysis and performance assessment of TCGA and GEO. (A) Forest map showing six signature genes identified by a multivariate Cox regression analysis. (B, E) Kaplan-Meier survival curves and log-rank tests of six prognosis genes in TCGA and GEO. (C, F) Distribution of risk score, survival status and gene expression among patients in TCGA and GEO. (D, G) AUC values for 1-, 3-, and 5-year OS in TCGA and GEO.

Low-risk scores were concentrated in other races, left site, and stage I/II; all of these subgroups were associated with higher OS. However, the white race, right site, and stage III/IV had an unfavorable prognosis and accumulated significant high-risk scores (Supplementary Figure 4). Since IRG risk score is highly correlated with the above clinical characteristics, we attempted to clarify whether the IRG signature has prognostic value independent of these clinical characteristics. Stratified survival analysis showed that the log-rank tests for the survival probability of both the high- and low-risk groups were significant in different subgroups (Figure 4). IRG risk score was an effective predictor of survival in subgroups of patients with different clinical characteristics.

Figure 4. Kaplan-Meier survival subgroup analysis according to IRG signature stratified by clinical characteristics. (A) Age <65 years and age ≥65 years. (B) Female and Male. (C) White race and other race. (D) Left site and right site. (E) Adenocarcinoma and other histological type. (F) Stage I/II and stage III/IV.

Verification of expression profiles and immunohistochemistry of six immune-related prognostic genes

In the GSE39582 dataset, we calculated the correlation between the immune-related independent prognostic genes and risk scores. The expression levels of AXIN2, CCL22, and CLEC10A were significantly higher in the low-risk group, and the expression levels of CRIP2, RUNX3, and TRPM5 were significantly positively correlated with the risk scores (Figure 5A). This is consistent with TCGA cohort analysis.

Figure 5. Differences in protein expression induced by six genes were verified in human tissue samples. (A) Pearson correlation of expression and risk score in GEO. (B) Representative immunohistochemical staining images of four genes in normal colon tissue and colon cancer specimen.

The HPA database was used to validate the protein expression levels of six genes. CCL22 and TRPM5 were not available on the website, so we compared the other four genes’ IHC in normal tissue and tumor tissue. We found that AXIN2 staining was medium in normal tissues but low in cancer tissues, and the number of normal tissues was higher than that of cancer tissues. RUNX3 was weakly positive in tumor tissue but negative in normal tissue. CLEC10A and CRIP2 were negative in both normal and tumor tissues (Figure 5B).

Correlations between somatic mutation, immune cell infiltration, immune checkpoints expression, and IRG signature

By analyzing MuTect2 mutation annotation files, the total TMB and mutation distribution from the TCGA cohort (Supplementary Figure 5) were obtained. All patients with somatic mutation information were divided into either a high-risk group (n = 177) or a low-risk group (n = 179) according to the above grouping rules. The frequency of mutation of the first 20 genes in the two groups was mostly similar (Figure 6A, 6B). The two groups of mutant genes with significantly different mutation frequencies are shown in the Figure 6C. The frequencies of all mutated genes were higher in the high-risk group than in the low-risk group, indicating that the frequency of somatic mutation was positively correlated with our IRG risk score (Figure 6C). The relationship between the IRG risk score and the first eight mutation pathways is shown in Supplementary Table 2. The mutation rate of RTK-RAS, NOTCH, MYC, Cell-Cycle, and TCF-Beta pathways in the high-risk group was significantly higher than that in the low-risk group. Taken together, these results suggested that colon cancer tumor cells in the high- and low-risk groups may have different mutation driver genes and pathways.

Figure 6. Mutant landscape of high-risk and low-risk groups in TCGA. (A, B) Waterfall plot representing the mutant landscape of the top 20 most frequently mutated genes in the high-risk group and low-risk group. (C) Forest plot representing the top 18 genes with significant differences in mutation rates between high- and low-risk groups.

The “CIBERSORT” algorithm was used to estimate the difference in immune infiltration between low-risk and high-risk colon tumors in 22 immune cell subsets. The composition of immune cells in the high-risk and low-risk samples is shown in Supplementary Figure 6A. The difference in the proportion of immune infiltrating cells between the high-risk and low-risk samples was exhibited on the heat map (Supplementary Figure 6B). The immune cells with significantly higher infiltration in the high-risk samples were resting NK cells, and regulatory T cells (P < 0.05). The infiltration of CD8 T cells, plasma cells, memory-activated CD4 T cells, resting dendritic cells, and activated dendritic cells were significantly higher in the low-risk samples than in the high risk samples (P < 0.05) (Figure 7). Therefore, different immune infiltrators in colon cancer patients might be used as prognostic indicators and targets of immunotherapy.

Figure 7. Difference analysis of 22 immune cells infiltration between high- and low-risk groups.

The correlation between the risk scores and expression of five common immune checkpoints was shown on the circular plot (Figure 8A). The results showed that the risk scores were significantly positively correlated with the expression of PD-1, PD-L1, LAG3, and VSIR. The expression of CTLA-4 in high and low risk groups had no significant difference (Figure 8B). The expression of PD-1, PD-L1, LAG3 and VSIR in the high-risk group was significantly higher than that in the low-risk group (Figure 8C–8F).

Figure 8. Correlation of risk scores with expression of five prominent immune checkpoints. (A) Circular plot visualizing correlation coefficient of risk scores with expression of five common immune checkpoints. Box plots showing comparison of the expression of (B) CTLA4, (C) PDL1, (D) VSIR, (E) PD1, and (F) LAG3 between high- and low-risk groups.

Independent prognostic validation of IRG risk score and construction of a nomogram

To explore whether IRG risk score is an independent predictor of colon cancer, age, gender, histological classification, pathological stage, tumor invasion, lymph node, metastasis, and tumor location were incorporated in the univariate analysis. The analysis indicated that age, metastasis, lymph node, tumor stage, and risk score were related to prognosis (Figure 9A). Risk score stratified and the meaningful clinical factors selected by univariate analysis were combined into a multivariate Cox regression analysis. There was collinearity among metastasis, lymph node, and tumor invasion and tumor stage. The multivariate Cox regression model did not include meaningful factors such as metastasis and lymph node. Old age (≥ 65 years: HR 2.09; 95% CI, 1.29–3.41), high pathological stage (stage III: HR 3.96; 95% CI, 1.19–13.14; stage IV: HR 13.02; 95% CI, 3.96–42.83), and high-risk score (HR 3.01; 95% CI, 1.83–4.95) were independent prognostic factors (Figure 9B).

Figure 9. Independent analysis of IRG and construction of nomogram and performance assessment. (A, B) Univariate and multivariate Cox regression analysis of prognostic factors in TCGA. (C) Nomogram based on clinical factors and risk grads in TCGA. (D) AUC values for 1-, 3-, and 5-year survival rates in a nomogram. (E) Calibration for the possibility of 1-, 3-, and 5-year survival in a nomogram.

A nomogram including risk score, age, and pathological stage was constructed that could visualize prognostic risk factors and provide a quantitative method for predicting the survival probability of colon cancer patients (Figure 9C). The AUC values of the nomogram for 1-, 3-, and 5-years were 0.799, 0.791, and 0.738 respectively (Figure 9D). The calibration of the nomogram for the possibility of 1-, 3-, and 5-year survival suggested strong coherence between the prediction and actual observations (Figure 9E).

Comparison of the IRG signature using the ESTIMATE algorithm with models using other methods

To determine whether the immune-related gene signature obtained with the ESTIMATE algorithm was superior to other models, our model was used to compare it with the CIBERSORT and IMMPORT algorithms. Nine immune-related gene signature based on the CIBERSORT algorithm showed that the AUC values of 3-year and 5-year OS were 0.676 and 0.661, respectively. Five immune-related genes signature based on the IMMPORT algorithm showed that the AUC values of 3-year and 5-year OS were 0.663 and 0.713, respectively. Although the other two algorithms do not calculate the 1-year AUC, our model has higher 3-year and 5-year AUC than other models (Table 1). We can preliminarily conclude that the immune-related gene signature we obtained through the ESTIMATE algorithm has higher accuracy.

Data download and preparation

Fragments per kilobase million (FPKM) data on RNA-Seq and corresponding clinical characteristics of colon cancer were downloaded from The Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/, up to June 06, 2020) using the Genomic Data Commons (GDC) tool. The original data included 482 tumor tissues and 42 adjacent tissues. FPKM data was translated into transcripts per million (TPM) expression data. Based on the requirement of data integrality, patients that met the following criteria were excluded from the subsequent analysis: (1) repeated patient records; (2) missing follow-up time, T/N/M staging, or pathological stage information; (3) follow-up time was 0 days. Finally, our study identified 415 colon cancer samples as analytical datasets.

This study selected the GSE39582 dataset from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) for external validation. The GSE39582 dataset was based on a GPL570 platform (Affymetrix Human Genome U133 Plus 2.0 Array) and included 558 colon cancer patients and 19 normal subjects. Our study eventually included 519 patients after deleting 38 samples that had a follow-up time of 0. The raw data were processed by RMA background correction, log2 transformation and normalization using an “affy” package. All data were from publicly available databases, and there were no ethical issues involved.

Identification of immune/stromal related-DEGs and functional enrichment analysis

The ESTIMATE algorithm (https://sourceforge.net/projects/estimateproject/) was used to calculate immune scores and stromal scores [14]. For samples that had a normal distribution, an independent t-test was performed to compare immune and stromal scores between two groups of clinical characteristics; one-way analysis of variance was used for comparison among three or more groups. Otherwise, Wilcoxon and Kruskal-Wallis analysis were conducted for two and three or more groups, respectively. We used the "maxstat" statistic from the "survminer" package to identify the optimal cut-off for continuous variables [56]. Immune and stromal scores were classified into high and low score groups, respectively, according to the optimal cut-off values. Overall survival (OS) was estimated by Kaplan-Meier, and a log-rank test was employed to compare survival differences between the two groups.

Differential expression analysis of high and low score groups was performed with a "limma" package [57]. The P-value was adjusted with the false discovery rate (FDR) [58]. The up- and down-regulated immune and stromal genes were obtained based on the criteria of fold change ≥1.5 and adjusted FDR <0.05. The intersection between immune-related differentially expressed genes (DEGs) and stromal-related DEGs was identified by using the VENNY website (https://bioinfogp.cnb.csic.es/tools/venny/). Heat maps and volcano plots of DEGs were constructed using the “pheatmap” and “ggplot2” packages.

The Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of DEGs were implemented by the Database for Annotation, Visualization, and Integrated Discovery (DAVID) (https://david-d.ncifcrf.gov/) [59–61]. GO analysis included three main parts: biological process (BP), cellular component (CC), and molecular function (MF). We selected the top ten gene ontology terms in three parts to draw the histogram. The top fifteen KEGG analysis terms were exhibited in the bubble chart (where an adjusted P-value <0.05 was considered statistically significant).

Construction, validation, and subgroup analysis of prognostic risk score model

First, a univariate analysis of differential genes was performed to screen out significant genes. Subsequently, the least absolute shrinkage and selection operator (LASSO) Cox regression model was utilized for further screening of prognostic genes to avoid overfitting the model. Finally, a Cox proportional hazard regression was used to determine the optimal prognostic genes of the model [62]. The formula of the gene signature is as follows: Risk score = Ʃ (βi * Expi) (i = the number of prognostic genes, βi represents the gene coefficient, and Expi represents gene expression). Colon cancer patients were divided into high-risk and low-risk groups based on the median risk score. The “survival”, “glmnet”, “survminer”, and “forestplot” packages were used to conduct the above analysis.

The risk score model was validated by 519 colon cancer patients from the GSE39582 dataset. The risk scores of the GSE39582 dataset were calculated based on the above formula. The GSE39582 dataset was also divided into high-risk and low-risk groups based on the median risk score. The performance of the risk score model in the training set and validation set was assessed based on a time-dependent receiver operating characteristic curve (ROC) analysis realized by “timeROC” [63]. An area under time-dependent receiver operating characteristic curve (AUC) greater than 0.75 indicated clearly useful discrimination, while an AUC of 0.60 to 0.75 indicated possibly useful discrimination.

To explore the association between the IRG signature and clinical characteristics, box plots were used to show risk scores for differences in age groups, gender, race, pathological type, tumor location, and tumor stage. A t-test was used to compare whether there were significant differences in risk scores between different clinical features. A Kaplan-Meier method and log-rank test were used to compare the survival curves of different clinical features. To determine that the IRG signature was significant for survival prediction between different clinical characteristics, a survival analysis was performed on each subgroup.

Expression profile and immunohistochemistry (IHC) validation of immune-related genes

Scatter plots were used to show the distribution of gene expression profiles and risk scores. Pearson correlation coefficients were used to indicate the correlation between gene expression profiles and risk scores in the GEO cohort. Images of immunohistochemistry (IHC) staining of the selected prognosis-related genes in normal tissue and colon cancer tissue were retrieved from the Human Protein Atlas (HPA) database (http://www.proteinatlas.org).

Correlations between somatic mutation, immune cell infiltration, immune checkpoint expression, and IRG signature

CIBERSORT is a deconvolution algorithm for immune cell subtype expression based on linear support vector regression [17]. LM22 provides the annotated gene expression signatures for 22 immune cell subtypes, including plasma cells, myeloid subsets, seven T cell types, naive and memory B cells, and natural killer (NK) cells. Compared with traditional immunohistochemistry and flow cytometry methods, the CIBERSORT algorithm can comprehensively, quickly, and accurately infer the relative proportion of the 22 invasive immune cells in tumors [17]. Immune cell fractions from high- and low-risk groups were analyzed using the CIBERSORT algorithm in our research. Immune checkpoint expression was used to predict the immunotherapeutic benefits of various malignancies. Therefore, we assessed the association between the expression of five common immune checkpoints—PD-1, CTLA-4, lymphocyte-activation gene 3 (LAG3), programmed cell death 1(PD-1), and V-Set immunoregulatory receptor (VSIR), with the IRG signature. A circular plot showed the correlation between the risk score and the expression of the immune checkpoint, which was implemented with the "circlize" package. A box plot showed the differences in immune checkpoint expression between high- and low-risk groups.

Verification that IRG risk score is an independent prognostic factor and construction of a nomogram

To assess the independent prognostic effect of the gene risk model, a Cox proportional hazards regression was used to perform univariate and multivariate analysis on all possible prognostic factors. Histological classification included adenocarcinoma and non-adenocarcinoma. Tumor location referred to left-side (descending colon, sigmoid colon, splenic flexure of colon), right-side (ascending colon, cecum, hepatic flexure of colon, transverse colon), or unknown [64]. Age, sex, pathological stage, and T/N/M stage were considered categorical variables.

The nomogram was applied as a visual tool for predicting survival probability. All significant prognostic factors selected by multivariate Cox regression analysis were used to construct a nomogram to obtain a total score, thereby predicting 1-, 3-, and 5-year survival probabilities. Subsequently, time-dependent ROC and calibration were used to assess the performance of the nomogram. If probabilities approach the 45-degree angle line in a calibration plot, it means there is strong agreement between the nomogram prediction and the actual observations.