A prognostic fibroblast-related risk signature in colorectal cancer

Training set collection

Processed RNA-seq data of TCGA CRC and paired clinical information (version: September 8, 2017) were obtained from UCSC Xena webserver (https://xenabrowser.net/). After data quality control (QC) with excluding missing values, a total of 55 normal samples and 375 tumor samples data set were downloaded, and 375 tumor samples data set were collected as training set in this study. The log2 (Counts + 1) value and log2 (FPKM +1) value were both downloaded, if not explained, log2 (FPKM + 1) value means the gene expression in this study.

Testing set collection

Processed microarray data of CRC and related meta-data were downloaded from NCBI Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo) with accession code GSE39582. By performing the same QC as training set, we have a validation set with 566 samples.

Quantification of tumor immune and stromal content

MCP-counter is a wide used computational tool for quantification of tumor immune and stroma content from bulk RNA-seq data [12], such as analysis of immune cells associated with immunotherapy response [13], identification of different TME subtypes [14]. For the ten major cell types in TME, MCP-counter can efficiently derive the semi-quantitative scores which predict the enrichment of specific cell types in a sample based on well-defined marker genes [15]. Thus, we could obtain the quantification value of a specific cellular content and perform inter-sample comparisons in the next analysis.

In this study, we applied R package MCP-counter (https://github.com/ebecht/MCPcounter) to obtain the abundance scores of ten major stromal and immune cell types (endothelial cells, fibroblasts, CD3+ T cells, CD8+ T cells, cytotoxic lymphocytes, natural killer cells, B lymphocytes, monocytic lineage cells, myeloid dendritic cells and neutrophils). And the enrichment scores of TME were hierarchically clustered by R package pheatmap (Figure 1A) and grouped by the fibroblast enrichment (Figure 1B). Also, after the construction of risk model, we compared the quantification of different cell types across samples grouped into high risk and low risk. Another deconvolution-based software EPIC was also applied for the cellular composition analysis as previously described [16], and the implementation of MCP-counter was documented in detail elsewhere [12].

Figure 1. Analysis of TME in TCGA database and analysis of OS in high- and low-fibroblast groups. (A) In TCGA, all types of immune cells and mesenchymal cells in adjacent normal tissues are higher than that in tumor tissues. (B) We divided the 375 patients into a high-fibroblast group and a low-fibroblast group according to the median of the fibroblast content, the OS in the high-fibroblast group is shorter than in the low-fibroblast group. (C) The abundance of immune cells in the high-fibroblast group was significantly higher than that in the low-fibroblast group.

Differential expression analysis

For differentiated expressed genes (DEGs) of TCGA with high/low fraction of fibroblasts or high/low risk fibroblast score, we used raw count value of genes as input with R package DESeq2 [17]. Briefly, we performed deseq() function with default parameters, then we defined the DEGs which had higher absolute log2 (fold-change) value than 1.5 and tested significantly (p.adj value < 0.05). Out of 3845 DEGs, we chose 1720 genes which had expression in validation set as well. The same criteria was applied to DEG of high/low risk samples.

Development and validation of the prognostic signature

The OS was defined as the time from diagnosis to death or the last follow-up date. Univariate Cox was used to obtain 222 prognostic genes in training set through traversing all DEGs in fibroblast high/low condition, and LASSO-penalized Cox regression was used to construct optimal prognostic risk signatures in training group. The COX regression model with the LASSO penalty successfully achieved compression and selected 14 fibroblast-related genes simultaneously. The risk score formula was as follows: $\text{risk score = }{\displaystyle {\sum }_{i=1}^n\exp i\beta i}$ where exp represented the gene expression value while β represented the LASSO coefficient. This model was utilized to estimate the OS of each patient in the training dataset and testing dataset. The predictive potential of the signature was evaluated via area under the curve (AUC) value of the receiver operating characteristic (ROC) curve. The performance of proposed risk score was also tested by calibration curve and DCA curve in both datasets, which were implemented by R package rms and R package rmda, respectively.

Survival analysis

R package survival and survminer were applied to investigate the difference of prognosis between two groups by log-rank test. By using median expression as a cutoff, we obtained 222 prognostic genes in training set through traversing all DEGs in fibroblast high/low condition. As for prognosis analysis of high risk/low risk in training set, validation set or subgroup analysis of two datasets, we set the median risk score as a metric of risk.

Enrichment analysis

For pathway enrichment and functional annotation, R packages clusterprofiler and GSEA software were performed [18, 19]. Basically, we investigated upregulated genes (or downregulated genes) whether significantly enriched in particular pathways or terms refer to GO (http://geneontology.org) database. Next, we used a ranked genelist to GSEA algorithm with default parameters, in order to find out which significant pathway (q < 0.25) involved the specific condition (i.e., high risk versus low risk) refer to HALLMARK and C2 from MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb).

Genomic alteration analysis

The related public mutation and copy number alteration data were surveyed and analyzed in online webserver cBioPortal (https://www.cbioportal.org/) [20]. All analysis and plot were followed as the webpage UI guide. Network of gene interactions was obtained from GeneMANIA (https://genemania.org).

Statistical analysis

All statistical analyses in this study were performed in R version 4.04. P < 0.05 was considered as statistical significance. Significance in comparisons of gene expression and putative microenvironment components in figures was calculated by the Wilcoxon Rank-sum test. ROC curves were plotted by R package ROCR. Spearman correlation was the metrics between fibroblast contents and modeled risk scores.

Data availability statement

All data for this study are available from the corresponding public database.

Overall survival correlation

The quantification of tumor stroma is critical to unveil the multi-faceted role of the TME, which may be involved in affecting overall survival (OS) of CRC patients. From TCGA bulk transcriptomic data, we compared the quantification of cellular components of immune cells and mesenchymal cells between 375 cancerous tissues and 55 adjacent normal tissues using the MCP-counter, a tool to estimate cellular content by scoring marker gene expression, the result showed that the content of all types of immune cells and mesenchymal cells in normal tissues are higher than that in cancerous tissues (Figure 1A). Therefore, the TME must play a key role in the occurrence and development of tumors. In view of the crucial roles of the fibroblasts in tumor progression, then we divided the 375 patients into a high-fibroblast group and a low-fibroblast group according to the median of the fibroblast content. From analyzing the relationship between OS and the content of fibroblasts in CRC patients, we observed OS was significantly longer in the low-fibroblast group than in the high-fibroblast group (Figure 1B). Because the TME has an important impact on patients with colorectal cancer, we further investigated the abundance of tumor infiltrating lymphocytes between low- and high-fibroblast groups. We found the abundance of immune cells in the high-fibroblast group was significantly higher than that in low-fibroblast group (Figure 1C).

Development of a fibroblast-related risk signature

To further explore the differences between the high- and low-fibroblast groups, 3845 differentially expressed genes (DEGs) between the two groups were identified in TCGA (Figure 2A). Then, univariate COX analysis was used to identify 222 prognostic genes commonly expressed in TCGA and GSE39582, including 216 risky genes and 6 protective genes. The change in trajectory of each variable was plotted (Figure 2B). We utilized 10-fold cross-validation to construct the model, and showed the confidence interval under each lambda (Figure 2C). When lambda equaled 0.02296418, the model reached the optimal value, and 14 variables were selected (Table 1). Then 14 genes identified by LASSO-COX analysis were used to establish a fibroblast-related risk signature of CRC in TCGA. The formula of the fibroblast-related risk signature was as follows: Risk score = ∑ (βn × expression of gene n). Using the median risk score as a cut-off, CRC patients were divided into low-risk group and high-risk group in TCGA, the OS of the low-risk group was significantly longer than that of the high-risk group (P < 0.05) (Figure 2D). To demonstrate the universality of the fibroblast-related signature, the GSE39582 cohort was used to validate the signature. The results confirmed that the fibroblast-related risk signature was significantly associated with the OS of CRC patients in the testing group (Figure 2E). And the analyses of the two databases also demonstrated the validity of the signature (Figure 2F).

Figure 2. Prognostic signature based on 14 hub genes. (A) Differentially expressed genes in TCGA between low-fibroblast group and high-fibroblast group. (B, C) LASSO-COX analysis of prognostic genes. (D) The OS is shorter in the high-risk group than in the low-risk group in TCGA. (E) The OS is shorter in the high-risk group than in the low-risk group in GSE39582. (F) The OS is shorter in the high-risk group than in the low-risk group in TCGA and GSE39582.

Stratified survival assays

To verify the stability of the model and its ability to predict survival in different clinical subgroups, the fibroblast-related risk signature was used to perform stratified survival analysis of CRC patients with survival information in the training and testing groups after adjusting for age, gender, stage, and TNM stage. The clinical characteristics of the training and testing groups were also shown (Table 2). Patients in the different subgroups with missing grouping information were excluded. The results showed that the risk signature could predict the prognosis of CRC patients in most cases, and its prognostic stability was superior in patients with advanced non-metastatic cancer (Figure 3).

Figure 3. Stratified survival analysis adjusted to age, gender, stage, and TNM stage. All CRC patients in the training and testing groups were summarized in the stratified survival analysis. 68 years old was the median age of 937 CRC patients.

Validation of the fibroblast-related risk signature

To determine whether the fibroblast-related risk signature was a clinically independent prognostic factor, we performed univariate and multivariate analyses and found that age, stage, and risk score were independent prognostic factors in the training group (Figure 4A and 4C). Similar analyses in the testing group GSE39582 led to the same conclusion (Figure 4B and 4D). The results demonstrated that the fibroblast-related signature could be used to predict the survival of CRC patients.

Figure 4. The validity of prognostic signature and the relationship between the risk score and fibroblast content. (A) Univariate COX regression analysis in TCGA. (B) Univariate COX regression analysis in GSE39582. (C) Multivariate COX regression analysis in TCGA. (D) Multivariate COX regression analysis in GSE39582. (E) The receiver operating characteristic (ROC) curve and the areas under the curve verified the accuracy of prognostic signature in the training and testing groups. (F) The risk score is associated with the fibroblast content.

To evaluate the accuracy of this signature, a Receiver Operating Characteristics (ROC) curve was plotted and the area under the curve (AUC) was calculated. Compared to other existing clinicopathological factors, the results showed that the signature is an effective index to predict the OS of CRC patients in both the training group and the testing group (Figure 4E). Moreover, calibration plots indicated that in comparison with an ideal model, the signature had a similar performance (Supplementary Figure 2A). The results of DCA also demonstrated that our signature had high potential for clinical utility (Supplementary Figure 2B). A higher risk score was related to a higher content of fibroblasts, suggesting the risk score was positively correlated with the number of fibroblasts (Figure 4F).

Clinical relevance of the risk signature and fibroblast content

The relevance of the fibroblast content relative to clinical traits including age, gender, stage, and TNM status was assessed using the TCGA database. The fibroblast content was significantly increased in the advanced stage, advanced T stage, and positive lymph node metastasis groups (Figure 5A). Assessment of the relationship between risk score and clinical traits in TCGA showed that the score was significantly higher in the advanced stage, advanced T stage, positive lymph node metastasis, and distant metastasis groups (Figure 5B).

Figure 5. Clinical relevance of fibroblast content and clinical relevance of risk score. (A) The fibroblast content was significantly increased in the advanced stage, advanced T stage, and positive lymph node metastasis groups. (B) The risk score was significantly higher in the advanced stage, advanced T stage, positive lymph node metastasis, and distant metastasis groups.

Mutation and copy number alteration analysis and protein-protein interactions of hub genes

In addition to the analysis at the transcriptome level, we examined the role of marker genes at genome level. For this purpose, the cBioPortal was used to analyze the mutations of these genes in CRC. At the genome level, the mutation frequencies of the 14 genes were not significant; the main mutation types were missense mutations and amplifications. MMP19, DACT1, SCG2, and CHST3 were the most frequent CNAs with a 3% mutation rate among the 14 hub genes, whereas RETNLB was the least frequent CNA with a 0.6% mutation rate (Figure 6A).

Figure 6. Mutation and copy number alteration (CNA) analysis and protein-protein interactions (PPI) of hub genes. (A) Mutation and copy number alteration of each hub gene. (B) The 20 functional similar genes were located in the outer circle, while hub genes were located in the inner circle. The color of nodes was related to the protein function while line color represented the type of protein interaction.

To predict functionally similar hub genes, we used GeneMANIA to obtain the 20 most similar hub genes. The hub genes were located in the inner circle and the predicted genes were in the outer circle. The results showed that their functions were mostly related to cell chemotaxis and protein kinase binding which is essential for tumor progression and metastasis (Figure 6B).

Molecular characteristics and pathways of the fibroblast-related risk signature

Gene Set Enrichment Analysis (GSEA) of hallmark gene sets and pathway gene sets in the low-risk and high-risk groups from the training cohort was performed, and the four most enriched characteristics and pathways were selected. For the hallmark gene sets, HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION, HALLMARK_HYPOXIA, and HALLMARK_INFLAMMATORY_RESPONSE were significantly upregulated in the high-risk group, and HALLMARK_OXIDATIVE_PHOSPHORYLATION was significantly upregulated in the low-risk group (Figure 7A). This suggested that the high-risk group has a higher progression potential and is characterized by inflammation and hypoxia. For pathway gene sets, KEGG_ECM_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, KEGG_LEUKOCYTE_TRANSENDOTHELIAL_MIGRATION, and REACTOME_EXTRACELLULAR_MATRIX_ORGANIZATION were significantly upregulated in the high-risk group (Figure 7B). Gene ontology (GO) biological process enrichment, GO cellular component enrichment, and GO molecular function enrichment were performed in the TCGA cohort. Functions related to the ECM were most enriched (Figure 8A–8C). For GO cellular component enrichment, some of the upregulated DEGs in high-risk group were related to the most enriched collagen-containing ECM (Figure 8D). For GO biological process enrichment, the relationship of different enriched biological processes was also established (Figure 8E).

Figure 7. GSEA between the high- and low-risk groups. (A) HALLMARK_EPITHELIAL _MESENCHYMAL_TRANSITION, HALLMARK_HYPOXIA, and HALLMARK_INFLAMMATORY_RESPONSE were significantly upregulated in the high-risk group, and HALLMARK_OXIDATIVE_PHOSPHORYLATION was significantly upregulated in the low-risk group. (B) KEGG_ECM_RECEPTOR_INTERACTION, KEGG_FOCAL_ADHESION, KEGG_LEUKOCYTE_TRANSENDOTHELIAL_MIGRATION, and REACTOME_EXTRACELLULAR_MATRIX_ORGANIZATION were significantly upregulated in the high-risk group.

Figure 8. GO enrichment analysis. (A) GO biological process enrichment. (B) GO cellular component enrichment. (C) GO molecular function enrichment. (D) The up-regulated DEGs were related to collagen-containing extracellular matrix. (E) The relationship of different enriched biological processes.

Correlations with immune cells and immune responses

To determine whether the signature indicated the immune cell infiltration in the TME of CRC patients, we used MCP-counter to analyze the TME. The results showed that B lineage cells, T cells, and NK cells were negatively correlated with the risk score, whereas endothelial cells and monocytic lineage content were positively correlated with the risk score (Figure 9A). In the high-risk and low-risk groups, main immune cell lineages that suppress tumors including B lineage cells, T cells, and NK cells were higher in the low-risk group (P < 0.05), whereas endothelial and monocytic lineage cells were higher in the high-risk group (P < 0.05) (Figure 9B). EPIC was also used to confirm the correlation between TME and risk score, the result showed a higher correlation between TME and risk score (Supplementary Figure 3). CD8+T cells are important anti-tumor immune cells, and the content of CD8+T cells did not differ significantly between the two groups. Analysis of common immune checkpoint genes showed high expression levels of HAVCR2 in the high-risk group, whereas GZMA was expressed at high levels in the low-risk group. This suggested that the content of exhausted CD8+T cells was higher in the high-risk group (Figure 9C). Analysis of inflammatory cytokine genes showed that the high-risk group had high expression of TGFβ1, whereas the low-risk group had high expression of TNF, implying an immunosuppressive TME in the high-risk group (Figure 9D).

Figure 9. Correlations with immune cells and immune responses. (A) The B lineage cells, T cells, and NK cells were negatively correlated with the risk score, whereas endothelial cells and monocytic lineage content were positively correlated with the risk score. (B) The main immune cell lineages including B lineage, T cells and NK cells are higher in low-risk group (P < 0.05) while endothelial cells and monocytic lineage are higher in high-risk group. (C) The immune checkpoint genes (PDCD1, CTLA4, HAVCR2, GZMA) expression between high-risk group and low-risk group. (D) The inflammation cytokine genes (IL10, TGFβ1, IFNG, TNF) expression between high-risk group and low-risk group.