scRNA-seq data GSE231559 (Hsu et al., 2023) associated with colorectal cancer were retrieved from the GEO public database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE231559). This dataset included 26 samples, focusing exclusively on colorectal cancer-related tumors and normal control samples. Subsequently, we narrowed down the dataset to a total of 9 samples, consisting of 6 tumor samples and three normal control tissue samples based on the quality and completeness of samples. The detection platform utilized was GPL18573 Illumina NextSeq 500 (Homo sapiens).

Data quality control was conducted applying the “Seurat” (Yu et al., 2022) (version 4.3.0.1, https://cran.r-project.org/web/packages/Seurat/index.html) in R4.1.2. The standard pre-processing workflow included filtering out cells (such as red blood cells, low-quality cells, doublets, mitochondrial, and ribosomal cells) based on quality control (QC) metrics, removing batch effects, normalizing, and scaling the data, and selecting highly variable features. Subsequently, datasets were merged, and integration anchors were defined using the “FindIntegrationAnchors” function (with reduction = ‘rpca’). The integrated objects were then clustered differentially (using the top 40 principal components with resolution = 0.4) and visualized using the UMAP algorithm. Cell types were annotated using “SingleR” (Aran et al., 2019) (version 2.2.0, http://www.bioconductor.org/packages/release/bioc/html/SingleR.html) and canonical cell marker genes from the update database “CellMarker” (Hu et al., 2023) (http://bio-bigdata.hrbmu.edu.cn/CellMarker/). For each cluster, the top gene was selected to generate a distribution diagram illustrating gene expression within each cluster, along with expression plots displaying top marker genes in each cluster.

The FindAllMarkers function in R4.1.2 Seurat package was used to screen differentially expressed genes (DEGs, threshold: min pct = 0.1, |log fold-change|>1 and p value <0.05) in each cluster. Functional enrichment analysis was carried out by using “clusterProfile” version 4.10.0 (Yu et al., 2012) with a threshold of corrected false discovery rate (FDR) value less than 0.05.

Cellchat could infer cell-state-specific signal communication within scRNA-seq profiles through analyzing the expression patterns of ligand-receptor among different clusters (Jin et al., 2021). Here, we utilized “iTALK” (Wang et al., 2019), a computational tool for characterizing and visualizing intercellular communication signals. iTALK categorizes receptor-ligands into four main groups: cytokines, growth factors, immune checkpoints, and others.

Genomic expression profiles associated with colon and rectal cancer were retrieved from the Xena database (https://xenabrowser.net/datapages/). Normalized data was generated from the Illumina HiSeq 2000 RNA Sequencing platform, comprising 438 tumor samples and 41 normal samples, after matching with corresponding clinical information.

Additionally, the GSE39582 (He et al., 2022) dataset from NCBI GEO (https://www.ncbi.nlm.nih.gov/geo/) serving as a validated set, was processed using the GPL570 Affymetrix Human Genome U133 Plus 2.0 Array platform. GSE39582 includes 585 samples, with 519 CRC samples providing prognostic information.

We analyzed the immune cell ratio in each sample based on the whole genomic expression patterns using GSVA (Ye et al., 2019) (version 1.36.3, http://www.bioconductor.org/packages/release/bioc/html/GSVA.html). The Kruskal–Wallis test was employed to assess the immune cell distribution between the tumor group and normal tissue. Subsequently, focusing on fibroblasts related to cancer progression, we integrated the single-cell results.

Using the “limma” package (Ritchie et al., 2015) (version 3.34.7, https://bioconductor.org/packages/release/bioc/html/limma.html), we identified DEGs in the tumor vs normal group from the TCGA profiles, with thresholds set at FDR<0.05 and |log2FC|>1.

For these TCGA-derived DEGs, we adopt WGCNA package (Langfelder and Horvath, 2008) (version 1.61, https://cran.r-project.org/web/packages/WGCNA/index.html) in R4.1.2 to identify modules associated with CAFs. We selected an appropriate soft threshold to ensure that the constructed gene co-expression network exhibited scale-free properties. The cutHeight parameter was set to 0.995, and we employed a clustering algorithm to delineate distinct gene modules. Each module consists of genes that are highly co-expressed, with high connectivity among the genes within the same module. For each identified module, we computed its eigengene, which is the first principal component of the expression data for all genes in that module. Additionally, we calculated the correlation coefficient between the eigengenes of the modules and the CAFs to assess their association. Subsequently, we compared the WGCNA-screened DEGs with previously identified fibroblast-related DEGs from scRNA data. Genes overlapping between these two cohorts were deemed crucial for disease progression. Functional analysis of these overlapping genes was performed using DAVID version 6.8 (Huang da et al., 2009) online tool (https://david.ncifcrf.gov/) to identify GO processes and KEGG pathways.

The screened DEGs underwent univariate analysis using the survival package (version 2.41–1, http://bioconductor.org/packages/survivalr/) (Wang et al., 2016), with a p-value less than 0.05 as threshold. Moreover, the LASSO algorithm regression analysis was employed using the “lars” package (version 1.2, https://cran.r-project.org/web/packages/lars/index.html) (Goeman, 2010) to identify the optimal gene combination. The resulting prognostic coefficients from LASSO were utilized to construct the risk score (RS) model, integrating the target gene expression. The RS formula was defined as follows: Risk score = e sum each   gene ′ s   expression   levels * corresponding   coefficient e sum each   gene ′ s   mean   expression   levels * corresponding   coefficient

RS value was computed for each sample in both the TCGA set and GSE39582 validation set, followed by their classification. The association between risk grouping and actual clinical outcomes was evaluated using Kaplan-Meier curves generated.

Mutation Annotation Format (MAF) files pertaining to CRC tissue were acquired from TCGA, and the “maftools” (Zhang et al., 2019) (version 2.6.05, https://bioconductor.org/packages/release/bioc/html/maftools.html) package in R4.1.2 was utilized to analyze somatic mutation status of signature genes across two risk groups.

Meanwhile, the “estimate” package (Hu D. et al., 2019) (http://127.0.0.1:29606/library/estimate/html/estimateScore.html) was employed to compute the estimate, immune, stromal scores, and tumor purity for all TCGA samples. Subsequently, the Kruskal–Wallis test was conducted to evaluate the immune cell proportion differences among two risk group, and assess the distribution variance of estimate scores. Furthermore, the correlation of model gene expression and immune cell ratio was analyzed.

The CC cell lines (HCT116 and HT29 cells) and NCM460 cells were procured from the Cell Bank of the Chinese Academy of Sciences. The HCT116 and HT29 cells were maintained in an incubator set to 37°C and 5% CO₂, utilizing McCoy’s 5A medium (Servicebio Technology Co., Ltd., China, G4541-500 ML) containing 10% Fetal Bovine Serum (Gibco, United States, 16,000–044) for growth. Additionally, the NCM460 cells were maintained in RPMI-1640 medium (Servicebio Technology Co., Ltd., China, G4532-500 ML).

Total RNA was isolated from the cells using an RNA extraction kit (Servicebio Technology Co., Ltd., China, product number G3013). Subsequently, a cDNA synthesis kit (TransGen Biotech, China, catalog number AU341-02) was employed for the reverse transcription process. The relative mRNA expression levels were assessed using the 2^−ΔΔCq method, normalizing by the expression levels of GAPDH. The CAFs sequences are presented in Table 1.

Statistical analyses were conducted through R software (version 4.1.2) and GraphPad Prism 8.0 software. Multiple groups differences were assessed using the One-way ANOVA. P < 0.05 was considered significant.

The scRNA-seq data were extracted from nine samples, comprising 6 tumor samples and three normal tissue samples. Following data quality control (minGene = 500, maxGene = 6,000, pctRibo = 50, mitochondrial-encoded genes <30%), a total of 21,992 cells were obtained, including 5,358 normal cells and 16,634 tumor cells. All cells underwent classification using dimensionality reduction algorithms, with median count RNA and feature gene in a single cell observed to be 4,462 and 1,372, respectively, as detailed in Table 2.

Cell group classifications were annotated using specific gene markers from the CellMarker database, resulting in the identification of 20 different clusters, including naive T-cells, NK cells, epithelial cells, enterocyte progenitor cells, enteroendocrine cells (EECs), Paneth 1/2 cells, cycling B cells, myeloid cells, Treg cells, LGR5+ stem cells, T follicular helper (Tfh) cells, dendritic cells (DCs), enterocytes, plasma cells, fibroblasts, memory B cells, stromal cells, and plasmacytoid DCs (Figures 1A,B). These results indicated that the cancer cells were highly heterogeneous among patients.

A total of 3,306 DEGs were identified from tumor fibroblasts in comparison to normal. Functional enrichment analysis showed that upregulated cohorts were mainly associated with nuclease activity, replication, and repair pathways, while downregulated genes were predominantly involved in fatty acid metabolic processes and amino acid metabolism pathways (Figures 1C,D).

Intercellular communication among the 20 cluster groups was analyzed based on the expression of receptor-ligand pairs. Cellular interactions were evaluated across four major modules: growth factor, cytokine, immune checkpoint, and others. Numerous cellular communication signals were identified within these cluster groups, serving as crucial factors in intercellular communication (Figure 1E). The top 20 genes related to cellular communication in each cluster were depicted in Figure 1F. Of note, our focus lay on the communicational signal alterations in cancer-related fibroblasts. We observed that fibroblasts interacted with Tfh cells, myeloid cells, and EECs through ligand-receptor pairs of BMPR2-RGMA, BMPR2-BMP8B, and BMPR2-GDF7. Additionally, fibroblasts could interact with B cells involved in disease progression through the ligand-receptor pairs of CAMP-P2RX7 (Figure 1F).

In the TCGA data, we assessed the immune landscape of CRC samples, wherein the immune cell proportion in tumor tissues was evaluated according to gene expression level. A comprehensive screening process identified 23 immune cell types with notable distribution changes in tumor samples compared to the normal group (Figure 2A), such as, fibroblasts, memory B cells, CD56bright NK cells, activated CD4 T cells, NK cells, immature DCs, and others. Notably, fibroblast proportion was significantly increased in tumor tissues, suggesting its prominent association with disease progression (Figure 2B). Leveraging the TCGA expression profile dataset, we utilized the limma package to screen 1,499 DEGs meeting threshold conditions in the tumor vs normal group (Figure 2C).

To pinpoint genes within modules that correlate with CAFs, we analysed 1,499 DEGs using the WGCNA algorithm. To achieve a scale-free network distribution, we explored the parameter value of the adjacency matrix weight, selecting the value of weight power = 8 (Figure 2D). By calculating the dissimilarity coefficient among gene nodes, we obtained a systematic clustering tree (Figure 2E). With a pruning height set at cutHeight = 0.995, we identified seven different modules, different colors represent different modules. We then generated heatmaps to evaluate the correlation between these modules and disease traits, employing the Spearman correlation coefficient (Figure 2F). Heatmap colors indicate the strength of the correlation: red for positive, blue for negative, with deeper hues denoting stronger associations. Focusing on modules with significant associations to CAFs expression within tumors, we identified five modules (blue, brown, green, red, and turquoise) where the correlation values with CAFs exceeded 0.5 (blue module: r = −0.55, P = 3e−117; brown module: r = 0.63, P = 3e−163; green module: r = 0.67, P = 5e−194; red module: r = −0.68, P = 5e−203; and turquoise module: r = −0.52, P = 6e−105). This indicates that 950 genes within these modules are linked to CAFs and the progression of CRC. Notably, while the MEgrey and MEyellow modules exhibited strong correlations with tumor tissue, their correlations with CAFs were comparatively low. This suggests that the gene expression patterns in these modules may not be specific to CAFs but rather more generally associated with tumor tissue characteristics. Consequently, we chose to exclude these modules from further analysis.

Comparison of these modules’ screened DEGs with the previously identified scRNA fibroblast-related DEGs yielded 131 overlapping genes for subsequent analyses (Figure 2G). These candidate genes were associated with 14 GO terms and 7 KEGG pathways, encompassing functions such as B cell proliferation related to immune response, cell adhesion molecule, signal transduction, B cell receptor signal, and peroxisome (Figure 2H).

Univariate analysis results showed 18 genes exhibited prognostic significance values (Figure 3A). Further LASSO analysis revealed that an optimal combination comprising eight DEGs, namely, CD177, RAB36, ACSL6, PBX4, CLDN11, PLIN1, CCDC78, and KCNJ14, was determined (Figure 3B). Subsequently, the RS model was established based on the formula: RS = (−0.02774 ∗ ExpRAB36) + (−0.00025 ∗ ExpCD177) + (−0.00136 ∗ ExpACSL6) + (0.011975 ∗ ExpPBX4) + (0.014539 ∗ ExpCLDN11) + (0.031381 ∗ ExpPLIN1) + (0.039063 ∗ ExpCCDC78) + (0.078007 ∗ ExpKCNJ14).

Samples could be stratified into high- and low-expressed groups according to each gene expression. Kaplan-Meier curve analysis revealed a significant correlation between signature gene expression and patient survival times (Figure 3C).

To assess the predictive capacity of the RS model, samples from TCGA sets and GSE39582 datasets were classified into different risk groups. The association of RS status and survival times of patients were depicted (Figures 3D,E). In the TCGA training set, patients in the low-risk group exhibited a more favorable prognosis compared to individuals in the high-risk group (Figure 3D, left). The AUC values for the ROC curve at 1-, 3-, and 5-years were 0.883, 0.874, and 0.857, demonstrating the high accuracy and efficacy of the risk model (Figure 3D, right). Similarly, the performance of the RS model was verified in the GSE39582 dataset, and AUC values for 1-, 3-, and 5-years were 0.854, 0.803, and 0.805 (Figure 3E, right).

The somatic mutation profiles of tumor samples were obtained from TCGA-CRC cohorts, and the mutational signatures of eight CAFs were visualized among different groups. Our results demonstrated consistent mutation patterns across the eight genes in both groups (Figure 4). Missense mutations were predominant, with single nucleotide polymorphisms occurring more frequently than deletions (Figure 4). Additionally, C>T was identified as the most frequent single nucleotide variant (SNV) in all samples (Figures 4A,B). The mutation profiles of eight CAFs in CRC were further depicted, showing ranked percentages of mutations. Notably, ACSL6 mutations were most prevalent in the low-risk group (23%, Figure 4A), whereas RAB36 mutations accounted for the highest proportion in high-risk group samples (32%, Figure 4B).

The immune cell proportion was evaluated across different risk group samples using the TCGA dataset. Notably, 8 cell types exhibited significantly different ratios, including activated B cells, neutrophils, central memory CD8 T cells, activated CD4 T cells, regulatory T cells, immature dendritic cells, gamma delta-T cells, and eosinophils (P < 0.05, Figure 5A). Additionally, the association between the RS and tumor indicators were explored (Figure 5B). The analysis revealed significant differences among two groups in the tumor purity (P = 0.007812), estimate score (P = 0.036983), stromal score (P = 0.007996), and immune score (P = 0.006668).

Furthermore, there was a significant correlation between the most signature gene expression and tumor indicators (P < 0.05, Figure 5C). Specifically, activated B cell was positively related to PLIN1 and CLDN11 expression, but negatively correlated with ACSL6 expression. Regulatory T cell proportion was positively related to PLIN1, CLDN11, and CD177 expression, while negatively associated with PBX4 and KCNJ14 expression (P < 0.05). The proportion of activated CD4 T cells was negatively correlated with CD177, PLIN1, CLDN11, and KCNJ14 expression (P < 0.05).

Additionally, most signature genes showed correlations with fibroblast proportions in TCGA samples (Figure 5D). Furthermore, the model gene expression was validated in scRNA fibroblasts. Consistent with the previous expression analysis, aberrant expression levels of the eight signature genes were observed in CRC samples (Figure 5E).

qRT-PCR results indicated that the expression level of RAB36 was significantly downregulated in the HCT116 and HT29 cell lines compared to the NCM460 cells. Conversely, CD177, PBX4, and CCDC78 were upregulated in the HCT116 and HT29 cell lines, and ACSL6 and KCNJ14 only in HCT116 cells (P < 0.05) (Figure 6). The expression trend of CD177 and CCDC78 was consistent with the results in Figure 5E. Therefore, these prognostic CAFs may be potential therapeutic targets for CRC.