Clinical Significance of a Novel Tumor Progression-Associated Immune Signature in Colorectal Adenocarcinoma

Background Some colorectal adenocarcinoma (CRC) patients are susceptible to recurrence, and they rapidly progress to advanced cancer stages and have a poor prognosis. There is an urgent need for efficient screening criteria to identify patients who tend to relapse in order to treat them earlier and more systematically. Methods We identified two groups of patients with significantly different outcomes by unsupervised cluster analysis of GSE39582 based on 101 significantly differentially expressed immune genes. To develop an accurate and specific signature based on immune-related genes to predict the recurrence of CRC, a multivariate Cox risk regression model was constructed with a training cohort composed of 519 CRC samples. The model was then validated using 129, 292, and 446 samples in the real-time quantitative reverse transcription PCR (qRT-PCR), test, and validation cohorts, respectively. Results This classification system can also be used to predict the prognosis in clinical subgroups and patients with different mutation states. Four independent datasets, including qRT-PCR and The Cancer Genome Atlas (TCGA), demonstrated that they can also be used to accurately predict the overall survival of CRC patients. Further analysis suggested that high-risk patients were characterized by worse effects of chemotherapy and immunotherapy, as well as lower immune scores. Ultimately, the signature was identified as an independent prognostic factor. Conclusion The signature can accurately predict recurrence and overall survival in patients with CRC and may serve as a powerful prognostic tool to further optimize cancer immunotherapy.


INTRODUCTION
Colorectal adenocarcinoma is the second most common cancer in men and the third most common cancer in women worldwide (Thrumurthy et al., 2016;Miller et al., 2019). The 5-year relative survival rate of CRC patients is 65%. For patients diagnosed with stage I and stage II CRC, the 5-year relative survival rates are 91 and 82%, respectively. However, the 5-year survival rate of stage IV disease is only 12% (Miller et al., 2019;Siegel et al., 2020). Additionally, even in early-stage patients, there is a significant recurrence rate after surgical removal, and patients with relapse tend to have a poor prognosis (Thrumurthy et al., 2016;Lefèvre et al., 2019). Novel molecular biomarkers that can precisely indicate the stage of disease progression and predict clinical outcomes are urgently needed.
The immune response in the tumor microenvironment is considered to be an important factor in determining tumor aggressiveness, progression, and response to immunomodulators. The density and type of tumor-infiltrating immune cells, as well as their cytokine and immune-related gene (IRG) expression, have been extensively studied as prognostic biomarkers for CRC (Koch et al., 2006;Dupaul-Chicoine et al., 2015). In addition, previous studies have reported the important value of using genes or lncRNA signatures to predict recurrence and prognosis in patients with CRC (Tian et al., 2017;Mu et al., 2020). However, whether these IRG signatures can also be used as predictors of recurrence and prognosis of CRC remains to be explored.
Considering the increasing proportion of patients diagnosed with early-stage CRC and the poor prognosis of advanced-stage patients, there is an urgent need to develop an IRG-based recurrence signature (IGBRS) in patients with CRC. Comprehensive analysis of IRG and the tumor microenvironment in CRC can improve the stratification of the risk of CRC patients clinically and allow exploration of possible biotherapeutic targets. In the current study, we integrated 1,386 CRC cases with recurrence-free survival (RFS) data and 1,375 cases with OS data from eight independent cohorts; a dataset from TCGA; and GSE14333, GSE17538, GSE33113, GSE37892, GSE38832, GSE39582, and real-time quantitative reverse transcription PCR (qRT-PCR) cohorts to establish a novel robust IGBRS. Further, we studied whether IGBRS can predict clinical outcomes independent of the clinical and pathological characteristics and molecular subtype. We also explored the differences in the immune landscape and somatic mutation pattern between high-and low-risk groups.

Data Acquisition and Preprocessing
The entire analytical process of the study is presented in Figure 1.
The GSE14333, GSE17538, GSE33113, GSE37892, GSE38832, GSE35640, GSE63557, and GSE39582 datasets were downloaded from the GEO 1 database with log2 transformation and quantilenormalized matrices. In general, the expression of genes with multiple probes is based on the median. The mRNA expression profiles were in the fragments per kilobase of transcript per million mapped reads (FPKM) format. The data used in this study met the following criteria: (1) the expression level of each probe must be greater than 0 in ≥75% of the samples and (2) data with respect to survival time and survival status must be available for each patient. GSE39582 was developed as a training cohort including 519 patients with RFS and 562 patients with OS. We combined the GSE17538 and GSE38832 datasets as a test cohort because they collectively have RFS (292 patients) and disease-specific survival (DSS, 299 patients) data. In addition, the GSE14333, GSE33113, and GSE37892 datasets were combined as a validation cohort owing to their single RFS (446 patients) data. When combining the GSE datasets, we used the combat function of the R software package sva to remove the batch effect. Additionally, GSE17538 contained 232 patients with OS data.
Moreover, between July 2013 and August 2015, a total of 129 frozen surgically resected tumor tissues were obtained from patients with a pathological diagnosis of CRC at the Department of Colorectal Surgery, National Cancer Center/National Clinical Research Center for Cancer/Cancer Hospital, Chinese Academy of Medical Sciences and Peking Union Medical College.
The TCGA-COAD and TCGA-SKCM datasets were downloaded from TCGA 2 , including mRNA expression profiles of 452 CRC specimens and the corresponding clinical follow-up data. The basic information of the dataset included in this study is shown in Supplementary Table 1.

Identification of Differentially Expressed IRGs in CRC and Adjacent Normal Tissues
We used the IRGs obtained from the ImmPort database 3 to intersect with the mRNA data in the TCGA and GEO database matrices. Then, the edgeR package in R software (Robinson et al., 2010) was used to analyze the difference in expression between CRC and adjacent normal tissues of the TCGA and GSE39582 datasets of common IRGs. Significantly differentially expressed RNAs were identified by setting adjusted P-values < 0.05 and | log 2 [fold change (FC)]| > 1. We chose the IRGs that were significantly differentially expressed in both databases at the same time. We plotted a volcano map using the ggplot2 package.

Unsupervised Clustering for Differentially Expressed IRGs
Unsupervised cluster analysis was used to classify the patients according to the expression of 101 differentially expressed IRGs identified by differential analysis. The number of clusters and their stability were determined by the consensus clustering algorithm (McLachlan et al., 2017). We performed these steps based on the ConsensusClusterPlus package (Wilkerson and Hayes, 2010) and repeated it 1,000 times to ensure the stability of the classification.

Prognostic Evaluation Using the IGBRS
The prognostic value of each IRG was first calculated by univariate Cox analysis with the R/survival package, and IRGs with P < 0.05 were selected as seed IRGs for Cox LASSO regression. Next, multivariate Cox regression was applied to identify prognostic signatures with the R packages glment, survminer, and survival. The risk scores for each patient in the training group were calculated based on the following formula: risk score = {expGene(n) × βGene(n)}, where exp is the prognostic gene expression level and β is the multivariate Cox regression model regression coefficient (Yang et al., 2019). All samples were randomly divided into high-and low-risk score sets, with the median risk score as the cutoff value. The Kaplan-Meier method was used to generate survival curves for each group, and the log-rank test was used to determine the statistical significance of differences. The R packages survivalROC and timeROC were used to plot and visualize ROC curves to calculate the AUC and confidence intervals (CIs) to evaluate the diagnostic accuracy of the IGBRS.

Real-Time Quantitative Reverse Transcription PCR
RNA was extracted from the frozen samples using TRIzol reagent (Takara, #9109) according to standard protocols. Total RNA was reverse transcribed into cDNA with random primers using the Transcriptor First Strand cDNA Synthesis Kit (Roche, Penzberg, Germany) following the manufacturer's instructions. RNA expression levels were measured by real-time quantitative reverse transcription PCR (qRT-PCR) using FastStart Essential DNA Green Master mix (Roche, Penzberg, Germany) on a Roche LightCycler 480 (Roche, Penzberg, Germany) with the following PCR conditions: pre-denaturation at 95 • C for 2 min, followed by 45 cycles of denaturation at 94 • C for 15 s, annealing at 55 • C for 15 s, and extension at 68 • C for 30 s. The melting curve was prepared at 65-95 • C. The mRNA expression of each of the seven genes was normalized to GAPDH levels. All quantitative PCR analyses were conducted in triplicate and the average value was calculated. Relative gene expression was analyzed by the 2 − Ct method. We verified the specificity of the PCR primers using BLAST. A single peak in the melting curve indicated that the PCR products were specific. The primers used in the study are presented in Supplementary Table 2.

Validation of the Prognostic Value of the IGBRS
The IGBRS was validated in different clinical subgroups and histopathological subtypes. A similar analysis process was also applied to the qRT-PCR, test, validation, and TCGA cohorts. Univariate and multivariate Cox regression analyses of the IGBRS and other clinicopathological factors were performed to evaluate whether the IGBRS is an independent risk factor for CRC.

Exploration of the Relationships Between the IGBRS and Immunity
The ESTIMATE algorithm was used to determine the stromal score, the ESTIMATE score, and immune scores of each sample with R software, and the differences in the degree of immune cell infiltration between the high-and lowrisk groups were further compared using the Wilcoxon test (Runa et al., 2017). The CIBERSORT.R package was used to assess the proportions of 22 immune cell subtypes based on the expression profile (Newman et al., 2015). CRC samples with P < 0.05 in the CIBERSORT analysis results were used in further analyses. The Mann-Whitney U test was used to compare differences in immune cell subtypes in the high-risk and low-risk groups. We also compared the differences in the expression of 49 immune checkpoints {including the B7-CD28 family [CD274 (PD-L1), B7-H3, CTLA4, ICOSLG, PD-L2, TMIGD2, ICOS, PD-1, and HHLA2] (Janakiram et al., 2015b;Zhang et al., 2018), the TNF superfamily [TNFRSF14, TNFRSF18, TNFRSF25, TNFRSF4, TNFRSF8, TNFRSF9, TNFSF14, TNFSF15, TNFSF18, TNFSF4, TNFSF9, CD40, BTLA, CD27, CD40LG, and CD70] (Ward- Kavanagh et al., 2016), and several other immune checkpoint members [ADORA2A, BTNL2, CD160, CD200, CD200R1, CD244, CD28, CD44, CD48, CD80, CD86, ENTPD1, FGL1, HAVCR2, IDO1, IDO2, KIR3DL1, LAG3, LAIR1, LGALS9, NRP1, NCR3, and TIGIT] (Chrétien et al., 2019;Wang et al., 2019a,b)} between the high-and low-risk groups to analyze the landscape of genetic variation (R package maftools). Moreover, we calculated the Pearson correlation coefficients between seven IRGs and differentially expressed immune checkpoint genes based on the gene expression. In addition, we calculated the Pearson correlation coefficients between seven IRGs and immune cells based on the gene expression of these IRGs and immune cell infiltration in each sample. P < 0.05 was considered to be statistically significant. We corrected P-values by the Benjamini-Hochberg method (Hochberg and Benjamini, 1990).

Gene Set Enrichment Analysis
Gene set enrichment analysis (GSEA) was performed using the software GSEA v4.0.3 4 . We divided all samples into highand low-risk score groups according to the median cutoff value of the risk score. We input the profiles of the adjusted expression data for all transcripts, groups of high-and low-risk score samples, and gene set files (c2. cp. kegg. v6. 1. symbols. gmt). Enrichment P-values were based on 10,000 permutations and were subsequently adjusted for multiple testing using the Benjamini-Hochberg procedure to control the false discovery rate (FDR) (Hochberg and Benjamini, 1990).

Tumor Somatic Mutation Analysis
The waterfall function of the maftools package was used to present the mutation landscapes in patients with high-and low-risk score subtypes in the TCGA-COAD cohort. The total number of somatic mutations was adopted to assess the mutation burden, which is convenient and significantly correlated with the number of non-synonymous mutations. Missense, nonsense, non-stop, silent, and frameshift/in-frame insertions and deletions were counted and summed, and germline mutations without somatic mutations were excluded . According to the median tumor mutation burden (TMB), all CRC samples with somatic mutations in the TCGA dataset were divided into the high-TMB group and the low-TMB group.

Prognostic Meta-Analysis
In order to clarify the comprehensive prognostic value of IGBRS in four different groups, prognostic meta-analysis was carried out using STATA software (version 12.0). Then, the random effect model was used to calculate the combined hazard ratio (HR) value.

Identification of Differentially Expressed IRGs and Unsupervised Cluster Analysis
First, we combined the two datasets GSE17538 and GSE38832 as the test cohort and merged GSE14333, GSE33113, and GSE37892 as the validation cohort and removed the batch effect (Supplementary Figures 1A,B). Then, IRGs coincident in the TCGA, GSE39582, test, validation, and ImmPort databases were screened (Supplementary Figure 1C), and difference analysis was performed in CRC and adjacent normal tissues of TCGA and GSE39582. Next, we obtained 328 and 124 differentially expressed IRGs from the two respective datasets. Ultimately, we identified a total of 101 overlapping IRGs from two cohorts (Supplementary Figure 1D). The log 2 (FC) values of seven IRGs in the IGBRS are presented in Supplementary Figure 1E. Based on 101 differentially expressed IRGs, GSE39582 was divided into two groups by unsupervised clustering (Figures 2A-C). The prognostic analysis of the two groups showed that cluster-1 had significant advantages in survival and no recurrence (Figures 2D,E).

Construction of the IGBRS for CRC With the Training Cohort
A total of 101 IRGs were evaluated by univariate Cox survival analysis, and 25 IRGs with P < 0.05 were filtered out and included in subsequent analyses ( Figure 3A). As shown in Figure 3B, LASSO regression analysis identified seven IRGs (based on lambda.lse criteria) that were subjected to multivariate Cox regression analysis ( Figure 3C). The results are presented in Supplementary  Table 3. Ultimately, we identified seven IRGs predicting CRC patient recurrence, namely, BMP4, CXCL3, GZMB, IL1R2, LGR5, PLAU, and PTGDR. The risk score of the IGBRS was calculated based on the following formula: risk The patients in the training set were divided into a high-risk group and a low-risk group, with the median risk score as the cutoff point (Supplementary Figure 3A). The patients in the high-risk group had significantly worse RFS than those in the low-risk group (log-rank test P-value < 0.0001, Figure 4B). In addition, ROC analysis was implemented to determine whether survival predictions made with the IGBRS were accurate ( Figure 4A). The AUC values were assessed for 3-year (AUC = 0.748) and 5-year (AUC = 0.731) RFS, and these values were higher than the AUC of the traditional TNM pathological staging system (Figures 4C,D). The results suggest that the IGBRS can be used to effectively evaluate the recurrence of CRC patients. Furthermore, patients in the high-risk group had worse RFS than those in the low-risk group for pathological stage II (P = 0.0058), pathological stage III (P < 0.0001), and pathological stage IV (P = 0.04) (Figures 4E-G), which means that the IGBRS can create a good risk stratification in different pathological stages of CRC. Similarly, to explore whether the IGBRS is equally valid for predicting the OS of CRC patients, we conducted a similar analysis process on the patients with OS data. Comparisons of OS showed that the mortality rate in the high-risk group was significantly higher than that in the low-risk group ( Figure 4I, P < 0.0001). The AUC for 5-year OS using the predictive nomogram reached 0.702 ( Figure 4H).

Verification of the Prognostic Classifier With the qRT-PCR Cohort
To determine the robustness of this signature, this verification process was also performed for the qRT-PCR cohort. By using the same risk score formula in the validation set, we classified patients into high-risk score (n = 64) and low-risk score (n = 65) groups, with the median score as the cutoff point ( Figure 5A). The results showed that the Kaplan-Meier curves for RFS suggested that the patients with high-risk scores had significantly worse RFS than those with low-risk scores (log-rank test P-value < 0.0001, Figure 5B). The AUC for 6-year RFS reached 0.884 ( Figure 5C). We conducted a similar analysis process on the patients with OS data. Comparisons of OS showed that the mortality rate in the high-risk group was significantly higher than that in the low-risk group ( Figure 5E, P < 0.0001). The AUC for 6-year OS using the IGBRS reached 0.890 ( Figure 5F).
Additionally, in all CRC patients who underwent chemotherapy, the patients in the high-risk groups had shorter RFS and OS times than those in the low-risk groups (Figures 5D,G; P < 0.0001). These results indicate that the IGBRS-based risk score is a promising prognosis predictor whether or not the alluded chemotherapy drugs were administered ( Figures 5H-K).

Verification of the Prognostic Classifier With the Test Cohort
In agreement with the abovementioned findings, we also validated the robustness of the IGBRS with the test dataset (n = 292). We classified the patients into high-risk score (n = 146) and low-risk score (n = 146) groups, with the median score as the cutoff point (Supplementary Figure 3B). The AUC for 3-year recurrence prediction by the IGBRS reached 0.754 (Supplementary Figure 4A). The Kaplan-Meier curves suggested that the patients with high-risk scores had significantly worse RFS than those with low-risk scores (log-rank test P-value < 0.0001, Supplementary Figure 4B). Furthermore, patients in the high-risk group had worse RFS than those in the low-risk group for pathological stage II (P = 0.0012), pathological stage III (P = 0.039), and pathological stage IV (P = 0.025) (Supplementary Figures 4C-E), which demonstrated that the IGBRS can create a good risk stratification in different pathological stages of CRC. To explore whether the IGBRS is equally valid for predicting the DSS of CRC patients, we conducted a similar analysis process on the patients with DSS data. The AUC for 3-year DSS using the predictive nomogram reached 0.752 (Supplementary Figure 4F). Comparisons of DSS showed that the patients in the high-risk group had a significantly higher mortality rate than those in the low-risk group (Supplementary Figure 4G, P < 0.0001). Furthermore, we validated the ability of the IGBRS to predict OS for the 232 patients with OS data in the GSE17538 set. The results show that the AUC for 3-year OS reached 0.721 (Supplementary Figure 4H), and the patients with a high-risk score had a significantly higher mortality rate than those with a low-risk score (Supplementary Figure 4I, P < 0.0001).

Verification of the Prognostic Classifier With the Validation and TCGA Cohorts
To validate the robustness of the IGBRS, this verification process was also performed for the validation and TCGA cohorts. Using the same risk score formula in the validation set, we classified patients into high-risk score (n = 223) and low-risk score (n = 223) groups, with the median score as the cutoff point (Supplementary Figure 5A). The results showed that the AUC for 5-year RFS reached 0.660 (Supplementary Figure 5B). The Kaplan-Meier curves for RFS suggested that the patients with high-risk scores had significantly worse RFS than those with low-risk scores (log-rank test P-value = 0.00015, Supplementary Figure 5C). In addition, we observed that the AUC values for the assessment of 5-year RFS (AUC = 0.624) by IGBRS were higher than those of the Duke staging system (AUC = 0.487) in the GSE14333 set (Supplementary Figure 5D). In the GSE33113 and GSE37892 cohorts, the AUC values for the assessment of 5-year RFS (AUC = 0.700) by the IGBRS were higher than those of the traditional TNM pathological staging system (AUC = 0.673) (Supplementary Figure 5E). It is important that, in the pathological stage II and III subgroups, the patients in the high-risk group had worse RFS than those in the low-risk group (Supplementary Figures 5F,G; P < 0.05).
Using the same signature score formula in the TCGA set (n = 452), we classified the patients into groups with a high-risk score (n = 226) and a low-risk score (n = 226) by taking the median score as the cutoff point (Supplementary Figure 6A). The AUC exhibited by the IGBRS for 5-year survival reached 0.686 (Supplementary Figure 6B). The Kaplan-Meier OS curves suggested that patients with high-risk scores had significantly worse OS than those with low-risk scores (log-rank test P-value = 0.037, Supplementary Figure 6C).

Validation of the IGBRS in Different Clinical Subgroups
Although the pathological stage is the most significant factor that influences CRC patient survival, other factors such as sex, age, and location are also important (Loupakis et al., 2015;Thrumurthy et al., 2016;Guo et al., 2019;Siegel et al., 2020); therefore, we stratified CRC patients in the training cohort by four clinical characteristics. The results suggested that in all subgroups including males and females, older (age ≥65 years) and younger (age <65 years) patients, and distal and proximal, the patients in the high-risk groups had shorter RFS times than those in the low-risk groups (Supplementary Figure 7, P < 0.05). More convincingly, similar results were obtained in the validation dataset (Supplementary Figure 8). Additionally, in all CRC patients who underwent chemotherapy, the patients in the highrisk groups had shorter RFS and OS times than those in the low-risk groups (Supplementary Figures 7G,H; P < 0.0002). These results suggest that the effects of chemotherapy in high-risk patients were significantly worse than those in low-risk patients. Furthermore, among the subtypes of chemotherapy, including 5FU, FOLFOX, and FUFOL, the RFS and OS of chemotherapy in high-risk patients were significantly worse than those in low-risk patients (Supplementary Figure 9).

IGBRS Can Predict Clinical Outcomes Independent of TP53, KRAS, or BRAF Mutation Status
In view of the fact that TP53, KRAS, and BRAF are commonly mutated genes in CRC (Russo et al., 2005;Sinicrope et al., 2015;Guo et al., 2019), we analyzed the performance of the IGBRS among patients with different TP53, KRAS, and BRAF mutation statuses. First, we analyzed the proportion of high-risk and low-risk patients with different mutation statuses in the training cohort. The results showed that IGBRS can predict RFS and OS independent of the TP53, KRAS, or BRAF mutation status (Supplementary Figures 10, 11, P < 0.01).

IGBRS Can Predict Clinical Outcomes Independent of CpG Island Methylation Phenotype, Mismatch Repair, and Chromosomal Instability Status
A CpG island methylation phenotype (CIMP) occurs when promoter CpG island methylation causes epigenetic silencing, which is an epigenetic phenotype of CRC, called CIMP cancer (Lao and Grady, 2011). Chromosomal instability (CIN) is a special phenotype of CRC, which is found in the majority of CRC and leads to a different pattern of gene alterations that contribute to tumor formation (Grady, 2004). We analyzed the performance of the IGBRS in patients with Moreover, compared with the CIN-positive group, the CIN-negative group had a higher risk score (Supplementary Figure 12H, P < 0.01), while compared with the high-risk group, CIN-positive patients were obviously concentrated in the low-risk group (Supplementary Figure 12G, 31% vs. 15%, P < 0.01), which may have something to do with the role of CIN in triggering the immune response. Among the MMRproficient (pMMR) and MMR-deficient (dMMR) subgroups, patients in the dMMR subgroup had a lower risk score than those in the pMMR subgroup (Supplementary Figure 13H, P < 0.0001). Furthermore, the higher risk score was obviously concentrated in the CIMP-positive group compared with the CIMP-negative group. Besides, compared with the low-risk

Immune Cell Infiltration and Immune Checkpoint Correlates of the IGBRS
Tumor microenvironments contain a variety of cell types, including immune cells, interstitial cells, endothelial cells, and inflammatory mediators, and extracellular matrix (ECM) molecules (Runa et al., 2017). To explore the potential mechanisms underlying the association between the IGBRS and CRC recurrence, we used a series of analysis methods related to the immune profile. In the training, test, and TCGA cohorts, the immune score of the high-risk group was significantly lower than that of the low-risk group (Figures 6A-C). Kaplan-Meier analysis of the CRC data in the three cohorts showed that RFS was significantly shorter in the low-immune score groups (Figures 6D-F).
Since risk scores are closely related to immune infiltration scores, we analyzed specific immunophenotypic differences between high-and low-risk groups, including immune cell infiltration and immune checkpoints. The CIBERSORT package was used to assess the proportions of 22 immune cell subtypes based on the expression profile. CRC samples with P < 0.05 in the CIBERSORT analysis results were used for further analyses. There were large differences in the proportions of infiltrating immune cells between the high-FIGURE 6 | The relationship between immune infiltration scores and risk scores and recurrence. (A,B) Differences in the stromal, immune, and ESTIMATE scores (the ESTIMATE score can be obtained by adding the immune score and the stromal score, which can be used to estimate tumor purity) between the high-and low-risk groups in the training, test, and TCGA cohorts. (C-F) Impact of the immune score on RFS of CRC patients based on Kaplan-Meier analysis in the training, test, and TCGA cohorts. The ordinate represents the values of the three scores. *FDR < 0.05, **FDR < 0.01, ***FDR < 0.001, and ****FDR < 0.0001. TCGA, The Cancer Genome Atlas; CRC, colorectal adenocarcinoma; RFS, recurrence-free survival; FDR, false discovery rate. and low-risk score groups in the training and test cohorts, which were predominantly M1 macrophages, M2 macrophages, memory CD4 T cells, regulatory T cells, and CD8 T cells (Supplementary Figures 14A,B). Moreover, we calculated the Pearson correlation coefficients of seven IRGs with various immune cells in the training and test cohorts (Supplementary  Figures 14C,D). Significant correlations were found in both matrices, especially for GZMB and CXCL3, which were significantly associated with CD8 T cells (r = 0.18, P < 0.01; r = −0.17, P < 0.01). When further restricted to immune cells for the TCGA cohort, we obtained similar results (Supplementary Figure 15).
Next, we extended the analysis to 49 immune checkpoint molecules, and the results for the training and test cohorts are shown in Supplementary Figures 16A,B, respectively. In total, we found that the expression of 13 immune checkpoints [CD274 (PD-L1), CTLA4, ICOS, BTLA, CD27, TNFRSF14, TNFRSF18, TNFRSF9, CD28, CD80, IDP1, LAG3, and TIGIT] was significantly downregulated in patients with high-risk scores in both cohorts (Supplementary Figures 16C,D; FDR < 0.05). In addition, we calculated the Pearson correlation coefficients of seven IRGs with various differentially expressed immune checkpoints (Supplementary Figure 17). Good correlations were found in both matrices, especially for GZMB, LGR5, and PLAU, which were significantly associated with many differentially expressed immune checkpoints. Surprisingly, CD274 was significantly positively correlated with GZMB, PLAU, CXCL3, and IL1R2, while it was significantly negatively correlated with LGR5 in the training and test cohorts. When further restricted to immune checkpoints for the TCGA cohort, we obtained similar results (Supplementary Figure 18).
In terms of guidance for immunotherapy, we studied the application of the signature in cutaneous melanoma (TCGA-SKCM). The results showed that IGBRS can also be used to predict the prognosis of patients with cutaneous melanoma (Supplementary Figures 19A,B). The effects of immunotherapy in cutaneous melanoma patients with high-risk score were weaker than those in patients with low-risk scores (5-year OS: 57% vs. 88%, Supplementary Figure 19C), and ROC analysis showed that our signature was used to predict the 4-year survival rate of immunotherapy patients with an AUC value as high as 0.700 (Supplementary Figure 19D). In addition, when our signature was applied to the cohort of anti-CTLA4 immunosuppressant therapy (GSE63557) and the anti-MAGE-A3 immunosuppressant therapy cohort (GSE35640), we found that the risk score of immunotherapy responders was significantly lower than that of non-responders (Supplementary  Figures 19E,F).

Tumor Somatic Mutation Correlates of the IGBRS
We analyzed the differences with respect to the distribution of somatic mutations between low-and high-risk score groups in the TCGA-COAD cohort using the maftools package. Supplementary Figure 20A shows the overall gene mutations in the TCGA-COAD cohort. Somatic mutations appeared in 390 (97.74%) of 399 samples. Figures 7A,B shows the top 30 gene mutations in the high-and low-risk score groups, respectively. The commonly mutated genes TP53, BRAF, and KRAS were altered in 311 (77.94%) of 399 samples (Supplementary Figure 20B), especially TP53 (Figures 7C,D). The Kaplan-Meier curves for OS suggested that there was no statistical difference in survival between the high-and low-TMB groups (log-rank test P-value = 0.56, Figure 7E). As shown in Figure 7F, the low-risk score group presented more extensive TMB than the high risk-score group (P = 0.0134). By contrast, compared with the high-TMB group, the low-TMB group had a higher risk score (Figure 7G, P = 0.0037). Moreover, we explored the mutation landscape of the seven IGRs in the IGBRS (Supplementary Figure 20C). Their mutation rates were very low and were well suited as diagnostic or prognostic biomarkers.

Identification of IGBRS as an Independent Risk Factor for CRC Patients
With the training cohort, we performed univariate and multivariate Cox regression analyses for each clinical factor ( Table 1) and screened factors with P < 0.05, which included the TNM stage, KRAS status, and the IGBRS. The multivariate modeling results showed that several patient and disease variables were significantly associated with survival, including a KRAS mutation (HR = 1.54, P = 0.0245), stage III (HR = 10.17, P = 0.0239), stage IV (HR = 12.78, P = 0.0178), and increasing risk score (HR = 2.71, P < 0.0001). The IGBRS was a significant independent predictor of RFS in CRC patients. Additionally, the results for the test (Supplementary Table 4, HR = 3.84, P < 0.0001) and validation (Supplementary Table 5, HR = 1.58, P = 0.0009) cohorts confirmed the independent predictive value of the IGBRS for RFS in CRC patients. In addition, as shown in Table 1, the IGBRS was an independent factor for OS (HR = 2.42, P < 0.0001). As expected, in the GSE17538 (Supplementary Table 6, HR = 2.18, P = 0.0011) and TCGA (Supplementary Table 7, HR = 2.00, P = 0.0065) cohorts, the IGBRS was still found to be an independent factor for OS after multivariate Cox regression analysis. In addition, the results for the test cohort showed that the IGBRS was an independent factor for DSS (Supplementary Table 4, HR = 3.09, P < 0.0001). Moreover, we performed a prognostic meta-analysis to investigate the comprehensive prognostic value across all groups. The results indicated that the IGBRS was a significant risk factor for RFS (n = 1,386, combined HR = 2.85, 95% CI = 1.77-4.58, P < 0.01) and OS (n = 1,375, combined HR = 2.48, 95% CI = 1.72-3.58, P < 0.01) in CRC patients (Supplementary  Figures 21A,B).

Gene Set Enrichment Analysis (GSEA)
We also identified pathways that were up-and downregulated between the high-and low-risk score groups by running a GSEA of the adjusted expression data for all transcripts. GSEA identified that, compared with the high-risk score group, the genes highly expressed in the low-risk score group were significantly enriched in 15 pathways including the cell cycle, MMR, nucleotide excision repair, cytokinereceptor interaction, ECM-receptor interaction, cell adhesion molecules (CAMs), and DNA replication (Supplementary  Figures 22A,B,D). Besides, compared with the patients with lowrisk scores, the genes with significantly higher expression in the high-risk score group were mainly concentrated in pathways involved in cancer, the MAPK signaling pathway, and the P53 signaling pathway (Supplementary Figure 22C).

DISCUSSION
Several lines of evidence have demonstrated that IRGs have indispensable roles in inflammation and innate immunity and antitumor effects. Considering the close relationships between the immune system and the occurrence and development of CRC (Ferrone and Dranoff, 2010), we believe that it is imperative to develop biomarkers related to RFS and OS for CRC. To this end, we analyzed the relationships between IGBRS and CRC patient recurrence for the first time. Using large-scale datasets from multiple centers, we identified an IGBRS based on seven IRGs that were significantly associated with RFS and OS in CRC patients. In addition, the prognostic value of the IGBRS in different molecular and clinical subtypes was verified, and this feature was related to different immune and somatic mutation landscapes. Ultimately, we concluded that the signature is an independent risk factor in CRC patients.
In the present study, we confirmed that IRGs were strongly correlated with recurrence based on unsupervised cluster analysis of 566 CRC patients from the GSE39582 dataset. Seven IRGs (BMP4, CXCL3, IL1R2, LGR5, GZMB, PTGDR, and PLAU) were applied to construct a recurrence signature for CRC. To avoid false positives in the sequencing data, the 867 CRC patients in the qRT-PCR, test, and validation cohorts were used to validate the stability of the IGBRS. In addition, through analysis of the training, qRT-PCR, validation, and TCGA cohorts, we found that the IRG combination accurately predicted OS and DSS in CRC patients. Finally, we conducted a meta-analysis based on four cohorts with RFS and OS data to fully understand the value of IGBRS in assessing the prognosis of CRC.
We have developed and tested a new prognostic IRG signature, which can carry out risk stratification and predict the OS and recurrence of patients with CRC more accurately than the AJCC stage. The effect of chemotherapy in high-risk patients (3-year OS: 35% in the training cohort and 5-year OS: 43% in the qRT-PCR cohort) was significantly weaker than that in the low-risk group (3-year OS: 70% in the training cohort and 5-year OS: 85% in the qRT-PCR cohort). Among different chemotherapy subtypes, the RFS and OS after single use of capecitabine, 5FU, FOLFOX, FUFOL, or XELOX (combined capecitabine with oxaliplatin) scheme in high-risk patients were significantly worse than those in low-risk patients. In terms of guidance for immunotherapy, the results showed that the effects of immunotherapy in cutaneous melanoma patients with highrisk score were weaker than those in patients with low-risk scores (5-year OS: 57% vs. 88%). In addition, we found that the risk score of immunotherapy responders was significantly lower than that of non-responders in the cohort of anti-CTLA4 immunosuppressant therapy (GSE63557) and the anti-MAGE-A3 immunosuppressant therapy cohort (GSE35640). This has great significance for the guidance of clinical stratified treatment. Patients with low-risk scores can appropriately consider adjuvant chemotherapy and immunotherapy before and after operation. Besides, patients with high-risk scores need to consider total surgical resection and radiotherapy more actively and recheck more frequently to monitor the recurrence for further early treatment.
Pathological staging is currently the most important factor in the clinical evaluation of CRC prognosis. Different stages of CRC have different immune statuses and respond differently to immunotherapy (Brahmer et al., 2012). We explored the applicability of our signature in different staging subgroups. As expected, the signature performed well in both early and advanced stages. Age, sex, and location were also important factors affecting the prognosis of patients with CRC (Loupakis et al., 2015;Thrumurthy et al., 2016;Guo et al., 2019;Siegel et al., 2020). Patients with CRC who are aged younger than 50 years have higher 5-year relative survival rates than their older counterparts at every stage of diagnosis (Siegel et al., 2020). Besides, the signature performed convincingly well in all subgroups including males and females, older (age ≥65 years) and younger (age <65 years) patients, and distal and proximal metastases. These findings indicate that our signature may help to identify high-risk CRC patients independent of other clinical factors affecting prognosis and may better guide clinical treatment.
TP53, KRAS, and BRAF are commonly mutated oncogenes in CRC (Russo et al., 2005;Sinicrope et al., 2015;Guo et al., 2019). We analyzed the performance of IGBRS among patients with different TP53, KRAS, and BRAF mutation statuses. The results suggest that the risk stratification effect of IGBRS is independent of these oncogenes. Then, we analyzed the differences with respect to the distribution of somatic mutations between the lowand high-risk groups in the TCGA-COAD cohort. The results indicated that the low-risk group presented more extensive TMB than the high-risk group.
CIN, CIMP, and microsatellite-unstable (MSI) are the three molecular types of CRC. MMR refers to the repair mechanism to restore the normal nucleotide sequence in DNA molecules containing mismatch bases, which is mainly used to correct the mismatched base pairs in the double helix of DNA. MSI refers to the phenomenon that when the function of MMR is abnormal, the replication errors of microsatellites cannot be corrected and accumulate, changing the length or composition of microsatellite sequences. We analyzed the performance of the IGBRS among patients with different CIN, MMR, and CIMP statuses. The results showed that compared with the CIN-positive group, the CIN-negative group had a higher risk score and a higher proportion of high-risk patients and showed better RFS and OS. Among the pMMR and dMMR subgroups, the patients in the high-risk group had shorter RFS and OS times than those in the low-risk group. Patients in the dMMR subgroup had lower risk scores than those in the pMMR subgroup. In addition, higher risk scores and high-risk patients were more concentrated in the CIMP-positive group than in the CIMP-negative group, which was associated with worse RFS and OS. A number of studies showed that lower 5-year survival rates were observed in microsatellite-stable cancer patients with CIMP-low or CIMPhigh status than in patients with no CIMP (Barault et al., 2008;Sinicrope et al., 2015), which strongly confirmed the ability of the IGBRS to conduct risk stratification of CRC patients.
The present study found that this signature was a good predictor of the outcomes of several cohorts and subgroups; therefore, we examined possible underlying mechanisms. The degree of immune infiltration significantly affects the prognosis of CRC (Mao et al., 2018). In this study, the high-risk group had significantly lower immune scores than the low-risk group, and the patients with high immune scores tended to have better prognoses. The increased accumulation of CD8 T cells, regulatory T cells, and proinflammatory macrophages (M1) and the reduced accumulation of immunosuppressive macrophages (M2) indicated better prognosis for CRC (Funada et al., 2003;Correale et al., 2012). High-risk patients were characterized by high proportions of M2 macrophages and low proportions of activated memory CD4 T cells, CD8 T cells, regulatory T cells, and M1 macrophages, suggesting that the IRGs included in our combination may affect prognosis by interacting with infiltrating immune cells. Therefore, we also calculated the Pearson correlation coefficients of seven IRGs with various immune cells. Good correlations were found in both matrices, especially for GZMB, which is significantly associated with a variety of immune cells.
We further found that the two groups of patients had different characteristics of immune checkpoints. The CD274 (PD-L1), CTLA4, ICOS, BTLA, CD27, TNFRSF14, TNFRSF18, TNFRSF9, CD28, CD80, IDP1, LAG3, and TIGIT levels were significantly downregulated in patients with high-risk scores. The PD-1 plus CTLA4 blockade is highly effective in advancedstage, dMMR CRC, yet not in pMMR tumors (Chalabi et al., 2020). Responses to immune checkpoint inhibitors correlated with PD-L1 expression (Salem et al., 2018). Studies have demonstrated that relatively high PD-L1 expression in cancer cells was associated with a good prognosis in CRC patients (Wyss et al., 2019;Noh et al., 2020). The high expression of PD-L1 and CTLA4 in the low-risk group not only demonstrated that the IGBRS is an efficient classifier of risk stratification in CRC and closely related to tumor immunity but also suggested that the IGBRS may be a reference for the classification of CRC patients treated with immune checkpoint inhibitors. To further verify this possibility, we analyzed the correlations of seven IRGs and immune checkpoints for differential expression in each cohort. CD274 was significantly positively correlated with GZMB, PLAU, CXCL3, and IL1R2, while it was significantly negatively correlated with LGR5 in the training and test cohorts.
All of the genes in the IGBRS were found to be involved in CRC progression (Kalmár et al., 2015;Mar et al., 2015;D'Eliseo et al., 2016;Morgan et al., 2018;Zhou et al., 2018;Liao et al., 2019;Lin et al., 2019). Bone morphogenetic protein 4 (BMP4) encodes a secreted ligand of the transforming growth factor-beta (TGFβ) superfamily of proteins. Deng et al. (2009) and Zhou et al. (2018) found that BMP4 knockdown could ameliorate CRC cell migration, and invasion and overexpression of BMP4 enhance the invasiveness of Smad4-deficient human CRC cells. The findings of Yokoyama et al. (2017) suggest inhibition of autocrine BMP4 as a candidate treatment strategy for CRC. In our study, BMP4 expression (HR = 1.213, P = 0.02) was an independent risk factor for the prognosis of patients with CRC and was positively correlated with the infiltration abundance of M2 macrophages in tumor tissue (r = 0.28, P < 0.0001). The above findings suggest that BMP4 may be secreted by M2 macrophages and participate in the malignant biological behavior of CRC cells.
The plasminogen activator urokinase (PLAU) encodes a secreted serine protease that converts plasminogen to plasmin and acts as one of the TGF-β downstream factors. Lin et al. (2019) reported that inhibition of PLAU can suppress CRC cell proliferation and progression. In the present study, PLAU (HR = 1.357, P = 0.001) also acted as an independent risk factor for the prognosis of patients with CRC and was negatively correlated with the abundance of infiltrating regulatory T cells in tumor tissue (r = −0.26, P < 0.0001). Moreover, PLAU was positively correlated with the expression of many immune checkpoint genes, including CD274 and CALT4, and may act as a key molecule for tumor cells to escape immune surveillance.
Granzyme B (GZMB) (HR = 0.742, P < 0.001), whose encoded preproprotein is secreted by natural killer cells and cytotoxic T-lymphocytes, also acted as an independent protective factor for the prognosis of patients with CRC. Koelzer et al. (2017) found that the immune-activated phenotype was associated with high counts of intratumoral CD8 cytotoxic T-lymphocytes (P = 0.007) and the expression of the immune effector molecule GZMB (P < 0.001). The increased infiltration of cytotoxic T-lymphocytes is the key to tumor immune rejection. Our study suggests a significant positive correlation between the expression of GZMB and the abundance of infiltrating cytotoxic T-lymphocytes in each cohort (r = 0.21, P < 0.001). The significant downregulation of GZMB expression in CRC patients with low risk scores suggests that GZMB may kill tumor cells and inhibit their malignant biological behavior to improve the prognosis of CRC patients by activating tumor immune rejection.
The antimicrobial gene C-X-C motif chemokine ligand 3 (CXCL3) encodes a member of the CXC subfamily of chemokines. The encoded protein is a secreted growth factor that signals through the G protein-coupled receptor CXC receptor 2 and plays a role in inflammation and as a chemoattractant for neutrophils. Liao et al. (2019) demonstrated that KRAS *mediated repression of IRF2 results in high expression of CXCL3, which binds to CXCR2 on myeloid-derived suppressor cells and promotes their migration to the tumor microenvironment, which drives immune suppression and immune therapy resistanc in CRC. Our results suggest a positive correlation (r = 0.31, P < 0.0001) between the mRNA expression levels of CXCL3 and CD274, which encodes PD-L1, and a negative correlation (r = −0.17, P < 0.0001) between the mRNA expression levels and the abundance of infiltrating cytotoxic T-lymphocytes. These results suggest that CXCL3 may transmit inhibitory signals through the CXCL3-CXCR2 axis, reduce the proliferation of CD8 T cells in lymph nodes, and promote immune escape. CXCL3 may act as one of the targets to enhance immune efficacy.
The protein encoded by LGR5 is a leucine-rich repeat containing receptor (LGR) and member of the G proteincoupled, seven-transmembrane receptor (GPCR) superfamily. The encoded protein is a receptor for R-spondins and is involved in the canonical Wnt signaling pathway (Glinka et al., 2011). It has been demonstrated that LGR5-positive cancer cells functionally act as stem cells in human CRCs (Fumagalli et al., 2020). Jang et al. (2018) found that LGR5 overexpression attenuates proliferation, migration, and colony formation in CRC cells. Besides, LGR5 functions as a tumor suppressor in the late stages of CRC progression and is an independent prognostic marker for better clinical outcomes in CRCs. In the present study, LRG5 (HR = 0.855, P = 0.003) also acted as an independent protective factor for the prognosis of patients with CRC.
However, the contributions of prostaglandin D2 receptor (PTGDR) and interleukin 1 receptor type 2 (IL1R2) to CRC immune microenvironment remodeling remain unknown. In the present study, these genes showed strong correlations with tumor immune cell infiltration and immune checkpoints, but these correlations require further exploration. Although the prospects for IRG signatures are promising, they also have certain limitations. On the one hand, all cohorts were retrospective, and this risk scoring system still needs to be prospectively verified. On the other hand, because of the high spatial heterogeneity of the tumor immune microenvironment, the relationships between the IGBRS and immune cell infiltration and immune checkpoints are based on estimates of tumor characteristics, which may lead to errors. Further research is needed to verify our findings.

CONCLUSION
In general, we developed and tested a new recurrence immunerelated gene signature for CRC. Our research provides new insights into the link between immunotherapy and CRC. This IGBRS may help clinicians develop personalized treatment plans, especially when choosing which patients will benefit from immunotherapy, and may improve the survival of CRC patients.

DATA AVAILABILITY STATEMENT
The original contributions generated for this study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Ethics Committee/Institutional Review Board of the Cancer Institute/Hospital, Peking Union Medical College and Chinese Academy of Medical Sciences (approval no. NCC2013RE-025). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
RM, TZ, and YL conceived the project and designed the experiments. ZW, FY, and RM carried out the experiments. RM and FY contributed equally to this work. YL, TZ, and RM wrote the manuscript. RM, CX, and QL carried out the statistical analysis and gave assistance in collecting tissue samples. TZ contributed to manuscript revision. All authors provided suggestions during manuscript preparation and read the final version.

FUNDING
This work was supported by grants from the National Natural Science Foundation of China (81502075) and the Foundation of Science and Technology of Sichuan Province (2019YJ0635). The funders had no role in study design and implementation.