The RNA-sequencing profile data in the format of Fragments Per Kilobase Million (FPKM) of 612 CRC samples (including 44 normal samples and 568 tumor samples) were obtained from The Cancer Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga/). Meanwhile, the clinical information of 548 CRC cases was also downloaded from TCGA. The samples with incomplete survival information were excluded, and finally, 472 samples with complete follow-up data were included in the next analysis. In addition, RNA-sequencing profile data and corresponding clinical information of 122 samples from Gene Expression Omnibus (GEO, GSE38832, https://www.ncbi.nlm.nih.gov/geo) database were downloaded as a validation cohort. IRGs were obtained using the Immunology Database and Analysis Portal (ImmPort, http://www.immport.org) database, a powerful online platform that provides a mass of immune-related molecular profiles in different tumor types (Bhattacharya et al., 2014).

We used the “limma” package in R to identify differentially expressed genes (DEGs) between tumor and normal samples in TCGA (Ritchie et al., 2015). The threshold for the differentially expressed gene was set as adjusted p value < 0.05 combined with | log2(fold change) | ≥1. Then, based on the ImmPort database, we obtained the expression profile of differentially expressed IRGs in CRC using the “VennDiagram” package in R (Chen and Boutros, 2011). To better understand the biological function of these differentially expressed IRGs in human CRC, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses on these genes via the “clusterProfiler” packages in R (Yu et al., 2012). Only the enrich terms with p < 0.05 were considered statistically significant.

First, all the CRC samples obtained from TCGA were randomly assigned into training and internal validation cohorts at a 1:1 ratio. Then, we used a univariate Cox regression model to identify the association between immune-related mRNAs in relation to the prognosis of CRC patients. IRGs with p < 0.01 in univariate Cox analysis were incorporated into the least absolute shrinkage and selection operator (LASSO) model to prevent overfitting, and the results of LASSO analysis were included into a multivariate Cox model to construct a risk signature. Ultimately, we constructed a prognostic IRG signature, and the risk score of each patient was calculated as follows: = Ʃ (βi × EXPi), where βi (coefficients) represented the weight of the respective signature and EXPi represented the expression level. Based on the median of the risk score, CRC patients in the training cohort were divided into the low- and high-risk groups. In addition, the internal validation cohort and entire cohort were also classified into low- and high-risk groups according to the same cutoff value.

In order to evaluate the prognostic value of the IRG signature, the survival difference between the low- and high-risk groups was illustrated by Kaplan–Meier analysis with the log-rank test. The area under the curve (AUC) in the receiver operating characteristic (ROC) curve was drawn to verify the diagnostic efficacies of the IRG signature based on the “survivalROC” package in R (Heagerty et al., 2000). In addition, the predictive value of the IRG signature was also assessed in the internal validation and entire cohorts. Furthermore, we also performed stratified survival analysis to further determine the predictive value of the IRG signature in distinct subgroups in the entire cohort.

To determine whether the IRG signature had a similar predictive value in different populations, its prognostic capability was further validated using the GEO cohort. The same prognostic scoring system developed in the training cohort was used to calculate the risk score for each patient in the GEO cohort, and patients were also divided into low- and high-risk groups based on the median of risk score. Then, survival curve, risk curve, and ROC curve analysis were performed.

To evaluate the clinical relevance of the developed IRG signature, we first obtained the clinicopathological information from TCGA, including age, gender, American Joint Committee on Cancer (AJCC) tumor stage, T classification, N classification, and M classification. Then, the differences of these clinicopathological features between low- and high-risk groups were evaluated using the Chi-square test. In addition, we performed the Wilcoxon test to assess the risk score differences across different subgroups based on the clinicopathological characteristics.

We constructed a nomogram consisting of clinical parameters and risk score based on multivariate Cox regression analysis using “survival,” “rms,” and “foreign” packages in R software. The calibration curve was plotted to compare the association between the actual outcomes and the predicted probabilities. The AUC of the ROC curve was plotted to compare the predictive abilities of the nomogram with other prognostic factors. Furthermore, we used decision curve analysis (DCA) to evaluate the net benefit of the nomogram in a clinical context.

To determine whether our IRG signature is better than previously identified models in CRC, we compared the predictive ability of our IRG signature with other signatures, including a six-mRNA signature (Dai et al., 2020), a five-mRNA signature (Zhu et al., 2020), and two eight-mRNA signature (Wang et al., 2020; Wu et al., 2020). We obtained the IRGs in these signatures from the corresponding published literature, and then calculated the C-index and AUC of 1-, 3-, and 5-year ROC curves for each signature.

To explore the correlation between the risk score and TIIC characteristics, we applied the currently developed deconvolution algorithm, including TIMER (Li et al., 2017), QUANTISEQ (Plattner et al., 2020), MCPcounter (Dienstmann et al., 2019), EPIC (Racle et al., 2017), and CIBERSORT (Chen et al., 2018) to calculate the abundances of distinct TIICs in each sample from the TCGA project. We developed a lollipop diagram to show the correlation between the risk score and TIIC abundance based on Spearman correlation analysis, and the difference in TIIC abundance between low- and high-risk groups was analyzed using the Wilcoxon test and shown in a box plot. In addition, we compared the expression levels of immunotherapeutic biomarkers between low- and high-risk groups based on the Wilcoxon test.

All statistical analyses were performed using R software (version 4.0.2). The Student’s t-test was applied to identify statistically differentially expressed genes, and the Wilcox test was used to identify statistically differentially infiltrative immune cells. Fisher’s exact test was used to evaluate the significant GO terms and KEGG enrichment pathways. Cox regression and Kaplan–Meier curves with the log-rank test were used for evaluating the survival data. The ROC curve was used to determine the predictive accuracy of the risk signature, and the DCA and calibration curves were used to explore the reliability of the nomogram. The correlation between the risk score and TIIC abundance was evaluated using the Spearman correlation coefficient. p value < 0.05 was considered statistically significant.

In total, 6,482 DEGs including 4,518 upregulated and 1,964 downregulated genes were identified based on the threshold of adjusted p value < 0.05 and | log2(fold change) | ≥1. The heatmap of the DEGs was shown in Figure 1A. Meanwhile, a total of 1,811 IRGs were extracted from the ImmPort database. Then, 476 differentially expressed IRGs were identified using the intersect function of Venn diagram (Figure 1B). The heatmap of the differentially expressed IRGs was shown in Figure 1C. We then performed GO and KEGG enrichment analyses to further explore the biological function of these IRGs. As is shown in Figure 1D, GO enrichment analysis found that these IRGs are mainly enriched in immune-related crosstalk, including complement activation, humoral immune response, and leukocyte migration. KEGG pathway enrichment analysis showed that these genes are mainly enriched in several areas of tumor-related pathway, including the Rap 1 signaling pathway, MAPK signaling pathway, and NF-κB signaling pathway (Figure 1E).

First, the entire TCGA cohort (n = 472) was randomly divided into the training cohort (n = 236) and internal validation cohort (n = 236), and the clinical parameters of the CRC patients in each cohort are shown in Table 1. Then, we used the expression profiles of the differentially expressed IRGs in the training cohort to construct a survival-related signature. In total, 17 IRGs were identified to be associated with the prognosis of patients with CRC using univariate Cox regression analysis (p < 0.01, Figure 2A) and 14 of which were filtered out using a LASSO regression analysis with a 10-fold cross validation (Figures 2B,C). Finally, 10 IRGs were identified as independent prognostic factors and were selected for constructing the IRG signature (Figure 2D; Table 2). The risk score for each sample was calculated based on the expression level and Cox regression coefficient of these 10 genes. Risk score = (−5.96137 × EXP_CD1B) + (0.19671 × EXP_TFR2) + (0.56983× EXP_FGF2) + (0.17961 × EXP_SEMA3G) + (0.12713 × EXP_PLXNA3) + (−4.98459 × EXP_NRG1) + (0.22476 × EXP_UCN) + (0.16455 × EXP_IL1RL2) + (0.43311 × EXP_PTH1R) + (0.12939 × EXP_TRDC).

Based on the cutoff value of 1.763 for the risk score, patients were divided into the low- (n = 118) and high-risk (n = 118) groups. The Kaplan–Meier analysis with the log-rank test found that patients in the high-risk group presented worse prognosis than those in the low-risk group (Figure 3A). CRC patients appeared to have an increased mortality rate with the increase of risk scores according to the risk plot (Figure 3B). The time-ROC curve showed that the AUC of 1-, 3-, and 5-year survival were 0.822, 0.834, and 0.888, respectively (Figure 3C). Then, we further validated the immune signature in the internal validation and the entire TCGA cohorts. The prognosis of CRC patients in the high-risk group was significantly worse than that in the low-risk group (Figure 3D). The patients with CRC have an increased mortality rate with the increase of risk scores according to the risk plot (Figure 3E). The AUC values of 1-, 3-, and 5-year survival in the internal validation cohort were 0.745, 0.620, and 0.607, respectively (Figure 3F). Similarly, in the entire TCGA cohort, patients in the high-risk group presented worse prognosis than those in the low-risk group (Figure 3G). The mortality rate of patients was increased with the increase of risk scores according to the risk plot (Figure 3H), and the AUC for 1-, 3-, and 5-year survival was 0.778, 0.737, and 0.757, respectively (Figure 3I). These aforementioned results indicated that the IRG signature was competent for predicting the prognosis of human CRC.

We performed stratification analysis to further determine the widespread utility of the IRG signature based on the clinicopathological variables. As shown in Figure 4, we found that the IRG signature has predictive significance for different subgroups, and the prognosis of patients in the low-risk group was significantly superior to that of patients in the high-risk group. Taken all above, these findings indicate that the IRG signature exerts vital roles in predicting the prognosis of CRC patients.

Finally, we further validated the prognostic value of the TRG signature in an independent cohort from the GEO project. Similarly, survival analysis in 122 CRC patients also found that the high-risk score indicated poor survival outcomes (Figures 5A,B). Correspondingly, the 1-, 3-, and 5-year AUC in predicting prognosis was 0.577, 0.636, and 0.647, respectively (Figure 5C).

We performed a series of Shi-square tests and Wilcoxon tests to determine the correlations between the risk score and the clinicopathological parameters in patients with CRC. As shown in the strip chart (Figure 6A) and scatter diagrams (Figure 6B), The risk scores were found to be significantly associated with patients’ clinicopathologic parameters, including T classification, N classification, M classification, and AJCC tumor stage. These aforementioned results indicated that the high-risk score was significantly correlated with more aggressive clinicopathological characteristics.

The age, gender, stage, and risk score were integrated to construct a nomogram to predict the survival probability of CRC patients at 1, 3, and 5 years (Figure 7A). The calibration plot showed that the nomogram performed well compared against the performance of an ideal model (Figure 7B). The ROC curve showed that the AUC values corresponding to nomogram, risk score, stage, gender, and age were 0.800, 0.781, 0.737, 0.436, and 0.622, respectively (Figure 7C). In addition, the DCA curve proved that the nomogram provided better net benefits than other prognostic factors (Figure 7D). These results indicated that the nomogram signified satisfied accuracy in predicting the prognosis of CRC patients.

We compared our IRG signature with previously identified models to determine whether the predictive efficiency of our signature is higher than others. As shown in Figure 8A, the C-index of our signature is 0.818, which is higher than that in other models. The Kaplan–Meier survival curve showed that the prognosis of the low- and high-risk groups significantly differed across these signatures (Figures 8B–F). The ROC curve showed that the AUC values of the 1-, 3-, and 5-year survival for our signature are 0.822, 0.834, and 0.888, respectively (Figure 8B). The AUC values of the 1-, 3-, and 5-year survival for Dai et al. signature are 0.779, 0.799, and 0.748, respectively (Figure 8C). The AUC values of the 1-, 3-, and 5-year survival for Zhu et al. signature are 0.593, 0.627, and 0.759, respectively (Figure 8D). Wang et al. signature showed that the AUC values of the 1-, 3-, and 5-year survival are 0.793, 0.763, and 0.781, respectively (Figure 8E). Wu et al. signature showed that the AUC values of the 1-, 3-, and 5-year survival are 0.779, 0.714, and 0.733, respectively (Figure 8F). By comparing with previously constructed signatures, we found that the predictive efficiency of our IRG signature is higher than that of other models.

In order to explore whether the IRG signature was associated with the immune characteristics of CRC, we applied five algorithms to estimate the abundance of TIICs. As shown in the lollipop diagram (Figure 9A), Spearman correlation analysis revealed that the risk score was significantly negatively correlated with the infiltration levels of multiple TIICs, including CD4⁺ T cells, CD8⁺ T cells, NK cells, macrophages, B cells, and neutrophils. The box plot showed that the abundance of these TIICs was significantly higher in the low-risk group than that in the high-risk group (Figure 9B). In addition, we further explored the expression of immunotherapeutic biomarkers in low- and high-risk groups. As shown in Figure 9C, compared to the high-risk group, CRC patients in the low-risk group express higher levels of immunotherapeutic biomarkers, including PD-1, PD-L1, TIM-3, and CTLA-4. These results indicated that CRC patients in the low-risk group had an immune “hot” phenotype and were more likely to benefit from immunotherapy. That is, to say, the IRG signature can predict the immune characteristics of human CRC.