Comprehensive analysis of prognostic immune-related genes in the tumor microenvironment of colorectal cancer

ESTIMATE analysis shows relationship between tumor purity, stromal scores, and immune scores with CRC prognosis

We downloaded OS and clinicopathological data for 1635 eligible CRC patients from the TCGA and GEO databases with a mean age of 68 years at diagnosis. Among these, 53.9% of the CRC patients were males. ESTIMATE results showed that immune scores ranged from -899.57 and 3202.84, stromal scores ranged from -2232.54 to 2193.08, and tumor purity ranged from 0.27 to 0.98. We then classified the CRC patients into low- and high-score groups based on their immune, stromal, and tumor purity scores, and investigated their correlation with OS rates. We observed positive correlation between stromal scores and OS of CRC patients (p=0.035; Figure 1A). However, immune scores did not show any significant correlation with OS of CRC patients (p=0.381; Figure 1B). We also observed negative correlation between tumor purity and OS (p=0.03; Figure 1C). We then assessed the relationship between immune scores, stromal scores and tumor purity with the pathological stages of CRC. We observed positive correlation between stromal scores and pathologic stage (p = 4.815e-10; Figure 1D). Immune scores were not associated with the pathological stage of CRC patients (p = 0.593; Figure 1E). Tumor purity showed inverse correlation with the pathological stages of CRC patients (p = 1.411e-04; Figure 1F).

Identification of differentially expressed immune- related genes in CRC tissues

As shown by the volcano plots, we identified several differentially expressed genes (upregulated and downregulated) between the low- and high-score groups in the TCGA-CRC cohort (n=611; Figure 2A–2C). We identified 9 downregulated and 1493 upregulated genes in the high stromal score group compared to the low stromal score group (Figure 2A), 61 downregulated and 1235 upregulated genes in the high immune score group compared to the low immune score group (Figure 2B), and 1830 upregulated genes and 69 downregulated genes in the low tumor purity group compared to the high tumor purity score group (Figure 2C). Overall, we identified 4 downregulated and 937 upregulated IRGs that were common among all the three groups (Figure 2D–2E).

Functional enrichment analysis of the differentially expressed IRGs

We performed functional enrichment analysis using the clusterProfiler R package and identified 1454 significantly enriched GO terms and 54 significantly enriched KEGG pathways (FDR<0.05) associated with the 941 differentially expressed IRGs. The top 20 GO terms and KEGG pathways are shown in Figure 3 and Supplementary Table 1. The top GO terms related to biological processes (BF) included regulation of leukocyte activation (GO:0002694), leukocyte migration (GO:0050900), and T cell activation (GO:0042110). The top GO terms related to cellular components (CC) were membrane (GO:0098552) and extracellular matrix (GO:0031012). The top GO terms related to molecular functions (MF) included receptor regulator activity (GO:0030545) and receptor ligand activity (GO:0048018). The top enriched KEGG pathways were cytokine-cytokine receptor interactions (hsa04060), chemokine signaling pathway (hsa04062), PI3K-Akt signaling pathway (hsa04151), and phagosomes (hsa04145).

Construction and validation of the IRG-based prognostic signature

We used univariate COX regression and Lasso regression analysis to identify IRGs that predict the prognosis of CRC patients. We identified 172 IRGs based on the LASSO regression analysis results that were significantly correlated with the OS of CRC patients (Figure 4A, 4B). Subsequently, based on the co-efficient values, we generated an IRG-based prognostic risk signature with ten IRGs. The risk score was calculated with the following formula: (0.1235 × expression of MAP2) + (0.0873 × expression of NKAIN4) + (0.2936 × expression of VAX2 status) + (0.0321 × expression of CALB2) + (0.3958 × expression of FCRL2) + (0.0471 × expression of HAND1) + (0.0986 × expression of A2ML1) + (0.002 × expression of IDO1) + (0.0134 × expression of COL22A1) + (-0.4143 × expression of CD1B). CRC patients from the TCGA database (n=611) were then classified into low- and high-risk groups according to the cutoff risk score of 1.099 (Figure 4C, 4D). The high-risk group CRC patients showed significantly lower OS than those in the low-risk group (p < 0.0001; Figure 4E). The area under the ROC curve (AUC) values for 1-, 3-, and 5-year OS were 0.708, 0.716, and 0.680, respectively (Figure 4F). In the internal validation cohort, patients from the high-risk group demonstrated significantly worse survival outcomes compared to the low-risk group (P<0.001; Supplementary Figure 1). The AUC values for the internal validation cohort were 0.722, 0.712, and 0.738, for 1-, 3- and 5-year OS, respectively (Supplementary Figure 1). Next, we verified the prognostic model in the GSE39582 cohort (n=531) by first classifying the CRC patients into low- and high-risk groups according to the cutoff risk score of 0.989 (Figure 4G, 4H). The high-risk CRC patients from the GSE39582 cohort showed significantly lower OS than the corresponding low-risk patients (p = 0.006; Figure 4I). The AUC values for 1-, 3-, and 5-year OS for CRC patients from the GSE39582 cohort were 0.711, 0.632, and 0.613, respectively (Figure 4J).

The ten-IRG signature is an independent prognostic predictor for CRC patients

We performed univariate and multivariate Cox regression analyses using the TCGA and GEO CRC patient datasets after adjusting for clinicopathological parameters such as M stage, N stage, T stage, tumor stage, age and gender to determine if the risk score based on the ten-IRG signature accurately predicted prognosis of CRC patients. The results showed that the ten-IRG signature-based risk score was an independent prognosis factor for determining the OS of CRC patients from the two datasets (Figure 5A, 5B, 5D, 5E and Supplementary Figure 1). Next, we analyzed the correlation between the ten-IRG signature and other clinicopathological parameters using the Chi-square test. In the TCGA cohort, the higher risk score based on the ten-IRG signature was associated with tumor stage (p<0.001), M stage (p<0.05), N stage (p<0.001), and T stage (p<0.001) (Figure 5C), but was not related to gender and age of CRC patients. We obtained similar results with the GSE39582 validation cohort of CRC patients (Figure 5F).

Nomogram construction and validation

Next, we analyzed if the ten-IRG prognostic signature would enhance the accuracy of predicting OS of CRC patients. Towards this, we constructed a nomogram with the ten-IRG signature-based risk score, M stage, N stage, T stage, and age to determine the 1-, 3-, and 5-year OS of CRC patients (Figure 6A). The calibration plots showed that the nomogram accurately predicted the 3-year OS of CRC patients (Figure 6B). Overall, the nomogram significantly improved the prediction of the 1-, 3-, and 5-year OS of CRC patients (Figure 6C). The OS prediction of the nomogram was significantly higher compared to the OS prediction based on TNM stages alone or age plus the ten-IRG signature-based risk score (Figure 6D). The C-index value for the nomogram was 0.769 (95% CI, 0.717-0.821). The AUC values for the 1-, 3-, and 5-year OS based on the nomogram were 0.788, 0.782 and 0.789, respectively (Figure 6F). We then verified the predictive value of the nomogram in the TCGA and GSE39582 validation cohorts. The C-index value for the internal validation cohort (TCGA) was 0.781 (95% CI, 0.709-0.853), and the AUC values for 1-, 3-, and 5-year OS were 0.834, 0.780, and 0.752, respectively (Supplementary Figure 2). The C-index value for the GSE39582 cohort was 0.732 (95% CI, 0.691-0.773) and the AUC values for 1-, 3-, and 5-year OS were 0.834, 0.780, and 0.752, respectively (Figure 6H). Furthermore, we stratified the CRC patients into high- and low-risk groups based on the median risk score cut off value of the nomogram. Kaplan-Meier survival curves and time-dependent ROC curves showed significantly worse survival outcomes for patients from the high-risk group compared to those from the low-risk group in both the TCGA cohort (P<0.001; Figure 6E, 6F) and GSE39582 cohort (P<0.001; Figure 6G, 6H and Supplementary Figure 2).

Survival and immune cells infiltration analysis of ten IRGs in signature

We performed Kaplan-Meier survival curve analysis of the TCGA-CRC patients based on the expression of the ten prognostic IRGs to determine their independent prognostic values. High expression of A2ML1 (Figure 7A), CALB2 (Figure 7B), COL22A1 (Figure 7D), FCRL2 (Figure 7E), HAND1 (Figure 7F), IDO1 (cut off value: 34, Figure 7G), MAP2 (Figure 7H), NKAIN4 (Figure 7I) and VAX2 (Figure 7J) correlated with shorter OS in CRC patients, whereas, high expression of CD1B (Figure 7C) significantly correlated with longer OS. We observed similar results in the GSE39582 cohort (Supplementary Figure 3A–3J). Interestingly, IDO1 may also have another tumor suppressor gene identity in the GSE39582 cohort (Supplementary Figure 3G) and in TCGA database when expression was between 1.06 and 34 (Supplementary Figure 4).

We performed TIMER database analysis to determine the correlation between the ten-IRG prognostic signature and tumor infiltration of six types of immune cells, namely B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells in the TME. The ten-IRG signature-related risk score showed positive correlation with the levels of CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells in the TME (all P < 0.05), but was not significantly associated with the levels of B cells (P > 0.05; Figure 8). The levels of CD1B expression correlated with the infiltration levels of neutrophils (r= 0.468, p= 3.31E-23) and dendritic cells (r = 0.505, P = 2.05E-27) in colon adenocarcinoma (COAD) samples (Figure 9A and Table 1). Moreover, FCRL2 expression levels showed significant correlation with the infiltration levels of B cells (r= 0.484, p= 3.74E-25), CD4⁺ T Cells (r= 0.447, p= 3.63E-21), and dendritic cells (r = 0.434, P = 6.73E-20) (Figure 9B and Table 1). IDO1 expression levels were significantly associated with the infiltration levels of CD8⁺ T cells (r = 0.4, P = 4.66E-17), neutrophils (r= 0.638, p= 3.97E-47), and dendritic cells (r = 0.564, P = 3.47E-35) (Table 1). MAP2 expression levels were significantly associated with the infiltration levels of CD4⁺ T Cells (r= 0.441, p= 1.38E-20), macrophages (r= 0.535, p= 3.03E-31), and dendritic cells (r = 0.465, P = 5.68E-23) (Figure 9C and Table 1).

Data collection

We searched the GSE and TGCA databases to identify CRC datasets with sample sizes exceeding 50 subjects. The selection criteria included: (1) RNA-Seq or microarray data from transcriptome studies; and (2) clinical data including data regarding clinicopathological parameters such as TNM stages and OS data. Based on these parameters, we selected six datasets, namely GSE103479, GSE72970, GSE41258, GSE39582, GSE12945, and TCGA [32–37] for this study. This included 1635 primary CRC patients that were staged according to the American Joint Committee on Cancer (AJCC) staging system. The clinicopathological data was downloaded from the GEO (https://www.ncbi.nlm.nih.gov/geo/) and TCGA (https://tcga-data.nci.nih.gov/tcga/) databases. The ESTIMATE algorithm was used to evaluate tumor purity as well as stromal and immune scores [16]. The entire TCGA dataset was used for genomic analysis and 60% of the patients from this dataset were randomly selected as an internal validation cohort. The results were also validated using the largest CRC dataset GSE39582 as the external validation cohort. This study complied with the approved publication guidelines for the TCGA and GEO databases. Since we used the data from these public databases, we did not require the approval of the ethics committee from our university or the consent of the patients included in this study.

Identification of immune-related genes (IRGs)

The ESTIMATE algorithm was used to calculate tumor purity, stromal scores, and immune scores of the CRC patient samples, which were then classified into low and high groups using X-tile [38]. X-tile is a novel tool for assessing the biological relationships between potential biomarker genes and prognosis, and identifies optimal cut-points based on the gene expression values [38]. Then, we used the limma R package [39] to identify differentially expressed genes using the cutoff values set at false discovery rate (FDR) < 0.05 and log₂ | fold change | > 1. Finally, we identified the prognostic immune-related genes (IRGs) by comparing their stromal scores, immune scores, and tumor purity.

Functional enrichment analysis of IRGs

We used the clusterProfiler R package [40] to perform Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis and identified the most enriched biological pathways and functions related to the IRGs. We used the false discovery rate (FDR) to determine the adjusted P values and set the statistical significance as FDR < 0.05 to identify the most enriched GO terms for the biological processes (BP), molecular functions (MF), and cellular components (CC) as well as the KEGG pathways related to the IRGs.

Establishment of the IRG prognostic signature for CRC

We investigated the association between IRGs and prognosis of CRC patients by performing univariate Cox regression analysis. P < 0.05 was considered statistically significant. Based on this analysis, we selected a panel of IRGs, which were then evaluated using the LASSO Cox regression analysis with the R package “glmnet”. We then set up a multi-gene signature to evaluate the prognosis of the CRC patients after cross-validating 1000 times with standard error within one standard deviation of the minimum. We then established the most simplified (smallest perimeter) immune gene expression signature model and calculated the risk scores for all patients based on the sum of the corresponding coefficients and expression value of each gene in the prognostic IRG signature model.

TIMER database analysis

We used the deconvolution algorithm from the TIMER database (https://cistrome.shinyapps.io/timer/) to analyze the association between the prognostic risk-related IRGs and the tumor-infiltration levels of immune cells, such as B cells, CD4⁺ T cells, CD8⁺ T cells, macrophages, neutrophils, and dendritic cells [41].

Statistical analysis

The median risk score was used as the cutoff value to classify CRC patients from the GSE39582 and TCGA datasets into low- and high-risk groups. Kaplan-Meier survival curves and log-rank test were used to evaluate the differences in OS rates between the low- and high-risk groups. Univariate and multivariate analysis of clinicopathological parameters and the ten-IRG prognostic signature was performed to determine the prognostic factors for the CRC patients from the TCGA dataset. A nomogram was constructed with the prognostic clinicopathological factors to determine the OS of the CRC patients in the TCGA cohort. The performance of the nomogram was evaluated using C-index and calibration plots with the “rms” R package. The efficiency of the nomogram was validated with the internal validation TCGA cohort and the GSE39582 cohort. Two-tailed t-tests were used to determine statistical significance, which was set at P < 0.05. R software version 3.6.0 was used to perform all statistical analyses.