Immune infiltration and immune gene signature predict the response to fluoropyrimidine-based chemotherapy in colorectal cancer patients.

Fluoropyrimidine-based chemotherapy is an essential component of systemic chemotherapy for colorectal cancer (CRC). The immune response is implicated in chemotherapy-induced cytotoxicity. Here, we reported an immune risk (Imm-R) model for prognostic prediction in patients receiving fluoropyrimidine-based chemotherapy. Gene expression profiles and corresponding clinical information were collected from four data sets and divided into training set (n = 183) and validation set (validation set1: n = 34; validation set2: n = 99). The composition of 22 tumor-infiltrating immune cells (TIICs) populations was characterized with the CIBERSORT deconvolution algorithm. A prognostic Imm-R model for predicting overall survival was established by performing least absolute shrinkage and selection operator (LASSO) penalized COX regression analysis. T follicular helper cells and M0 macrophages were associated with better survival, while eosinophils were associated with worse survival. TIICs signature was constructed based on the above three immune cell types. Furthermore, a Imm-R model was created by integrating TIICs signature with immune-related genes (IRGs), which effectively in distinguishing CRC patients with poorer prognosis. The Imm-R model was associated with activation of the TGF-beta signaling and suppression of DNA damage. Results of this research provide new insights into the role of immunity for in fluoropyrimidine-based chemotherapy as well as a useful tools to predict the outcome of CRC patients receiving fluoropyrimidine-based chemotherapy.

Introduction

Colorectal cancer (CRC) is the third most common cancer and a leading contributor to cancer-related mortality worldwide.[1] Although significant advances have been made in screening, diagnosis, and treatment of CRC, clinical outcome of CRC patients remains unsatisfactory. Radical resection is currently the only curative strategy for CRC. However, approximately 50% of CRC patients who underwent curative resection developed recurrence within two years after surgery and half of these recurrences proved to be fatal.[2] To prevent recurrence in and prolong survival, fluoropyrimidine-based chemotherapy, is recommended for all suitable patients with stage III disease or stage II disease with high-risk features[3] after curative resection. Fluoropyrimidine-based combination chemotherapy, such as FOLFOX, CapeOX, or FOLFIRI, with/without biologics, is also recommended as the first-line treatment.[4] However, the therapeutic benefit to patients is frequently compromised by the development of chemoresistance.[5] Adjuvant chemotherapy improves survival in only 3% of patients with stage II CRC, and up to 15%–20% among those with stage III CRC.[6] Therefore, identifying high-risk CRC patients, those who received fluoropyrimidine-based therapy but with poor responses, may contribute to developing individualized treatment.[7]

Tumor-infiltrating immune cells (TIICs) are of various consist of the tumor immune microenvironment[8] TIICs were composed of various types of immune cells, like T cell, macrophages, neutrophils etc. Different subtype of TIICs may exhibit either immune-promoting or immune-suppressing ability under various situation.[9] The composition of the TIICs is tightly regulated by the immune interactions between host and tumor,[10] and reflects the intrinsic characteristics of tumor.[11] Increasing evidence indicates that composition and abundance of TIICs are potential predictors of patient survival.[12-14] For example, high levels of infiltrating CD8 + T cell is associated with better prognosis, while regulatory T (Treg) cells predict poor clinical outcomes in CRC patients.[15] A CRC classification based on TIICs has been developed, and five subgroups of CRC with distinct survival patterns have been identified.[16] Therefore, the comprehensive identification of prognostic TIICs may complement current strategies used for the prediction of CRC prognosis.

Most previous studies have used immunohistochemical approaches to assess the landscape of TIICs. Recently, CIBERSORT algorithm was developed to enumerate immune cell subtypes. It generates the proportion of TIICs based on changes of immune-related and other genes in expression profiles.[17] Due to the excellent performance of CIBERSORT, its application in studying cell heterogeneity has aroused increasing attention.[18-20] For example, Xiong et al analyzed the clinical significant of TIICs in CRC by using CIBERSORT[16] and Narayanan et al performed CIBERSORT algorithm to identify TIICs associated with OS of microsatellite unstable CRC patients.[21] However, The impact of TIICs on the prognosis of patients treated with fluoropyrimidine is not yet conclusive.

The immune response is regulated by a complex immune molecular network that includes chemokines, cytokines, major histocompatibility complex (MHC) molecules, as well as co-stimulatory, immunomodulatory, and cytotoxicity-mediating molecules.[22] The genes encode those molecules are referred as immune-related genes (IRGs).[23] Besides the composition and abundance of TIICs, the expression of IRGs is associated with prognosis in CRC patients as well.[24] For example, higher expression of PD-1 and PD-L1 is correlated with better prognosis.[25] Furthermore, expression of genes encoded MHC class II complex could predict metastasis in CRC patients.[26] An increasing number of studies also implicate the immune response in chemotherapy-induced cytotoxicity.[27,28]

In the current study, we aimed at illustrating the landscape of TIICs and screen TIICs associated with prognosis in CRC patients receiving fluoropyrimidine-based chemotherapy. We also construct an immune risk (Imm-R) model by integrating TIICs- and IRGs to identify potential high-risk CRC patients for fluoropyrimidine-based chemotherapy.

Results

Cohort characteristics

Firstly, 865 patients from GSE39582 (585 cases), GSE103479 (156 cases), GSE72968 (68 cases), and GSE72969 (56 cases) were selected. The follow-up information and chemotherapy regimens of each case was collected. After filtering out patients who received preoperative treatment and whose microarry data did not pass the CIBERSORT quality control step (p value ≥ 0.05),[29] 316 CRC patients receiving fluoropyrimidine-based chemotherapy, included 182 male patients and 134 female patients, were finally included in the study (Table 1 and Figure 1). Among them, 199 cases were in stage II–III (62.97%) and 117 cases in stage IV (37.03%). The median follow-up time was 46 months. 172 patients (54.43%) experienced disease recurrence or progression during follow-up. The training set consisted of 183 cases from GSE39583. Validation set1 included 34 cases from GSE103479, who were with stage II–III disease and received 5-fluorocrail (5-FU) based chemotherapy. Validation set2 was the combination of 61 cases from GSE72968 and 38 cases from GSE72969, who were stage IV and received fluoropyrimidine-based combination treatments (Table 1).

Table 1.

Clinical features of patients included in the study.

	Training set	Validation set1	Validation set2		All dataset
Features	GSE39582(n = 183)	GSE103479(n = 34)	GSE72968(n = 61)	GSE72969(n = 38)	Total number(n = 316)
Age
<60 years	70 (38.25)	11(32.35)	27(44.26)	18(47.37)	126(39.87)
≥60 years	113(61.75)	23(67.65)	34(55.74)	20(52.63)	190(60.13)
Gender
Female	82(44.81)	16(47.06)	22(36.07)	14(36.84)	134(42.41)
Male	101(55.19)	18(52.94)	39(63.93)	24(63.16)	182(57.59)
Pathological T category
T1-T2	5(2.73)	2(5.88)	2(3.28)	2(5.26)	11(3.48)
T3-T4	168(91.80)	32(94.12)	38(62.30)	31(81.58)	269(85.13)
Unknown	10(5.46)	0	21(34.43)	5(13.16)	36(11.39)
Pathological N category
N0	52(28.42)	16(47.06)	5(8.20)	6(15.79)	79(25.00)
N1-3	131(71.58)	18(52.94)	35(57.38)	27(71.05)	211(66.77)
Unknown	0	0	21(34.43)	5(13.16)	26(8.23)
M category
M0	165(90.16)	34(100)	0	0	199(62.97)
M1	18(9.84)	0	61(100)	38(100)	117(37.03)
TNM stage
II–III	165(90.16)	34(100)	0	0	199(62.97)
IV	18(9.84)	0	61(100)	38(100)	117(37.03)
Chemotherapy regimens
5FU	59(32.24)	34(100)*	0	0	93(29.43)
5FU+LV	63(34.43)	0	0	0	63(19.94)
FOLFIRI	6(3.28)	0	39(63.93)	28(73.68)	73(23.10)
FOLFOX	12(6.56)	0	22(36.07)	6(15.79)	40(12.66)
FUFOL	41(22.40)	0	0	0	41(12.97)
other	2(1.09)	0	0	4(10.53)	6(1.90)
Tumor location
Proximal	65(35.52)	16(47.06)	12(19.67)	12(31.58)	105(33.23)
Distal	118(64.48)	18(52.94)	49(80.33)	26(68.42)	211(66.77)
Recurrence or progression
No	109(59.56)	25(73.53)	6(9.84)	4(10.53)	144(45.57)
Yes	74(40.44)	9(26.47)	55(90.16)	34(89.47)	172(54.43)
Median follow-up (sd), months	56.00(37.93)	62.09(27.63)	22.05(17.85)	26.04(19.65)	46.15(36.03)

* : Fluorouracil based

Figure 1.

The flow diagram of this study. In brief, four colorectal cancer (CRC) microarray datasets in the GEO database were included in the study at first. After filtering out patients who received preoperative treatment and whose microarray data did not pass the CIBERSORT quality control step, 183 cases in GSE39582 (Training set), 34 cases in GSE103479 (Validatiaon set1), 61 cases from GSE72968, and 38 cases from GSE72969 (combined as Validation set2) were selected. The prognostic effects of the 22 subpopulations of tumor-infiltrating immune cells (TIICs) were analyzed and three TIICs were selected to create the TIICs signature. Further, TIICs and immune-related genes were integrated to create the immune risk (Imm-R) model.

Correlation between the immune infiltration landscape and clinical features in CRC

TIICs in the tumor microenvironment modulate the antitumor response and represent attractive therapeutic targets.[30] Using the CIBERSORT algorithm, we systematically evaluated the abundance of 22 subpopulations of TIICs in the included patients of training set. Significant variation in the proportion of TIICs was observed among individuals, reflecting the varied intrinsic individual immune characteristics (Figure S1A). Correlation analysis showed that subpopulations of TIICs showed a weak to moderate correlation with each other (Figure S1B). Next, we analyzed the the corelation of each TIICs subgroup and clinical features, likes TNM stages, TP53 status, KRAF/BRAF status, RFS and OS, of all patients in the training set (Figure S1C–E). As a result, four TIICs were identified to be associated with TNM stages (Figure 2a). CD8 + T cell was down-regulated with increasing TNM stage (p = .036), while Monocytes (p = .045), Eosinophils (p = .0061), and activated Mast cells (p = .032) were upregulated with increasing TNM stage (Figure 2a). Correlation analysis also showed that proportion of CD8 + T cell was negatively correlated to the three other cells (Figure 2b)

Figure 2.

Correlation between immune infiltration landscape and clinical features in CRC. (a) The abundance of four tumor-infiltrating immune cell subpopulations (TIICs) which were significantly associated with TNM stage in the training set. CD8 + T cells (T cell CD8) was downregulated while Monocytes, Eosinophils and Mast cells activated were upregulated with increasing TNM stages (p < .05). (b) Correlation matrix of the CD8 + T cells, Monocytes, Eosinophils, and Mast cells activated. (c) Six TIICs (T cells CD4 memory resting, T cells follicular helper, T cells gamma delta, NK cells activated, Monocytes, and Macrophages M1) were associated with MMR status significantly (p < .05). dMMR, mismatch repair deficient; pMMR mismatch repair proficient. (d and e) T cells CD4 memory activated (d) and Macrophages M0 (e) were associated with KRAS mutation status (p < .05). WT-KRAS, wild-type KRAS; MT-KRAS, mutant KRAS.

We comprehensively explored the relationship between TIICs and clinical features in CRC. Mismatch-repair-deficient (dMMR) CRC represents a strong indication for immunotherapy.[31] We compared the abundance of the 22 TIICs between dMMR and mismatch-repair-proficient (pMMR) CRC (Figure S1D). We found that the abundance of T follicular helper cells (Tfh), gamma delta T cells, NK cell activated and M1 macrophages were significantly increased in the dMMR CRC group compared with pMMR (p< .05), while that of resting memory CD4 + T cells and Monocytes were increased in the pMMR CRC group (Figure 2c). Recent evidence indicates that the KRAS mutation might also shape tumor immunity.[32,33] In this study, we found that the abundance of activated memory CD4 + T cells and M0 macrophages were lower in KRAS-mutated CRC than in wild-type CRC (Figure S1E, Figure 2d and e). These results suggested that TIICs represented some of the intrinsic characteristics of CRC patients with different clinical features.

Development and validation of the tumor-infiltrating immune cell (TIICs) signature

To further explore the potential of using TIICs as an indicator of the therapeutic response to cytotoxic regimens, we firstly performed Cox regression analysis in the training set to screen TIICs associated with overall survival (OS). 13 of 22 TIICs were associated with OS (Figure S2 and Figure 3a). Next, we performed multivariate Cox regression analysis to identify independent prognostic TIICs and the results showed that Tfh, plasma cells, activated mast cells, M0 macrophages, eosinophils, and memory B cells were potential independent prognostic TIICs for chemotherapy (Figure S3A and B). The association between the above six TIICs and OS were validated in both validation set1 and validation set2 and confirmed that Tfh and M0 macrophages were associated with better OS while eosinophil was associated with worse OS in the three datasets (Figure 3(b and c)). Therefore, collection of Tfh, M0 macrophages and eosinophil were defined as the TIICs signature to predict CRC patients response to fluoropyrimidine-based chemotherapy. We also performed a subgroup analysis to assess the application of TIICs signatures in stage II–III patients in the training set (Figure S4). Tfh and M0 macrophages were associated with better survival, while eosinophils were associated with worse survival in the stage II–III patients of training set (Figure S4A).

Figure 3.

Development and validation of the tumor-infiltrating immune cell signature. (a)–(c) High abundance of T cells follicular helper (Tfh)and M0 macrophages (macrophages M0) were associated with better OS, while eosinophils was associated with worse OS in the training set (a), validation set1 (b) and validation set2 (c). (d)–(f) Kaplan–Meier curves for RFS/PFS of patients with high- and low TIICs risk scores in the training set (d), validation set1 (e), and validation set2 (f). (g)–(i) Kaplan–Meier curves for OS of patients with high- and low TIICs risk scores in the training set (g), validation set1 (h), and validation set2 (i). (j)–(l) ROC curve for measuring the predictive value of the TIICs signature for OS and RFS/PFS in the training set (j), validation set1 (k), and validation set2 (l).

Based on the TIICs signature, we constructed a method to calculate scores of TIICs signature (TIICs score) to predict the prognosis of patients receiving fluoropyrimidine-based chemotherapy. The TIICs score was calculated as follows: (−0.990 × abundance of Tfh)+(−0.731 × abundance of Macrophages M0) + (0.655 × abundance of Eosinophils). Using the median of TIICs scores as cutoff value, the TIICs scores successfully distinguished CRC patients into high-risk or low-risk groups. Patients with high risk had significantly poorer relapse-free survival (RFS) and OS compared with those with low risk in the training set (RFS, p=.013 Figure 3g; OS, p< .001; Figure 3d). With threshold set as in training set, TIICs scores in validation set1 and validation set2 yielded similar results. Compared with those with low risk, high-risk patients had poorer RFS and OS in the validation set1 (RFS, p=.054 Figure 3e; OS, p = .043, Figure 3h), or poorer progression-free survival (PFS) and OS in the validation set2 (PFS, p=.043, figure 3f; OS, p= 5.9 × 10−4, Figure 3i). The AUC under ROC curve of TIICs scores for predicting OS in the training set, validation set1 and 2 were 0.691, 0.695, and 0.755 respectively (Figure 3j–l). The AUC of the TIICs scores for predicting RFS in training set and the validation set1 was 0.602 and 0.707; for PFS in validation set2 was 0.670 (Figure 3j–l). Application of the TIICs score successfully generated similar results in stage II–III subgroup of training set (Figure S4B-D).

Development and validation of the immune risk (Imm-R) model by intergrating TIICs and IRGs

Although the TIICs scores seemed to distinguish patients with high-risk well, the AUC of ROC curve did not reach 0.8, which indicated limited precision. Immune-related genes (IRGs) were reported to orchestrate tumor-associated immune responses.[34] The integration of the TIICs signature and IRGs signature may enable more comprehensive assessment of immune status and more precise prognostic prediction. By integrating two immune gene lists from ImmPort[35] and TISIDB,[36] we obtained a list of 2404 IRGs. Univariate Cox regression analysis was performed for preliminary screening of survival-associated IRGs in the training set. Based on the threshold value of p value < 0.05, we identified 61 IRGs significantly related with overall survival. Then, these 61 IRGs and the TIICs scores were input into the LASSO-Cox regression model for feature selection. Under penalizing conditions, 13 IRGs and TIICs scores which with nonzero coefficients were selected (Figure 4(a,b)). Next, we performed stepwise multivariate Cox regression analysis based on the features selected by the LASSO-Cox regression model in order to construct the Imm-R model (Table S1 and Figure 4c). Ultimately, a Imm-R model consisting of seven IRGs (PSMD3, CCL22, FGF19, LGR5, SKAP1, TACSTD2, S100A4) and TIICs scores, was constructed. Hazard ratios (HR), 95% confidence intervals of HR, and p-values in univariate and multivariate Cox regression of the seven IRGs and TIIC scores for Imm-R model were summarized in Table S1 and Figure 4c. In the Imm-R model, an immune risk score (Imm-R score) was generated using the formula:

Figure 4.

Development and validation of the immune risk (Imm-R) model. (a) Tuning parameter (λ) selection in the LASSO-Cox regression model was performed using 10-fold cross-validation via the 1 standard error for the minimum criteria. The black vertical lines were plotted at the optimal λ based on the minimum criteria and 1 standard error for the minimum criteria. (b) The LASSO coefficient profiles of the 61 immune-related genes and TIICs signature. A coefficient profile plot was produced versus the log (λ). (c) Parameters of Imm-R model. p values of features were indicated by the color scale presented on the side. Horizontal bars represent 95% confidence intervals. (d) and (e) Kaplan–Meier curves for RFS (d) and OS (e) of patients with high- and low Imm-R scores in the training set. (f) and (g) Kaplan–Meier curves for RFS (f) and OS (g) of patients with high- and low Imm-R scores in the validation set1. (h) and (i) Kaplan–Meier curves for PFS (h) and OS (i) of patients with high- and low Imm-R scores in the validation set2. (j)–(l) ROC curves for measuring the predictive value of the Imm-R model in the training set (j), validation set1 (k), and validation set2 (l). (m) and (n) Construction and analysis of risk scores in the training set and validation set2. Top panels: the risk scores of individual patients. Middle panels: the survival status and survival times of the patients distributed by risk score. Bottom panel: heatmap of the levels for the eight predictive factors distributed by risk score.

Imm-R score = (−0.939× PSMD3)+(−1.304× CCL22)+(0.319× FGF19)+(−0.228× LGR5)+(0.665× SKAP1)+(0.182× TACSTD2)+(0.329× S100A4)+(0.809× TIICs scores).

In the training set, CRC patients could be divided into high-risk and low-risk groups by the median value of Imm-R scores. Survival analysis demonstrated that the Imm-R model was powerful to distinguish CRC patients with good or bad prognosis (Figure 4(d,e)). Patients with high-risk had significantly shorter RFS and OS compared with those with low risk in the training set (RFS, p = .0019; OS, p< .0001). Patients with high risk had shorter RFS in the validation set1 with a marginally significant p-value of 0.069 Figure 4f and shorter OS (p = .018, Figure 4g). In validation sets2, the PFS and OS were significantly shorter in the patients with high risk (PFS: p = .008, Figure 4h; OS: p< .0001 Figure 4i).

Time-dependent ROC analysis indicated that the AUC of Imm-R scores for predicting OS in the training set was 0.826, better than that of TNM stage (0.643) (Figure 4j). In the validation set1 and validation 2, the AUC of the Imm-R scores for predicting OS was 0.789 and 0.807, all significantly greater than TNM stage (Figure 4k and l). The Imm-R scores (upper panel) and corresponded OS (middle panel) of individuals in training set (Figure 4m) and validation set2 (Figure 4n) showed a trend that the OS roughly decreased along with the increasing of Imm-R scores. The expression levels of the seven IRGs and the TIICs scores were shown in the bottom panels (Figure 4m and n).

The efficiency of using Imm-R scores for predicting OS and RFS in stage II–III patients subgroup of training set was conincident with validation set2 (Figure S5A and B). The RFS and OS were significantly shortened in high-risk patients compared with low-risk patients (RFS p= .0062; OS, p< .0001). The ROC curve also indicated that Imm-R model had good discrimination in prognosis estimation, with AUC of 0.804 for OS (Figure S5C). These results indicated that the Imm-R model showed good discrimination ability in both training and validation sets.

We compared the prognostic effects of Imm-R model with common clinical parameters such as age, gender, BRAF status, KRAS status, MMR (Mismatch Repair) status, tumor location and TNM stage in the training set. Univariate Cox regression analysis indicated that KRAS status, TNM stage and Imm-R model were associated with OS (KRAS: p = .04; TNM stage: p = .018; Imm-R model p = 6.00 × 10−20; Table S2). Multivariate analysis confirmed that the Imm-R model was an independent prognostic factor (p = 4.61 × 10−11) for OS as shown in Table S2.

Nomogram construction and evaluation of calibration

The nomogram is a user‑friendly graphical regression model with excellent applicability in clinical settings. Recently, nomograms have been widely used in the diagnosis and prognostic prediction of a variety of cancers. To improve the usability of the Imm-R model, we constructed a nomogram to depict the Imm-R model better (Figure 5a). The nomogram included above eight features, and a point for each feature was assigned based on the scale on the top. The total score was defined as the sum of the points of the eight variables. By drawing a perpendicular line from the total point axis to the two-outcome axis, estimated three- and five-year OS probabilities could be obtained (Figure 5a). To assess the goodness-of-fit of the nomogram, we compared the predicted three- and five-year survival probabilities to the actual three- and five-year survival probabilities using calibration plots (Figure 5b–d). The calibration curve revealed good concordance between the predicted and observed probabilities in both of the training set, validation set1 and validation set2. These results proved that the Imm-R model nomogram had very appropriate calibration.

Figure 5.

Nomogram for predicting OS probabilities of CRC patients receiving fluoropyrimidine-based chemotherapy. (a) The nomogram for predicting the three- and five-year OS probabilities of CRC patients receiving fluoropyrimidine-based chemotherapy. Points are assigned for eight features. The score for each feature was calculated by drawing a line upward to the 'Points' line, and the sum of the eight scores was 'Total Points'. The total points on the bottom scales correspond to the predicted three- and five-year survival. (b–d) Calibration curves for predicting three- and five-year survival OS in training set (b), validation set1 (c), and validation set2 (d). X-axis: predicted survival produced by the nomogram; Y-axis: actual survival. Dashed lines represented an identical calibration model in which predicted OS approximates to actual OS. (e) Kaplan–Meier curves for OS of patients receiving 5-fluorouracil (5-FU); 5-FU and folinic acid (FUFOL); 5-FU, leucovorin, and irinotecan (FOLFIRI); and 5-FU, leucovorin, and oxaliplatin (FOLFOX), respectively.

The effect of Imm-R model in patients receiving different fluoropyrimidine-based chemotherapy

It was not uncommon that clinical outcomes in patients who of the same stage and received the same therapeutic regimens varied, suggesting the limitation of the current prognostic model. We explored whether the Imm-R model could distinguish between high- and low-risk patients receiving different fluoropyrimidine-based chemotherapy. We observed that this model could well stratify CRC patients who treated with 5-FU, FUFOL and FOLFIRI (p<.01). The Imm-R model could also distinguish high-risk patients who received FOLFOX treatment but only with marginal significance (p= .057, Figure 5e).

Biological annotation of Imm-R model

To explore the underlying biological mechanisms of the Imm-R model, we performed Gene Set Enrichment Analysis (GSEA) based on the microarray datasets from the training set. Patients were divided into two groups according to the median value of their Imm-R scores. The results observed that several pathways related to tumor growth and proliferation were significantly activated in those high-risk patients, included the TGF-beta signaling pathway, PI3K-Akt signaling pathway, Ras signaling pathway, and Wnt signaling pathway (Figure 6a), reflecting the active tumor proliferation in high-risk patients. Inversely, several pathways related to DNA damage and repair, such as the mismatch repair, base excision repair, and DNA replication pathways were significantly upregulated in low-risk patients (Figure 6b).

Figure 6.

Gene set enrichment analysis indicated significant pathways associated with the immune risk (Imm-R) model. (a) Upregulated KEGG pathways in high-risk patients. (b) Downregulated KEGG pathways in high-risk patients. (c) Volcano plots showed the differential expressed genes (DEGs) between HCT8/5-FU (5-FU resistant) and HCT8/WT (Wild type). Red dots represented the significantly upregulated DEGs (UP) in HCT8/5-FU group. Blue dots represented the significantly downregulated DEGs (DOWN) in HCT8/5-FU group. Black dots represented non DEGs (NS). (d) Enrichment diagram of TGF-beta signaling pathway. (e) Heatmap showed the expression of hub genes in the TGF-beta signaling pathway in HCT8/5-FU and HCT8/WT group. (f) The expression of eight hub genes in the TGF-beta signaling pathway in high-risk and low-risk patients in training set.

In order to further clarify the connection of Imm-R and the mechanism of 5-FU resistance, we performed differential expression analysis based on the GSE81005 datasets derived from wild type and 5-FU-induced resistant HCT8 human CRC cell lines.[37] In total, 801 genes were significantly upregulated and 399 genes were significantly down-regulated in the 5-FU resistant HCT8 cells (HCT8/5-FU) compared with wild type HCT8 cells (HCT8/WT) (Figure 6c). Pathway enrichment analysis indicated that TGF-beta signaling pathway was significantly upregulated in the HCT8/5-FU cells and the expression of hub genes in the TGF-beta signaling pathway were significantly increased in the HCT8/5-FU cells (Figure 6(d and e)).

We investigated the expression of hub genes in the TGF-beta signaling pathway between the high risk and low risk patients in the training set. As predicted, seven out of eight hub genes (TGFB2/3, SMAD1/2/5/7/9, and TGFBR1) were significantly increased in the high risk group (figure 6f). These results indicated that activation of TGF-beta signaling pathway might play a pivotal role in 5-FU resistance.

Immune checkpoints represent an immunosuppressive mechanism that allows cancer cells to escape anti-tumor immunity. Some immune checkpoint molecules have been identified as potential immunotherapeutic targets,[38] such as PD-1, PD-L1, CTLA-4, LAG-3, and TIM3.[21] We analyzed the association between Imm-R scores and the expression of immune checkpoint molecules. In patients with low risk, the expression of CTLA-4 was increased compared with patients with high risk (p = .043, Figure S6). However, there was no evidence of significant association between Imm-R scores and expression of PD-1, PD-L1, LAG-3, and TIM3 (Figure S6).

Correlations between the Imm-R scores and the clinical prognostic factors

We analyzed the correlations between the Imm-R scores and clinical prognostic factors, including gender, T stage, N stage, M stage, tumor location, KRAS mutation status, microsatellite state, BRAF mutation status, p53 mutation status, CpG island methylator phenotype, and chromosomal instability. Patients were stratified into different groups by each clinical factor and the Imm-R scores of different subgroups were compared. As shown in Figure S7, the Imm-R scores were significantly lower in patients with wild-type KRAS genotype and dMMR, suggesting that Imm-R score was intrinsically related to KRAS status and dMMR. However, no significant correlations were found between Imm-R scores and gender, T stage, N stage, M stage, tumor location, BRAF mutation status, p53 mutation status, CpG island methylator phenotype, or chromosomal instability.

Discussion

Fluoropyrimidine is an important component of systemic chemotherapeutic regimens for CRC patients in palliative and adjuvant settings.[5] However, patient response rates to fluoropyrimidine-based chemotherapy remain low, and a significant proportion of patients experience substantial side effects, without benefiting from chemotherapy. The prediction of response to chemotherapy is critical for effective treatment of CRC. An increasing number of studies indicate that the immune response is implicated in chemotherapy-induced cytotoxicity.[27] Here, we developed and validated a TIICs- and IRGs-based Imm-R model for prognostic prediction of CRC patients receiving fluoropyrimidine-based chemotherapy. The Imm-R model achieved favorable discrimination and calibration in both the training and the validation sets for identification of high-risk patients. Biological pathways involved in the risk model include activation of cell growth, proliferation and immune-related pathways, as well as defects in DNA damage and repair pathways. Our findings provide new insights into the role of immune responses in modulating patient response to fluoropyrimidine-based chemotherapy, with potential implications for personalized follow-up and decision-making regarding individualized chemotherapy in CRC patients.

The balance between the TIICs and the cancer cell population in the tumor microenvironment plays an important role in tumor progression.[39] The presence of TIICs within the tumor microenvironment is considered to reflect the immune interaction between the host and tumor.[40] Complex biological interactions occur between immune cells and cancer cells in the tumor microenvironment, which have substantial prognostic relevance because the immune system plays both tumor-promoting and -suppressing roles.[41] In the present study, we developed TIICs-based immune infiltration risk scores comprising Tfh, M0 macrophage, and eosinophil populations. Tfh and M0 macrophages were associated with favorable prognosis while eosinophils was an indicator of poor prognosis. Tfh is a specialized subset of CD4+ helper T cells present in the lymphoid organs and in peripheral blood.[42] Inducible costimulator (ICOS), C-X-C chemokine receptor type 5 (CXC-R5), and PD-1 are classic cell surface markers of Tfh cells,[43] which play a key role in the development of antigen-specific B cells.[44] Numerous studies suggest that Tfh cells exhibit antitumoral immunity[45] by promoting the effector functions of CD8 + T cells via the secretion of IL-21, which is downregulated due to PD-1/PD-L1-mediated suppression.[46] We found that Tfh cells were associated with favorable prognosis in patients receiving fluoropyrimidine-based chemotherapy, which indicated that Tfh cells might be a potential prognostic marker in the context of fluoropyrimidine treatment. Among various TIICs, macrophages of the M0 and M1 subtypes show the highest infiltration in TME compared with that in normal tissues.[47] Indeed, recent studies support the view that fully polarized macrophages (M1 and M2) are the extremes of a continuum of macrophage polarization.[48] M0 is the unactivated subtype, while M1 activates the inflammatory response and participates in the antitumoral response. M2 inhibits T cell proliferation and promotes tumor proliferation. A recent study found that the abundance of M0 and M1 subtypes was increased in TME of CRC tissues compared with that in normal tissues, indicating that these subtypes were implicated in CRC development and progression.[49] There was evidence to show that high densities of tumor-associated macrophages strongly predict the benefitting from 5-fluorouracil adjuvant therapy in stage III CRC patients.[50] Besides, in-vitro study indicated that exposure of M0 macrophages to 5-FU increased the expression of the M1-marker, which significantly enhanced the cytotoxicity. Eosinophils are evolutionarily conserved cells and are mainly studied in the context of allergies. Recent data indicate that eosinophils function as immune effectors and immunomodulators in the tumor microenvironment.[51] Tumor-infiltrating eosinophils is associated with improved survival of CRC patients.[52,53] However, in the present study we found that infiltration of eosinophil was associated with a poor prognosis in CRC patients receiving 5-FU chemotherapy. Whether there is an interaction between eosinophils and 5-FU needs to be further investigated. In addition to the tumor immune microenvironment, which can affect the response to chemotherapy, chemotherapeutic agents can also affect the tumor immune microenvironment.[54] For example, oxaliplatin-induced tumor cell death invoked tumor-targeting immune responses, as shed tumor antigens captured and presented by dendritic cells causing activation of cytotoxic T lymphocytes.[55,56] Interestingly, we found that no significant difference in OS between the Imm-R model-defined high- and low-risk patients when receiving FOLFOX (5-FU, leucovorin, and oxaliplatin). A possible explanation is that oxaliplatin promotes immunogenic death in high-risk patients.

Significant breakthroughs have been made in cancer immunotherapy; in particular, IRGs offer potential for the identification of new molecular targets for cancer immunotherapy.[57] In the present study, we identified seven IRGs with potential prognostic significance in patients receiving fluoropyrimidine-based chemotherapy. CCL22 plays an important role in the regulation of the tumor microenvironment.[58] A recent study indicated that CCL22 could be used as a biological factor to evaluate tumor chemotherapy response and progression in CRC.[59] These previous findings along with the current results highlight the potential application of CCL22 for enabling improved management of adjuvant chemotherapeutic regimens in patients with CRC during preoperative and postoperative evaluation. FGF19 belongs to the fibroblast growth factor (FGF) family, which exhibits broad mitogenic and cell survival activities. FGF19 expression is increased in CRC samples, and FGF19 blockade has been shown to inhibit the growth of colon tumors.[60] LGR5 is a promising marker of cancer stem-like cells in colorectal cancer.[61] Overexpression of LGR5 is associated with poor survival and recurrence in CRC.[62] Furthermore, LGR5 is a predictive marker for fluoropyrimidine-based adjuvant chemotherapeutic response in CRC, which was in line with our results.[63] TACSTD2 is up-regulated in CRC and associated with response to short-term 5-fluoropyrimidine treatment in P53-depleted HCT116 colon cancer.[64] IRGs in our established Imm-R scores can be regarded as individual biomarkers, and their immune features and prognostic significance may enable better management of patients receiving fluoropyrimidine-based chemotherapy.

Although the present findings provide new insights into the impact of immune responses on outcomes of fluoropyrimidine-based chemotherapy, there are limitations to our work. First, the data on which the prediction model was built were obtained from public databases, and it was not possible to obtain all information needed for each patient which may affect the accuracy of results. Second, there was significant heterogeneity in the patients in the study. For example, GSE39582 contained patients with stages II–IV, and GSE103479 derived from CRC patients with stage II–III, while GSE72968 and GSE72969 included patients with metastatic CRC. Although the risk model was validated in different populations, this study might overlooked factors that had a potential impact on the risk model. Third, although the expression profiling data included in this study were all derived from gene microarray, different sources of gene microarrays, different analysis procedures may affect the accuracy of this study. Fourth, the statistical tests widely used in this study and the resulting false positives required attention.

In summary, we systematically identified prognostic TIICs and IRGs and developed and validated a TIICs- and IRGs-based Imm-R model for identifying CRC patients receiving fluoropyrimidine-based chemotherapy with poor prognosis. These results have potential implications for personalized follow-up and decision-making regarding individualized chemotherapy in CRC patients.

Materials and Methods

CRC datasets and preprocessing

We systematically searched publicly available CRC gene expression data sets on Gene Expression Omnibus (GEO) database fulfilled the following criteria: (1) tumor samples derived from CRC patients receiving fluoropyrimidine-based chemotherapy, (2) gene expression was detected by microarray or high-throughput sequencing, (3) chemotherapy and prognosis was available. Exclusion criteria included: (1) samples derived from in vitro or in vivo experiments or human xenografts, (2) data set cannot be processed. In total, we collected four datasets (GSE39582,[65] GSE103479,[66] GSE72968,[67] and GSE72969[67]). Dataset GSE39582 contained patients with stage II–IV, and GSE103479 derived from CRC patients with stage II–III, while GSE72968 and GSE72969 included patients with metastatic CRC. GSE39582 served as the training set and GSE103479 served as the validation set1. GSE72968 and GSE72969 were generated using the same microarray platform by the same research group in France, so we combined them into validation set2. GSE39582, GSE72968 and GSE72969 were generated using the Affymetrix Human Genome U133 Plus 2.0 Array and GSE103479 was generated by Almac Diagnostics Custom Xcel array. The raw data of all microarray datasets were downloaded from GEO and subjected to background adjustment using the RMA algorithm.[68] The microarray probe set was then mapped to the gene symbol according to the platform annotation file. Corresponding clinical characteristics of CRC samples were also obtained. Overall survival is defined as the length of time from the date of diagnosis to death. Recurrence-free survival is defined as the length of time after primary treatment ends that the patient survives without any signs or symptoms of CRC.

In order to investigate the potential mechanisms of 5-FU resistance, we downloaded a microarray dataset (GSE81005) derived from wild type and 5-FU-induced resistant HCT-8 human colorectal cancer cells. GSE81005 was detected by Affymetrix Human Gene Expression Array. The raw data of the microarray datasets were downloaded from GEO. Datasets were normalized and annotated as described above.

Inference of infiltrating cells in the tumor microenvironment

To evaluate the abundance of TIICs in the CRC samples, we employed the LM22 gene signature and the CIBERSORT algorithm, which could sensitively and specifically discriminate 22 human immune cell phenotypes (B cells, T cells, NK cells, macrophages, DCs, and myeloid subsets).[69] CIBERSORT is a deconvolution algorithm based on support vector regression, which uses a set of reference gene-expression values corresponding to a minimal representation for each cell type to infer cell type proportions in data from bulk tumor samples with mixed cell types. The CIBERSORT algorithm can, in particular, be used to derive the proportion of cells in complex microarray data.[17] Using Monte Carlo sampling, CIBERSORT calculates the empirical p value of the deconvolution to indicate the accuracy of the results, while a p value of <0.05 indicates that the inferred cell composition is highly reliable.[29] Therefore, we only retained CRC samples with CIBERSORT p values <0.05 for subsequent analysis.

Construction and validation of tumor-infiltrating immune cell signature

In the training set, we first applied the univariable Cox proportional hazards regression model to define the prognostic value of the TIICs for overall survival (OS). The optimal cutoff values of the proportions of different TIICs in the training cohort were calculated on the basis of the prognostic significance using X-Tile software.[70] TIICs with significant prognostic value were first validated in the validation set. TIICs with significant prognostic value in both training set and validation set were analyzed using the multivariate Cox proportional hazards regression model, and those with a p value < 0.05 were used to construct the TIICs signature. The individual TIICs scores were calculated based on the immune cell signature. A formula for the TIICs scores was established to predict patient survival: TIICs scores = (−0.990 × abundance of Tfh)+(−0.731 × abundance of Macrophages M0) + (0.655 × abundance of Eosinophils). The prognostic performance of the TIICs signature was evaluated using the receiver operating characteristic (ROC) curve and area under the curve (AUC).

Construction and validation of the immune risk (Imm-R) model

We integrated the IRGs and TIICs signature to develop the Imm-R model. A functionally annotated list of IRGs (n = 1811), using Gene Ontology, was downloaded from ImmPort, and a list of genes related to antitumor immunity as reported in literatures was obtained from TISIDB (n = 988). Merging the two lists and removing duplicate genes resulted in 2404 unique IRGs. We first performed univariate Cox regression analysis for preliminary screening of survival-associated IRGs in the training set. Genes entered the model as continuous variables. Based on a significance threshold of p value < 0.05, the significative IRGs determined by univariate Cox regression analysis and TIICs risk scores were subjected to LASSO-penalized Cox regression analysis. In this regression analysis, the penalty parameter lambda.1se for prevention of overfitting was selected using 10-fold cross validation.[71] The features selected by LASSO-penalized Cox regression were subsequently entered into the multivariate Cox proportional hazards regression model. Finally, features with p values < 0.05 in the multivariate Cox proportional hazards regression model were used for construction of the Imm-R model. A formula for the Imm-R score was established to predict patient survival:

Imm-R score = (−0.939× PSMD3)+(−1.304× CCL22)+(0.319× FGF19)+(−0.228× LGR5)+(0.665× SKAP1)+(0.182× TACSTD2)+(0.329× S100A4)+(0.809× TIICs scores)

The individual Imm-R score in both training set and validation set was calculated according to the formula accordingly.

A nomogram was used to visualize the Imm-R model and predict the three-year and five-year survival rate. The ROC curve and AUC were used to evaluate the discrimination, and the calibration curve was used to evaluate the goodness of fit of the Imm-R model.

Biological significance of the Imm-R model

To explore the biological implication of the Imm-R model, we performed a gene set enrichment analysis (GSEA) using the clusterProfiler R package.[72,73] The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene ontology (GO) biological process were evaluated. Differential expression analysis was performed using the “limma” package with False Discovery Rates (FDR) <0.05 and |logFC| >1 as the threshold.

Statistical analysis

In this study, all statistical analyses were performed using R software (version 3.5.1; https://www.Rproject.org/). Continuous variables were compared using Student’s t-tests or Wilcoxon rank-sum tests. Survival analyses were conducted using the Kaplan–Meier method with log-rank test by ‘survival’ package. Feature selection and model construction were conducted with Univariate and Multivariate Cox regression. Time-dependent ROC curve and LASSO Cox regression analysis with 10-fold cross-validation and nomogram were estimated using R packages[74] (survivalROC, glmnet, rms). Unless otherwise stated, statistical significance was defined at p values < 0.05. The planned sample size was calculated with α of 0.05, β of 0.2 and HR of the risk score in the training set.

Click here for additional data file.