Abstract
Objective. Stromal cells and immune cells have important clinical significance in the microenvironment of colorectal cancer (CRC). This study is aimed at developing a CRC gene signature on the basis of stromal and immune scores. Methods. A cohort of CRC patients () were adopted from The Cancer Genome Atlas (TCGA) database. Stromal/immune scores were calculated by the ESTIMATE algorithm. Correlation between prognosis/clinical characteristics and stromal/immune scores was assessed. Differentially expressed stromal and immune genes were identified. Their potential functions were annotated by functional enrichment analysis. Cox regression analysis was used to develop an eight-gene risk score model. Its predictive efficacies for 3 years, 5 years, overall survival (OS), and progression-free survival interval (PFI) were evaluated using time-dependent receiver operating characteristic (ROC) curves. The correlation between the risk score and the infiltering levels of six immune cells was analyzed using TIMER. The risk score was validated using an independent dataset. Results. Immune score was in a significant association with prognosis and clinical characteristics of CRC. 736 upregulated and two downregulated stromal and immune genes were identified, which were mainly enriched into immune-related biological processes and pathways. An-eight gene prognostic risk score model was conducted, consisting of CCL22, CD36, CPA3, CPT1C, KCNE4, NFATC1, RASGRP2, and SLC2A3. High risk score indicated a poor prognosis of patients. The area under the ROC curves (AUC) s of the model for 3 years, 5 years, OS, and PFI were 0.71, 0.70, 0.73, and 0.66, respectively. Thus, the model possessed well performance for prediction of patients’ prognosis, which was confirmed by an external dataset. Moreover, the risk score was significantly correlated with immune cell infiltration. Conclusion. Our study conducted an immune-related prognostic risk score model, which could provide novel targets for immunotherapy of CRC.
1. Introduction
CRC, as a heterogeneous disease, is a common cause of cancer-related deaths worldwide [1]. TNM staging is usually considered to be one of the main tools for CRC prognosis [2]. However, the prognosis varies greatly among CRC patients with the same TNM stage, suggesting that the current TNM stage does not well provide complete prognostic information for CRC patients. Therefore, it is necessary to adopt a new strategy to increase the predictive efficiency of prognosis and survival outcomes of CRC patients.
Due to the considerable heterogeneity between CRCs, determination of the optimal treatment strategy at the individual level faces the large challenges. Thus, it is an urgent need to conduct robust models to identify high-risk CRC patients and to find novel molecular targets. In the tumor microenvironment (TME), stromal and immune cells are involved in the development of CRC [3, 4]. Increasing evidence suggests that stromal and immune cells possess critical clinical significance for CRC. It has been reported that stromal cells can contribute to transcriptome and clinical features of CRC subtype [5]. Furthermore, stromal gene expression can more robustly predict the prognosis of CRC subtypes compared to epithelial tumor cells [6]. In a large cohort of CRC patients, infiltrating immune cell data could better predict patients’ survival than histopathological methods [7]. Growing studies have found that infiltrating immune cells are involved in chemoresistance [8] and metastasis [9]. Thus, it is essential to further analyze the biological characteristics of stromal and immune genes and to determine their prognostic value for CRC patients. However, there is a lack of stromal and immune scores that can predict CRC patients’ prognosis based on multiple clinical factors. Moreover, robust prognostic models on the basis of stromal and immune scores are also lacking.
In this study, we established a reliable prognostic immune-related risk score for CRC. Our results could offer novel insights for prediction of CRC patients’ prognosis and development of individualized immunity therapy strategies.
2. Materials and Methods
2.1. CRC Datasets
TCGA RNA-seq data (including Counts and FPKM data) of GDC CRC (including 469 CRC tissue samples and 41 adjacent normal tissue samples) were downloaded from the xenabrowser website (https://xenabrowser.net/). Among all CRC samples, 433 samples contained complete clinical information, including gender, age, TNM stage, tumor grade, microsatellite instability (MSI), and mismatch repair (MMR). The clinical information of 433 CRC patients is listed in Table 1. Survival information including OS status, OS time, progression-free survival (PFS) status, and PFS time was derived from the pan-cancer on the GDC website (https://gdc.cancer.gov/about-data/publications/PanCan-Clinical-2018). Furthermore, mutation data (including BRAF, KRAS, and TP53) were from CRC MuTect. An overview of the workflow is shown in Figure 1.
2.2. Estimation of Stromal/Immune Scores
ESTIMATE algorithm was used to calculate the stromal/immune scores on the basis of unique expression profiles of stromal/immune cells by the ESTIMATE package in R (https://bioinformatics.mdanderson.org/estimate/) [10].
2.3. Kaplan-Meier Survival Analysis
According to the optimal cutoff of stromal/immune scores, CRC samples were classified into high and low stromal/immune score groups. Kaplan-Meier plot of overall survival between the two groups was analyzed, and the results were evaluated by log-rank test.
2.4. Correlation between Clinical Characteristics and Stromal/Immune Scores
To probe into the clinical significance of stromal/immune scores, we analyzed the correlation between clinical characteristics (including pathologic T, pathologic N, pathologic M, and tumor stage) and stromal/immune scores.
2.5. Differential Expression Analysis
Differential expression analysis between high and low stromal/immune score groups was carried out using the edgeR package in R, following the screening criteria of and . Then, up- or downregulated stromal/immune genes were intersected by the VennDiagram package in R, respectively.
2.6. Functional Enrichment Analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of differentially expressed stromal and immune genes were carried out through the clusterProfiler package in R [11]. GO analysis contains three terms, cellular component (CC), molecular function (MF), and biological process (BP). value after adjustment < 0.05 was significantly enriched.
2.7. Protein-Protein Interaction (PPI) Analysis
PPI analyses of differentially expressed stromal and immune genes were carried out via The Search Tool for the Retrieval of Interacting Genes (STRING, https://string-db.org/; version 11) [12]. Then, the PPI network was visualized through Cytoscape (version 3.7.2) [13].
2.8. Univariate and Multivariate Cox Regression Analyses
Univariate cox regression analysis of differentially expressed stromal and immune genes was carried out via the survival package in R. Genes with value < 0.05 were screened for multivariate cox regression analysis. To validate the sensitivity and accuracy of the risk score for prediction of CRC, an ROC curve was drawn to evaluate the predictive performance of the risk core for 3 years, 5 years, OS, and PFI using the “tdROC” package in R. The results were visualized with the “ggplot2” package in R. The AUC was then calculated. The GSE39582 dataset from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov) was used to validate the prognostic value of the risk score. The dataset was composed of 566 CRC samples.
2.9. Immune Infiltration Analysis
The tumor-immune infiltration cells including B cells, CD4+T cells, CD8+T cells, macrophages, neutrophils, and dendritic cells were estimated via the TIMER (https://cistrome.shinyapps.io/timer/) [14]. Spearman’s correlation between the risk score and the infiltrating levels of immune cells was evaluated through the psych package in R. Furthermore, we also assessed the correlation between the genes in the risk score and marker genes of immune cells. The strength of correlation followed the criteria: suggested a high correlation, suggested a moderate correlation, and suggested a weak correlation [15].
3. Results
3.1. Immune Score Is in Significant Association with Prognosis and Clinical Features of CRC Patients
According to the optimal cutoff of stromal/immune scores, the CRC patients were divided into two groups. Kaplan-Meier OS analysis results showed that patients with high stromal score had shorter OS time than those with low stromal score; however, it was not statistically significant (Figure 2(a); value = 0.057). As depicted in Figure 2(b), we found that patients with low immune score implied a poor prognosis ( value = 0.0094). Furthermore, we analyzed the correlation between stromal/immune scores and clinical features. As depicted in the results, stromal score was not significantly associated with pathologic T (Figure 3(a); value = 0.61), pathologic N (Figure 3(b); value = 0.28), pathologic M (Figure 3(c); value = 0.63), tumor stage (Figure 3(d); value = 0.68), and age (Figure 3(e); value = 0.76). Similarly, we also found that there was no statistical significance between immune score and pathologic T (Figure 3(f); value = 0.88) and pathologic N (Figure 3(g); value = 0.17). As expected, immune score was in significant association with pathologic M (Figure 3(h); value = 0.0045) and tumor stage (Figure 3(i); value = 0.0093). However, no significant correlation between immune score and age was found in Figure 3(j) ( value = 0.29). Furthermore, ESTIMATE scores were not correlated with pathologic T (Figure 3(k); value = 0.98), pathologic N (Figure 3(l); value = 0.73), pathologic M (Figure 3(m); value = 0.095), tumor stage (Figure 3(n); value = 0.28), and age (Figure 3(o); value = 0.74). These findings indicated that immune score was in significant association with CRC patients’ prognosis and clinical features.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
3.2. Identification of Differentially Expressed Stromal and Immune Genes for CRC
We analyzed differentially expressed genes (DEGs) with and between the high and low stromal/immune score groups. As volcano plots, there were 1197 up- and 28 downregulated stromal genes in the high stromal score group (Figure 4(a)). Furthermore, 899 immune genes were upregulated and eight immune genes were downregulated in the high immune score group (Figure 4(b)). Hierarchical clustering analysis results showed that both stromal DEGs and immune DEGs could distinguish high stromal/immune score from low stromal/immune score (Figures 4(c) and 4(d)). 736 genes were upregulated both in high stromal and immune scores (Figure 4(e)). Moreover, among eight downregulated immune genes, two genes were also downregulated in the high stromal score group (Figure 4(f)). We further performed functional enrichment analysis of these common stromal and immune genes. These genes were significantly correlated with inflammatory biological processes like regulation of inflammatory response and pathways such as cytokine-cytokine receptor interaction and chemokine signaling pathway (Figures 5(a)–5(d)). As shown in the PPI network, COL6A2, COL6A1, COL5A2, C1S, and C1R were the top five genes, which were considered hub genes (Figure 5(e)).
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
3.3. Construction of an Eight-Gene Prognostic Signature for CRC
Among 738 differentially expressed stromal and immune genes, 23 genes were significantly associated with CRC patients’ prognosis according to univariate Cox regression analysis results. Of them, 20 genes were risk factors, and the remaining three (CCL22, CPA3, and MMP1) were protective factors (Table 2). These genes were used for multivariate Cox regression analysis. Finally, an eight-gene signature was constructed for CRC. The risk score was calculated on the basis of the coefficients and expression values of the eight genes. All CRC patients were divided into two groups in accordance with the median value of risk score (Figure 6(a)). Heat maps depicted the difference in expression patterns of the eight genes (including CD36, KCNE4, CPT1C, SLC2A3, RASGRP2, NFATC1, CCL22, and CPA3) between the high and low risk scores (Figure 6(b)). As shown in Figure 6(c), the risk score was capable of predicting CRC patients’ prognosis. High risk score implied a poor prognosis ( value < 0.0001). Among the eight genes, KCNE4 and CCL22 were protective factors of CRC, while CD36, CPT1C, SLC2A3, RASGRP2, NFATC1, and CPA3 were risk factors of CRC, as shown in the forest diagram (Figure 6(d)). We further validated the sensitivity and accuracy of the model. AUCs of the model for 3 years, 5 years, OS, and PFI were 0.71, 0.70, 0.73, and 0.66, respectively (Figures 6(e) and 6(f)). Thus, the risk score model could well predict CRC patients’ prognosis, with high sensitivity and accuracy. As shown in the multivariate Cox regression analysis, the model and MMR could become independent factors for CRC prognosis after adjustment of other clinical characteristics (Table 3).
(a)
(b)
(c)
(d)
(e)
(f)
3.4. Eight Genes in the Risk Score Model Are Significantly Associated with CRC Patients’ Prognosis
Box plot depicted the difference in expression patterns of CCL22 (Figure 7(a)), CD36 (Figure 7(b)), CPA3 (Figure 7(c)), CPT1C (Figure 7(d)), KCNE4 (Figure 7(e)), NFATC1 (Figure 7(f)), RASGRP2 (Figure 7(g)), and SLC2A3 (Figure 7(h)) between the high risk score and low risk score. Among them, CCL22 ( value < 2.22-16), CPA3 ( value < 2.22-16), CPT1C ( value = 0.00078), KCNE4 ( value = 0.023), NFATC1 ( value = 0.00062), and SLC2A3 ( value = 0.00081) were differentially expressed between the high and low risk scores. Furthermore, the expression levels of these genes between CRC samples and normal samples were visualized (Figures 8(a)–8(g)). CD36 ( value < 2.22-16), CPA3 ( value < 2.22-16), NFATC1 ( value = 9.1-08), and RASGRP2 ( value < 2.22-16) were highly expressed and SLC2A3 ( value = 0.0015) was lowly expressed in tumor samples. As shown in Figures 9(a)–9(h), low expression of CCL22 ( value = 0.0047) and CPA3 ( value = 0.035) indicated shorter OS time than high expression. Moreover, we found that highly expressed CPT1C ( value = 0.0017), KCNE4 ( value = 0.002), and SLC2A3 ( value = 0.0048) was significantly correlated with poor PFS (Figures 9(i)–9(p)).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
(j)
(k)
(l)
(m)
(n)
(o)
(p)
3.5. The Eight-Gene Model Is in Significant Correlation with Immune Cell Infiltration
The correlation between the model and the infiltrating levels of six immune cells was analyzed. We found that the model was in significant association with the infiltrating levels of six immune cells, including B cell (Figure 10(a); , value = 0.0064) and CD4+T cell (Figure 10(b); , value =4.3-06). However, no significant correlation between the model and CD8+T cell was found in Figure 10(c) (, value = 0.34). Furthermore, there were distinct correlations between the model and dendritic cell (Figure 10(d); , value = 0.0072), macrophage (Figure 10(e); , value = 3.3-05), neutrophil (Figure 10(f); , value = 9.4-05). We also found that the expression levels of the eight genes in the model were significantly correlated with the infiltrating levels of six immune cells, including CCL22 (Figure 11(a)), CD36 (Figure 11(b)), CPA3 (Figure 11(c)), CPT1C (Figure 11(d)), KCNE4 (Figure 11(e)), NFATC1 (Figure 11(f)), RASGRP2 (Figure 11(g)), and SLC2A3 (Figure 11(h)). Moreover, the eight genes were in significant association with markers of immune cells (Supplementary Table 1). These results suggested that the model could be in association with immune cell infiltration.
(a)
(b)
(c)
(d)
(e)
(f)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
3.6. Validation of the Risk Score Using an Independent Dataset
Based on 566 CRC samples from the GSE39582 dataset, the prognostic value of the risk score was validated. The risk score distribution and survival status of CRC patients are shown in Figure 12(a). Heat maps showed the expression differences of CD36, KCNE4, CPT1C, SLC2A3, RASGRP2, NFATC1, CCL22, and CPA3 between the high and low risk scores (Figure 12(b)). As expected, CRC patients with high risk score had a poorer prognosis than those with low risk score (Figure 12(c)). Among the eight genes, CD36, NFATC1, and CCL22 were significantly associated with prognosis of CRC patients (Figure 12(d)). AUCs of the model for 3 years and 5 years were 0.65 and 0.66, respectively (Figure 12(e)), indicating that the risk score could well predict CRC patients’ prognosis. The expression levels of CCL22 (Figure 13(a)), CD36 (Figure 13(b)), CPA3 (Figure 13(c)), CPT1C (Figure 13(d)), KCNE4 (Figure 13(e)), NFATC1 (Figure 13(f)), RASGRP2 (Figure 13(g)), and SLC2A3 (Figure 13(h)) between the high risk score and low risk score were validated based on the 566 CRC samples. Univariate Cox regression analysis results showed that age, KRAS mutation, pathologic T, pathologic N, pathologic M, tumor stage, and risk score were notably associated with CRC patients’ prognosis. After multivariate Cox regression analysis, we found that age, KRAS mutation, pathologic M, and risk score could be independent prognostic factors for CRC (Table 4).
(a)
(b)
(c)
(d)
(e)
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
3.7. The Eight Genes in the Risk Score Are Distinctly Correlated with Molecular Markers of CRC Prognosis
In Figure 14(a), CCL22 was significantly correlated with BRAF mutation ( value = 0.014) and KRAS mutation ( value = 0.041). For CD36, there was a distinct correlation between its expression and KRAS mutation ( value = 0.00034) and MMR ( value = 0.025) in Figure 14(b). CPA3 ( value = 0.0066; Figure 14(c)) and CPT1C ( value = 0.005; Figure 14(d)) had higher expression levels in KRAS mutation. As shown in Figure 14(e), KCNE4 expression was in significant correlation with BRAF mutation ( value = 0.0014), KRAS mutation ( value = 0.049), and MSI ( value = 0.05). NFATC1 expression was prominently correlated with BRAF mutation ( value = 2.2-11), KRAS mutation ( value = 0.00051) and MSI ( value = 1.1-13) in Figure 14(f). In Figure 14(g), RASGRP2 expression was significantly decreased in KRAS mutation. For SLC2A3, we found that there was a distinct correlation between its expression and BRAF mutation ( value = 0.0011) and MSI ( value = 0.00013) in Figure 14(h).
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
4. Discussion
In TME, stromal and immune cells are involved in the development of CRC. In this study; using the ESTIMATE algorithm, stromal and immune scores were calculated. A significant correlation between the immune score and CRC patients’ prognosis was observed. Both the stromal score and immune score were in significant correlation with clinical characteristics of CRC patients. Furthermore, we identified differentially expressed stromal and immune genes for CRC. Functional enrichment analysis results suggested that these genes were positively related with immune-related pathways like cytokine-cytokine receptor interaction [16, 17] and chemokine signaling pathway [18, 19].
Individual prognosis for CRC patients varies widely. Individual genes often cannot accurately predict the prognosis of patients with CRC. Genes in most prognostic risk scores are screened via differential expression analyses [20–22]. Yet, there are few prognostic models associated with CRC immune infiltration. Therefore, in this study, we selected eight differentially expressed stromal and immune genes related to prognosis for constructing a risk score model. However, focusing only on CRC-related immune-related genes may limit its clinical value. For this reason, through multivariate regression analysis, after adjustment of the clinical characteristics of CRC, we assessed the association between the risk score and CRC prognosis. The results showed that the model may become an independent prognostic factor for CRC. Our risk score exhibited well efficiency in predicting CRC patients’ prognosis. Therefore, the risk score model possessed potential prognostic value, which was confirmed using an independent dataset. Among the eight genes in the model, both in the discovery and independent datasets, CCL22 was a protective factor of CRC, while CD36 and NFATC1 were two risk factors of CRC. However, other genes exhibited inconsistent results in the two datasets. This is partly due to the heterogeneity of the samples in the two datasets. Patients in the same pathological stage have different prognosis. Both in the discovery and independent datasets, CCL22 and CPA3 were lowly expressed and KCNE4, NFATC1, and SLC2A3 were highly expressed in the high-risk samples compared to the low-risk samples. However, there were inconsistent results about other genes between the high- and low-risk samples in the two datasets, partly due to the heterogeneity of the samples, different sequencing platforms, different background correction and normalization methods and so on. Thus, it is unreliable to predict CRC patients’ prognosis by an individual gene. However, our risk score composed of these genes may accurately suggest the patient’s prognosis.
As described in a previous study, high CCL22 expression was found in CRC tissues [23]. Recent study has found that CCL22 secreted by M2 macrophages could mediate CRC 5-FU-mediated chemoresistance [24]. Furthermore, it has been reported that CCL22 was in significant correlation with the infiltrating levels of different T cell subsets for CRC [25]. Our results showed that CD36 was significantly downregulated in CRC tissues compared to normal tissues, which was validated in vitro and in vivo [26]. Genome-wide DNA methylation analysis revealed that hypermethylation of CD36 could contribute to its low expression [27]. Fang et al. found that CD36 expression gradually decreased from adenoma to cancer and CD36 loss implied a poor prognosis in patients with CRC [28]. NFATC1 was deregulated in CRC tissues, which was consistent with previous findings [29]. In vitro, its overexpression significantly promoted CRC cell invasion and metastasis [30]. Kumar et al. reported that NFATC1 indicated poor survival outcomes of CRC patients [31]. High SLC2A3 expression was observed in CRC tissues and its high expression indicated a poor prognosis, consistently with previous research [32, 33]. Furthermore, downregulated CPA3 and RASGRP2 and upregulated CPT1C and KCNE4 were found in CRC tissues, which implied poor prognosis.
As for immune cell infiltration, we found that the eight genes in the risk score model were moderately correlated with the infiltering levels of CD4+T cell, dendritic cell, macrophage, and neutrophil. It has been confirmed that TME affects the efficacy of immunotherapy, and immune cells in TME possess predictive value for immunotherapy treatment [34–36]. Increasing genes have been shown to participate in the regulation of immune cells [37–40]. Therefore, our risk score model could possess potential value to predict CRC patients’ prognosis, and the eight genes could become promising immunotherapeutic targets, which deserve further study.
Our correlation analysis results confirmed that the eight genes in the risk score were distinctly correlated with molecular markers of CRC prognosis. However, our study has several limitations. First, our retrospective study limited the application of this risk score. Second, the heterogeneity of the immune microenvironment would inevitably contribute to result bias. Therefore, it is necessary to validate our findings in a prospective clinical study.
5. Conclusion
In this study, we conducted an immune-related prognostic model for CRC on the basis of stromal and immune scores. The model had well predictive efficacy for CRC patients’ prognosis. Our findings could provide novel biomarkers for predicting the prognosis of CRC patients and developing individualized immunity therapy strategies.
Abbreviations
| CRC: | Colorectal cancer |
| TCGA: | The Cancer Genome Atlas |
| ROC: | Receiver operating characteristic |
| OS: | Overall survival |
| PFI: | Progression-free survival interval |
| AUC: | Area under the ROC |
| TME: | Tumor microenvironment |
| FC: | Fold change |
| GO: | Gene Ontology |
| KEGG: | Kyoto Encyclopedia of Genes and Genomes |
| CC: | Cellular component |
| MF: | Molecular function |
| BP: | Biological process |
| PPI: | Protein-protein interaction |
| STRING: | The Search Tool for the Retrieval of Interacting Genes. |
Data Availability
The (data type) data used to support the findings of this study are included within the supplementary information file(s).
Conflicts of Interest
The authors declare no conflicts of interest.
Supplementary Materials
Supplementary Table 1: the correlation between the eight genes in the risk score model and markers of immune cells. (Supplementary Materials)