A 12-immune cell signature to predict relapse and guide chemotherapy for stage II colorectal cancer

The management of stage II colorectal cancer is still difficult. We aimed to construct a new immune cell-associated signature for prognostic evaluation and guiding chemotherapy in stage II colorectal cancer. We used the “Cell Type Identification by Estimating Relative Subsets of RNA Transcripts” (CIBERSORT) method to estimate the fraction of 22 immune cells by analyzing bulk tumor transcriptomes and a LASSO Cox regression model to select the prognostic immune cells. A 12-immune cell prognostic classifier, ISCRC, was built, which could successfully discriminate the high-risk patients in the training cohort (GSE39582: HR = 3.16, 95% CI: 1.85–5.40, P < 0.0001) and another independent cohorts (GSE14333: HR = 3.47, 95% CI: 1.18–10.15, P =0.0167). The receiver operating characteristic analysis revealed that the AUC of the ISCRC model was significantly greater than that of oncotypeDX model (0.7111 versus 0.5647, p=0.0152). We introduced the propensity score matching analysis to eliminate the selection bias; survival analysis showed relatively poor prognosis after chemotherapy in stage II CRC patients. Furthermore, a nomogram was built for clinicians and did well in the calibration plots. In conclusion, this immune cell-based signature could improve prognostic prediction and may help guide chemotherapy in stage II colorectal cancer patients.


INTRODUCTION
Colorectal cancer (CRC) is still a major health burden worldwide, it ranks as the third leading cause of cancer death, with an estimated 881,000 deaths in 2018 [1]. According to the SEER cancer statistics in 2017, the stage I/II CRC patients (Localized stage) account for about 40%, and about 35% for stage III patients (Regional stage) [2]. Surgical intervention is the basis of curative treatment for CRC, however, about 13.6% of stage II CRC patients and 21.5% of stage III CRC patients develop relapse after surgery [3]. Chemotherapy may be given after surgical resection to eradicate the remaining cancer cells. For stage III CRC, post-surgical chemotherapy is now the standard treatment [4], but the benefit of post-surgical chemotherapy remains controversial in stage II CRC. QUASAR trial reveals that chemotherapy with fluorouracil and folinic acid could improve outcomes in stage II CRC, but the absolute improvements are small, indicating that the decision to provide post-surgical chemotherapy in stage II CRC needs to be more cautious [5]. With the early screening of tumors worldwide, cancer patients are becoming younger, and more CRC patients are diagnosed in the AGING early tumor stage [2], exacerbating the dilemma for treatment of stage II CRC patients. Thus, there is an urgent need to construct new prognostic markers in stage II CRC to discriminate the patients who may benefit from post-surgical chemotherapy.
The tumor is highly heterogeneous, as is stage II CRC. The disparities in CRC survival and the benefit of postsurgical chemotherapy may be related to the complex mechanism of tumorigenesis and development. Chromosomal instability and somatic mutation are critical genetic factors implicated in tumorigenesis [6][7][8]; and mis-match repair (MMR) has been extensively investigated in CRC. Studies show that CRC patients with microsatellite instability (MSI) tumors are associated with favorable outcomes, compared with the patients with microsatellite stable (MSS) tumors [9,10]. However, MMR status cannot predict the benefits of chemotherapy in stage II CRC, nor can KRAS or BRAF mutations [11]. Recently, gene expression profiles have shown great promise in prognosis assessment. Several gene signatures have been built to evaluate patients' prognosis, as well as predict the benefit of chemotherapy [12][13][14][15][16][17]. However, few have been applied clinically and the robustness and reliability still need further evaluation. The most widely used gene signature is the OncotypeDX colon cancer assay [18,19], which has been commercialized since 2010 and is mainly used to predict the recurrence risk in stage II CRC [20,21]. Thus, there is of great clinical significance to identify novel molecular markers from a new perspective.
The role of tumor-infiltrating immune cells is currently getting increasing attention in cancer research. As one of the pivotal components of the tumor microenvironment, tumor-infiltrating immune cells exert their pathophysiological functions through reciprocal communication with neoplastic cells [22,23]. Tumor-infiltrating immune cells are consist of both mononuclear and polymorphonuclear immune cells, such as macrophages, T cells, B cells, natural killer cells, etc. Studies reveal that the abundance of tumor-infiltrating immune cells is closely related to tumor stage and has significant tumor specificity [24][25][26]. The prognostic value of tumorinfiltrating immune cells has also been investigated in various cancers, and the intra-tumoral γδ T-cell signatures have emerged as the most significant favorable cancerwide prognostic populations [23,25,27,28]. Besides, immune cells are also implicated in chemotherapeutic response in cancer, and macrophages are found to reduce chemotherapy sensitivity [29,30].
Comprehensive quantification of the immune cell infiltrates in tumors can be accessed by multiple computational analyses of the gene expression profiles, which are relatively easy to obtain from the accumulating public microarray data and RNA sequencing data of human tumors [31]. The algorithm of "Cell Type Identification by Estimating Relative Subsets of RNA Transcripts" (CIBERSORT) is the current most accurate and extensively used computational algorithm available for enumeration of various immune cell types [32,33]. In our study, CIBERSORT was used to assess the distribution of 22 immune cells from tumor RNA transcripts, and then the least absolute shrinkage and selection operator method (LASSO) method was introduced to construct a 12-immune cell signature to predict disease-free survival (DFS) and assist in evaluating the benefit of chemotherapy for stage II CRC patients.

Development and validation of the ISCRC model
The LASSO Cox regression was introduced to interrogate the relevance between 22 immune cells and the patients' survival in the training cohort-GSE39582, (Supplementary file 1 and Supplementary Figure 1), and a 12-immune cell model was identified, which was significantly associated with DFS of stage II CRC. The 12 immune cells and associated coefficients generated through LASSO analysis were shown in Supplementary file 2 and Supplementary Table 1. The immune scores for CRC patients, namely ISCRC, were calculated based on the fractions of 12 immune cells in each sample and the associated coefficients, and the formula of calculation was shown in Supplementary files 3 and Supplementary Table 2. A dichotomous ISCRC was adopted in survival analysis, based on the optimum cut-off value of ISCRC (0.7473), patients were divided into low-and high-ISCRC groups. The Kaplan-Meier survival analysis revealed that patients in the high-ISCRC group had worse outcomes compared with the low-ISCRC group [hazard ratio (HR) =3.16, 95% confidence interval (CI) =1.85-5.40, P < 0.0001; Figure 1A]. The efficacy of the ISCRC model for the prognosis prediction in stage II CRC patients was further validated in another independent dataset GSE14333. Patients were also classified into two subgroups using the same cut-off value (0.7473), and it generated consistent results (HR=3.47, 95% CI=1.18-10.15, P = 0.0167, Figure 1B). We also validated the ISCRC model in the combined cohort of GSE39852 and GSE14333, and a significantly different prognosis can be seen between low-and high-ISCRC groups (HR=3.35, 95% CI=2.09-3.93, P < 0.0001, Figure 1C).
In the univariate Cox regression model, the ISCRC classifier was found to be a strong variable correlated with CRC recurrence in both training and validation cohorts (GSE39582: HR=3.1630, P <0.0001; GSE14333: HR=3.4460, P = 0.0234, Figure 2A). After adjustment for the common clinical covariates, including age, sex, and chemotherapy, multivariate Cox regression analysis demonstrated that the ISCRC classifier remained an independent prognostic factor for DFS in the training dataset (GSE39582: HR=3.51, 95%CI=2.03-6.06, P < 0.0001; Figure 2B) and the validation dataset (GSE14333: HR=3.05, 95%CI=1.02-9.16, P = 0.0468; Figure 2C). To investigate the sensitivity and specificity of survival prediction, the receiver operating characteristic (ROC) analysis was also performed to calculate the area under ROC curve (AUC) of the ISCRC model. Supplementary file 4 and Supplementary  Figure 2 showed that our ISCRC model owned considerable predicted power of prognostic evaluation for stage II CRC patients in the training cohort (GSE39582: AUC=0.7111) and the validation dataset (GSE14333: AUC= 0.7041).

Comparison ISCRC with OncotypeDX colon
To further evaluate the prognostic value of the ISCRC model, comparison analysis was performed between our AGING ISCRC model and other known gene makers. We did not mean to make a comprehensive review of all prognostic biomarkers in CRC, thus the widely used gene signature, oncotypeDX colon cancer assay, was selected as the representative. The prognostic indexes of the oncotypeDX model were calculated according to the associated formula (See Supplementary file 3  and Supplementary Table 2). We first performed the univariable Cox regression analysis in GSE39582 to assess the prognostic value, where the prognostic index was used as a continuous variable. As shown in Figure  3A, the ISCRC and oncotypeDX models were all AGING significantly associated with DFS in stage II CRC, but the HR of our ISCRC model was significantly larger, with an even lower p-value. ROC analysis was also performed in GSE39582 to assess the sensitivity and specificity of survival prediction. As shown in Figure  3B, the AUC of the ISCRC model was significantly greater than that of the oncotypeDX model (0.7111 versus 0.5647, p=0.0152), indicating well prognostic accuracy of our ISCRC model.

ISCRC and the benefit of chemotherapy
Studies have revealed the close association between chemotherapy and tumor microenvironment. One of the crucial components of the tumor microenvironment is the immune cell, and macrophages have been reported to promote tumor angiogenesis and drive chemotherapy resistance [29,30,34]. In our study, the ISCRC model demonstrated well prognostic value for stage II CRC patients, so we speculated that our immune cell-derived prognostic biomarker might be associated with chemosensitivity. All the two datasets in our study provided information on chemotherapy, so we intended to examine the benefit of chemotherapy in stage II patients. As a random assignment of chemotherapy to samples is not feasible in the retrospective study, the presence of an imbalance in baseline characteristics between the control and treatment groups can lead to a biased estimation [35]. Thus, the propensity score matching analysis was performed to eliminate the selection bias. After propensity score matching analysis, 82 matched samples were generated in GSE39582, but only 16 matched samples left in GSE14333, so the cohort GSE14333 was not suitable for further analysis (see supplementary file 5 and Supplementary Table 3). The matched patients in GSE39582 were classified into two subgroups based on the status of chemotherapy, as shown in Figure 4A, Kaplan-Meier survival analysis revealed that the prognosis of the patients in chemotherapy group was significantly worse than in non-chemotherapy group (HR=5.34,95%CI=1.70-16.81,p=0.0013). We further investigated the effect of chemotherapy in low-and high-ISCRC groups, respectively. Figure 4B demonstrated the similar result in the low-ISCRC group (HR=7.64, 95%CI=1.62-36.11, P=0.0025), however, no significant difference was seen in the high-ISCRC group (HR=4.23, 95%CI=0.68-26.22, p=0.0941) ( Figure 4C), these results suggested that stage II patients in the low-ISCRC group were not suitable for chemotherapy.

Distribution of ISCRC and clinical characteristics
Our ISCRC model was composed of 12 immune cells, and then the distribution of 12 immune cells and the common clinical parameters in the low-and high-ISCRC groups was illustrated in Figure 5. Our study revealed that the low-ISCRC group was characterized by high expression of M1 Macrophages and B memory cells but low expression of M2 Macrophages, while the distribution of these three immune cells in the high-ISCRC group was opposite in training cohort GSE39582 ( Figure 5A). A similar distribution characteristic can also be observed in GSE14333 ( Figure 5B). The DFS status was significantly different between the low-and high-ISCRC groups in both GSE38582 and GSE14333, and there were more recurred patients in the high-ISCRC group. However, the distribution of the common clinical factors between different ISCRC groups did not vary significantly.

Construction of nomogram based on ISCRC
To develop a quantitative recurrence prediction method for clinical application, a nomogram was built, which integrated the ISCRC model, age, sex, and AGING chemotherapy status ( Figure 6A). The nomogram revealed that age had the largest impact on the patients' prognosis, followed by the ISCRC model and adjuvant chemotherapy. The total points for each patient were calculated based on these variables to estimate the possibility of recurrence in the following 3-, 5-and 7-year for each stage II CRC patient. The performance of the nomogram was evaluated by calibration plots ( Figure 6B). The line-segment was close to 45-degree, suggesting that the nomogram did quite well. DCA showed a well clinical utility of the nomogram ( Figure 6C).

DISCUSSION
Recently, the study of the immune contexture has receiving accumulating attention in cancer research [24,36]. Several immune score models have been developed from analyzing the fractions of immune cells to predict prognosis in various types of tumors, including CRC [28,37,38]. However, the immune score model, especially for stage II CRC, has not been reported yet. In our study, we have constructed a novel immune cellderived prognostic marker, ISCRC model, to predict the relapse of stage II CRC. The prognostic value of the AGING ISCRC model was confirmed in an independent validation cohort and the combined cohort, indicating good reproducibility. When compared with the representative known gene signature-oncotypeDX colon, the ISCRC model exhibited a better predictive capability for DFS, so we concluded that our ISCRC model could improve the recurrence prediction in stage II CRC. Furthermore, our study demonstrated the adverse effect of chemotherapy in stage II CRC, especially for patients with low-ISCRC. This finding might shed new light on the eligibility of adjuvant chemotherapy for stage II CRC patients.
Immune cells are implicated in the prognosis of cancer patients, and the underlying mechanism is partly due to the immune-suppression effect induced by tumorassociated dendritic cells, neutrophil, and M2 macrophages, etc. [25,[39][40][41]. However, it is very difficult to measure the immune infiltrates in tumors by traditional methods, such as IHC or flow cytometry; they are not qualified for the task of discriminating varieties of immune cells and quantifying them simultaneously. It is well known that solid tumors are highly heterogeneous and lack distinctive markers for various cells. Therefore, a new method is urgently needed to comprehensively study the immune contexture in cancers. In our study, a newly developed computer-based analytical algorithm, CIBERSORT, was introduced to evaluate the fractions of immune cells, which allows estimating the abundance of immune cells simultaneously in a large patient cohort with high accuracy. Then the LASSO analysis was performed to construct the immune cell-derived prognostic model-ISCRC; it is a reliable instrument in selecting the prognostic markers and has been extensively used in high-dimensional microarray data [42][43][44]. We validated the ISCRC model in another independent cohort-GSE14333. Furthermore, comparison analysis showed that the ISCRC model had better predictive power compared with the representative gene signature-oncotypeDX, which has been used clinically and validated in a large prospective study [20]. These findings indicated that our ISCRC model could improve the prognostic prediction in stage II CRC with fair reliability.
Interestingly, we found that the distribution of common clinical factors was not significantly different between low-and high-ISCRC groups in stage II CRC patients, so the ISCRC model might serve as a good supplement to the clinical-pathological features for stage II CRC. Further analysis revealed that the low-ISCRC group was characterized by high expression of M1 Macrophages and low expression of M2 Macrophages, while the high-ISCRC group presented the opposite distribution feature.
Macrophages play a pivotal role in innate immunity and are closely associated with inflammation and host defense [45]. The patterns of macrophage activation can roughly divide into two main types: classically activated macrophages (M1) and alternatively activated macrophages (M2). The M1-M2 macrophages are characterized by distinctive secretome profile and biological functions [46]. M1 Macrophages can promote Th1 response and show strong antineoplastic activity [47]. M2 Macrophages are implicated in immunosuppression, chronic inflammation, etc. [41]; M2 macrophage polarization can promote tumor progression [48,49]. In our study, the high-ISCRC group showed poor prognosis compared with the low-ISCRC group. The potential mechanism might partly attribute to the immune characteristics of our ISCRC model; the high-ISCRC group was found to be characterized by high expression of M2 macrophages and low expression of M1 macrophages, thus resulting in a poor outcome.
Adjuvant chemotherapy is of great significance in cancer treatment. However, the rationality of chemotherapy in stage II CRC remains controversial. Up to date, there is still some conflicting evidence in clinical trials [50,51]. The different responses to chemotherapy in stage II patients might attribute to the inherent heterogeneity of the tumor. The significant difference of immune contexture can be seen between the low-and high-ISCRC groups, although with the same TNM stage. Our study demonstrated that adjuvant chemotherapy might exert adverse effects on stage II patients, especially for patients with low-ISCRC, so chemotherapy in patients with stage II CRC needs to be more cautious. Notably, when we interrogated the prognosis between the low-and high-ISCRC groups, the propensity score matching method was introduced to remove the selection bias, which has been widely used for retrospective analysis [52,53]. In our study, the potential confounding factors, such as patients' age, etc., were balanced in the two subgroups to accurately weigh the risk-benefit ratio of chemotherapy in stage II CRC. These findings indicated that the ISCRC model might be a useful biomarker to predict the patients' prognosis and guide chemotherapy in stage II CRC.
There are some limitations to our study. First, our study is with limited sample size, and the ratio of patients with chemotherapy to non-chemotherapy patients is relatively low in the two datasets. After propensity score matching, the sample size of GSE14333 becomes even smaller, which is not suitable for further analysis. Second, it is a retrospective analysis using publicly available datasets. The detailed clinical information for each patient is unavailable; the selection bias is hard to avoid, although propensity score matching has been AGING performed. Thus, the reliability of our analysis requires to be further verified in a prospective study.

CONCLUSIONS
In our study, a new immune cell-related marker was built through investigating the whole immune landscape to predict prognosis and evaluate the benefit of post-surgical chemotherapy in stage II CRC. The prognostic value of the ISCRC model was confirmed in another independent cohort and the combined cohort. Moreover, our study demonstrated the adverse effect of chemotherapy in stage II CRC patients, especially for patients with low-ISCRC. Thus, this analysis might have a beneficial effect on the personalized management of stage II CRC patients.

CRC datasets preparation
CRC RNA transcripts data and corresponding clinical information in our study were downloaded from ArrayExpress (http://www.ebi.ac.uk/arrayexpress/) and Gene Expression Omnibus (GEO) (https://www.ncbi.nlm.nih.gov/geo/). The raw files of microarray data were normalized using a robust multiarray averaging method [54], and the process of which was performed with "affy" and "affycoretools" packages of R software. A total of 339 CRC stage II patients were selected from two datasets for this study, that were GSE39582 (n=253) and GSE14333 (n=86). Dukes' B CRC patients in the GSE14333 cohort were used for our analyses, which are corresponding to the American Joint Committee on Cancer (AJCC) stage II patients based on AJCC Colon and Rectum Cancer staging 7th Edition.

Estimation of the immune cell type fractions through CIBERSORT
For estimation of the abundance of immune cells in CRC cancer samples, the gene expression data were processed through the CIBERSORT web portal (http://cibersort.stanford.edu/). The immune cells are distinguished with a leukocyte gene signature matrix, termed LM22. The 22 human hematopoietic cells are grouped into seven T cell types, naïve and memory B cells, plasma cells, NK cells, and myeloid subsets, and the cell subsets can be further divided into 11 major leukocyte types [32]. The CIBERSORT method can discriminate the human immune cells with high sensitivity and specificity and has been well validated on microarray data. The estimated fractions of immune cell subsets by this method are considered to be accurate with the threshold of P <0.05 [33,55]. The estimated CIBERSORT results were normalized and the sum of proportions of 22 immune cells in each sample was equal to 1.

Construction and validation of the ISCRC model
GSE39582 dataset was used as the training cohort to perform the LASSO Cox regression analysis with the "glmnet" package of R software (version 3.5.1) (R Foundation for Statistical Computing, Vienna, and Austria). The penalized Cox regression model with LASSO penalty can achieve shrinkage and variable selection simultaneously with 10-times cross-validations [44,56]. Based on the optimal penalty parameter lambda, the most useful prognostic markers with associated coefficients were screened out from the 22 candidate immune cells. The immune score for each sample in CRC was calculated based on the fraction of each immune cell and its associated coefficient, then they were normalized and the "survminer" package was used to determine the best cut-off value. Thus, a multiimmune cell classifier, ISCRC, was constructed to predict the recurrence probability for CRC patients. All samples were split into low-ISCRC and high-ISCRC groups, and the prognostic value of ISCRC was further validated in GSE14333. The workflow of this study was depicted in Supplementary file 6 and Supplementary Figure 3.

Statistical analysis
DFS was from the time of surgical intervention to the first confirmed relapse. The DFS differences were estimated by the Kaplan-Meier estimator and the logrank test. To investigate the effect of chemotherapy on prognosis, the propensity matching analysis was conducted with the "MatchIt" package to remove the selection bias and the "exact" method was used to match samples between the chemotherapy group and the non-chemotherapy group [57].
The correlation between our ISCRC model and clinicopathological features was assessed using the chisquare test and one-way ANOVA test. Univariate and multivariate Cox proportional hazards regression models were performed to estimate the prognostic value of ISCRC and clinical characteristics.
To assess the sensitivity and specificity of the survival prediction, ROC analysis was performed to calculate the AUC with the "pROC" package, which served as a measure to assess the accuracy of diagnosis [58]; the method "delong" was used to investigate the differences between ROC curves. For ROC analysis, the samples were excluded who still did not relapse at last follow-up and whose follow-up times were relatively short (less than the median DFS) [59].

AGING
The nomogram was constructed using Cox regression coefficients of ISCRC and other clinical factors and the nomogram plots were generated with the "rms" R package. The performance characteristics of the nomogram were assessed by calibration plots, where the 45-degree line stood for the ideal prediction model. Decision curve analysis (DCA) was introduced to explore the clinical utility of the nomogram [60,61]. All the statistical analyses in our study were performed with R software (3.5.1) and a p-value less than 0.05 was regarded to be statistically significant.

AUTHOR CONTRIBUTIONS
XT and XZ and WM contribute equally to the work. XT and XZ drafted the manuscript; WM, XT and XZ performed data analysis; SB, MS, SX, and CZ prepared all the figures and tables. XT and YW designed this project. All the authors reviewed and approved the final manuscript.

CONFLICTS OF INTEREST
The authors declare that they have no conflicts of interest.

Supplementary Tables
Supplementary Table 1. Immune cells significantly associated with the disease-free survival in the training cohort (N=253).

CALCULATION OF PROGNOSTIC INDEX FOR 12-IMMUNE CELL SIGNATURE (ISCRC)
The calculation of the immune risk score for individual patients was based on the multivariate model including all twelve immune cells as shown in table S1. Firstly, through LASSO Cox regression analysis, we have identified 12 prognostic immune cells and their associated coefficients. Then, multiply the fraction of each immune cell by its associated coefficient to generate the value for each immune cell. At last, sum all the values of 12 immune cells to get the immune risk score for each patient. The formula was as follows: Immune risk score = (

ALGORITHMS TO CALCULATE RISK SCORES FOR ONCOTYEDX
Recurrence Risk score is calculated using the prespecified genes and algorithm [1,2]: The O'Connell Recurrence Risk (RS) score is composed of 12 genes among which 5 reference genes and 7 genes associated to recurrence.
For the reference genes, when several probe set were possible, the less variant one was selected.
For the other genes, data were median gene centered and aggregated by mean if several probe sets were available.