A Comprehensive Multi-Omics Analysis Identifies a Robust Scoring System for Cancer-Associated Fibroblasts and Intervention Targets in Colorectal Cancer

doi:10.21203/rs.3.rs-3341136/v1

Download PDF

Research Article

A Comprehensive Multi-Omics Analysis Identifies a Robust Scoring System for Cancer-Associated Fibroblasts and Intervention Targets in Colorectal Cancer

https://doi.org/10.21203/rs.3.rs-3341136/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 13 Mar, 2024

Read the published version in Journal of Cancer Research and Clinical Oncology →

You are reading this latest preprint version

Background

Cancer-associated fibroblasts (CAF) play a critical role in promoting tumor growth, metastasis, and immune evasion. While numerous studies have investigated CAF, there remains a paucity of research on their clinical application in colorectal cancer (CRC).

Methods

In this study, we collected differentially expressed genes between CAF and normal fibroblasts (NF) from previous CRC studies, and utilized machine learning analysis to differentiate two distinct subtypes of CAF in CRC. To enable practical application, a CAF-related genes (CAFGs) scoring system was developed based on multivariate Cox regression. We then conducted functional enrichment analysis, Kaplan-Meier plot, consensus molecular subtypes (CMS) classification, and Tumor Immune Dysfunction and Exclusion (TIDE) algorithm to investigate the relationship between the CAFGs scoring system and various biological mechanisms, prognostic value, tumor microenvironment, and response to immune checkpoint blockade (ICB) therapy. Moreover, single-cell transcriptomics and proteomics analyses have been employed to validate the significance of scoring system related molecules in the identity and function of CAF.

Results

We unveiled significant distinctions in tumor immune status and prognosis not only between the CAF clusters, but also across high and low CAFGs groups. Specifically, patients in CAF cluster 2 or with high CAFGs scores exhibited higher CAF markers and were enriched for CAF-related biological pathways such as such as epithelial-mesenchymal transition (EMT) and angiogenesis. Additionally, CAFGs score was identified as a risk index and correlated with poor overall survival (OS), progression-free survival (PFS), disease-free survival (DFS), and recurrence-free survival (RFS). High CAFGs scores were observed in patients with advanced stages, CMS4, as well as lymphatic invasion. Furthermore, elevated CAFG scores in patients signified a suppressive tumor microenvironment characterized by the upregulation of programmed death-ligand 1 (PD-L1), T cell dysfunction, exclusion, and TIDE score. And high CAFGs scores can differentiate patients with lower response rates and poor prognosis under ICB therapy. Notably, single-cell transcriptomics and proteomics analyses identified several molecules related to CAF identity and function, such as FSTL1, IGFBP7, and FBN1.

Conclusion

We constructed a robust CAFGs score system with clinical significance using multiple CRC cohorts. Additionally, we identified several molecules related to CAF identity and function that could be potential intervention targets for CRC patients.

CAF

Colorectal cancer

Prognosis

ICB response biomarkers

According to the 2020 global cancer statistics, colorectal cancer (CRC) is one of the most common malignant tumors in the world, ranking third and second in incidence and mortality, respectively¹. At present, there are many treatments for CRC, including radiotherapy and chemotherapy, immunotherapy, surgical treatment, which greatly improve the survival of patients. However, the response of patients to treatment is different, resulting in different outcomes. In view of this phenomenon, CRC consensus molecular subtypes (CMS)², a recently established classification based on the transcriptome data may provide an explanation, in which the CMS4 (Mesenchymal), characterized by prominent transforming growth factor β (TGF-β) activation, stromal invasion, and angiogenesis, are associated with poor prognosis, and suggest the importance of cancer-associated fibroblasts (CAF), the major component of the stroma, to the prognosis of CRC patients .

Cancer-associated fibroblasts (CAF), one of the plastic cells types in TME, has different origins, including resident fibroblasts, pericytes, endothelial cells, adipocytes and so on. For the origin of CAF in CRC, most of CAF have been proved to be produced by the proliferation of intestinal pericryptal leptin receptor (Lepr) + cells³. In the progression of CRC, CAF may promote or inhibit tumor, but the prevailing idea is considered to promote tumor. CAF contribute significantly to tumorigenesis^{4, 5}, angiogenesis^{6, 7}, immunosuppression⁸, metastasis, and drug resistance⁹ in CRC. CAF has several subtypes, among which inflammatory CAF (iCAF) and myofibroblast CAF (myCAF) have been studied. In pancreatic cancer, iCAF can secrete high levels of inflammatory cytokines far away from tumor cells, while myCAF can produce matrix contractile phenotype, and are adjacent to tumor cells¹⁰. In CRC, myCAF and iCAF are induced by high and low levels of Wnt activity, respectively. And iCAF promote EMT phenotype, while myCAF reverse the phenotype¹¹. In preoperative radiotherapy and chemotherapy for rectal cancer, the inflammatory polarization of CAF leads to the resistance of radiotherapy and chemotherapy, and promotes tumor progression¹². After cetuximab treatment, CAF can make neighboring cancer cells resistant to cetuximab in CRC¹³. Additionally, it has been reported that some molecules expressing on CAF, including WNT2¹⁴, WNT5a¹⁵, CLEC3B⁵, IFNAR1¹⁶, IL-34¹⁷, miR-1246¹⁸, and FAP¹⁹, play an important role in the development of CRC, and are expected to become a target for anticancer therapy. And there are many models have been established to study CAF in CRC, such as 3D model of tumor tissue in vitro to simulate the physiological function of cells in vivo²⁰, in vitro co-culture model of patient-derived organ-like organ (PDO) and patient-derived CAF^{21, 22}, and mouse xenotransplantation model co-injected with CAF and CRC cell line²³. These models are helpful to investigate the functions of various subtypes of CAF, and find the therapy strategies targeting CAF in CRC. Although there have been many in vitro and in vivo experiments focusing on the functional characteristics of CAF, applying the research findings in clinical practice is still an urgent issue.

Previous studies have shown that the heterogeneity of CAF significantly correlates with the efficacy of ICB therapy. For example, Wang et al. has used single-cell RNA-seq to analyze the heterogeneity of CAF, and identified a novel fibroblast subtype, independent of iCAF and myCAF, which was termed meCAF characterized by highly active glycolysis, and associated with better response to anti-PD-1 therapy in pancreatic ductal adenocarcinoma²⁴. However, compared with bulk RNA-seq, single-cell RNA-seq is dramatically expensive, and not suitable for wide clinical application.

In this study, we performed unsupervised clustering in a large-scale CRC cohort based on CAF-related genes, and identified two groups of patients with distinct biological functions related to CAF (named CAF cluster 1 and 2). Furthermore, we constructed a novel CAF scoring system composed of 15 genes, which were associated with poor overall, disease-free, recurrence-free, and disease-specific survival, (OS, DFS, RFS and DSS), and had the potential to guide ICB treatment. Moreover, single-cell sequencing and proteomics data suggest that these 15 genes might be linked to CAF identity and function, thus rendering them potential therapeutic targets for CAF intervention in CRC.

2.1 Construction of CAF clusters with different prognosis and immune states

The study flow diagram is presented in Fig. 1. Firstly, we collected 596 differentially expressed genes (named CAF genes) between CAF and normal fibroblasts (NF) in CRC from a previous study²⁵ (Supplementary Data 1). To find prognostic genes, we then performed univariate Cox regression analysis in GEO combined cohort and recognized 115 potential prognostic genes among above genes (Supplementary Data 2). Furthermore, we conducted unsupervised clustering in GEO combined cohort based on 115 prognostic CAF genes using the ConsesusClusterPlus R package. As shown in Supplementary Fig. 1, the clustering results were most stable when patients were divided into two groups (defined as CAF cluster 1 and 2). The PCA plot shows significant differences in gene expression profiles between the two clusters (Fig. 2A). Remarkably, we found several previously identified CAF-related markers^{26, 27}, including ACTA2, FAP, FOXL1, MCAM, and PDGFRA, were substantially upregulated in CAF cluster 2 relative to CAF cluster 1 (Fig. 2B), suggesting that the distinct status of CAF is correlated to group classification. In addition, MCPcounter analysis showed that fibroblast scores of patients in CAF cluster 2 were markedly higher than patients in CAF cluster 1 (Fig. 2C), while KM survival analysis illustrated that the OS was notably better for patients in CAF cluster 1 than those in CAF cluster 2 (Fig. 2D, log rank test, p = 0.0024). These results imply that the varied CAF status significantly impacts the survival of patients with CRC.

To dissect the underlying biological functions between the CAF clusters, we collected 11 tumorigenesis-related pathways from prior research. Our findings showed that the angiogenesis, TGF β, and F-TBRS signature scores were dramatically elevated in CAF cluster 2 compared to CAF cluster 1 (Fig. 2E). Additionally, the genes expression of APM, CD8 + Teff, and ICI pathways were dramatically increased in CAF cluster 2 (Fig. 2E, F), suggesting potentially distinct biological functions between the two groups. To further understand the immune status between the two groups, we analyzed the expression patterns of 122 immunomodulators (including MHC, receptors, chemokines, and immunostimulants) between CAF cluster 1 and 2, most of which were highly expressed in CAF cluster 2 (Supplementary Fig. 2A). The abundance of most tumor infiltrating lymphocytes inferred by MCPcounter analysis were also significantly higher in CAF cluster 2 than CAF cluster 1, such as T cells, cytotoxic lymphocytes, and neutrophils (Fig. 2C). Besides, ssGSEA-inferred adaptive and innate immunity scores were also significantly increased in CAF cluster 2 (Supplementary Fig. 2B). Whereas, the expression of immune checkpoint molecules was significantly higher in CAF cluster 2 than in CAF cluster 1. These results demonstrate that the CAF clusters exhibit distinct immune microenvironments, with CAF cluster 2 exhibiting an inhibitory immune microenvironment.

2.2 Identification the key genes affecting the prognosis of patients in different CAF clusters

In order to identify the key genes affecting the survival of patients in CAF clusters, WGCNA analysis was carried out, and CAF clusters 1 and 2 were used as the traits. As shown in Fig. 3A, B, the soft threshold power of β was set as 4 when scale-free topology model-fit R = 0.9. Then, we identified 16 modules, except the grey module (Fig. 3C-E). Module-trait heatmap shows that the blue module was the most closely related to CAF clusters (Fig. 3F, G). And the biological functions of genes in the blue modules were explored using GO and KEGG analysis. When the adjusted P-value was less than 0.05, 27 and 12 items were identified by GO and KEGG analyses, respectively (Supplementary Data 3, 4). The top 10 enrichment items in GO and KEGG analyses included extracellular matrix structural constituent, collagen binding, ECM-receptor interaction, PI3K-Akt, and TGF − beta signaling pathway (Fig. 3H, I), which were in consistent with the functions of CAF. Therefore, the blue module containing 1089 genes was identified as the key module, among which 100 genes meeting GS > 0.2 and MM > 0.8 were considered as the critical genes related to CAF clusters (CAFGs). In addition, univariate cox regression analysis identified 43 of the 100 CAFGs was prognostic genes with a p value less than 0.01. Then, these genes were again used for unsupervised clustering in GEO combined cohort (Supplementary Data 5). The detail processes of unsupervised clustering were shown in the Supplementary Fig. 3. Surprisingly, when the patients were again divided into two groups, defined as CAFGs clusters 1 and 2, the clustering results were the most stable. The PCA plot shows that there were significant differences in gene expression profiles between the CAFGs clusters (Fig. 4A). The markers associated with CAF and the fibroblast score inferred by MCPcounter analysis were significantly higher in CAFGs cluster 1 than CAFGs cluster 2 (Fig. 4B, F). Consistently, the OS of CAFGs cluster 1 was significantly worse than that of CAFGs cluster 2 (Fig. 4C, log rank test, p 0.04). Additionally, the heatmap of 11 tumorigenesis-related pathways shows that the angiogenesis, TGF β, and F-TBRS signature score of CAFGs cluster 1 were higher than that of CAFGs cluster 2, and the levels of APM, CD8 + Teff, and ICI signatures were also increased in CAFGs cluster 1 (Supplementary Fig. 4A, B). Besides, the expression 122 immunomodulators (Fig. 4D), adaptive and innate immunity (Fig. 4E), and most tumor infiltrating lymphocytes (Fig. 4F, Supplementary Fig. 4C) were higher in CAFGs cluster1 than CAFGs cluster 2. Furthermore, we observed most of the patients consisting of CAFGs cluster 1 were from CAF cluster 2 (Fig. 4G). Therefore, these findings demonstrated the crucial roles of CAFGs, which can reproduce the biological category of CAF clusters.

2.3 Colorectal cancer patients with high CAFGs score have poor outcomes in multiple colorectal cohorts

Gene model plays an important role in clinical application. In order to construct a scoring model for clinical application, the CAFGs meeting a p value < 0.2 in univariate analysis were included in multivariate unicox regression analysis. Finally, 15 genes with a p value < 0.05 were obtained, including FNDC1, FRMD6, FBN1, RAB31, GLT8D2, COL1A2, GLIS2, COL8A1, GPC6, COL3A1, PRICKLE1, FSTL1, HLX, IGFBP7, and EFS. These genes were considered as the important prognostic factors, and again incorporated in the cox model. Next, using the expression values of these 15 genes and their corresponding regression coefficients, a scoring model, named CAFGs scoring system, was constructed (Supplementary Data 6). We then included TCGA COAD and GSE39582 cohorts as external and internal validation sets. According to the best cutoff value determined by survminer R package, patients in these cohort were divided into high and low CAFGs score groups (Supplementary Data 7). We found that the expression levels of CAF markers were significantly higher in the high CAFGs score group than in the low CAFGs score group (Fig. 5A: GEO combined; Fig. 6A: TCGA COAD; Supplementary Fig. 5A: GSE39582). Then, we observed that patients with high CAFGs scores showed worse OS in GEO combined cohort, which were also verified in the internal and external cohorts (Fig. 5B: GEO-combined; Fig. 6B: TCGA COAD cohort; Supplementary Fig. 5B: GSE39582). In additional, we analyzed the DFS, RFS, and DSS in GSE39582, GSE17536, and GSE17537 cohorts. The results show that the RFS, DFS, and DSS of patients with high CAFGs scores were also significantly worse than that of patients with low CAFGs scores (Fig. 5C-F, GSE39582 RFS, GSE17537 DFS, GSE17536 DFS, GSE17536 DSS). Further analysis demonstrated that high CAFGs scores also represented poor OS in both early (Stag I&II) and advanced (stage III&IV) patients (Fig. 5G, H; Fig. 6C, D; Supplementary Fig. 5C, D). And the CAFGs scores of patients in stage III&IV were significantly higher than that of patients in Stag I&II (Fig. 5I; Fig. 6E; Supplementary Fig. 5E). CMS classification, a widely used classification system in CRC, has strong prognostic implications, and includes four subtypes, such as CMS1 (MSI Immune), CMS2 (Canonical), CMS3 (Metabolic), and CMS4 (Mesenchymal). To our surprise, CAFGs scores of patients in CMS subtype 4 were significantly higher than those of patients in CMS subtypes 1–3 in the combined GEO cohort, as well as TCGA COAD and GSE39582 cohorts (Fig. 5J; Fig. 6F; Supplementary Fig. 5F; Supplementary Data 8). CMS subtype 4 is mesenchymal subtype characterized by prominent transforming growth factor β (TGF-β) activation, stromal invasion, and angiogenesis; and has been reported to be associated with poor prognosis.

Next, we conducted a comprehensive investigation of the CAFGs scoring system in the TCGA COAD cohort and examined its association with various clinicopathologic features, such as T, N, and M stages, as well as venous and lymphatic invasion. Our observations revealed that patients in the T3&4 stage, N + stage, and those with lymphatic invasion had significantly higher CAFGs scores than patients in the T1&2 stage, N0 stage, and without lymphatic invasion (Fig. 6G-I; Supplementary Fig. 6). This may underscore the importance of early detection and monitoring of lymphatic invasion and metastasis in patient with high CAFGs score. Moreover, we also performed GSVA by using hallmark gene sets from the MSigDB website, and revealed that several signaling pathways, including EMT, were significantly upregulated in patients with high CAFG scores (Supplementary Fig. 7A).

2.4 Patients with a high CAF score develop resistance to immunotherapy

It’s well known that the therapeutic effect of ICB treatment is closely related to the tumor immune microenvironment, including the abundance of CAF. Therefore, we used TIDE analysis to explore the relationship between CAFGs score and therapeutic response to ICB, which focused on two mechanisms of tumor immune evasion, namely, T cell dysfunction and exclusion. Through TIDE analysis, we found a significant correlation between CAFGs score and T cell disfunction score, T cell exclusion score, and TIDE score (Fig. 7A). In addition, the expression levels of ICI genes were markedly higher in the high CAFG score group compared to the low CAFG score group (Supplementary Fig. 7B). These findings indicated that patients with high CAFGs score may suffer from immune evasion and resistance to immunotherapy. Furthermore, we collected two patient cohorts receiving anti–PD-1 or anti–PD-L1 therapy (GSE78220, and IMvigor210) to analyze the relationship between CAFGs score and ICB efficacy. Indeed, in the high CAFGs score group, we found a higher proportion of no responders (Fig. 7B, C). Consistently, patients with high CAFGs score had significantly worse overall survival (OS) after receiving ICB therapy (Fig. 7D, E). These findings indicated that CAFGs score may have the potential to identify CRC patients who were sensitive to ICB therapy.

2.5 Exploring potential molecules associated with the identity and function of CAF

Several of these 15 genes belong to CAFGs scoring system were previously reported to be specifically associated with the identity and functionality of CAF, including COL1A2, COL3A1, COL8A1, and FBN1, implying the remaining genes may also correlate with CAF. Further analysis of the relationship between these 15 genes and CAF levels in TIDE website revealed that almost all of the genes were significantly correlated with CAF levels (Fig. 8A). Single-cell transcriptomic data from patients with CRC also indicated that nearly all of the genes were significantly expressed in stromal cell types, such as IGFBP7, GLT8D2, FSTL1, GPC6, and FRMD6 (Fig. 8B-D, Supplementary Fig. 8). The proteomic data obtained from CAF derived from AOM/DSS-induced CRC mice and normal mice-derived NF indicates that in addition to commonly known CAF related proteins (such as COL3A1, COL1A2, and FBN1), IGFBP7 and FSTL1 were significantly upregulated in the CAF conditioned medium (Fig. 8E). And transcriptional data reveals that compared to the adjacent non-tumorous samples, expression of COL3A1, COL1A2, FBN1, IGFBP7, and FSTL1 was significantly elevated in the tumor tissue (Fig. 8F). These results suggest that the 15 marker models we identified might play crucial roles in CAF-mediated CRC progression.

Cancer-associated fibroblasts (CAF), an important component of the tumor microenvironment, are involved in tumor initiation, progression and metastasis^28–30. Although many promising achievements have been made in the basic research of CAF, applying the research findings in clinical practice is an urgent issue.

In our study, we have established two types of classification based on gene expression related to CAF, named CAF clusters and CAFGs scoring system. Of note, the levels of CAF-related biological terms, such as TGF-b, F-TBRS, angiogenesis, were highly expressed in CAF cluster 2 and high CAFGs score group. Among these groups, CAF markers collected from previous researches significantly increased, including ACTA2, FAP, MCAM, and PDGFRA^{26, 31}. These results indicated that there were more CAF infiltrating in CAF Cluster 2 and high CAFGs score group. And the high CAF scoring not only correlated with poor OS, but also with poor DFS, PFS, and DSS. Additionally, patients in the CMS4 group had the highest CAFGs scores. CMS is a wide used CRC classification, among which CMS4 is characterized by prominent TGF-β activation, stromal invasion, and angiogenesis; and associated with prognosis, indicating the CAFGs scoring system may have important application potential.

It encourages us to explore the association between CAF scoring and ICB therapy efficacy. Notably, in our study, patients in CAF cluster 2 and high CAFGs score group also had higher expression of inhibitory immune checkpoints (such as, CD274 [PD-L1], CTLA4, and TIGIT), suggesting that the immune system is in a suppressed state. Also, we found a positive correlation between CAFGs score and T cell dysfunction, T cell exclusion, and TIDE scores. Strikingly, there was a higher proportion of non-responders within the high CAFGs score group, and patients exhibiting high CAFGs score had a significantly reduced overall survival (OS) rate following their receipt of ICB therapy.

Moreover, single-cell transcriptomics and proteomics analyses have been employed to validate the significance of scoring system related molecules in the identity and function of CAF, including FNDC1, FRMD6, FBN1, RAB31, GLT8D2, COL1A2, GLIS2, COL8A1, GPC6, COL3A1, PRICKLE1, FSTL1, HLX, IGFBP7, and EFS. Consistent with the proteomic data, a previous study established that FSTL1, secreted by activated fibroblasts, promotes hepatocellular carcinoma metastasis and stemness³². CAF expressing IGFBP7 induce colony formation when co-culturing with CRC cells through paracrine tumor-stromal interaction³³. These results suggest that the 15 marker models we identified might play crucial roles in CAF-mediated CRC progression.

Overall, these findings signify the critical role that CAFs play in tumor immune phenotypes and response to ICB therapy, thereby evidencing the potential value of assessing CAFGs as a prognostic tool for those undergoing cancer immunotherapy.

4.1 Data source and process

The combined GEO cohort (1175 samples) used in this study were integrated by GSE39582, GSE14333, GSE17536, GSE17537, and GSE72968 cohorts. The transcriptome data of above five cohorts was the microarray data from GPL570 platform. The method of merging multiple data sets and the procedure of removing batch effects were carried out as reported in our previous study³⁴.The transcriptome data (FPKM) and clinical information of the TCGA COAD cohort were downloaded from the UCSC website³⁵. Then FPKM was converted to transcripts per kilobase million (TPM) and further log-2 transformed in next analysis.

4.2 Machine Learning

The R package ConsensusClusterPlus was applied for cluster the combined GEO cohort based on the input genes³⁶. In order to make the clustering result robust, we set the following parameters: 80% item resampling (pItem), 100% gene resampling (pFeature), a maximum evaluated k of 9 (maxK), 1000 resamplings (reps), and pam clustering algorithm (clusterAlg) upon spearman distances (distance).

4.3 Evaluation of immunological characteristics

R package MCPcounter³⁷ were applied to infer the abundance of immune cells infiltrating in the TME using the transcriptome data. Additionally, the adaptive and innate immune scores of patients were also calculated though ssGSEA algorithm in GSVA package, and the parameters were set as follows: method='ssgsea', KCDF ='Gaussian'. And 122 immunomodulators (Supplementary Data 9), including major histocompatibility complex (MHC), receptors, chemokines, and immunostimulants; and several common immune checkpoints with therapeutic potential were collected from previous studies^38–40.

4.4 Weighted correlation network analysis (WGCNA)

Weighted correlation network analysis (WGCNA) analysis enables to identify gene modules most associated with traits ⁴¹. In this study, the CAF clusters 1 and 2 were used as the traits. An appropriate soft threshold β (β = 4 in this study) was used to meet the criteria for the scale-free network. In next steps, WGCNA analysis was performed with default parameters.

4.5 Functional enrichment analyses

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) ⁴² analyses were applied to explore the biological functions of the blue modules in WGCNA using the R package “clusterprofiler”⁴³. Additionally, GSEA and GSVA analyses were performed using the Hallmark gene sets from MSigDB website with default parameters. An adjusted P-value of less than 0.05 was regarded as a statistically significant difference.

4.6 Construction and Validation of CAFGs scoring system

First, 100 CAFGs meeting GS > 0.2 and MM > 0.8 in WGCNA analysis was applied to univariate Cox regression. Then, the genes with P-value less than 0.2 in univariate Cox regression were considered as candidates, and inputted to multivariate Cox regression, which finally identified 15 genes (FNDC1, FRMD6, FBN1, RAB31, GLT8D2, COL1A2, GLIS2, COL8A1, GPC6, COL3A1, PRICKLE1, FSTL1, HLX, IGFBP7, and EFS) with a P-value less than 0.05 in multivariate Cox regression. Next, the CAFGs scoring system was constructed based on the 15 genes and corresponding regression coefficients in multivariate Cox regression, as follows:

$$\varvec{C}\varvec{A}\varvec{F}\varvec{G}\varvec{s} \varvec{s}\varvec{c}\varvec{o}\varvec{r}\varvec{e}=\sum _{\varvec{i}}\varvec{C}\varvec{o}\varvec{e}\varvec{f}\varvec{f}\varvec{i}\varvec{c}\varvec{i}\varvec{e}\varvec{n}\varvec{t} \varvec{o}\varvec{f} \left(\varvec{i}\right)\times \varvec{E}\varvec{x}\varvec{p}\varvec{r}\varvec{e}\varvec{s}\varvec{s}\varvec{i}\varvec{o}\varvec{n} \varvec{o}\varvec{f} \varvec{g}\varvec{e}\varvec{n}\varvec{e} \left(\varvec{i}\right)$$

The regression coefficient of the gene was designated (i) in the multivariate Cox proportional hazards regression.

4.7 Survival analysis

A total of 864 samples in the combined GEO cohort, 435 samples in the TCGA COAD cohort, have overall survival (OS) data (Supplementary Data 11). In addition, the recurrence-free survival (RFS) data in GSE39582; disease-free survival (DFS) data in GSE17536 and GSE17537; and disease-specific survival (DSS) data in GSE17536 were summarized in Supplementary Data 11, which were used to validate prognostic power of the CAFGs scoring system. The survival time was converted to months format, and patients with survival time less than one month were removed during survival analysis. Based on the optimal cutoff value identified by the survminer package, the patients were divided into high and low CAFGs score groups. Log-rank test was used to evaluate statistical differences. Kaplan–Meier (KM) plots were visualized using the survminer package.

4.8 Inferring the Consensus Molecular Subtypes (CMS) classification

The Consensus Molecular Subtypes (CMS), a widely used classification system currently available for CRC, has strong prognostic implications in clinical application². There are four subtypes of CMS, including CMS1 (MSI Immune), CMS2 (Canonical), CMS3 (Metabolic), and CMS4 (Mesenchymal). Among them, CMS4 is characterized by prominent transforming growth factor β (TGF-β) activation, stromal invasion, angiogenesis, poor OS and RFS. In this study, we inferred CMS classification based on transcriptome data using R package CMScaller with the default parameter⁴⁴.

4.9 ICB response prediction

A predictive algorithm known as the Tumor Immune Dysfunction and Exclusion (TIDE) algorithm was utilized to forecast the response to immune checkpoint blockade (ICB) by analyzing the gene expression profiles related to T cell dysfunction (Dysfunction) and exclusion (Exclusion). A lower TIDE score indicates a more favorable immunotherapy response. The scores of T cell Dysfunction, T cell Exclusion, and TIDE were obtained from the TIDE website. The IMvigor210 cohort, a vast population of patients with metastatic urothelial cancer receiving anti–PD-L1 therapy (atezolizumab), was downloaded from the Creative Commons 3.0 license. GSE78220 is a cohort of pre-treatment melanomas receiving anti–PD-1 therapy.

4.10 Statistical analysis

All analyses were performed in R 4.0.3. The Wilcox test was used to test the difference between two groups. The log-rank and Pearson test were used in KM survival and correlation analyses, respectively. In present study, heatmaps were visualized with the ComplexHeatmap package ⁴⁵. The ggplot2 package was used to visualize boxplots, scatter plots, and sankey plots. *, **, ***, and **** represent a P-value less than 0.05, 0.01, 0.001, and 0.0001, respectively.

Data Availability Statement

The original data presented in the study can be downloaded from GEO and TCGA websites.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a conflict of interest.

Consent to participate

All authors have read and approved the final manuscript for publication. We confirm that the work described has not been published previously, it is not under consideration for publication elsewhere, and publication has been approved by all authors and relevant authorities at the institution(s) where the work was carried out. All authors agree to participate as co-authors and are accountable for the authorship, accuracy and integrity of the work.

Author Contributions

FW, and ZLL conceived, designed, and guided the study, and provided financial support. XHW TLX designed the study, analyzed the bioinformatic data, and wrote the manuscript draft. ZLL and FW assisted in generating the figures and tables and article structure design. QZ, SJL and TYM assisted in bioinformatics analysis and article structure design. FW, ZLL, and XHW assisted in revising the manuscript and figures. QZ and SJL assisted in collecting and collating data of public CRC cohorts.

Funding

This work was supported by the Beijing Medical Administration Bureau Cultivation Program Project [No. PX2022041].

Acknowledgments

We would like to thank Prof. Xiaohui Wang, Department of General Surgery, Xuanwu Hospital, Capital Medical University, Beijing, China, who helped us complete this study.

Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA: a cancer journal for clinicians 2021; 71 (3): 209-249. doi: 10.3322/caac.21660.
Guinney J, Dienstmann R, Wang X, de Reyniès A, Schlicker A, Soneson C, et al. The consensus molecular subtypes of colorectal cancer. Nature medicine 2015; 21 (11): 1350-1356. doi: 10.1038/nm.3967.
Kobayashi H, Gieniec KA, Lannagan TRM, Wang T, Asai N, Mizutani Y, et al. The Origin and Contribution of Cancer-Associated Fibroblasts in Colorectal Carcinogenesis. Gastroenterology 2022; 162 (3): 890-906. doi: 10.1053/j.gastro.2021.11.037.
Kasashima H, Duran A, Martinez-Ordoñez A, Nakanishi Y, Kinoshita H, Linares JF, et al. Stromal SOX2 Upregulation Promotes Tumorigenesis through the Generation of a SFRP1/2-Expressing Cancer-Associated Fibroblast Population. Developmental cell 2021; 56 (1): 95-110.e110. doi: 10.1016/j.devcel.2020.10.014.
Zhu HF, Zhang XH, Gu CS, Zhong Y, Long T, Ma YD, et al. Cancer-associated fibroblasts promote colorectal cancer progression by secreting CLEC3B. Cancer biology & therapy 2019; 20 (7): 967-978. doi: 10.1080/15384047.2019.1591122.
Unterleuthner D, Neuhold P, Schwarz K, Janker L, Neuditschko B, Nivarthi H, et al. Cancer-associated fibroblast-derived WNT2 increases tumor angiogenesis in colon cancer. Angiogenesis 2020; 23 (2): 159-177. doi: 10.1007/s10456-019-09688-8.
Pape J, Magdeldin T, Stamati K, Nyga A, Loizidou M, Emberton M, et al. Cancer-associated fibroblasts mediate cancer progression and remodel the tumouroid stroma. Br J Cancer 2020; 123 (7): 1178-1190. doi: 10.1038/s41416-020-0973-9.
Li Z, Zhou J, Zhang J, Li S, Wang H, Du J. Cancer-associated fibroblasts promote PD-L1 expression in mice cancer cells via secreting CXCL5. International journal of cancer 2019; 145 (7): 1946-1957. doi: 10.1002/ijc.32278.
Hu JL, Wang W, Lan XL, Zeng ZC, Liang YS, Yan YR, et al. CAFs secreted exosomes promote metastasis and chemotherapy resistance by enhancing cell stemness and epithelial-mesenchymal transition in colorectal cancer. Molecular cancer 2019; 18 (1): 91. doi: 10.1186/s12943-019-1019-x.
Öhlund D, Handly-Santana A, Biffi G, Elyada E, Almeida AS, Ponz-Sarvise M, et al. Distinct populations of inflammatory fibroblasts and myofibroblasts in pancreatic cancer. The Journal of experimental medicine 2017; 214 (3): 579-596. doi: 10.1084/jem.20162024.
Mosa MH, Michels BE, Menche C, Nicolas AM, Darvishi T, Greten FR, et al. A Wnt-Induced Phenotypic Switch in Cancer-Associated Fibroblasts Inhibits EMT in Colorectal Cancer. Cancer research 2020; 80 (24): 5569-5582. doi: 10.1158/0008-5472.Can-20-0263.
Nicolas AM, Pesic M, Engel E, Ziegler PK, Diefenhardt M, Kennel KB, et al. Inflammatory fibroblasts mediate resistance to neoadjuvant therapy in rectal cancer. Cancer cell 2022; 40 (2): 168-184.e113. doi: 10.1016/j.ccell.2022.01.004.
Garvey CM, Lau R, Sanchez A, Sun RX, Fong EJ, Doche ME, et al. Anti-EGFR Therapy Induces EGF Secretion by Cancer-Associated Fibroblasts to Confer Colorectal Cancer Chemoresistance. Cancers (Basel) 2020; 12 (6): doi: 10.3390/cancers12061393.
Huang TX, Tan XY, Huang HS, Li YT, Liu BL, Liu KS, et al. Targeting cancer-associated fibroblast-secreted WNT2 restores dendritic cell-mediated antitumour immunity. Gut 2022; 71 (2): 333-344. doi: 10.1136/gutjnl-2020-322924.
Hirashima T, Karasawa H, Aizawa T, Suzuki T, Yamamura A, Suzuki H, et al. Wnt5a in cancer-associated fibroblasts promotes colorectal cancer progression. Biochemical and biophysical research communications 2021; 568 37-42. doi: 10.1016/j.bbrc.2021.06.062.
Cho C, Mukherjee R, Peck AR, Sun Y, McBrearty N, Katlinski KV, et al. Cancer-associated fibroblasts downregulate type I interferon receptor to stimulate intratumoral stromagenesis. Oncogene 2020; 39 (38): 6129-6137. doi: 10.1038/s41388-020-01424-7.
Franzè E, Di Grazia A, Sica GS, Biancone L, Laudisi F, Monteleone G. Interleukin-34 Enhances the Tumor Promoting Function of Colorectal Cancer-Associated Fibroblasts. Cancers (Basel) 2020; 12 (12): doi: 10.3390/cancers12123537.
Si G, Li S, Zheng Q, Zhu S, Zhou C. miR-1246 shuttling from fibroblasts promotes colorectal cancer cell migration. Neoplasma 2021; 68 (2): 317-324. doi: 10.4149/neo_2020_200924N1018.
Yuan Z, Hu H, Zhu Y, Zhang W, Fang Q, Qiao T, et al. Colorectal cancer cell intrinsic fibroblast activation protein alpha binds to Enolase1 and activates NF-κB pathway to promote metastasis. Cell death & disease 2021; 12 (6): 543. doi: 10.1038/s41419-021-03823-4.
Chen H, Cheng Y, Wang X, Wang J, Shi X, Li X, et al. 3D printed in vitro tumor tissue model of colorectal cancer. Theranostics 2020; 10 (26): 12127-12143. doi: 10.7150/thno.52450.
Luo X, Fong ELS, Zhu C, Lin QXX, Xiong M, Li A, et al. Hydrogel-based colorectal cancer organoid co-culture models. Acta biomaterialia 2021; 132 461-472. doi: 10.1016/j.actbio.2020.12.037.
Naruse M, Ochiai M, Sekine S, Taniguchi H, Yoshida T, Ichikawa H, et al. Re-expression of REG family and DUOXs genes in CRC organoids by co-culturing with CAFs. Scientific reports 2021; 11 (1): 2077. doi: 10.1038/s41598-021-81475-2.
Fernando-Macías E, Fernández-García MT, García-Pérez E, Porrero Guerrero B, López-Arévalo C, Rodríguez-Uría R, et al. A new aggressive xenograft model of human colon cancer using cancer-associated fibroblasts. PeerJ 2020; 8 e9045. doi: 10.7717/peerj.9045.
Wang Y, Liang Y, Xu H, Zhang X, Mao T, Cui J, et al. Single-cell analysis of pancreatic ductal adenocarcinoma identifies a novel fibroblast subtype associated with poor prognosis but better immunotherapy response. Cell discovery 2021; 7 (1): 36. doi: 10.1038/s41421-021-00271-4.
Herrera M, Berral-González A, López-Cade I, Galindo-Pumariño C, Bueno-Fortes S, Martín-Merino M, et al. Cancer-associated fibroblast-derived gene signatures determine prognosis in colon cancer patients. Molecular cancer 2021; 20 (1): 73. doi: 10.1186/s12943-021-01367-x.
Han C, Liu T, Yin R. Biomarkers for cancer-associated fibroblasts. Biomarker research 2020; 8 (1): 64. doi: 10.1186/s40364-020-00245-w.
Gascard P, Tlsty TD. Carcinoma-associated fibroblasts: orchestrating the composition of malignancy. Genes & development 2016; 30 (9): 1002-1019. doi: 10.1101/gad.279737.116.
Kalluri R. The biology and function of fibroblasts in cancer. Nature reviews Cancer 2016; 16 (9): 582-598. doi: 10.1038/nrc.2016.73.
Öhlund D, Elyada E, Tuveson D. Fibroblast heterogeneity in the cancer wound. The Journal of experimental medicine 2014; 211 (8): 1503-1523. doi: 10.1084/jem.20140692.
Sahai E, Astsaturov I, Cukierman E, DeNardo DG, Egeblad M, Evans RM, et al. A framework for advancing our understanding of cancer-associated fibroblasts. Nature reviews Cancer 2020; 20 (3): 174-186. doi: 10.1038/s41568-019-0238-1.
Togo S, Polanska UM, Horimoto Y, Orimo A. Carcinoma-associated fibroblasts are a promising therapeutic target. Cancers (Basel) 2013; 5 (1): 149-169. doi: 10.3390/cancers5010149.
Loh JJ, Li TW, Zhou L, Wong TL, Liu X, Ma VWS, et al. FSTL1 Secreted by Activated Fibroblasts Promotes Hepatocellular Carcinoma Metastasis and Stemness. Cancer research 2021; 81 (22): 5692-5705. doi: 10.1158/0008-5472.Can-20-4226.
Rupp C, Scherzer M, Rudisch A, Unger C, Haslinger C, Schweifer N, et al. IGFBP7, a novel tumor stroma marker, with growth-promoting effects in colon cancer through a paracrine tumor-stroma interaction. Oncogene 2015; 34 (7): 815-825. doi: 10.1038/onc.2014.18.
Wang H, Li Z, Ou S, Song Y, Luo K, Guan Z, et al. Tumor Microenvironment Heterogeneity-Based Score System Predicts Clinical Prognosis and Response to Immune Checkpoint Blockade in Multiple Colorectal Cancer Cohorts. Front Mol Biosci 2022; 9 884839. doi: 10.3389/fmolb.2022.884839.
Navarro Gonzalez J, Zweig AS, Speir ML, Schmelter D, Rosenbloom KR, Raney BJ, et al. The UCSC Genome Browser database: 2021 update. Nucleic acids research 2021; 49 (D1): D1046-d1057. doi: 10.1093/nar/gkaa1070.
Wilkerson MD, Hayes DN. ConsensusClusterPlus: a class discovery tool with confidence assessments and item tracking. Bioinformatics (Oxford, England) 2010; 26 (12): 1572-1573. doi: 10.1093/bioinformatics/btq170.
Becht E, Giraldo NA, Lacroix L, Buttard B, Elarouci N, Petitprez F, et al. Estimating the population abundance of tissue-infiltrating immune and stromal cell populations using gene expression. Genome biology 2016; 17 (1): 218. doi: 10.1186/s13059-016-1070-5.
Charoentong P, Finotello F, Angelova M, Mayer C, Efremova M, Rieder D, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell reports 2017; 18 (1): 248-262. doi: 10.1016/j.celrep.2016.12.019.
Auslander N, Zhang G, Lee JS, Frederick DT, Miao B, Moll T, et al. Robust prediction of response to immune checkpoint blockade therapy in metastatic melanoma. Nature medicine 2018; 24 (10): 1545-1549. doi: 10.1038/s41591-018-0157-9.
Wang H, Luo K, Guan Z, Li Z, Xiang J, Ou S, et al. Identification of the Crucial Role of CCL22 in F. nucleatum-Related Colorectal Tumorigenesis that Correlates With Tumor Microenvironment and Immune Checkpoint Therapy. Frontiers in genetics 2022; 13 811900. doi: 10.3389/fgene.2022.811900.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC bioinformatics 2008; 9 559. doi: 10.1186/1471-2105-9-559.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic acids research 2000; 28 (1): 27-30. doi: 10.1093/nar/28.1.27.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics : a journal of integrative biology 2012; 16 (5): 284-287. doi: 10.1089/omi.2011.0118.
Eide PW, Bruun J, Lothe RA, Sveen A. CMScaller: an R package for consensus molecular subtyping of colorectal cancer pre-clinical models. Scientific reports 2017; 7 (1): 16618. doi: 10.1038/s41598-017-16747-x.
Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics (Oxford, England) 2016; 32 (18): 2847-2849. doi: 10.1093/bioinformatics/btw313.

No competing interests reported.

SupplementaryFigure.zip
SupplementaryData1596CAFgenes.csv
Supplementary Data1: 596 differentially expressed genesbetween CAF and NF in CRC from a previous study.
SupplementaryData2115prognosticgenes.csv
Supplementary Data2: 115 potential prognostic CAF genesused as the input genes of unsupervised clustering.
SupplementaryData3GOresultsofbluemodule.csv
Supplementary Data3: GO terms enriched in blue module.
SupplementaryData4KEGGresultsofbluemodule.csv
Supplementary Data4: KEGG terms enriched in blue module.
SupplementaryData543genesmeetingunicox0.01.csv
Supplementary Data5: 43 genes meeting unicox p value less than 0.01 used as the input genes of unsupervised clustering.
SupplementaryData6modelgenesandcoefficients.csv
Supplementary Data6: The expression values of these CAFGs model genes and their corresponding regression coefficients.
SupplementaryData7Patientsgroupinmultiplecohorts.xlsx
Supplementary Data7: High and low CAFGs score groups in multiple CRC cohorts.
SupplementaryData8CMSofmultiplecohorts.xlsx
Supplementary Data8: CMS subtypes in multiple CRC cohorts.
SupplementaryData9122immunemolecules.xlsx
Supplementary Data9: 122 immunomodulators collected from previous studies.
SupplementaryData10survivaldata.xlsx
Supplementary Data10: Survival data of multiple CRC cohorts.
SupplementalFiguresandlegends.docx

Download PDF

Journal Publication

published 13 Mar, 2024

Read the published version in Journal of Cancer Research and Clinical Oncology →

Editorial decision: Major revision
17 Oct, 2023
Reviews received at journal
08 Oct, 2023
Reviewers agreed at journal
04 Oct, 2023
Reviewers invited by journal
03 Oct, 2023
Submission checks completed at journal
20 Sep, 2023
Editor assigned by journal
20 Sep, 2023
First submitted to journal
09 Sep, 2023

You are reading this latest preprint version

A Comprehensive Multi-Omics Analysis Identifies a Robust Scoring System for Cancer-Associated Fibroblasts and Intervention Targets in Colorectal Cancer

Status:

Journal Publication

Version 1

Abstract

Background

Methods

Results

Conclusion

Figures

1 Introduction

2 RESULTS

2.1 Construction of CAF clusters with different prognosis and immune states

2.2 Identification the key genes affecting the prognosis of patients in different CAF clusters

2.3 Colorectal cancer patients with high CAFGs score have poor outcomes in multiple colorectal cohorts

2.4 Patients with a high CAF score develop resistance to immunotherapy

2.5 Exploring potential molecules associated with the identity and function of CAF

3 DISCUSSION

4 MATERIALS AND METHODS

4.1 Data source and process

4.2 Machine Learning

4.3 Evaluation of immunological characteristics

4.4 Weighted correlation network analysis (WGCNA)

4.5 Functional enrichment analyses

4.6 Construction and Validation of CAFGs scoring system

4.7 Survival analysis

4.8 Inferring the Consensus Molecular Subtypes (CMS) classification

4.9 ICB response prediction

4.10 Statistical analysis

Declarations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1