Constructing a prognostic model for colon cancer patients on basis of coagulation genes enriched in cancer-associated fibroblasts to guide personalized immunotherapy

Data collection of bulk transcriptomes

The clinical data and RNA-seq data originated from the TCGA dataset. The TCGA-COAD dataset included 476 tumor tissues and 41 normal tissues. Additionally, to validate our model’s accuracy, we obtained expression data from GSE39582 [11] and corresponding platform annotation files from the Gene Expression Omnibus database, comprising 585 tumor samples. In the analysis of colon cancer patient prognosis, individuals with incomplete clinical data were excluded, leading to 452 patients in the TCGA-COAD dataset and 578 patients in the GSE39582 dataset. Detailed clinicopathological characteristics are offered in Table S1 of the Supplementary materials. The immunohistochemical pictures used in this study were provided by the online database The Human Protein Atlas (HPA) (https://www.proteinatlas.org/ENSG00000102265-TIMP1).

Construction and validation of coagulation-related signatures

The coagulation pathway data were offered by the KEGG database (https://www.genome.jp/kegg/), encompassing hsa04610 and hsa04611 [12], as well as HALLMARK_COAGULATION from the MsigDB database (https://www.gsea-msigdb.org/gsea/msigdb). There were 290 coagulation-associated genes (CRGs) included. CRGs related to prognosis were recognized with univariate Cox regression analysis. Later, a prognostic model on basis of LASSO regression was built with the R package “glmnet.” The coagulation-related risk score (CRRS) was derived from the model and adopted to fall patients into high- and low-risk groups. The production of Kaplan–Meier survival curves were made with the “gsurvplot” function in the R software package “survminer”. The TCGA-COAD cohort acted as the training set for internal validation, while the GSE39582 cohort was employed for external validation. The accuracy and predictive capability of CRRS were assessed via receiver operating characteristic curves (ROCs) using the “survivalROC” package. This study used 4.3.0 R software.

Building and evaluation of nomogram for colon patients

The “Survival” and “RMS” packages were employed to build a predictive nomogram for colon patients, utilizing variables from the TCGA-COAD training set. The prognostic efficiency of the predictive model was gauged with Harrell’s Concordance Index.

Functional enrichment analysis of DEGs

The “limma” R package was adopted to analyze differentially expressed genes between high- and low-risk groups. Gene set enrichment analysis on basis of DEGs and their logFCs was made with the R package “clusterProfiler” [13]. The consequent enrichment maps were visualized by applying R software. Additionally, hidden variations in pathway activity in every sample were assessed with GSVA, analyzing enrichment scores for pathways across all samples with the “GSVA” package in R [14].

Exploration of the tumor immune microenvironment

The calculation of immune scores and estimate scores for colon patients was made with the “estimate” R package [15]. We obtained the TCGA pan-cancer immune cell infiltration file (“infiltration_estimation_for_TCGA.csv”) from the TIMER 2.0 database (http://timer.cistrome.org/), enabling correlation analysis between immune cell infiltration and risk values.

Estimation of CAFs infiltration

Four methods were adopted to forecast the abundance of CAFs, such as the Estimation of Ratio of Immune and Cancer Cells (ERIC) algorithm [16], the xCell algorithm [17], the Microenvironmental Cell Population Counter (MCPcounter) algorithm [18], and the Tumor Immune Dysfunction and Rejection algorithm [19]. Spearman correlations were then analyzed between CRRS and stromal scores, and different CAFs infiltration level scores (EPIC, xCell, MCPcounter, and TIDE), validating the correlation between the model and CAFs infiltration.

Role of CRRS in predicting immunotherapy

To predict immunotherapy outcomes using immune checkpoint blockers (ICBs), we computed predictive scores such as TIDE scores and Immunophenotype scores (IPS). The gene expression profile of the TCGA-COAD cohort was uploaded to the Cancer Immunome Atlas, and IPS for each patient was downloaded [20]. The TIDE score, a new way to evaluate the efficacy of ICBs, is available on the TIDE website [19].

TIDE through two main mechanisms simulation of tumor immune escape, in cytotoxic T lymphocytes (CTL) to induce T cell dysfunction in high levels of tumor, as well as the prevention of T cells in the low level of CTL tumor infiltrating [19]. Tumor Inflammatory Signature includes genes related to antigen presentation, chemokine expression, cytotoxic activity, and adaptive immune resistance and has been retrospectively studied in clinical experiments [21]. We compared the prognostic power of CRRS, TIDE, and TIS by analyzing ROC curves for the three models.

Discussion single-cell sequencing data

The single-cell dataset GSE166555 [22] originated from the Tumor Immune Single-cell Hub (TISCH) database (http://tisch.comp-genomics.org/home/), encompassing 13 tumor tissues and 12 corresponding para-cancer tissues. Quality control procedures involved filtering and normalizing the original gene expression matrix using the “Seurat” R software package [23]. Subsequently, data dimensions were reduced using the Unified Move Approximation and Projection (UMAP) method, and the results were clustered and visualized using T-distributed Random Neighborhood Embedding (t-SNE) projection. Cells were annotated using standard cell surface markers.

Cell culture

We bought HCT116 and DLD1 cells from American Types Culture Collection, no mycoplasma infection. HCT116 cells were cultured in DMEM medium completely (KGL1206, KeyGEN BioTECH), DLD1 cell culture in the 1640 complete medium (KGL1501, KeyGEN BioTECH). All the complete media contained 10 % fetal bovine serum (FBS), 80 U/mL penicillin and 0.08 mg/mL streptomycin. Cells were cultured in an incubator at 37 °C containing 5 % CO₂.

Cell counting kit-8 (CCK-8) assay

Cell viability was estimated applying CCK-8 (KGA9305, KeyGEN BioTECH, China) based on producer’s guidance. Transfected and untransfected HCT116 and DLD1 cells were inoculated in 96-well plates with a density of 1×10⁴/well. Then CCK-8 reagent was added at different times (0 h, 24 h, 48 h, 72 h) to measure the OD value under the wavelength of 450 nm (Infinite 200 PRO, Tecan, Switzerland) and evaluate cell survival rate.

Cell transfection

The cloning of TIMP1 cDNAs into Lentiviral pLVX plasmid (SyngenTech, Beijing, China) was made, next to verification by sequencing. Then, Lipofectamine 3000 (L3000015, Invitrogen, USA) was adopted to transfect the pLVX-TIMP1 plasmids into HCT116 and DLD1 cells. CRC cells overexpressing TIMP1 were selected after 7 days of treatment with 10 μg/mL purinomycin (ST551, Beyotime, China).

Colony formation assay

The seeding of 1,000 cells per well into 6-well plates was made. The medium was changed every three days, next to 2-week incubation. With 4 % paraformaldehyde fixed colonies for 10 min, 15 min and then use 0.1 % crystal violet stain (C0121 Beyotime, China), in the end, rinsed in running water clean, air-dried under send pictures on the board.

Transwell assay

Transwell chamber was adopted to make the migration and invasion assays. For migration assay, the seeding of transfected and untransfected cells into the upper chamber (3472, Corning, USA) with serum-free medium (KeyGEN BioTECH, China) and the bottom of the chamber containing medium containing 20 % serum was made. For the invasion assay, Matrigel (356234, Corning, USA) was used to coat the chamber, and the following steps were like the migration assay. After 48 h incubation, the fixture of cells was made with 4 % paraformaldehyde, next to staining with 0.1 % crystal violet, and photographing.

Western blot

Proteins of colon cancer cells with and without TIMP1 overexpression were extracted using RIPA lysis buffer (P0013B, Beyotime, China). Equivalent amounts of proteins were isolated by sodium dodecyl sulfo-polyacrylamide gel electrophoresis in a 12 % gel, next to a PVDF membrane. Blocking was performed with 5 % skim milk powder (P0216, Beyotime, China) for 2 h, followed by incubation of membranes with TIMP1 antibody (1:1,500, 16644-1-AP, Proteintech, USA) overnight at 4 °C, next to incubation with secondary antibody (BL003A, Biosharp, China) the following day after three washes with TBST (G0001-2L, Servicebio, China). Resistance to GAPDH antibodies (1:10000, 10494-1-AP, Proteintech, USA) as a control. Finally, the stripes were visualized using a chemiluminescence instrument (Tanon 5200, China).

Statistical analysis

R software was adopted to make all statistical analyses. Results were of statistical significance in case of p-values of below 0.05 (two-sided). Kaplan–Meier analysis assessed and compared survival between high- and low-risk groups. Univariate Cox regression analysis screened genes related to prognosis, and independent risk elements influencing colon patients were identified after factors with p values <0.05 underwent multivariate Cox regression analysis. Cell counts were performed using ImageJ (1.54g), and the corresponding statistical plots were counted and plotted by GraphPad Prism 8. T-test or Wilcoxon rank sum test was adopted for comparing the two groups. Correlation analysis of continuous variables was made with Spearman’s test.

Construction of coagulation-related signatures for colon cancer prognosis prediction

We initiated our analysis with a univariate Cox regression assessment to identify CRGs within the TCGA-COAD cohort. This scrutiny led to the selection of 23 prognostic genes exhibiting p-values <0.05 (Figure 1A). Later, LASSO Cox regression analysis was made on these 23 genes to construct a multi-gene prognosis model that minimized overfitting risk. The optimal value of λ yielded the inclusion of eight genes in our model (Figure 1B). The CRRS was then computed using the formula: CRRS=(0.447229484560575×PLCB1 expression)+(0.362128614897294×RASGRP2 expression)+(0.0457263309301581×PLCG2 expression)+(0.350092455079463×SERPINA5 expression)+(−0.443300975035417×F2RL2 expression)+(0.283241161078688×DUSP14 expression)+(−0.167780688007362×MMP10 expression)+(0.225325634783633×TIMP1 expression). With the median value as a limit, the classification of patients into high-risk and low-risk groups was made. Figure 1C shows risk curves, survival status, and heatmap of individual gene expression.

Figure 1:

Construction of coagulation-related signatures for colon cancer prognosis prediction. (A) The forest plots of prognostic genes with p <0.05 were recognized by univariate Cox regression analysis. (B) The core genes were screened with LASSO Cox regression analysis. (C) The allocation of the risk score for every patient (top panel), survival situation of patients (middle panel), heatmaps for eight-gene signature between the high-risk group and low-risk group (bottom panel).

Validation of the model in training and verification sets

We proceeded to assess the prognostic capacity of our model in both the TCGA-COAD training set and the verification set (GSE39582). Survival analysis comparing low-risk and high-risk patients in the TCGA-COAD training set revealed greatly reduced overall survival (OS) among colon cancer patients in the high-risk group, a trend consistent with the validation set (Figure 2A and B). Receiver operating characteristic curves were adopted to gauge the sensitivity and specificity of our coagulation-related signatures in predicting colon cancer patient prognosis (Figure 2C and D). According to the outcomes, the areas under the curve for 1, 3, and 5-year OS predictions in the TCGA-COAD training set were 0.708, 0.713, and 0.701, while in the verification set, these AUCs were 0.575, 0.586, and 0.570. In summary, our 8-gene risk model related to coagulation effectively distinguishes high-risk from low-risk patients and demonstrates robust prognostic capabilities.

Figure 2:

Verification of the model in training and verification sets. (A and B) Kaplan–Meier survival curve of OS between high-risk scores and low-risk scores in the TCGA training cohort and GSE39582 verification cohort. (C and D) The ROC curves at 1-, 3- and 5-year in the training and verification cohorts.

Construction of nomogram for colon cancer patients

To ascertain if CRRS qualifies as an independent clinical prognostic factor for colon cancer patients, univariate and multivariate Cox regression analyses were made. Age, sex, stage, TNM stage, and CRRS were evaluated as covariates, revealing that CRRS emerged as an independent predicting agent of colon cancer patient prognosis (Figure 3A). To improve the prediction of colon cancer patient survival, a nomogram amalgamating CRRS and clinical determinants were crafted based on the TCGA-COAD cohort (Figure 3B). After assigning scores to each prognostic parameter for patients, we calculated the cumulative survival rate, with higher total scores indicating a worse prognosis. Calibration plots portrayed 1-, 3-, and 5-year survival forecasts for colon cancer patients in the training set (Figure 3C). Furthermore, CRRS demonstrated significant associations with tumor stage and TNM stage, as confirmed by correlation analysis in both the training and verification sets (Figure 3D and E).

Figure 3:

Construction of nomogram for colon cancer patients. (A) Univariate analysis and multivariate analysis including CRRS and clinical elements. (B) Nomogram integrating risk scores and clinical factors. (C) Calibration of the nomogram at 1-, 3- and 5-year survival rate in the training cohort. (D) Heat map of risk scores and clinical relevance in the TCGA-COAD set. (E) Heat map of risk scores and clinical relevance in the GSE39582 set (^*p<0.05; ^**p<0.01; ^***p<0.001).

Enrichment analysis of risk score-related functions

To delve deeper into the underlying mechanisms associated with CRRS, we scrutinized enriched KEGG gene sets within the high-risk and low-risk groups. As displayed in Figure 4A, genes in the high-risk group showed enrichment in pathways such as cell adhesion molecules (cams), dilated cardiomyopathy, ECM receptor interaction, focal adhesion, and systemic lupus erythematosus. Conversely, the enrichment of genes in the low-risk group was made in pathways like aminoacyl tRNA biosynthesis, butanoate metabolism, citrate cycle (TCA cycle), peroxisome, and retinol metabolism (Figure 4B). Additionally, GSVA analysis highlighted significant correlations between CRRS and pathways including the TGF-β signaling pathway, MAPK signaling pathway, and calcium signaling pathway (Figure 4C).

Figure 4:

Enrichment analysis of risk score-related functions. (A) Enrichment pathways in the KEGG dataset of the high-risk group. (B) Enrichment pathways in the KEGG dataset of low-risk group. (C) Heat map of GSVA analysis CRRS and signaling pathway correlation (^*p<0.05; ^**p<0.01; ^***p<0.001).

Tumor microenvironment analysis in high- and low-risk groups

A contrast of the tumor microenvironment (TME) features between high- and low-risk groups was made. The ESTIMATE algorithm provided stromal marks, immune marks, and estimate marks to assess cancer-related stroma, immune cell infiltration, and tumor purity, respectively. The results indicated that the high-risk group exhibited greatly higher stromal marks, immune marks, and estimate marks by comparing with the low-risk group (Figure 5A). To ensure reliability, we also analyzed the correlation between immune cells and CRRS applying various algorithms, including xCell, TIMER, QUANTISEQ, MCPcounter, EPIC, CIBERSORT-ABS, and CIBERSORT. The findings demonstrated a strong positive correlation between CRRS and CAFs predicted by the MCPcounter and EPIC algorithms (Figure 5B and C). Spearman correlation analysis confirmed the robustness of the CRRS model as a CAFs infiltration predictor. CRRS was positively correlated with stromal scores and CAFs infiltration scores estimated by four algorithms (xCell, MCPcounter, EPIC, and TIDE) in both the TCGA-COAD and GSE39582 cohorts (Figure 5D). These results collectively affirmed the reliability of the risk score as an index of CAFs infiltration. Furthermore, the MCPcounter and EPIC algorithms were employed to score CAFs in high- and low-risk patients from the TCGA-COAD and GSE39582 cohorts, underscoring stronger CAFs infiltration in the high-risk group (Figure 5E and F).

Figure 5:

Tumor microenvironment analysis in high- and low-risk groups. (A) Stromal mark, immune mark and estimate mark between high- and low-risk groups in TCGA-COAD and GSE39582. (B) Bubble chart of correlation between immune cells and CRRS. (C) Scatterplot of MCPcounter and EPIC algorithms to assess the correlation between CAFs infiltration and CRRS. (D) According to Spearman’s correlation analysis, CRRS was positively related to stromal scores and multiple estimates of CAFs infiltration in the TCGA-COAD and GSE39582 cohorts. (E) EPIC algorithm and MCPcounter algorithm score high- and low-risk groups in the TCGA-COAD for CAFs. (F) EPIC algorithm and MCPcounter algorithm score high- and low-risk groups in the GSE39582 for CAFs (^*p<0.05; ^***p<0.001).

Comparison of immunotherapy in high- and low-risk groups

Immunotherapy has gained prominence as a cornerstone in the treatment of different solid cancers, such as colon cancer. Notably, the FDA accelerated the approval of pembrolizumab and nivolumab, two programmed cell death 1 (PD1) blocking antibodies, for mCRC patients with mismatch repair flaws and high microsatellite unsteadiness (dMMR-MSI-H), demonstrating promising efficacy [24]. To investigate whether CRRS correlates significantly with immunotherapy outcomes, we employed the TIDE method. In the TCGA-COAD cohort, a higher percentage of high-risk patients (67 %) displayed unresponsiveness to treatment compared to low-risk patients (43 %). A similar tendency was discovered in the GSE39582 cohort, where a greater number of high-risk patients (60 %) did not respond to treatment as opposed to low-risk patients (51 %). High-risk patients also exhibited higher TIDE scores (Figure 6A), collectively indicating a greater potential for immune evasion. Additionally, we evaluated the likelihood of patients benefiting from immune checkpoint inhibitors (ICIs) by analyzing Immunophenotype scores (IPS) in high- and low-risk groups within the TCGA-COAD cohort. The low-risk group appeared to possess a more immunogenic phenotype (Figure 6B). Furthermore, we assessed the association between risk scores and immune checkpoints in both the TCGA-COAD and GSE39582 cohorts, revealing positive associations between risk marks and the majority of immune checkpoints (Figure 6C). While TIDE and TIS have been employed to predict patient responses to immunotherapy, the area under the ROC curve of CRRS exceeded that of TIDE and TIS, signifying superior predictive value and establishing CRRS as a promising biomarker for immune-related prognosis (Figure 6D).

Figure 6:

Comparison of immunotherapy in high- and low-risk groups. (A) TIDE scoring and percentage of responders and non-responders in the high- and low-risk groups in the TCGA-COAD and GSE39582. (B) IPS scores for high- and low-risk groups in the TCGA-COAD. (C) Heatmap of CRRS and immune checkpoint correlation for the TCGA-COAD and GSE39582 cohorts. (D) ROC analysis of CRRS, TIDE, and TIS on OS at 1-, 2-, and 3-year follow-up in the TCGA-COAD cohort (^*p<0.05; ^**p<0.01; ^***p<0.001).

Identification of TIMP1⁺ CAF as a key prognostic component through scRNA transcriptomic analysis

Given the close connection between CAFs and CRRS within the immune microenvironment, our focus turned to express model genes in CAF cells applying scRNA transcriptome analysis. After dimensionality reduction and cell clustering of GSE166555 single-cell data, 33 distinct cell clusters emerged, subsequently annotated into 13 cell types, including fibroblasts (Figure 7A–C). Importantly, the HALLMARK_COAGULATION gene set was predominantly expressed in fibroblasts (Figure 7D). Further correlation analysis between model gene expression and CAFs infiltration using the EPIC and MCPcounter algorithms identified PLCG2, RASGRP2, and TIMP1 as genes positively correlated with CAFs infiltration (Figure 7E). Notably, TIMP1 is abundantly expressed in fibroblasts (Figure 7F), and was notably overexpressed in CAFs of tumor tissues compared to normal fibroblasts of normal tissues (p<0.001; Figure 7G). Additionally, analysis of TIMP1 protein expression in relation to CAFs using IHC images from the HPA online database [25] further confirmed TIMP1 as a specific gene marker for CAFs, showing strong staining in the stromal regions of colon cancer tissues and no expression in normal tissues (Figure 7H).

Figure 7:

Identification of TIMP1⁺ CAF as a key prognostic component through scRNA transcriptomic analysis. (A and B) UMAP plot showing clusters of cells defined according to marker genes. (C) Bubble plot showing expression of labeled genes in all cell clusters. Dot size suggests the percentage of expressed genes and color suggests the level of expression intensity. (D) UMAP plot showing coagulation gene pathways enriched in fibroblast clusters. (E) Scatterplot of correlation between model gene expression and degree of CAFs infiltration. (F) The expression of RASGRP2, PLCG2 and TIMP1 on all cells. (G) Expression of TIMP1 in fibroblasts between normal and tumor tissues. (H) Protein expression of TIMP1 in normal colon tissue (left) and colon cancer specimens (right) from the Human Protein Atlas database. Scale bar=200 μm.

TIMP1 promotes the proliferation and metastasis of CRC cells

To research the role of TIMP1 in CRC, we upregulated TIMP1 in HCT116 and DLD1 cells (Figure 8A) and found that TIMP1 overexpression promoted CRC cell proliferation (Figure 8B), colony formation (Figure 8C), migration and invasion (Figure 8D). Collectively, these data showed that TIMP1 serves as an oncogene in CRC.

Figure 8:

TIMP1 drives the propagation and metastasis of CRC cells. (A) TIMP1 overexpression efficiency was determined by Western blot at the protein level in HCT116 and DLD1 cells. (B) CCK8 reagent was used to detect the propagation of HCT116 cells and DLD1 cells after TIMP1 overexpression. (C) The changes in the propagation ability of HCT116 cells and DLD1 cells were explored with colony formation assay after TIMP1 overexpression. (D) Transwell assay investigated the changes of migration and invasion ability of HCT116 cells and DLD1 cells after TIMP1 overexpression, scale bar=100 μm (^*p<0.05; ^**p<0.01; ^***p<0.001).