We searched the NCBI gene expression omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) using the keywords “colon cancer fibroblast,” “cancer-associated fibroblast,” “CAFs” “stroma cells,” “colon tumor stroma cells” and “tumor stroma,” and identified two cancer-associated fibroblasts gene expression datasets GSE46824 (n = 34) [18] and GSE70468 (n = 14) [19] of colonic CAFs. In GSE70468, we selected the 7 CAFs samples and 7 colon normal fibroblasts samples and removed the other samples from this dataset. Besides, we downloaded the TCGA COAD cohort and normalized the data into log₂-base transformation (https://gdc-portal.nci.nih.gov/) [20]. For analyzing the survival differences in the TCGA COAD dataset (https://gdc-portal.nci.nih.gov/), we used the gene expression profiling interactive analysis (GEPIA) [21] tool (http://gepia2.cancer-pku.cn/#index). Furthermore, we used PrognoScan-based [22] colon cancer data, the GSE17536 (n = 177) [23], for validating the survival-associated genes found in TCGA data. Finally, we used the colon tumor-stromal dataset GSE35602 [24] for verifying the expression levels of key genes. The summary of the GEO datasets is presented in Table 1.

The differential expression was screened by using GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r), which is an interactive web tool to identify the DEGs. GEO2R tool identified the significant DEGs by utilizing the GEOquery [25] and limma R packages [26] from the Bioconductor project (http://www.bioconductor.org/). The R package “limma” was employed for identifying the significant DEGs between CAFs and normal colonic fibroblasts (NCFs) [26]. We identified the DEGs with a threshold absolute value of fold change (LogFC) > 0.585 and a p-value ≤0.05 [27, 28]. Finally, we identified the common DEGs between both datasets. We utilized the online tool “Calculate and draw custom Venn diagrams” (http://bioinformatics.psb.ugent.be/webtools/Venn/) for identifying common DEGs in both datasets.

We performed a gene-set enrichment analysis of the DEGs by using the gene-set enrichment analysis (GSEA) [29]. We inputted all significant commonly found upregulated and downregulated DEGs into the GSEA tool for identifying deregulated pathways. In the GSEA, we selected the KEGG pathways panel to identify the significant pathways. The KEGG [30] pathways significantly associated with the upregulated DEGs and the downregulated DEGs were identified, respectively. The p-value <0.05 was considered significant when separating the pathways.

To better know the relationship among these identified DEGs, the PPI network was established using the STRING-based analysis [31]. To identify the rank of hub genes, we used Cytoscape plug-in tool cytoHubba [32]. Hub genes were identified based on the degree of interactions with neighbour genes. We set the minimum required interaction score as 0.40 for identifying the PPI of DEGs. Hub genes were defined as a gene that was connected to a minimum of 10 other DEGs in the PPI. We visualize the PPI networks by utilizing the Cytoscape 3.6.1 software [33]. A Cytoscape plug-in molecular complex detection (MCODE) was employed to detect the modules from the PPI network [34]. This is a computational method for searching molecular complexes in large protein interaction networks. We identified the significant modules based on the MCODE score and node number. The threshold of the MCODE was node score cut-off: 0.2, haircut: true, k-core: 2, and maximum depth from Seed: 100.

We compared the overall survival (OS) and the DFS of colon cancer patients. Kaplan-Meier survival curves were used to show the survival differences between the high expression group and low expression groups. The survival significance of all DEGs in the TCGA COAD cohort was analysed using GEPIA [21] databases. GEPIA is a web-based tool to deliver fast and customizable functionalities based on TCGA data. Furthermore, we used PrognoScan-based [22] GSE17536 (n = 177) [23] for validating the survival-associated genes found in the TCGA COAD dataset. PrognoScan is a database for meta-analysis of the prognostic value of genes. Cox regression p-value <0.05 was considered as significant when comparing the survival between the two groups.

ESTIMATE is an algorithmic tool based on the R package for predicting tumor purity, immune score (predicting the infiltrations of immune cells), and stromal score (predicting the infiltrations of stromal cells) which uses the gene expression profiles of 141 immune genes and 141 stromal genes [35]. The presence of infiltrated immune cells and stromal cells in tumor tissues were calculated using related gene expression matrix data, represented by immune score and stromal score, respectively [35]. Then we calculated the correlations of key genes with immune scores and stromal scores. The threshold value of correlation is R > 0.30, and p-value is not less than 0.001 (Spearman’s correlation test).

One of the extension packages of GSEA, single-sample gene-set enrichment analysis (ssGSEA) was used to identify the enrichment scores of immune cells for each pairing of a sample and gene set in the tumor samples [36]. We collected the marker gene set for immune signatures and utilized each gene set to quantify the ssGSEA scores of specific immune signatures [37–40]. We identified the ssGSEA scores of CD8⁺ T cells, CD4⁺ regulatory T cells, NK cells, TAM, macrophages, M2 macrophages, Tregs, T cell exhaustion, myeloid-derived suppressor cells (MDSCs), and CAFs. Moreover, we used the metastasis-promoting gene set for identifying the ssGSEA scores of metastasis-promoting genes [40, 41]. All of the marker genes are displayed in Supplementary Table S1. Then we investigated Spearman’s correlation between the ssGSEA scores and specific survival-associated genes. The absolute value of this correlation was greater than 0.30 with a p-value less than 0.001.

To assess diagnostic values of the survival-associated genes, the receiver operating characteristic (ROC) curve was plotted and the area under the ROC curve (AUC) was calculated using the “pROC” R package [42] to evaluate the capability of distinguishing cancer-associated fibroblasts and normal tissues. The “pROC” R package can visualize, smooth, and compare ROC curves. The AUC can be compared with statistical tests based on U-statistics or bootstrap. A greater AUC value of individual genes indicated the differences between tumor and normal samples, and the key gene of AUC > 0.5 in the CAFs datasets was defined as a diagnostic efficiency of the gene [43]. To verify the survival-associated genes in a colorectal cancer stromal tissue, we selected GSE35602 [24] to distinguish the expression levels of key genes between non-tumor and tumor samples. We identified the ssGSEA scores (prediction of the content) of CAFs by using the marker genes (Supplementary Table S1) of CAFs. Then, we divided the TCGA COAD samples into the high-CAFs group versus the low-CAFs group based on the median value of ssGSEA scores of CAFs.

We used the R software version 4.0.1 for all statistical analysis. In the Log-rank test, p < 0.05 was considered statistically significant for survival analysis. To investigate the correlation of survival-associated genes, Spearman’s correlation between the ssGSEA scores and specific survival-associated genes was performed (p < 0.001). We used Welch’s t-test for identifying the significance of specific genes between the two groups of samples (p ≤ 0.05). We utilized the R package ggplot2 for the graphical presentation of the Heatmap and correlation graph.

We identified 246 commonly differentially expressed genes (DEGs) between the GSE46827 and GSE70468 datasets. The principal component analysis (PCA) plot of samples is shown in Figures 1A,B and the commonly upregulated and downregulated genes in CAFs and normal colonic fibroblasts (NCFs) are shown in Figures 1C,D respectively. We found 72 upregulated (Figure 1C; Supplementary Table S2) and 174 downregulated (Figure 1D; Supplementary Table S3) genes in the CAFs when compared with NCFs. CDH13, DSP, EFNB2, IL7R, KIAA1217, KRT19, NTM, PDLIM3, SRGN, SYNPO2L, TGFB2, and TNFSF4 were the most significantly upregulated genes in colon CAFs (Supplementary Table S2), while A2M, AOC3, ADH1B, ASPA, CHRDL1, COL14A1, CYP24A1, FBN2, FBLN1, GREM2, IL33, NOVA1, PRKCH, SEPP1, and SMOC2 are the most downregulated genes over the selected datasets (Supplementary Table S3).

The enriched upregulated and downregulated biological pathways were identified by using the GSEA tool (Figures 2, 3). GSEA identified 28 KEGG pathways significantly associated with common upregulated 72 DEGs (Supplementary Table S4). The pathways are mainly involved with cancer or malignant tumor (colorectal cancer, pancreatic cancer, pathways in cancer, melanoma, chronic myeloid leukemia, small cell lung cancer, and prostate cancer) and cellular signaling and developments (regulation of actin cytoskeleton, TGF-beta signaling pathway, focal adhesion, p53 signaling pathway, adherens junction, gap junction, MAPK signaling pathway, cell cycle, cell adhesion molecules (CAMs), and Wnt signaling pathway) (Figure 2).

Besides, we identified the 54 KEGG pathways that are significantly linked with commonly found downregulated 174 DEGs (Supplementary Table S5). The pathways are involved with immune regulation (such as complement and coagulation cascades, chemokine signaling pathway, and cytokine-cytokine receptor interaction) and cellular metabolism (glycolysis/gluconeogenesis, fatty acid metabolism, tyrosine metabolism, purine metabolism, beta-Alanine metabolism, metabolism of xenobiotics by cytochrome P450, histidine metabolism, propanoate metabolism, pyruvate metabolism, limonene and pinene degradation, retinol metabolism, drug metabolism—cytochrome P450, ascorbate and aldarate metabolism, butanoate metabolism, tryptophan metabolism, lysine degradation, valine, leucine and isoleucine degradation, glycerolipid metabolism, arginine and proline metabolism, and arachidonic acid metabolism) (Figure 3).

We investigated the PPI of all significant CAFs-derived DEGs. The PPI information of STRING is inputted into the Cytoscape for identifying and visualizing the hub genes and significant clusters. We identified 15 hub genes (minimum degree of interaction is 10 with other DEGs) which included 3 upregulated hub genes (PLAUR, RAC2, and TGFB2) and 12 downregulated hub genes (A2M, ADCY3, ADCY4, ADCY9, BIN1, CFD, CLU, CXCL12, LDLR, PPARG, PTK2B, and VTN) (Table 2).

We found that 183 DEGs were involved in PPI (Supplementary Table S6). The other DEGs were not involved in the PPI network. We investigated the significant cluster-based analysis of this PPI network. The MCODE-based analysis identified 9 clusters from the original PPI networks. The description of MCODE-derived clusters is illustrated in Table 3. The top significant cluster 1 included 7 nodes and 21 edges (Table 3). We identified the significantly enriched KEGG pathways for the identified 9 clusters by using the GSEA. Interestingly, we found that seven clusters are associated with the enrichment of KEGG pathways (p < 0.05). Gene set of cluster 1 is associated with the enrichment of complement and coagulation cascades (p = 1.67e-7). Gene set of cluster 2 is mainly involved with immune regulation and cellular signaling (Table 3). In addition, the gene set of cluster 3 is mainly associated with the regulation of metabolism (Table 3). Also, Cluster 4, Cluster 6, Cluster 8, and Cluster 9 are associated with the enrichment of KEGG pathways (Table 3).

We investigated the survival significance of the CAFs-derived all significant common DEGs (72 upregulated and 174 downregulated genes) in TCGA COAD data. Our analysis revealed that the upregulated group of CAFs-derived upregulated genes included CTHRC1, EMB, FOXS1, LYPD1, NTM, PDGFC, PDLIM3, SLC16A3, SYNPO2L, and TNFSF4 (Figure 4A) and the lower expression group of downregulated genes included ABCA5, FBN2, IMPA2, and TIMP4 (Figure 4B) are significantly correlated with shorter survival time of colon cancer patients (Figure 4). Moreover, we verified this result in a GEO dataset GSE17536 (n = 177) of colon cancer. We found that upregulated group of CTHRC1, NTM, PDGFC, PDLIM3, and SLC16A3 genes are consistently correlated with the shorter survival of colon cancer patients in GSE17536 (Figure 4C). In contrast, the downregulated group of FBN2 gene is consistently correlated with the shorter survival of colon cancer patients in GSE17536 (Figure 4C), indicating the tumor-suppressive role of FBN2 in colon cancer.

Since the upregulated group of CTHRC1, NTM, PDGFC, PDLIM3, and SLC16A3 and the downregulated group of FBN2 genes are consistently correlated with the shorter survival time in TCGA and GEO datasets, we investigated the association of these genes with immune score and stromal scores in colon cancer. Interestingly, we found that the expression levels (log2 transformation) of CTHRC1, NTM, PDGFC, and PDLIM3 are moderately correlated with immune scores (Spearman’s correlation test, p < 0.001) (Figure 5A). In addition, the expression levels (log2 transformation) of CTHRC1, NTM, PDGFC, PDLIM3, and FBN2 are strongly correlated with the stromal score (Spearman’s correlation test, p < 0.001) (Figure 5B). In contrast, SLC16A3 is not correlated with immune scores and stromal scores in the TCGA-COAD cohort (R > 0.30 and p < 0.01). This indicated that the expression levels of CTHRC1, FBN2, NTM, PDGFC, and PDLIM3 genes are associated with the modulation of immune and stromal activity in the colon cancer tumor microenvironment.

Since survival time of patients is correlated with immunological responses in human cancers [44], we investigated the correlation of six survival-associated genes (CTHRC1, NTM, PDGFC, PDLIM3, SLC16A3, and FBN2) with the several immune stimulatory and inhibitory signatures including CD8⁺ T cells, CD4⁺ regulatory T cells, NK cells, TAM, macrophages, M2 macrophages, Tregs, T cell exhaustion, and MDSCs. We found that the expression levels of upregulated CTHRC1, NTM, PDGFC, and PDLIM3 are positively correlated with ssGSEA scores of TAMs, macrophages, M2 macrophages, Tregs, T cell exhaustion, and MDSCs (Spearman’s correlation test, p < 0.001) (Figure 6). Besides, SLC16A3 is positively correlated with the infiltration of MDSCs (Spearman’s correlation test, p < 0.001) (Figure 6). In addition, the expression level of FBN2 is correlated with ssGSEA scores of TAMs, macrophages, M2 macrophages, Tregs, and MDSCs (Spearman’s correlation test, p < 0.001). Interestingly, the CD8⁺ T cells, CD4⁺ regulatory T cells, and NK cells are not significantly correlated with the expression levels of CTHRC1, NTM, PDGFC, PDLIM3, SLC16A3, and FBN2.

We analysed the expression level differences of five shorter survival-associated genes (CTHRC1, NTM, PDGFC, PDLIM3, and SLC16A3) and one longer survival-associated gene (FBN2) among the normal fibroblast (n = 9), primary tumor fibroblast (n = 14), and metastatic tumor fibroblast (n = 11) samples of GSE46824 [18]. Interestingly, we found that CTHRC1, NTM, PDGFC, PDLIM3, and SLC16A3 are gradually upregulated from normal to a metastatic tumor, indicating the association of these genes in tumor progression (Welch’s t-test, p < 0.05) (Figure 7). Besides, the expression of FBN2 is gradually downregulated from normal to the metastatic tumor, indicating that the downregulation of this gene is associated with tumor progression (Welch’s t-test, p < 0.05) (Figure 7). Since the expression levels of CTHRC1, FBN2, NTM, PDGFC, PDLIM3, and SLC16A3 are significantly differentiated among the normal samples, primary tumor, and metastatic tumor, we anticipated that these genes are associated with metastasis in colon cancer.

To prove this hypothesis, we investigated the correlation of these genes with the metastasis-promoting gene set scores. We identified the ssGSEA scores of metastasis-promoting genes and identified the correlation with the expression levels of these six gene signatures (Spearman’s correlation test, p < 0.001). We found that the expression levels of CTHRC1, FBN2, NTM, PDGFC, and PDLIM3 are positively correlated with the ssGSEA scores of metastasis-promoting genes in the TCGA COAD cohort (Figure 8).

We speculate that these six genes (CTHRC1, FBN2, NTM, PDGFC, PDLIM3, and SLC16A3) have diagnostic value in colonic CAFs. We used the GSE46824 dataset to validate our hypothesis, and the results showed that the ROC curve of the expression levels of these six genes (CTHRC1, FBN2, NTM, PDGFC, PDLIM3, and SLC16A3) showed excellent diagnostic value for colonic CAFs and normal colonic fibroblast (Figure 9A). Since CAFs are major regulatory components of tumor stroma [45], we validated the expression levels of survival-associated key genes between colon tumor stroma and normal colon stroma. Interestingly, we found that three upregulated survival-associated genes were also upregulated in colon tumor stroma versus normal colon stroma (p < 0.05) (Figure 9B). In addition, downregulated FBN2 is consistently downregulated in colon tumor stroma (Figure 9B). Moreover, we divided the TCGA COAD samples into the higher-CFAs group versus the lower-CAFs group to identify the expression level of CTHRC1, FBN2, NTM, PDGFC, and PDLIM3. Interestingly, we revealed that the expression of CTHRC1, FBN2, NTM, PDGFC, and PDLIM3 was highly expressed in the higher-CAFs group (Figure 9C). However, SLC16A3 is not significantly altered between the higher-CFAs group versus the lower-CAFs group. It indicates that the aberrant expression of CAFs-derived CTHRC1, NTM, and PDGFC is critically involved in the pathogenesis of CRC.