Abstract
The role of cancer-associated fibroblasts (CAFs) in modulating the tumor microenvironment (TME) is gaining attention, yet their impact on prognosis and therapeutic response in colon cancer remains unclear. Here, we identified genes associated with CAF infiltration via weighted gene co-expression network analysis (WGCNA) utilizing data from The Cancer Genome Atlas (TCGA) and GSE39582 cohorts. Univariate Cox regression and least absolute shrinkage and selection operator (LASSO) regression analyses were used to construct CAF molecular signatures (CAFscore). Patients were categorized into high and low CAFscore groups to analyze clinicopathological traits, somatic mutations, immune evasion, and treatment responses. In this study, a total of 244 genes were correlated with CAF infiltration, with 11 linked to overall survival. Notably, FSTL3, CRIP2, and SLC2A3 were selected for the CAFscore. A higher CAFscore was associated with poorer prognoses, increased malignancy, and therapeutic resistance, particularly among patients with high tumor mutation burden and microsatellite instability. Furthermore, elevated FSTL3 expression was associated with reduced CD8+ T cell infiltration, indicating a suppressive TME. Mechanistically, CAFs may promote immune evasion via NAMPT ligand-receptor interactions based on single-cell RNA sequencing data. Thus, the CAFscore is crucial for personalizing treatment strategies and identifying patients who require more aggressive management.
Similar content being viewed by others
Introduction
Colon cancer is one of the most common malignant gastrointestinal tumors1. However, the heterogeneity among colon cancer patients limits the accuracy of pathological staging in predicting clinical outcomes2,3. The “one gene, one drug” paradigm is no longer applicable to colon cancer, as diagnostic and treatment approaches require more tailored strategies. Consensus molecular subtypes (CMS) have emerged as one of the most robust molecular classification systems for colon cancer4. Among these, the CMS4 subtype (mesenchymal) is characterized by immunologically deserted regions and desmoplastic infiltration, including the activation of transforming growth factor-β (TGF-β), epithelial-mesenchymal transition (EMT) pathways and vasculogenesis, which correlate with the poorest prognosis5. However, the classification of patients using CMS is costly, and current research lacks sufficient evidence linking CMS subtypes to targeted therapies or immunotherapy, thereby limiting their utility in personalized treatment4,5. Therefore, there is an urgent need for a novel molecular signature for colon cancer that can predict prognosis and treatment response, thereby mitigating avoidable costs and complications.
Recent research has shown that targeting TGF-β alone is insufficient to effectively reduce tumor burden or restore immune function, even though it can diminish desmoplastic reactions6,7,8. Consequently, some studies have shifted focus toward targeting the extracellular matrix (ECM) within the tumor microenvironment (TME)8. Cancer-associated fibroblasts (CAFs), crucial stromal components of the TME, are gaining attention for their potential role in evaluating colon cancer outcomes9,10. CAFs significantly influence disease progression, therapy resistance, and the modulation of the immune response9,11,12,13. For instance, CAFs promote the accumulation of tumor-associated macrophages and regulatory T cells while inhibiting the function of natural killer cells12. Despite ongoing efforts to develop predictive models for clinical effectiveness, such as the CAF-derived exosome-related prognostic signature9, complete stratification of prognosis remains elusive. The advent of transcriptomics presents a unique opportunity to reassess the most effective biomarker combinations using a ‘multi-molecular’ approach to predict prognosis and evaluate responses to anticancer therapies.
Immune checkpoint blockade (ICB) is a fundamental treatment option for colon cancer. However, the effectiveness of ICB is limited by challenges such as the lack of accurate screening methods to identify potential beneficiaries and the presence of various immune escape mechanisms14,15,16. Patients harboring POLE/POLD1 mutations exhibit an elevated tumor mutation burden (TMB) and enhanced responsiveness to ICB, particularly within the microsatellite stability (MSS) subgroup17. This finding expands the eligibility for ICB treatment and challenges the notion that microsatellite instability (MSI) patients with colon cancer respond more favorably to ICB therapy18,19. Consequently, advancements in the molecular classification of colon cancer have optimized existing ICB biomarkers. This study aimed to establish a comprehensive molecular signature based on CAFs to identify high-risk individuals, predict treatment outcomes, and explore the underlying mechanisms of immune evasion using bulk RNA sequencing (RNA-seq) and single-cell RNA-seq (scRNA-seq).
Weighted gene co-expression network analysis (WGCNA) is a systematic bioinformatics algorithm that integrates highly co-expressed genes into the same gene module20 by studying the biological relationships between the genes in the module and clinical traits. In this study, WGCNA was concurrently utilized in two transcriptome cohorts comprising colon cancer samples from The Cancer Genome Atlas (TCGA) and the Gene Expression Omnibus (GEO), allowing for the identification of hub modules most commonly associated with CAF infiltration. Our findings provide a systematic investigation of the clinical significance of a CAF-related molecular signature (CAFscore) in predicting colon cancer prognosis and identifying beneficiaries of treatment.
Materials and methods
Patients and data collection
Patients from TCGA and GEO were screened as the discovery and validation cohorts, respectively. The TCGA GDC portal (https://www.tcga.org) was used to download bulk RNA-seq data for the TCGA-COAD samples (n = 446) with complete clinical data. Fragments per kilobase per million (FPKM) were converted to transcripts per million (TPM) for training purposes. Validation data were retrieved from GSE39582 (n = 562) cohort via GEO (https://www.ncbi.nlm.nih.gov/geo/). To correct for batch effects, the “ComBat” function of the “sva” R package (version 3.50.0) was employed21, following the transformation of all bulk RNA-seq data through log2 (TPM + 1) normalization. The enrollment criteria were as follows: (1) availability of RNA-seq data for CAFscore calculation; (2) clinical information, including overall survival (OS), age, gender, and stage; (3) an explicit pathological diagnosis of colon cancer; and (4) no prior history of anticancer medications.
A colon cancer cohort from Zhongda Hospital Southeast University (ZHSU) was used to validate the roles of specific genes. Tissue samples embedded in wax blocks were collected retrospectively from colon cancer patients (n = 157) who underwent radical surgery between January 1, 2014, and August 31, 2021. The last follow-up date was September 1, 2023. The enrollment criteria for these patients were as follows: (1) pathological diagnosis of colon cancer, (2) completely preserved tissue wax blocks, and (3) no history of anticancer therapy. The exclusion criteria included: (1) patients with severe cardiopulmonary insufficiency, severe infections, blood disorders, autoimmune diseases, hepatitis, or other tumors, (2) cases missing essential clinical or follow-up data. It contained 83 pairs of tumors and paracancerous tissues. The American Joint Committee on Cancer Staging Manual, 8th edition, was utilized for patient staging. This study involved a retrospective analysis of anonymized patient data, and was conducted in accordance with the Declaration of Helsinki. And the study has been approved by the Clinical Research Ethics Committee of Zhongda Hospital, Affiliated to Southeast University. Clinical characteristics of above three cohorts are detailed in Table 1.
Additionally, scRNA-seq data for primary colon cancer samples (n = 8) were retrieved from the GSE110009 cohort22. Furthermore, scRNA-seq data from tissue samples obtained prior to therapy from colorectal cancer patients (n = 20) who underwent ICB were included from the GSE236581 cohort23. The IMvigor210 cohort (n = 298), which comprised urothelial carcinoma patients treated with anti-PD-L1, was also accessed6 (http://research-pub.gene.com/IMvigor210CoreBiologies).
Estimation of fibroblast infiltration in the TME
To evaluate the infiltration of fibroblasts in the TME, several computational methods were employed, including the Estimate of Proportions of Immune and Cancer Cells (EPIC)24, Microenvironment Cell Populations-counter (MCP-counter)25, xCell26, and ESTIMATE27. These methods were implemented using the respective “EPIC”, “MCPcounter”, “xCell”, and “estimate” R packages, following protocols established in the literature28,29. Additionally, based on established transcriptomic biomarkers, the Tumor Immune Dysfunction and Exclusion (TIDE) tool30,31 (http://tide.dfci.harvard.edu/) was utilized to classify individuals as responders or non-responders to ICB therapy, and to assess CAF infiltration via Pearson correlation analysis between expression profile and FAP+ CAF signature. The results pertaining to CAF infiltration and Stromalscore were extracted from these five infiltration methods based on the mRNA expression profile of TCGA-COAD and GSE39582 cohorts. The results of five infiltration methods are presented in Table S1-2. To analysis the relationship between infiltration fraction of immune cells and expression levels of marker genes, infiltration score was calculated through RNA-seq data based on the xCell algorithm26.
WGCNA
The identification of CAF hub genes strongly associated with CAF infiltration was performed using the “WGCNA” R package (version 1.72-5)32. The analysis began with the calculation of Pearson’s correlation coefficients to assess the direct correlation between any two genes. The soft-thresholding power (β) was then determined based on the scale-free topology criterion. An adjacency matrix was generated and clustered using a topological overlap measure, enabling the calculation of both direct and indirect correlations between module genes. To identify the hub gene module with the strongest correlation to CAF infiltration and Stromalscore, epigene significance was analyzed using the “corPvalueStudent” function. Subsequently, a scatter plot was created to illustrate the relationship between module membership (MM) and gene significance (GS). Module genes with MM > 0.6 and GS > 0.5 were defined as CAF hub genes. Ultimately, 244 CAF hub genes were identified as overlapping between the TCGA-COAD and GSE39582 cohorts.
Construction and validation of the CAFscore
The identification of 244 hub genes was achieved through univariate Cox regression analysis, selecting genes with p < 0.05, and hazard ratios (HR) > 1. This analysis highlighted 11 OS-related CAF hub genes within the TCGA-COAD cohort. These 11 genes were then subjected to least absolute shrinkage and selection operator (LASSO) regression analysis, executed with 1,000 iterations using the “glmnet” R package (version 4.1-8)33. A 10-fold cross-validation was performed to determine the optimal value of the regularization parameter λ, ensuring robust model generalization. From the LASSO regression, three genes emerged as significantly associated with patient survival, each assigned non-zero regression coefficients. These coefficients indicate the substantial contribution of genes to prognosis prediction. Thus, the CAFscore was established as follows: CAFscore = sum (mRNA expression of each gene × corresponding coefficient). The optimal cut-off value for categorizing patients into low and high CAFscore subgroups was determined using the “Survminer” R package (version 0.4.9).
Construction and evaluation of a Nomogram
The nomogram was constructed utilizing the “rms” R package (version 6.7-1). To assess the prognostic utility of the nomogram, we employed several statistical methods, including Kaplan-Meier analysis, univariate and multivariate Cox regression, as well as time-dependent receiver operating characteristic curves (tROC). Harrell’s concordance index (C-index) and decision curve analysis (DCA) curves were calculated using the R packages “survcomp” (version 1.52.0) and “ggDCA” (version 1.1), respectively.
Functional enrichment and gene set enrichment analysis
Gene set enrichment analysis (GSEA)34 was conducted using reference files obtained from the Molecular Signatures Database (MsigDB) (https://www.gsea-msigdb.org/gsea/msigdb). Enrichment analyses focused on Gene Ontology (GO, c5.go.Hs.symbols) and Kyoto Encyclopedia of Genes and Genomes (KEGG, c2.cp.kegg.v7.3.symbols) to elucidate biological functions, which included biological processes (BP), cellular components (CC), and molecular functions (MF), along with relevant pathways. These analyses were performed using the “clusterProfiler” R package (version 3.14.3)35.
Somatic alteration data collection and analysis
Nucleotide variation data in the MAF format were obtained from VarScan2 through the TCGA GDC portal. The “maftools” R package (version 2.20.0) was utilized to analyze the frequency and types of DNA mutations, including POLE/POLD1. TMB for each sample was calculated by determining the number of mutant bases per million bases. The genomic data were then visualized in circular genomic coordinates using the “Rcircos” R package (version 1.2.2). Additionally, MSI status clinical data were acquired from The Cancer Immunome Atlas (TCIA), which provides insights from comprehensive immunogenomic analyses of next-generation sequencing data (https://tcia.at/)36.
Copy number data collection and analysis
Copy number variation (CNV) data were downloaded from the Firehose Broad GDAC (https://gdac.broadinstitute.org/) and preprocessed using the GISTIC 2.0 web tool in Genepattern (https://cloud.genepattern.org/). The results were visualized using the “ggplot2” R package (version 3.5.1).
Drug susceptibility predictions
The IMvigor210 cohort was utilized to validate the predictive capabilities of the CAFscore in relation to ICB responses. The immunophenscore (IPS) was downloaded from TCIA36. Subsequently, data on first-line medications and potentially effective small-molecule drugs for colon cancer were retrieved from the Cancer Therapeutics Response Portal (CTRP2.0) and PRISM databases. A ridge regression model was constructed to evaluate drug sensitivity across low and high CAFscore subgroups. Additionally, the Genomics of Drug Sensitivity in Cancer (GDSC) database (https://www.cancerrxgene.org/) and the “pRRophetic” R package (https://github.com/paulgeeleher/pRRophetic) were employed to identify potentially efficacious small-molecule drugs. The half-maximal inhibitory concentration (IC50) was calculated using ridge regression to assess therapeutic resistance, and ten-fold cross-validation was applied to evaluate the accuracy of the model.
Immunohistochemistry (IHC) assay
For the IHC assay, an anti-FSTL3 antibody (Abcepta, #10272, 1:100, Suzhou, China) was used to analyze 83 paired samples following a two-step protocol37. The IHC score was calculated by multiplying the categorical grades for the number of positively stained cells and their staining intensity38. The quantification process involved assessing the number of positively stained cells out of 100 cells in each of five selected high-power fields (400x magnification) to determine the percentage of positively stained cells per section. The percentages were categorized as 0% (grade 0), 10–25% (grade 1), 26–50% (grade 2), 51–75% (grade 3), and 76–100% (grade 4). Staining intensity was further classified into grades 0, 1, 2, and 3, corresponding to no stain, light brown, brown, and dark brown, respectively. The IHC score was reviewed by two pathology experts. High expression of FSTL3 was defined as an IHC score ≥ 1. In addition, the Human Protein Atlas (HPA) (https://www.proteinatlas.org/)39 was utilized to assess protein expression in colon cancer tissues, with basic annotation parameters that categorized staining levels as undetectable, low, medium, or high.
Immunofluorescence (IF) staining
IF staining was performed on paraffin-embedded samples using a two-plex staining protocol. The order of staining was as follows: anti-CD8 (Abcam, #ab1714, 1:50, USA) stained with Alexa Fluor 647 (Beyotime, #A0473, 1:500, China), anti-FSTL3 (Abcepta, #10272, 1:100, Suzhou, China) stained with CoraLite488 (Proteintech, #SA00013-2, 1:100, China). All primary antibody solutions were incubated overnight at 37 °C, followed by DAPI staining of the nuclei prior to mounting. The stained slides were then scanned, and images were analyzed using Virtual Microscope Software Slide Viewer 2.5 3D Histech®.
Analysis of scRNA-seq data
Tumor immune single-cell hub 2 (TISCH2) was used to explore the expression levels of model genes across various cell types within the EMTAB8107, GSE146771, GSE166555, and GSE179784 cohorts (http://tisch.comp-genomics.org/). To analyze the GSE236581 dataset, the “AverageExpression” function from the “Seurat” R package (version 4.3.0) was utilized to calculate the average gene expression for each sample. Subsequently, the expression of model genes was extracted to compute the CAFscore. Following the Response Evaluation Criteria in Solid Tumors (RECIST) version 1.1, patients with a partial response or complete response were classified as responders, while all other outcomes were designated as non-responders.
In addition, the “Seurat” R package was primarily utilized to analyze the scRNA-seq data from GSE110009 cohort based on a standardized procedure (https://github.com/satijalab/seurat)40. Genes expressed in fewer than three cells, cells with fewer than 50 detected genes, and cells with mitochondrial percentages exceeding 5% were excluded from the analysis. The count matrix was normalized with default settings, and top 1,500 variable genes from each dataset were scaled for principal components analysis (PCA). The first 20 principal components were computed, and cell clusters were identified at a resolution of 0.5. To annotate cell types based on marker genes, the “SingleR” R package (version 1.0.0) was employed41. Furthermore, cell-cell communication was investigated using the “CellChat” R package (version 2.1.2).
Cellular level validation
For cellular level validation, TPM-normalized RNA-seq transcriptomic data of the genes in colorectal cancer cells (n = 56) and CAFs (n = 2) were downloaded from the Cancer Cell Line Encyclopedia (CCLE) database (https://portals.broadinstitute.org/ccle).
Statistical analysis
Statistical analyses and data visualization was performed using R software version 4.1.2 (https://www.r-project.org/). Differences between two subgroups were analyzed using Wilcoxon test, while paired sample t-test was used to compare the differences of IHC score in FSTL3 and IF staining intensity in FSTL3 and CD8 between 83 pairs of tumors and paracancerous tissues. The Chi-squared test was employed to analyze the differences in ICB response between the high and low CAFscore groups. Spearman’s correlation analysis was utilized to determine the relationship between two indicators. The diagnostic utility of the model was evaluated using ROC analysis to calculate the area under the curve (AUC). The “DeLong” method of MedCalc (version 23.0.1) was used to conduct analyses of the ROC curves. Survival time was compared using the Kaplan-Meier method and the log-rank test. HR and 95% confidence intervals (CI) for the prognostic factors were calculated using Cox proportional hazard models. PCA was performed using the “stats” R package (version 4.1.3). A two-sided p < 0.05 was considered statistically significant.
Results
Co-expression network of gene modules and CAF infiltration
To identify hub genes associated with CAF infiltration, we performed correlation analyses between CAF traits and gene modules exhibiting similar expression patterns. In both the TCGA-COAD and GSE39582 cohorts, OS was consistently found to be significantly shorter in the high CAF infiltrating group, as determined by CAF_EPIC, CAF_MCPcount, and Stromalscore, compared to the low CAF infiltration group (p < 0.05, Fig. 1A,B). Therefore, these three infiltration scores were identified as hub traits related to CAFs.
CAF Hub Genes were identified through the analysis of co-expression networks of CAF infiltration using WGCNA. (A,B) Kaplan–Meier analysis showed that OS was generally shorter in the high CAF infiltration group, as determined by algorithms including CAF_EPIC, CAF_MCPcount, CAF_TIDE, CAF_xCell, and Stromalscore, compared to the low infiltration group in patients from the TCGA-COAD (A) and GSE39582 (B) cohorts. (C,D) Clustering dendrograms indicated that genes with similar expression patterns were grouped into co-expression modules based on a height value below 0.2. (E,F) The heat map illustrated the correlations between each gene module and CAF-related phenotypes. Red represented a positive correlation, while blue represented a negative correlation; the deeper the color, the stronger the correlation. The value inside the parentheses represented the p-value of the correlation analysis.
For the network analysis, we established a soft-thresholding power (β) for the adjacency function, where the square of the model R approached 1. This value was determined to be 3 for the TCGA-COAD cohort and 4 for the GSE39582 cohort (Fig. S1A,B). Following this, a distributed network was constructed to illustrate the relationships between gene modules (Fig. 1C,D). Notably, the relationship between the genes in the turquoise module and the three aforementioned infiltration scores exhibited the highest significance (Fig. 1E,F).
Construction of a molecular signature
Hub genes were identified based on MM values exceeding 0.6 and GS values exceeding 0.5 in the turquoise module of the TCGA-COAD and GSE39582 cohorts (Fig. 2A,B and Table S3). A total of 244 intersectional CAF hub genes were found to co-occur between the turquoise modules of both cohorts (Fig. 2C). GO and KEGG enrichment analyses42,43,44 were conducted on these 244 genes, revealing significant enrichment in biological processes such as extracellular matrix (ECM) organization, cell-substrate adhesion, and the TGF-β signaling pathway (Fig. 2D,E). Subsequently, we performed univariate Cox regression analysis in the TCGA-COAD cohort, identifying 11 OS-related CAF hub genes (HR > 1, p < 0.05, Fig. 2F). Loss of CNV was much more frequent in SFRP2, FSTL3, ZNF532, ARL4C, CRIP2, KCNE4, and KCNJ8, whereas CNV gain was more prevalent in SLC2A3, SERPINE1, S1PR3, and BGN (Fig. 2G). Then, the LASSO regression model was exploited to extract three genes, and the CAFscore was calculated in both the TCGA-COAD and GSE39582 cohorts using the following formula: CAFscore = (the mRNA expression of CRIP2) * 0.0189 + (the mRNA expression of FSTL3) * 0.1223 + (the mRNA expression of SLC2A3) * 0.1043 (Fig. 2H,I and Tables S4 and S5).
Functional enrichment analysis of CAF hub genes and construction the CAFscore. (A,B) Scatter plots were used to visualize the Module Membership (MM) and Gene Significance (GS) of genes within the turquoise module. (C) Venn diagram showing 244 intersectional CAF hub genes. (D,E) The bubble diagram showed the GO enriched analyzes of the biological process (BP), cellular component (CC), and molecular function (MF) terms (D), and KEGG pathway analysis (E) of these 244 genes. (F) Univariate Cox analysis was conducted to identify 11 OS-related genes in the TCGA-COAD. (G) CNV of the above 11 genes. (H) Coefficient profiles of LASSO regression analysis. (I) The adjustment parameter (lambda) was determined by calculating the partial likelihood of deviance using ten-fold cross-validation.
High CAFscore correlated with poor prognosis and served as an effective predictor
In the TCGA-COAD and GSE39582 cohorts, a high CAFscore was found to be positively correlated with several established infiltration of CAFs, including CAF_EPIC, CAF_MCPcount, CAF_TIDE, CAF_xCell, and Stromalscore (Fig. 3A and Fig. S2A). As the CAFscore increased, the expression levels of three model genes along with CAF-related marker genes also showed a significant upregulation (Fig. 3B and Fig. S2B). The CAFscore and mRNA expression of the model genes were positively correlated with the majority of CAF-related marker genes (Fig. 3C and Fig. S2C). These findings indicate that the CAFscore is closely associated with existing CAF-related markers, suggesting its potential importance in relation to clinical characteristics.
Clinical characteristics analysis based on CAFscore and nomogram construction in the TCGA-COAD. (A) CAFscore was positively correlated with CAF_EPIC, CAF_MCPcount, CAF_TIDE, CAF_xCell, and Stromalscore. The numbers inside the boxes represented the correlation coefficients in the analysis; the deeper the yellow color, the higher the correlation. (B) Heatmap showed that the expression of 3 model genes and CAF-related marker genes. (C) The chart displayed the correlation between CAFscore and mRNA expression of CAF-related marker genes. Red meant positive correlation. The larger the correlation coefficient, the darker the color. (D) The PCA diagram demonstrated the ability of the CAFscore to divide samples into two subgroups. (E) The riskplot displayed poor survival status in the high CAFscore subgroup. (F) The Kaplan–Meier analysis evaluated the prognosis. (G) Distribution of CAFscore in stage, classification of T, N, and M, status of survival, progression, and perineural invasion. Distribution of CAFscore in the TCGA-COAD stage II/III samples (n = 213) with tumor ___location data. (H) A nomogram based on the CAFscore and clinical indicators provided a reference for 1-, 3-, and 5-year survival predictions. (I) The calibration curve displayed the fitting of the prediction models to the actual events. (J) Time-dependent ROC curves were used to evaluate the predictive ability of the nomogram for 1-year, 3-year, and 5-year OS. (K) The C-index of the nomogram was totally over 0.65 during the ten years. (L) DCA curve showed the accuracy of the nomogram and the existing clinical indicators. (M) CAFscore was served as an independent prognostic factor. ***, p < 0.001; **, p < 0.01; *, p < 0.05.
PCA revealed that the CAFscore was helpful to segregate samples into two distinct subgroups (Fig. 3D and Fig. S2D), with the high CAFscore subgroup associated with unfavorable outcomes (Fig. 3E,F and Fig. S2E,F). Notably, the CAFscore was significantly higher in the dead and progression groups compared to the alive and no-progression groups (Fig. 3G and Fig. S2G). Additionally, the CAFscore showed a gradual increase with advancing tumor stage, T classification, and N classification (Fig. 3G and Fig. S2G). Although there was no significant difference in CAFscore between the M0 and M1 groups, the M1 group exhibited a higher score.
To further investigate the role of CAFscore in assessing recurrence risk in specific stages. The results showed that colon cancer patients with high CAFscore in stage II/III had a higher risk of disease progression or recurrence (HR > 1, p < 0.05 both in the TCGA-COAD and GSE39582 cohorts, Table S6). Interestingly, the CAFscore in the right-sided colon (RSC) and proximal samples of stage II/III tumors were significantly elevated compared to those in the left-sided colon and distal samples (Fig. 3G and Fig. S2G). This finding aligned with the commonly observed poor prognosis and treatment resistance associated with RSC tumors, further confirmed by other malignant clinicopathological features. The CAFscore for perineural invasion and positive Cimp_status generally increased (Fig. 3G and Fig. S2G). Moreover, in the GSE39582 cohort, the CAFscore was significantly elevated in the C4 (cancer stem cell) group, correlating with poor prognosis and showing a strong positive correlation with cancer stem cell-related genes (Fig. S2H–J). Consequently, the CAFscore serves as a valuable indicator of the aggressive biological characteristics of tumors.
To assess whether the CAFscore enhances prognostic accuracy compared to existing indicators, a nomogram was constructed incorporating the CAFscore along with common clinical indicators. This nomogram facilitated a more precise determination of patient prognosis risk and provided reference points for 1-, 3-, and 5-year survival predictions (Fig. 3H; Fig. S3A). Calibration curve indicated that the nomogram accurately predicted true prognosis (Fig. 3I and Fig. S3B), and tROC curves were employed to evaluate its predictive ability, achieving AUC values of 0.758, 0.837, and 0.803 at 1, 3, and 5 years, respectively, in the TCGA-COAD cohort (Fig. 3J and Fig. S3C). The C-index for the nomogram consistently exceeded 0.65 over a ten-year span (Fig. 3K and Fig. S3D). Importantly, DCA curve showed that the nomogram outperformed existing clinical indicators (Fig. 3L and Fig. S3E).
Furthermore, multivariate regression analysis confirmed that the CAFscore is an independent prognostic factor (Fig. 3M and Fig. S3F). Overall, a high CAFscore correlates with poor outcomes and serves as a practical classifier for colon cancer patients at high risk, validated across discovery and validation cohorts.
Correlation of CAFscore with ICB-related biomarkers and ICB efficacy
The CAFscore plays a crucial role in identifying high-risk colon cancer patients, and its assessment extends to evaluating treatment efficacy, which holds significant clinical value. According to TIDE tool, the proportion of responders to ICB was significantly lower in the high CAFscore subgroup compared to the low CAFscore subgroup (p < 0.001 for both the TCGA-COAD and GSE39582 cohorts, Fig. 4A). In the GSE236581 real-world cohort, patients were stratified based on the median CAFscore, revealing a striking 100% response rate to ICB in the low CAFscore subgroup (Fig. 4A and Table S7). ROC curves were utilized to evaluate the predictive power of the CAFscore for ICB responses. The AUC values for CAFscore exceeded those of ICB biomarkers, with values of 0.769, 0.798, and 0.909 observed in the TCGA-COAD, GSE39582, and GSE236581 cohorts, respectively (Fig. 4B). Additionally, the TIDE score of the high CAFscore subgroup was significantly higher than that of the low subgroup (Fig. 4C). Conversely, the IPS, which reflected the ICB response, was markedly lower in the high CAFscore subgroup compared to the low subgroup (Fig. 4D).
Relationship of CAFscore with ICB-related biomarkers and therapeutic efficacy. (A) The bar plot showed the proportion of responder within the high and low CAFscore subgroups. (B) ROC curves of the CAFscore in predicting ICB responses. (C,D) The violin diagram displayed the TIDE score and IPS in the high and low CAFscore subgroups. (E,F) The bar plot displayed the CAFscore in different MS and MMR states. (G,I) MSI and TMB values were shown in two CAFscore subgroups. (H,J) Comparison of correlation analysis of MSI and TMB with different CAF infiltration scores. The red font on behalf of p < 0.05. (K) TMB, MSI, and CAFscore values were compared between POLE/POLD1 mutated and wild-type patients. (L–N) Survival curve analysis revealed different prognoses among the subgroups in TCGA-COAD cohort. (O) The Kaplan-Meier curve showed that the high CAFscore subgroup had poor outcomes in the IMvigor210 cohort. MS mismatch state, MMR mismatch repair.
Furthermore, patients with microsatellite instability-high (MSI-H) and deficient mismatch repair (dMMR) exhibited higher CAFscore than those with microsatellite instability-low (MSI-L), microsatellite stable (MSS), and proficient mismatch repair (pMMR) (Fig. 4E,F). The high CAFscore subgroup also exhibited significantly elevated levels of MSI and TMB compared to the low CAFscore subgroup (Fig. 4G,I). A positive correlation was observed between the CAFscore and both MSI and TMB (Spearman’s correlation analysis, Fig. 4H,J).
In addition, TMB, MSI, and CAFscore were significantly elevated in the POLE/POLD1 mutation group (Fig. 4K). The CAFscore further distinguished the prognosis of wild-type status patients with POLE/POLD1 mutations indicating a poorer prognosis in the high CAFscore group (Fig. 4L). In a clinical subgroup analysis of stage II/III colon cancer within the TCGA-COAD cohort, patients in the high CAFscore subgroup of the MSI-H or TMB-H groups exhibited worse prognoses compared to those in the low subgroup, based on optimal cutoff values (MSI = 0.34, CAFscore = 1.44, Group 3 vs. Group 4, p = 0.016; TMB = 2.42, CAFscore = 1.44, Group 3 vs. Group 4, p = 0.044; Fig. 4M,N). This suggests that additional stratification based on CAFscore patterns may serve as a more effective prognostic indicator. Furthermore, the high CAFscore indicated potential drug resistance and poor prognosis, as the high CAFscore subgroup exhibited significantly shorter survival times compared to the low CAFscore subgroup in the IMvigor210 cohort (Fig. 4O). Therefore, the CAFscore is expected to become a predictive marker for resistance to ICB, facilitating the development of personalized treatment strategies.
CAFscore indicated immune evasion based on the NAMPT signaling pathway
Activation of the TGF-β signaling pathway and vasculogenesis are typical features of the CMS4 subtype of colon cancer. A positive correlation was observed between the CAFscore and both the TGF-β signaling pathway and vasculogenesis (Fig. 5A). TGF-β promotes tumor progression by creating a suppressive TME7,45,46 and reducing the efficacy of ICB6; however, the specific mechanism of immune evasion in colon cancer remains unclear. This study found that the mRNA expression levels of numerous immunosuppressive genes in the high CAFscore group were higher than those in the low CAFscore group (Fig. 5B). Additionally, in the TCGA-COAD cohort, GO analysis revealed significant enrichment of immune receptor activity in the high CAFscore subgroup (p < 0.001, Fig. 5C and Table S8).
The enrichment analysis of the CAFscore explored the mechanisms of immune evasion in colon cancer. (A) The correlation analysis of the TGF-β signaling pathway and vasculogenesis with CAFscore in the TCGA-COAD cohorts. (B) The bar plot showed the mRNA expression levels of nine immunosuppressive genes in high and low CAFscore groups of TCGA-COAD cohorts. (C) GO enrichment analysis was conducted for the high CAFscore subgroup. (D) Ligand-receptor pair interactions were analyzed using sc-RNA seq data from GSE110009. The x-axis denoted the interacting cell types, while the y-axis represented the ligand-receptor pairs. A higher intensity in red in the circle indicated a higher probability of communication. (E) The ligand-receptor networks of 7 types of cells, with line thickness representing both the number (top) and strength (bottom) of interactions. (F–G) Heatmap and violin plot illustrated the interactions within the NAMPT signaling network, focusing on the expression levels of ligand receptors. (H) The networks delineated the direction of interaction of the NAMPT pathway network. (I) Correlation analysis between CAFscore and NAMPT mRNA expression in TCGA-COAD. (J) The bar plot showed differential expression of NAMPT between the two groups.
To further elucidate the underlying mechanisms, sc-RNA seq was employed to analyze the ligand-receptor interaction networks associated with fibroblasts. Among these, Nicotinamide phosphoribosyl transferase (NAMPT)-related pathways, such as NAMPT-INSR and NAMPT-(ITGA5 + ITGB1), exhibited the highest communication probability (Fig. 5D). The intricate ligand-receptor networks displayed the highest number and strength of interactions among fibroblasts, monocytes (including macrophages), and endothelial cells (Fig. 5E). Specifically, the communication probabilities of ligand receptors associated with the NAMPT pathway were the highest in fibroblasts and monocytes (Fig. 5F). NAMPT was found to be highly expressed in monocytes, while its receptors, INSR, ITGA5, and ITBG1, were predominantly expressed in fibroblasts and endothelial cells (Fig. 5G). The ligand-receptor pairs of NAMPT-INSR and NAMPT- (ITGA5 + ITBG1) were of particular interest, with monocytes identified as the primary cell source of NAMPT (Fig. 5H). In addition, NAMPT was significantly positively correlated with the CAFscore, showing higher expression levels in the high CAFscore group (Fig. 5I,J). M2 polarization of tumor-associated macrophages in colon cancer, driven by NAMPT, results in the creation of an immunosuppressive TME47. Thus, the CAFscore has the potential to serve as a quantifiable marker for immune evasion in colon cancer, highlighting its relevance in understanding tumor biology and informing therapeutic strategies. The findings clarify the status of high CAFscore patients with an immunosuppressive TME.
Genomic variations related to CAFscore informed the selection of targeted therapy
Tumors sharing a pathological diagnosis can exhibit diverse molecular characteristics48. This study explored the relationship between the CAFscore and genomic mutation, given the potential of genomic variations to influence patient prognosis and serve as druggable targets. Analysis of CNV topography revealed that the TCGA-COAD samples displayed the highest gains in chromosomes 8 and 20 compared to the reference samples. Notably, the incidence of amplification or deletion was generally higher in the high CAFscore group than in the low CAFscore group (Fig. 6A).
Genomic mutation profile in different CAFscore groups. (A) The topography of CNV in the high and low CAFscore subgroups, including gistic score and frequency. The violin plots showed the amplification/deletion mutations in both the high and low CAFscore groups. (B) Waterfall plots were used to intuitively display the top 30 genes with the highest mutational frequencies in the low and high CAFscore subgroups. (C) 30 mutated genes were mapped to circular genome coordinates. (D) CNV of the top 30 mutational genes. (E) Mutations in the low and high CAFscore subgroups were counted, and the bar plot displayed the OR values in the TCGA-COAD. (F) The box plot showed that CAFscore in the BRAF mutant and wild-type group. (G,H) Survival curve analysis revealed different prognoses among subgroups in the TCGA-COAD. (I) Survival curve analysis revealed the prognoses of co-altered and non-altered patients in TCGA-COAD cohort. *, p < 0.05; **, p < 0.01; ***, p < 0.001.
The first 30 genes with the highest mutation frequencies were subsequently mapped to the circular genome coordinates (Fig. 6B). This analysis indicated that 13 genes demonstrated a higher frequency of loss events than of gain events, whereas 16 genes exhibited a higher frequency of gains than losses (Fig. 6C). Specifically, in chromosome 8, the CNV of CSMD3 and ZFHX4 showed a higher frequency of gain, and the CNV of CSMD1 displayed a higher frequency of loss (Fig. 6C-D). Further quantification of mutation proportions within the high and low CAFscore subgroups revealed that BRAF exhibited the highest odds ratio (OR) among the top 30 high-variability genes in TCGA-COAD (Fig. 6E). Importantly, the CAFscore was significantly higher in the BRAF-mutant group than in the wild-type group (Fig. 6F). However, the BRAF mutation status alone may not be sufficient to determine patient prognosis (Fig. 6G). When combined with the CAFscore, the prognostic stratification of patients with wild-type BRAF improved (Fig. 6H). Thus, the CAFscore serves as a valuable reference for prognostic assessment following genetic testing.
Recent studies indicate that patients with advanced colon cancer harboring TP53/RAS/BRAF co-mutations according to genetic testing exhibit a significantly poorer prognosis than those with either no mutations or a single mutation49. Within the stage III-IV colon cancer subgroup in TCGA-COAD cohort, individuals with TP53/RAS/BRAF co-mutations demonstrated reduced OS compared with other patients (Fig. 6I). In terms of potential therapeutic targets, a variety of targeted therapies predicted to be effective in colon cancer patients were identified using the GDSC database. A low CAFscore was associated with lower IC50 value and greater sensitivity to targeted therapeutics, including afatinib, erlotinib, gefitinib, lapatinib, and KRAS (G12C) inhibitors (Fig. 7A,B). In addition, the same conclusion was reached for standard targeted therapies, including regorafenib (Fig. 7C). Regarding chemotherapy agents, the low CAFscore subgroup exhibited increased sensitivity to fluorouracil based on analyses from the CTRP databases (Fig. 7D). Other first-line chemotherapeutic agents for colon cancer, including capecitabine and oxaliplatin, yielded consistent conclusions with previously reported findings (Fig. 7E,F). Hence, the subgroup with a low CAFscore demonstrates greater responsiveness to treatment, suggesting that CAFscore may serve as a potential tool for identifying individuals who are particularly responsive to pharmacological interventions in colon cancer.
Sensitivity of targeted therapy and chemotherapeutic agents in different CAFscore subgroups. (A,B) The bar chart showed the estimated IC50 values of targeted therapeutics within the TCGA-COAD (A) and GSE39582 (B) cohorts. (C,D) The bar chart displayed the estimated AUC values of regorafenib and fluorouracil in the CAFscore subgroups of the TCGA-COAD cohort using the CTRP database. (E,F) Sensitivity to the capecitabine and oxaliplatin was analyzed in the TCGA-COAD cohort using the Primer database.
Multidimensional validation of model genes
This study examined model gene expression in various cell subpopulations and their tissue distribution in colon cancer. After quality control, 2,496 cells from the GSE110009 cohort were selected (Fig. S4A–F). SLC2A3 was broadly expressed across subpopulations, CRIP2 in endothelial cells and fibroblasts, and FSTL3 in fibroblasts and epithelial cells (Fig. 8A). Subsequent sc-RNA seq analysis of colon cancer samples from the EMTAB8107, GSE146771, and GSE179784 cohorts, accessed through the TISCH online database, confirmed that CRIP2 and FSTL3 were predominantly expressed in stromal cells (Fig. 8B). And analysis of the CCLE database confirmed that FSTL3 expression was significantly higher in CAFs compared to colorectal cancer cells (p < 0.05, Fig. 8C and Table S9). According to data from HPA, IHC analysis revealed high expression of FSTL3 and SLC2A3 in colon cancer tissues, with FSTL3 displaying particularly strong staining (Fig. 8D). These findings suggest that FSTL3 may serve as a biomarker for CAFs in colon cancer.
Multidimensional validation of model genes. (A) Violin plot showed that the expression of three genes in the different cell types of GSE110009. (B) The heatmap displayed the expression of CRIP2 and FSTL3 in four different scRNA-seq cohorts from Tumor immune single-cell hub 2 (TISCH2) database. (C) The mRNA expression of FSTL3 in colorectal cancer cells (n = 56) and CAFs (n = 2) based on the CCLE database. (D) Protein expression of these three genes in colon cancer from HPA database. (E,F) IHC images and boxplots of matched colon cancer and paracancerous tissue samples were illustrated (n = 83). (G) ROC curves for the diagnosis of FSTL3 in the ZHSU, TCGA, and GSE39582 cohorts. (H) The expression of CD8 and FSTL3 were showed in the representative IF staining images of matched colon cancer and paracancerous tissue samples. (I) The heatmap visualized the correlation between FSTL3 expression levels and the infiltration of immune cells, as determined by the xCell algorithm. The numbers within the boxes represented the correlation coefficients. A deeper red color indicated a stronger positive correlation, while a deeper blue color reflected a stronger negative correlation. (J) Survival curves of the three cohorts based on FSTL3 expression.
To evaluate the impact of FSTL3 on colon cancer prognosis, we validated 83 paired samples from the local ZHSU cohort. The IHC scores for FSTL3 were significantly elevated in tumor tissues compared to matched normal tissues (Fig. 8E,F). Across three different cohorts (ZHSU, TCGA-COAD, and GSE39582), the corresponding AUC values for diagnosing colon cancer using FSTL3 were 0.833, 0.799, and 0.959, respectively (Fig. 8G). IF staining analysis indicated a significant reduction in CD8+ T cell populations within tumor tissues expressing high levels of FSTL3 compared to paired paracancerous tissues (n = 83, p < 0.01, Fig. 8H). A correlation analysis indicated that FSTL3 was negatively correlated with CD4 + memory T-cells, CD8 + naive T-cells, NK cells, and plasma cells (correlation coefficients < 0 both in TCGA-COAD and GSE39582 cohorts, Fig. 8I). In contrast, FSTL3 showed a significant positive correlation with M2 macrophages and CAFs (correlation coefficients > 0, p < 0.05, Fig. 8I). Furthermore, patients exhibiting high FSTL3 expression had substantially worse clinical outcomes than those with low expression levels (Fig. 8J). Thus, FSTL3 emerges as a model gene that can serve as both a diagnostic and prognostic indicator for colon cancer.
Discussion
Our understanding of cancer genome complexity and the multifaceted effects of therapeutic agents suggests a paradigm shift from a traditional “one gene, one drug” approach to a more integrated “multi-gene, multi-drug” model for making therapeutic decisions. While the CAFscore developed in this study is based on a limited number of genes, its robust efficacy has been clearly demonstrated, indicating its potential for clinical application.
In this study, we detected the expression of three CAF-related genes in tumor tissue and calculated an individualized CAFscore. This approach enabled us to effectively identify patients with high-risk factors for colon cancer. Notably, the CAFscore emerged as an independent risk factor for OS, underscoring its significance in prognostication. When integrated with existing clinical indicators, the nomogram incorporating the CAFscore outperformed traditional biomarkers, including age, gender, and stage, enhancing predictive accuracy. Simultaneously, the CAFscore proved clinically valuable in assessing the ability of tumor cells to evade the immune system and evaluating responses to ICB. By identifying patients with colon cancer who were likely to experience poor therapeutic outcomes early in their treatment journey, the CAFscore facilitated timely adjustments in therapeutic strategies, ultimately leading to improved prognoses. Additionally, the CAFscore served as a tool for assessing the sensitivity of colon cancer to targeted therapies and chemotherapeutic agents. This capability enhances treatment precision, minimizes adverse effects, and supports a more tailored approach to cancer treatment.
Our results demonstrate that the CAFscore, comprised of CRIP2, FSTL3, and SLC2A3, performs exceptionally well as a stand-alone predictor of prognosis and treatment efficacy across both training and external validation datasets. CRIP2 is known to be associated with the VEGFA-VEGFR2 signaling pathway, which plays a crucial role in angiogenesis and tumor progression50. Additionally, SLC2A3, which encodes glucose transporter 3, is linked to mTORC1-mediated glycolysis, significantly impacting tumor growth51. Significant overexpression of FSTL3 has been observed in various malignant tumors, including non-small cell lung cancer52, breast cancer53, thyroid cancer54, and colorectal cancer55. FSTL3 enhances tumor development via TGF-β pathway activation55,56 and promotes lymph node metastasis55. Furthermore, FSTL3 has emerged as a promising biomarker of ECM remodeling55 and is associated with the response to ICB in colorectal cancer patients57. Thus, these findings offer a solid theoretical foundation for utilizing the CAFscore as a prognostic tool for colon cancer patients.
Further research is warranted to identify this tool for predicting therapeutic efficacy and to explore the mechanisms underlying tumor exclusion. First, this study found that colon cancer patients with a high CAFscore had higher levels of TMB, MSI, and the mRNA expression levels of immunosuppressive genes, and the high CAFscore of colon cancer samples was tightly correlated with genomic instability. The CAFscore further stratified the prognosis of patients harboring POLE/POLD1 mutations or exhibiting wild-type status, revealing a significantly poorer prognosis in the cohort with elevated CAFscore. This study suggests that for patients who respond to ICB but have a high CAFscore, combined anti-tumor treatment may be required to reduce the incidence of metastasis and improve prognosis. CAFscore may also be used to evaluate the efficacy of chemotherapy and targeted therapy and assist in guiding the development of comprehensive anti-tumor treatment plans. Therefore, CAFscore may function as a more comprehensive combined classifier.
CAFs play a crucial role in regulating local innate and adaptive immune responses58. They can inhibit the proliferation of CD8+ T cells by secreting TGF-β59 and serve as central nodes for the synthesis and modification of ECM, coordinating regulatory signals among epithelial, immune, and endothelial cell populations60. In this study, we identified significant enrichment of immune receptor activity within the high CAFscore subgroup of colon cancer samples, revealing that a high CAFscore is associated with the formation of an immunosuppressive TME and activation of the NAMPT signaling pathway.
NAMPT expression is upregulated in cancer cells and promotes NAD+ and ATP metabolic reprogramming61,62,63,64. NAMPT-driven polarization of SPP1+ tumor-associated macrophages leaded to an immunosuppressive TME in colon cancer, thereby reducing sensitivity to ICB47,65. Our findings indicate a significant positive correlation between the CAFscore and NAMPT expression, with NAMPT being markedly upregulated in the high CAFscore group of colon cancer. These results align with the existing literature, suggesting that CAFs and NAMPT collaborate to create an inhibitory TME. Importantly, the CAFscore is anticipated to serve as a detectable marker for identifying NAMPT-induced immunosuppressive TME. Additionally, real-world studies are essential to elucidate the crosstalk between CAFs and TME in the context of NAMPT expression, which may provide valuable insights for improving therapeutic strategies targeting the immunosuppressive environment in colon cancer.
The development of pharmacological agents and clinical trials targeting mutational profiles is advancing rapidly. The ZFHX4 mutation has shown a strong correlation with enhanced tumor immunogenicity and an inflamed antitumor immune response, resulting in improved outcomes for ICB treatment across various immunotherapy cohorts66. CSMD1 expression is frequently lost in many epithelial cancers, and the loss of this tumor suppressor gene may contribute to the poor prognosis observed in colon cancer patients67. In contrast, the amplification of CSMD3 is associated with differentiation and lymphatic invasion in colon cancer68. According to tumor proteogenomics, KRAS mutations are emerging as promising candidates for public neoantigens and tumor-associated antigens, serving as viable targets for therapy69. The FDA recently approved a combination of Adagrasib and cetuximab for CRC patients harboring KRAS (G12C) mutations70. Combination therapy with KRAS inhibitors and EGFR inhibitors can delay adaptive resistance, with a better overall response rate of 34%, which is more than 15% higher than monotherapy71. However, this benefit comes with increased risks of adverse effects, such as hypertension and thrombosis. Based on these findings, high CAFscore in colon cancer patients may indicate resistance to KRAS inhibitors. Under the control of adverse reactions, assessing CAFscore can help guide these patients in undergoing combination therapy with KRAS inhibitors and EGFR inhibitors. Meanwhile, targeted therapy and chemotherapeutics in the low CAFscore group were more sensitive than those in the high CAFscore group, especially first-line chemotherapy agents such as fluorouracil, capecitabine, and oxaliplatin. These findings demonstrate the potential usefulness of the CAFscore in precision medicine for treating colon cancer.
Our study had some limitations that should be mentioned. Firstly, the impact of MSI status on prognosis is closely related to cancer staging. MSI-H is considered a protective factor in the early stages of the disease (Stage I-II)72,73; whereas as the disease progresses (Stage III-IV), MSI-H may gradually reverse to become a risk factor74,75. In this study, we preliminarily found that patients with MSI-H and dMMR exhibited higher CAFscore based on all stages (I-IV). However, the relationship between CAFscore and MSI status in different stages requires further investigation. Secondly, CAFscore was effective in identifying potentially high-risk individuals among Stage II/III MSI-H patients. An important consideration is whether there exists a modifiable transition state under the stress of TME in the biological behavior of MSI-H colon cancer, in which it shifts from being protective72,73 to a risk factor74,75 as the disease progresses. More importantly, whether this transition is linked to CAF requires further exploration. Thirdly, it was a retrospective study based on public and local cohorts. The prognostic prediction capability of CAF-related molecular signatures needs cross-validation using multicenter and prospective datasets. There has been a notable increase in clinical trials examining the combination of targeted CAF therapies with ICB in recent years76. Future research should leverage clinical data to validate the relationship between CAFscore and therapies in colon cancer patients. Despite these limitations, our findings provide a valuable foundation for the development of comprehensive CAF-based classifier for colon cancer.
Conclusion
Our study identified high-risk colon cancer patients by detecting three CAF-related genes in tumor tissue and calculating an individualized CAFscore, which was an independent risk factor for OS. Additionally, it proved clinically valuable for evaluating the immunosuppressive TME and the response to ICB, enabling the early identification of colon cancer patients at high risk for poor outcomes and facilitating timely treatment adjustments. Moreover, the CAFscore predicted the sensitivity to targeted therapies and chemotherapy, enhancing treatment accuracy while minimizing adverse events. This integrated analysis deepens our understanding of the molecular signatures related to CAFs and their significance in predicting prognosis and therapy sensitivity in colon cancer.
Data availability
The datasets generated and analyzed in the current study are accessible through the TCGA and GEO databases, including the TCGA-COAD, GSE39582, GSE236581, and GSE110009 cohorts. All other raw data are available upon request from the corresponding author.
References
Siegel, R. L., Giaquinto, A. N. & Jemal, A. Cancer statistics, 2024. CA Cancer J. Clin. 74 (1), 12–49 (2024).
Seshagiri, S. The burden of faulty proofreading in colon cancer. Nat. Genet. 45 (2), 121–122 (2013).
Wong, N. A. C. S. et al. Current dilemmas in the pathological staging of colorectal cancer: the results of a National survey. Histopathology 78 (4), 634–639 (2021).
Dienstmann, R. et al. Consensus molecular subtypes and the evolution of precision medicine in colorectal cancer. Nat. Rev. Cancer. 17 (2), 79–92 (2017).
Fontana, E., Eason, K., Cervantes, A., Salazar, R. & Sadanandam, A. Context matters-consensus molecular subtypes of colorectal cancer as biomarkers for clinical trials. Ann. Oncol. 30 (4), 520–527 (2019).
Mariathasan, S. et al. TGFβ attenuates tumour response to PD-L1 Blockade by contributing to exclusion of T cells. Nature 554 (7693), 544–548 (2018).
Tauriello, D. V. F. et al. TGFβ drives immune evasion in genetically reconstituted colon cancer metastasis. Nature 554 (7693), 538–543 (2018).
Martinez-Ordoñez, A. et al. Hyaluronan driven by epithelial aPKC deficiency remodels the microenvironment and creates a vulnerability in mesenchymal colorectal cancer. Cancer Cell. 41(2), 252–271.e9 (2023).
Herrera, M. et al. Cancer-associated fibroblast-derived gene signatures determine prognosis in colon cancer patients. Mol. Cancer. 20 (1), 73 (2021).
Zhang, L. et al. Cancer-associated fibroblast-related gene signatures predict survival and drug response in patients with colorectal cancer. Front. Genet. 13, 1054152 (2022).
Liu, T. et al. Cancer-associated fibroblasts: an emerging target of anti-cancer immunotherapy. J. Hematol. Oncol. 12 (1), 86 (2019).
Zhang, R. et al. Cancer-associated fibroblasts enhance tumor-associated macrophages enrichment and suppress NK cells function in colorectal cancer. Cell. Death Dis. 10 (4), 273 (2019).
Biffi, G. & Tuveson, D. A. Diversity and biology of Cancer-Associated fibroblasts. Physiol. Rev. 101 (1), 147–176 (2021).
Smith, H. G. et al. Challenges presented by complete response to immune checkpoint Blockade in patients with dMMR colorectal cancer: A case report. Int. J. Surg. Case Rep. 106, 108286 (2023).
Mandal, R. et al. Genetic diversity of tumors with mismatch repair deficiency influences anti-PD-1 immunotherapy response. Science 364 (6439), 485–491 (2019).
Lei, Y., Li, X., Huang, Q., Zheng, X. & Liu, M. Progress and challenges of predictive biomarkers for immune checkpoint Blockade. Front. Oncol. 11, 617335. https://doi.org/10.3389/fonc.2021.617335 (2021).
Ambrosini, M. et al. Immune checkpoint inhibitors for POLE or POLD1 proofreading-deficient metastatic colorectal cancer. Ann. Oncol. 35 (7), 643–655 (2024).
Overman, M. J. et al. Nivolumab in patients with metastatic DNA mismatch repair-deficient or microsatellite instability-high colorectal cancer (CheckMate 142): an open-label, multicentre, phase 2 study. Lancet Oncol. 18 (9), 1182–1191 (2017).
Fumet, J. D., Truntzer, C., Yarchoan, M. & Ghiringhelli, F. Tumour mutational burden as a biomarker for immunotherapy: current data and emerging concepts. Eur. J. Cancer. 131, 40–50 (2020).
Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 9, 559 (2008).
Leek, J. T., Johnson, W. E., Parker, H. S., Jaffe, A. E. & Storey, J. D. The Sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics 28 (6), 882–883 (2012).
Wang, R. et al. Single-cell genomic and transcriptomic landscapes of primary and metastatic colorectal cancer tumors. Genome Med. 14 (1), 93 (2022).
Chen, Y. et al. Spatiotemporal single-cell analysis decodes cellular dynamics underlying different responses to immunotherapy in colorectal cancer. Cancer Cell. 42 (7), 1268–85e7 (2024).
Racle, J., de Jonge, K., Baumgaertner, P., Speiser, D. E. & Gfeller, D. Simultaneous enumeration of cancer and immune cell types from bulk tumor gene expression data. Elife 6. https://doi.org/10.7554/eLife.26476 (2017).
Becht, E. et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome Biol. 17 (1), 218 (2016).
Aran, D., Hu, Z. & Butte, A. J. xCell: digitally portraying the tissue cellular heterogeneity landscape. Genome Biol. 18 (1), 220 (2017).
Yoshihara, K. et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat. Commun. 4, 2612. https://doi.org/10.1038/ncomms3612 (2013).
Zheng, H., Liu, H., Ge, Y. & Wang, X. Integrated single-cell and bulk RNA sequencing analysis identifies a cancer associated fibroblast-related signature for predicting prognosis and therapeutic responses in colorectal cancer. Cancer Cell. Int. 21 (1), 552 (2021).
Liu, Z. et al. Machine learning-based integration develops an immune-derived LncRNA signature for improving outcomes in colorectal cancer. Nat. Commun. 13 (1), 816 (2022).
Jiang, P. et al. Signatures of T cell dysfunction and exclusion predict cancer immunotherapy response. Nat. Med. 24 (10), 1550–1558 (2018).
Fu, J. et al. Large-scale public data reuse to model immunotherapy response and resistance. Genome Med. 12 (1), 21 (2020).
Zheng, H., Liu, H., Li, H., Dou, W. & Wang, X. Weighted gene Co-expression network analysis identifies a Cancer-Associated fibroblast signature for predicting prognosis and therapeutic responses in gastric Cancer. Front. Mol. Biosci. 8, 744677 (2021).
Simon, N., Friedman, J., Hastie, T. & Tibshirani, R. Regularization paths for Cox’s proportional hazards model via coordinate descent. J. Stat. Softw. 39 (5), 1–13 (2011).
Subramanian, A. et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc. Natl. Acad. Sci. U S A. 102 (43), 15545–15550 (2005).
Yu, G., Wang, L. G., Han, Y. & He, Q. Y. ClusterProfiler: an R package for comparing biological themes among gene clusters. OMICS 16 (5), 284–287 (2012).
Charoentong, P. et al. Pan-cancer Immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint Blockade. Cell. Rep. 18 (1), 248–262 (2017).
Ying, L. et al. An oncogene regulating chromatin favors response to immunotherapy: oncogene CHAF1A and immunotherapy outcomes. Oncoimmunology 13 (1), 2303195 (2024).
Li, J. C. et al. Up-regulation of Krüppel-like factor 8 promotes tumor invasion and indicates poor prognosis for hepatocellular carcinoma. Gastroenterology 139 (6), 2146–2157e12 (2010).
Uhlén, M. et al. Proteomics. Tissue-based map of the human proteome. Science 347 (6220), 1260419 (2015).
Stuart, T. et al. Comprehensive integration of Single-Cell data. Cell 177 (7), 1888–902e21 (2019).
Aran, D. et al. Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20 (2), 163–172 (2019).
Kanehisa, M., Furumichi, M., Sato, Y., Matsuura, Y. & Ishiguro-Watanabe, M. KEGG: biological systems database as a model of the real world. Nucleic Acids Res. 53 (D1), D672–D677 (2025).
Kanehisa, M. Toward Understanding the origin and evolution of cellular organisms. Protein Sci. 28 (11), 1947–1951 (2019).
Kanehisa, M. & Goto, S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28 (1), 27–30 (2000).
Huang, X., Zhang, G. & Liang, T. Cancer environmental immunotherapy: starving tumor cell to death by targeting TGFB on immune cell. J. Immunother Cancer. 9 (6), e002823 (2021).
Tharp, K. M. et al. Tumor-associated macrophages restrict CD8(+) T cell function through collagen deposition and metabolic reprogramming of the breast cancer microenvironment. Nat. Cancer. 5 (7), 1045–1062 (2024).
Hong, S. M. et al. NAMPT-driven M2 polarization of tumor-associated macrophages leads to an immunosuppressive microenvironment in colorectal cancer. Adv. Sci. (Weinh). 11 (14), e2303177 (2024).
Houlahan, K. E. et al. Germline-mediated immunoediting sculpts breast cancer subtypes and metastatic proclivity. Science 384 (6699), eadh8697 (2024).
Datta, J. et al. Coaltered Ras/B-raf and TP53 is associated with extremes of survivorship and distinct patterns of metastasis in patients with metastatic colorectal cancer. Clin. Cancer Res. 26 (5), 1077–1085 (2020).
Cheung, A. K. et al. Cysteine-rich intestinal protein 2 (CRIP2) acts as a repressor of NF-kappaB-mediated proangiogenic cytokine transcription to suppress tumorigenesis and angiogenesis. Proc. Natl. Acad. Sci. U S A. 108 (20), 8390–8395 (2011).
Zha, X. et al. NFκB up-regulation of glucose transporter 3 is essential for hyperactive mammalian target of rapamycin-induced aerobic Glycolysis and tumor growth. Cancer Lett. 359 (1), 97–106 (2015).
Gao, L., Chen, X., Wang, Y. & Zhang, J. Up-Regulation of FSTL3, regulated by LncRNA DSCAM-AS1/miR-122-5p axis, promotes proliferation and migration of non-small cell lung cancer cells. Onco Targets Ther. 13, 2725–2738 (2020).
Razanajaona, D. et al. Silencing of FLRG, an antagonist of activin, inhibits human breast tumor cell growth. Cancer Res. 67 (15), 7223–7229 (2007).
Panagiotou, G., Ghaly, W., Upadhyay, J., Pazaitou-Panayiotou, K. & Mantzoros, C. S. Serum follistatin is increased in thyroid cancer and is associated with adverse tumor characteristics in humans. J. Clin. Endocrinol. Metab. 106 (5), e2137–e2150 (2021).
Yang, C., Cao, F., Huang, S. & Zheng, Y. Follistatin-Like 3 correlates with lymph node metastasis and serves as a biomarker of extracellular matrix remodeling in colorectal cancer. Front. Immunol. 12, 717505. https://doi.org/10.3389/fimmu.2021.717505 (2021).
Li, C. H. et al. The cytoplasmic expression of FSTL3 correlates with colorectal cancer progression, metastasis status and prognosis. J. Cell. Mol. Med. 27 (5), 672–686 (2023).
Li, H. et al. FSTL3 promotes tumor immune evasion and attenuates response to anti-PD1 therapy by stabilizing c-Myc in colorectal cancer. Cell. Death Dis. 15 (2), 107 (2024).
Sahai, E. et al. A framework for advancing our Understanding of cancer-associated fibroblasts. Nat. Rev. Cancer. 20 (3), 174–186 (2020).
Takahashi, H. et al. Immunosuppressive activity of cancer-associated fibroblasts in head and neck squamous cell carcinoma. Cancer Immunol. Immunother. 64 (11), 1407–1417 (2015).
Parker, A. L. et al. Extracellular matrix profiles determine risk and prognosis of the squamous cell carcinoma subtype of non-small cell lung carcinoma. Genome Med. 14 (1), 126 (2022).
Menssen, A. et al. The c-MYC oncoprotein, the NAMPT enzyme, the SIRT1-inhibitor DBC1, and the SIRT1 deacetylase form a positive feedback loop. Proc. Natl. Acad. Sci. U S A. 109 (4), E187–E196 (2012).
Sampath, D., Zabka, T. S., Misner, D. L., O’Brien, T. & Dragovich, P. S. Inhibition of nicotinamide phosphoribosyltransferase (NAMPT) as a therapeutic strategy in cancer. Pharmacol. Ther. 151, 16–31 (2015).
Lee, J. et al. Selective cytotoxicity of the NAMPT inhibitor FK866 toward gastric cancer cells with markers of the epithelial-mesenchymal transition, due to loss of NAPRT. Gastroenterology 155 (3), 799–814e13 (2018).
Lucena-Cacace, A., Otero-Albiol, D., Jiménez-García, M. P., Muñoz-Galvan, S. & Carnero, A. NAMPT is a potent oncogene in colon cancer progression that modulates cancer stem cell properties and resistance to therapy through Sirt1 and PARP. Clin. Cancer Res. 24 (5), 1202–1215 (2018).
Travelli, C. et al. Nicotinamide phosphoribosyltransferase acts as a metabolic gate for mobilization of myeloid-derived suppressor cells. Cancer Res. 79 (8), 1938–1951 (2019).
Mao, Y. et al. The landscape of objective response rate of anti-PD-1/L1 monotherapy across 31 types of cancer: a system review and novel biomarker investigating. Cancer Immunol. Immunother. 72 (2), 2483–2498 (2023).
Sheffer, M. et al. Association of survival and disease progression with chromosomal instability: a genomic exploration of colorectal cancer. Proc. Natl. Acad. Sci. U S A. 106 (17), 7131–7136 (2009).
Zhang, R. & Song, C. Loss of CSMD1 or 2 May contribute to the poor prognosis of colorectal cancer patients. Tumour Biol. 35 (5), 4419–4423. https://doi.org/10.1007/s13277-013-1581-6 (2014).
Savage, S. R. et al. Pan-cancer proteogenomics expands the landscape of therapeutic targets. Cell 187 (16), 4389–407e15 (2024).
Yaeger, R. et al. Adagrasib with or without cetuximab in colorectal cancer with mutated KRAS G12C. N Engl. J. Med. 388 (1), 44–54 (2023).
Yaeger, R. et al. Efficacy and safety of adagrasib plus cetuximab in patients with KRAS G12C-mutated metastatic colorectal cancer. Cancer Discov. 14 (6), 982–993 (2024).
Benson, A. B. et al. Colon cancer, version 3.2024, NCCN clinical practice guidelines in oncology. J. Natl. Compr. Canc Netw. 22 (2 D), e240029 (2024).
Sargent, D. J. et al. Prognostic impact of deficient mismatch repair (dMMR) in 7,803 stage II/III colon cancer (CC) patients (pts): A pooled individual Pt data analysis of 17 adjuvant trials in the ACCENT database. J. Clin. Oncol. 32, 3507 (2014).
Venderbosch, S. et al. Mismatch repair status and BRAF mutation status in metastatic colorectal cancer patients: a pooled analysis of the CAIRO, CAIRO2, COIN, and FOCUS studies. Clin. Cancer Res. 20 (20), 5322–5330 (2014).
Mohan, H. M. et al. Microsatellite instability is associated with reduced disease specific survival in stage III colon cancer. Eur. J. Surg. Oncol. 42 (11), 1680–1686 (2016).
Zhang, H. et al. Define cancer-associated fibroblasts (CAFs) in the tumor microenvironment: new opportunities in cancer immunotherapy and advances in clinical trials. Mol. Cancer. 22 (1), 159 (2023).
Acknowledgements
We acknowledge the staff of the Biobank of Zhongda Hospital, affiliated with Southeast University, for their invaluable technical assistance.
Funding
This study was supported by the Bethune Charitable Foundation (BCF-XD-ZL-20221101-055), the National Natural Science Foundation of China (Grant numbers 82372023 and 82172010), and the National Key Research and Development Program of China (Grant number 2022YFE0116700).
Author information
Authors and Affiliations
Contributions
LY and LL were responsible for the study design and data analysis and contributed to drafting the manuscript. LZ, YC, CH, and JZ conducted the procedures and collected data. JX participated in drafting the manuscript. All authors have read and approved the final manuscript.
Corresponding author
Ethics declarations
Competing interests
The authors declare no competing interests.
Ethics approval and consent to participate
This study retrospectively gathered data from prior cases for anonymous analysis, ensuring that patient privacy information was not disclosed. This study was conducted in accordance with the principles of the Declaration of Helsinki. The study was approved by the Clinical Research Ethics Committee of Zhongda Hospital, Affiliated to Southeast University (#2022ZDSYLL273-P01). Informed consent was obtained from all study participants.
Additional information
Publisher’s note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.
About this article
Cite this article
Ying, L., Zhang, L., Chen, Y. et al. Predicting immunotherapy prognosis and targeted therapy sensitivity of colon cancer based on a CAF-related molecular signature. Sci Rep 15, 6387 (2025). https://doi.org/10.1038/s41598-025-90899-z
Received:
Accepted:
Published:
DOI: https://doi.org/10.1038/s41598-025-90899-z
