A CFH- and SPINT2-based prognostic signature for cholangiocarcinoma

Introduction

Cholangiocarcinoma (CCA), a highly aggressive malignancy originating from biliary epithelial cells, has an increasing global incidence, particularly in Asian regions (1,2). According to anatomical location, CCA is classified into intrahepatic, perihilar, and distal types, and it exhibits significant heterogeneity in etiology, molecular characteristics, clinical presentation, and prognosis (3,4). Due to a lack of early conspicuous symptoms and specific screening biomarkers, the majority of patients with CCA are diagnosed at an advanced stage, precluding curative surgical resection (5). Even after radical resection, the postoperative recurrence rate remains high, with a 5-year overall survival (OS) rate of less than 20% (6,7). Conventional chemotherapy, such as that consisting of gemcitabine plus cisplatin, provides only limited survival benefits and is associated with the development of chemoresistance (8,9). Therefore, the development of reliable prognostic biomarkers is of great significance for risk stratification and the personalized treatment of patients with CCA.

A prominent pathological feature of CCA is its dense desmoplastic stroma, characterized by an abundance of cancer-associated fibroblasts (CAFs), immune cells, and extracellular matrix components (10). The stromal components within the tumor microenvironment (TME) engage in complex bidirectional interactions with malignant cells, collectively driving tumor progression, immune evasion, and therapeutic resistance (11). Studies have demonstrated that CAFs play a central role in remodeling the TME and promoting tumor invasion and metastasis through the secretion of various cytokines, growth factors, and chemokines (12). In recent years, precision therapeutic strategies, mainly those with immune checkpoint inhibitors (ICIs), and targeted therapies have provided optimism for patients with CCA (13). However, due to the inherently immunosuppressive nature of the TME, the clinical efficacy of systemic therapies, including immunotherapies and targeted agents, remains limited in a substantial subset of patients with CCA (14). Therefore, the identification of reliable biomarkers for predicting therapeutic response and guiding treatment selection holds considerable clinical value for implementing precision medicine in patients with CCA (15).

Although studies on machine learning-based prognostic model have made progress in the field of CCA in recent years, they often have limited sample sizes, lack independent external validation, or rely on a single database (16,17). Large-scale, multicenter external validation is critical for ensuring model generalizability (18). To address these shortcomings, we analyzed multiple independent CCA datasets [GSE26566, GSE89749, E-MTAB-6389, and The Cancer Genome Atlas Cholangiocarcinoma (TCGA-CHOL) dataset], along with single-cell transcriptomic data (GSE138709) and a drug sensitivity database [Genomics of Drug Sensitivity in Cancer 2 database (GDSC2)], to construct and rigorously validate a dual-gene prognostic model. It is hoped the findings from this study can aid in the risk stratification and personalized treatment of patients with CCA. We present this article in accordance with the TRIPOD reporting checklist (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0798/rc).

Methods

This study used publicly available datasets from the Gene Expression Omnibus (GEO) database, ArrayExpress, and TCGA. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. Ethical approval and informed consent were obtained in the original studies; no additional institutional review board approval was required for this secondary analysis.

Data source and differential expression analysis

The CCA gene expression dataset GSE26566 (104 CCA samples and 59 adjacent normal samples; published in February 2012) was downloaded from GEO. Differential expression analysis was performed with the GEO2R online tool based on the limma package from Bioconductor in R (The R Foundation for Statistical Computing, Vienna, Austria). Differentially expressed genes (DEGs) with |log₂ fold change (FC)| >2 and adjusted P (adj. P) value <0.05 (Benjamini-Hochberg correction) were selected. The DEGs were intersected with the expression matrix of the GSE89749 training set to obtain candidate genes for subsequent analysis.

Training set data

The training set was the GSE89749 dataset (n=115; published July 2017), which contains complete gene expression profiles and survival information (OS time and survival status). Expression data were log₂-transformed and normalized. Survival status was defined as follows: 1= death and 0= censored. Sample size was determined according to data availability in the public repository; no a priori calculation was performed.

Univariate Cox regression screening

Univariate Cox regression analysis was performed for each candidate gene, and genes significantly associated with OS (P<0.05) were selected for subsequent machine learning screening.

Integrated screening and model construction

Three complementary algorithms were applied to the genes identified by univariate Cox regression using a consensus intersection strategy (equal contribution across algorithms; no weighted scoring).

• Least absolute shrinkage and selection operator (LASSO)-Cox regression: 10-fold cross-validation was used to select λ.min, and genes with non-zero coefficients were retained;
• Random forest: a model with 500 decision trees was constructed, and the top 10 genes were selected according to the mean decrease in the Gini index;
• Univariate Cox ranking: the 10 genes with the smallest P values were selected.

Only genes selected by all three methods were retained. A multivariate Cox regression model was then constructed using the overlapping genes to calculate the risk score. All gene expression values were treated as continuous variables. Patients were divided into high- and low-risk groups according to the median risk score.

External expression validation

Two independent datasets were used to validate the expression patterns of CFH and SPINT2: GSE76297 (90 CCA samples and 92 adjacent normal samples; published June 2017) and GSE107943 (30 CCA samples and 27 adjacent normal samples; published December 2018). Differential expression analysis was performed with the GEO2R online tool based on the limma package from Bioconductor in R. For validation, we assessed the direction of effect (up and down) and statistical significance (adj. P<0.05) without imposing an additional FC threshold.

Model evaluation

Kaplan-Meier survival curves and the log-rank test were used to assess survival differences between the high- and low-risk groups. Time-dependent receiver operating characteristic (timeROC) curves were used to evaluate the predictive accuracy of the model at 1, 2, and 3 years, and the area under the curve (AUC) was calculated. The concordance index (C-index) was also calculated to assess the discriminative ability of the model.

External validation

External validation was performed via Kaplan-Meier curves, log-rank tests, and timeROC curves with two independent cohorts: E-MTAB-6389 (n=75; published December 2019) and TCGA-CHOL (n=36; published 2016). Multivariate Cox regression adjusting for sex and vascular invasion was conducted in the E-MTAB-6389 cohort. To assess whether the combined prognostic signature outperformed its individual components, single-gene and combined models were compared in the training and validation cohorts using the C-index, Kaplan-Meier/log-rank test, and nested Cox regression (likelihood ratio test).

Functional enrichment analysis

Patients were divided into high- and low-risk groups based on the median risk score. DEGs between the two groups were screened with the limma package from Bioconductor in R using the following criteria: |log₂FC| >1 and adj. P<0.05. Hallmark gene set enrichment analysis (GSEA) was performed on the DEGs via the clusterProfiler package version 4.4.4 in R, and pathways with P<0.05 were considered significantly enriched.

Immune infiltration analysis

Single-sample GSEA (ssGSEA) was used to evaluate immune cell infiltration levels. Signature gene sets for 28 immune cell types, including B cells, CD4⁺ T cells, CD8⁺ T cells, dendritic cells, macrophages, and natural killer (NK) cells, were collected from a previous study (19). The ssGSEA algorithm implemented in the GSVA package from Bioconductor in R was used under default parameters to calculate enrichment scores for the 28 immune cell types in each sample. Patients were divided into high- and low-risk groups based on the median risk score, and the Wilcoxon rank-sum test was used to compare immune cell infiltration levels between the two groups. Spearman correlation analysis was used to assess the correlations of the model gene expression with the infiltration scores of each immune cell type.

Single-cell RNA-sequencing data analysis

The single-cell dataset GSE138709 was downloaded from the GEO database and analyzed via Seurat version 4.3.0 in R. Cells with a unique molecular identifier (UMI) count >2,000, a gene count of 500–6,000, and a mitochondrial gene percentage <20% were retained. Data were normalized via the LogNormalize function in the Seurat package, and the 2,000 most highly variable genes were selected for principal component analysis. t-distributed stochastic neighbor embedding and uniform manifold approximation and projection were applied based on the first 10 principal components, with clustering conducted at a resolution of 0.3. The following cell types were annotated with known marker genes: malignant cells (EPCAM, KRT19, and KRT7), fibroblasts (ACTA2 and COL1A2), T cells (CD2, CD3D, and CD3E), macrophages (CD14), B cells (MS4A1 and CD79A), endothelial cells (ENG and VWF), cholangiocytes (FXYD2, TM4SF4, and ANXA4), and hepatocytes (APOC3, FABP1, and APOA1). Fibroblasts were subsequently extracted from the annotated single-cell object for secondary subclustering, and CAF subtypes (myCAF, iCAF, and apCAF) were assigned by module scores of established marker genes. Pearson correlation was calculated between CFH expression and CAF subtype module scores at the single-cell level. Gene expression distribution [relative expression (arbitrary units)] was visualized with the FeaturePlot function in the Seurat package. Average expression levels were calculated according to cell type and tissue type. The Wilcoxon rank-sum test was used to compare expression differences between groups, with P<0.05 considered statistically significant.

Drug sensitivity analysis

Drug sensitivity analysis was performed with the GDSC2 database. The predicted half-maximal inhibitory concentration (IC₅₀) values of 198 compounds for each sample were estimated via the ridge regression model implemented in the oncoPredict version 0.2 package in R. The Wilcoxon rank-sum test was used to compare IC₅₀ values between high- and low-risk groups, and P values were adjusted via the Benjamini-Hochberg method, with an adj. P<0.05 [false-discovery rate (FDR) <0.05] being considered statistically significant. Observed IC₅₀ values of biliary tract cell lines were further retrieved from GDSC2 based on tissue annotation. The prognostic risk score was calculated for each cell line using the training-derived coefficients. Spearman correlation analysis was performed between the risk score and observed IC₅₀ across GDSC2 cell lines.

Statistical analysis

All statistical analyses were performed with R software version 4.2.3. LASSO-Cox analysis was performed with the glmnet package, random forest with the randomForest package, survival analysis with the survival and survminer packages, timeROC with the timeROC package, ssGSEA with the GSVA package, and functional enrichment analysis with the clusterProfiler package. A two-sided P<0.05 was considered statistically significant.

Results

Identification of DEGs

The overall study workflow is illustrated in Figure 1A. To identify CCA-related DEGs, we analyzed the GSE26566 dataset with the GEO2R online tool based on the R limma package from Bioconductor, comparing gene expression profiles between CCA tissues (n=104) and adjacent normal tissues (n=59). The criteria for screening DEGs were |log₂FC| >2 and adj. P<0.05 (Benjamini-Hochberg correction). A total of 265 significant DEGs were identified and are shown in the volcano plot in Figure 1B (red dots represent upregulated genes, and blue dots represent downregulated genes). These DEGs were used for subsequent prognostic model construction and functional analysis.

Figure 1 Study workflow and identification of differentially expressed genes in CCA. (A) Schematic diagram of the overall study workflow. (B) Volcano plot of DEGs from the GSE26566 dataset (104 CCA samples and 59 adjacent normal samples). Red dots: upregulated genes; blue dots: downregulated genes; gray dots: genes with no significant difference. Dashed lines indicate |log₂FC| =2 and adj. P=0.05. adj. P, adjusted P; CCA, cholangiocarcinoma; DEG, differentially expressed gene; FC, fold change; LASSO, least absolute shrinkage and selection operator; NS, not significant; ssGSEA, single-sample gene set enrichment analysis; TCGA-CHOL, The Cancer Genome Atlas Cholangiocarcinoma.

Construction and validation of the prognostic signature

The 265 DEGs were intersected with the expression matrix of the GSE89749 training set, yielding 236 candidate genes. Univariate Cox regression identified 21 genes significantly associated with OS, with the top 10 shown in Figure 2A. LASSO-Cox regression (λ.min =0.0348) further identified 8 genes (Figure 2B,2C), and random forest selected the top 10 genes (Figure 2D). The intersection of these three sets (LASSO-selected 8, top 10 from random forest, and top 10 univariate Cox regression) yielded two overlapping genes: CFH and SPINT2 (Figure 2E). External validation using independent datasets GSE76297 and GSE107943 confirmed that the expression trends of CFH and SPINT2 were consistent with those observed in GSE26566 (Figure 2F). A multivariate Cox regression model was then constructed based on these two genes, with the risk score calculated as follows: risk score =0.2594 × CFH + (−0.3215) × SPINT2. Patients were divided into high-risk (n=57) and low-risk (n=58) groups based on the median risk score. Kaplan-Meier survival curves showed that the high-risk group had significantly shorter OS (log-rank P=0.02; Figure 2G). The model achieved a C-index of 0.636 [95% confidence interval (CI): 0.539–0.733]. timeROC analysis yielded AUCs of 0.584 (95% CI: 0.392–0.773), 0.639 (95% CI: 0.477–0.807), and 0.693 (95% CI: 0.544–0.837) for 1-, 2-, and 3-year survival, respectively (Figure 2H). To evaluate whether both genes were required, we separately compared the prognostic performance of CFH alone, SPINT2 alone, and the combined CFH + SPINT2 model in the training cohort. The combined model showed the highest C-index (0.636) and remained significant in nested Cox regression (P=0.02), whereas SPINT2 alone showed limited discrimination (C-index =0.556, log-rank P=0.20) (Figure S1, Table S1).

Figure 2 Construction of the prognostic signature in the training set (GSE89749; n=115). (A) Univariate Cox regression forest plot showing the 10 genes most associated with OS. (B) LASSO coefficient paths for the 21 genes. Each colored curve represents the trajectory of a coefficient as λ increases. (C) Ten-fold cross-validation partial likelihood deviance as a function of log (λ). The vertical dashed line indicates λ.min (0.0348). (D) Random forest variable importance plot showing the top 10 genes ranked by the mean decrease in the Gini index. (E) Venn diagram showing the overlap among LASSO-Cox, random forest, and univariate Cox regression. Two overlapping genes, CFH and SPINT2, were identified. (F) External validation of expression patterns via GSE76297 and GSE107943. (G) Kaplan-Meier survival curves (log-rank P=0.02). (H) Time-dependent ROC curves, with 1-, 2-, and 3-year AUCs of 0.584, 0.639, and 0.693, respectively. AUC, area under the curve; CI, confidence interval; LASSO, least absolute shrinkage and selection operator; OS, overall survival; ROC, receiver operating characteristic.

External validation

To evaluate the generalizability of the model, we performed validation in an independent external cohort (E-MTAB-6389; n=75). The model achieved a C-index of 0.592 (95% CI: 0.509–0.675). Kaplan-Meier survival analysis showed that patients in the high-risk group had a significantly worse OS than did those in the low-risk group (log-rank P=0.006; Figure 3A). In the timeROC analysis, the 3-year AUC was 0.647 (Figure 3B). Multivariate Cox regression analysis including risk group, sex, and vascular invasion demonstrated that the risk group was an independent prognostic factor [hazard ratio (HR) =2.27, 95% CI: 1.18–4.37, P=0.01], whereas sex (P=0.69) and vascular invasion (P=0.39) were not statistically significant (Figure 3C). These findings confirm that the CFH-SPINT2 signature (CFH + SPINT2 signature) has independent clinical predictive value. Additional validation in the TCGA-CHOL cohort (n=36) further demonstrated acceptable discriminative ability (C-index =0.721, 95% CI: 0.620–0.821; 3-year AUC =0.692, 95% CI: 0.544–0.837; Figure S2). Baseline characteristics of patients stratified by risk group in the training and validation cohorts are summarized in Table S2. Similar single-gene versus combined comparisons were performed in E-MTAB-6389 (n=75) and TCGA-CHOL (n=36). CFH alone showed the strongest univariate discrimination in both cohorts, whereas nested Cox regression confirmed significantly better fit of the combined model (E-MTAB-6389, P=0.001; TCGA-CHOL, P=0.03), supporting independent prognostic contribution of SPINT2 beyond CFH.

Figure 3 Validation of the consensus model in the E-MTAB-6389 validation cohort. (A) Kaplan-Meier survival curves for the high- and low-risk groups (log-rank P=0.006). (B) Time-dependent ROC curves, with 1-, 2-, and 3-year AUCs of 0.548, 0.609, and 0.647, respectively. (C) Forest plot of multivariate Cox regression incorporating risk group, sex, and vascular invasion. The risk group was an independent prognostic factor (HR =2.27, 95% CI: 1.18–4.37; P=0.01), independent of sex and vascular invasion. AUC, area under the curve; CI, confidence interval; HR, hazard ratio; ROC, receiver operating characteristic.

Functional enrichment analysis

To elucidate the biological functions associated with the CFH + SPINT2 risk signature, we compared gene expression profiles between the high- and low-risk groups. A total of 109 DEGs were identified, including 74 upregulated and 35 downregulated genes in the high-risk group (|log₂FC| >1, adj. P<0.05). The distribution of these DEGs is visualized in the volcano plot in Figure 4A. The heatmap in Figure 4B shows the expression patterns of the top 30 DEGs, with red representing high expression and blue representing low expression, clearly distinguishing the two risk groups. Hallmark GSEA revealed that the high-risk group was characterized by significant enrichment of pathways associated with inflammation, immune modulation, and metastasis [e.g., epithelial-mesenchymal transition (EMT), inflammatory response, TNFA/NF-κB signaling, complement, and angiogenesis]. In contrast, pathways governing cell cycle progression and proliferation (e.g., MYC targets, E2F targets, and G2M checkpoint) were significantly suppressed in the high-risk group (Figure 4C). Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses of the upregulated genes confirmed the enrichment of complement-related and immune response pathways (Figure S3), while downregulated genes showed no significant enrichment.

Figure 4 Functional enrichment analysis of the high- and low-risk groups. (A) Volcano plot of DEGs. Red: upregulated (n=74); blue: downregulated (n=35); gray: not significant. Cutoff: |log₂FC| >1 and adj. P<0.05. (B) Heatmap of the top 30 DEGs. Rows: genes; columns: patients. (C) Hallmark GSEA comparing high- and low-risk groups. The NES is shown. Positive NES indicates activation in the high-risk group; negative NES indicates suppression. adj. P, adjusted P; DEG, differentially expressed gene; FC, fold change; GSEA, gene set enrichment analysis; NES, normalized enrichment score; NS, not significant.

Immune infiltration analysis

To investigate the relationship between the CFH + SPINT2 prognostic signature and the tumor immune microenvironment, we employed the ssGSEA algorithm to evaluate the infiltration levels of 28 immune cell types and analyzed the correlations of CFH/SPINT2 expression and risk score with immune cell infiltration. Correlation analysis revealed that CFH expression was positively correlated with most immune cell types, with the strongest correlations observed for macrophages (ρ=0.560; P<0.001), natural killer T (NKT) cells (ρ=0.496; P<0.001), and follicular helper T cells (ρ=0.478; P<0.001). In contrast, SPINT2 expression was negatively correlated with the majority of immune cell types and most strongly with NKT cells (ρ=−0.411; P<0.001), monocytes (ρ=−0.379; P<0.001), and NK cells (ρ=−0.351; P<0.001) (Figure 5A). Comparison of immune infiltration between the high- and low-risk groups revealed that the high-risk group had significantly higher infiltration of macrophages (log₂FC =3.66; P<0.001), myeloid-derived suppressor cells (MDSCs) (log₂FC =2.19; P<0.001), and NKT cells (log₂FC =1.62; P<0.001) (Figure 5B). These findings indicate that the CFH + SPINT2 risk score is closely associated with the TME. The high-risk group was characterized by increased infiltration of macrophages and MDSCs, with the latter being a well-established immunosuppressive cell population, suggesting a potentially immunosuppressive microenvironment.

Figure 5 Immune infiltration analysis. (A) Heatmap of correlations between CFH/SPINT2 expression and infiltration levels of 28 immune cell types. Red indicates positive correlation, and blue indicates negative correlation. ns, not significant; *, P<0.05; **, P<0.01; ***, P<0.001. (B) Boxplot comparing immune cell infiltration levels between the high- and low-risk groups, with immune cell types with significant differences between the two groups being shown (Wilcoxon test; *, P<0.05; **, P<0.01; ***, P<0.001). ssGSEA, single-sample gene set enrichment analysis.

Drug sensitivity analysis

To evaluate the potential clinical utility of the CFH + SPINT2 risk signature in guiding therapeutic strategies, we performed drug sensitivity analysis with the GDSC2 database. A total of 198 compounds were analyzed, and the predicted IC₅₀ values were compared between the high- and low-risk groups. Notably, 94 (47.5%) drugs showed significantly different predicted IC₅₀ values between the two risk groups (Wilcoxon rank-sum test P<0.05), of which 80 remained significant after FDR correction (FDR <0.05). The 10 drugs with the most significant differences included NU7441, AZD8055, entospletinib, SB216763, niraparib, GSK2606414, PCI-34051, dactolisib, IGF1R, and JQ1, all showing higher sensitivity in the high-risk group (Figure 6A). Encouragingly, the high-risk group exhibited significantly lower IC₅₀ values for the standard first-line chemotherapeutic agent gemcitabine (log₂FC =−1.33; P<0.001), as shown in Figure 6B and Figure S4A. Cisplatin also showed greater predicted sensitivity in the high-risk group (log₂FC =−0.54; P=0.004; Figure S4A). Additionally, the high-risk group showed increased sensitivity to several targeted therapies, including the poly(ADP-ribose) polymerase (PARP) inhibitors olaparib, niraparib, and talazoparib; the mammalian target of rapamycin (mTOR) inhibitors rapamycin, dactolisib, and AZD8055; and the BET inhibitor JQ1 (Figure 6B). In contrast, the low-risk group demonstrated greater sensitivity to the human epidermal growth factor receptor 2 (HER2)/epidermal growth factor receptor (EGFR) inhibitor lapatinib (log₂FC =0.31; P=0.007; Figure S4A,S4B). The heatmap in Figure 6C of the 30 most significant drugs shows distinct sensitivity patterns between the two risk groups. Five GDSC2 biliary cell lines (EGI-1, ETK-1, HuCCT1, TGBC1TKB, and TGBC24TKB) were identified, and their risk scores and observed IC₅₀ values are summarized in Table S3. We then extended this analysis to all available GDSC2 cell lines. The risk score was negatively correlated with observed IC₅₀ for gemcitabine (Spearman ρ=−0.197, P<0.001), cisplatin (ρ=−0.202, P<0.001), olaparib (ρ=−0.179, P<0.001), and palbociclib (ρ=−0.162, P<0.001), and positively correlated with IC₅₀ for lapatinib (ρ=0.142, P<0.001) (Table S4).

Figure 6 Drug sensitivity analysis of the CFH + SPINT2 risk signature. (A) Boxplots showing predicted IC₅₀ values of the 10 drugs with the greatest significant sensitivity differences between the high- and low-risk groups. P values were calculated with the Wilcoxon test (****, P<0.0001). (B) Bar plot showing the 20 drugs to which the high-risk group was most sensitive (log₂FC <0 and P<0.05). Negative log₂FC indicates a lower IC₅₀ (higher sensitivity) in the high-risk group. (C) Heatmap of drug sensitivity for the 30 most significant drugs. Rows: drugs. Columns: patients. Red: higher sensitivity (lower IC₅₀); blue: lower sensitivity (higher IC₅₀). Patients are ordered by risk group (high risk in red and low risk in blue). FC, fold change; IC₅₀, half-maximal inhibitory concentration.

Single-cell localization of CFH and SPINT2

To clarify the cellular origins of CFH and SPINT2 in the CCA microenvironment, we performed single-cell RNA-sequencing analysis on the GSE138709 dataset (30,729 cells). After quality control, dimensionality reduction, clustering, and annotation, eight major cell types were identified (Figure 7A). CFH was predominantly expressed in fibroblasts, whereas SPINT2 was predominantly expressed in cholangiocytes and malignant cells (Figure 7B,7C). A scatter plot further confirmed the cell-type-specific expression pattern of these two genes (Figure 7D). These results indicate that CFH is primarily derived from fibroblasts, while SPINT2 originates from malignant cells and cholangiocytes. To further characterize fibroblast heterogeneity, fibroblasts were subclustered into myCAF, iCAF, and apCAF (Figure S5A). Canonical CAF subtype markers showed expected subtype-specific enrichment (Figure S5B). CFH was preferentially enriched in apCAF and iCAF, with substantially lower expression in myCAF, whereas SPINT2 remained low across all fibroblast subtypes (Table S5). CFH expression was negatively correlated with the myCAF module score (Pearson r=−0.495).

Figure 7 Single-cell localization of CFH and SPINT2 in CCA. (A) t-SNE plot showing eight major cell types (N=30,729 cells). (B) t-SNE plots showing expression of CFH (red) and SPINT2 (purple). CFH was predominantly expressed in fibroblasts, and SPINT2 was expressed in cholangiocytes and malignant cells. (C) Average expression levels across cell types. CFH had the highest expression in fibroblasts, while SPINT2 had the highest expression in cholangiocytes and malignant cells. (D) Scatter plot comparing CFH and SPINT2 expression, demonstrating a cell-type-specific expression pattern. CCA, cholangiocarcinoma; t-SNE, t-distributed Stochastic Neighborhood Embedding.

Discussion

In this study, we successfully constructed and externally validated a CFH- and SPINT2-based dual-gene prognostic model for CCA using an integrated machine learning approach combining LASSO-Cox regression, random forest, and univariate Cox regression. Immune infiltration analysis revealed that the high-risk group exhibited an immunosuppressive microenvironment characterized by increased infiltration of macrophages and MDSCs. Drug sensitivity analysis demonstrated that the high-risk group was more sensitive to gemcitabine, cisplatin, PARP inhibitors, and mTOR inhibitors. Notably, single-cell RNA sequencing revealed the cell-type-specific expression of these two genes: CFH was predominantly expressed in fibroblasts, whereas SPINT2 was mainly expressed in malignant cells. These findings provide novel insights into the biological basis, immune landscape, and therapeutic implications of this prognostic model from the perspective of stroma-tumor interactions.

In our prognostic model, CFH served as a risk-associated gene, with a higher expression linked to worse outcomes, consistent with a recent pan-cancer analysis (20). Although bulk transcriptome analysis showed CFH was significantly downregulated in tumor tissues (log₂FC =−2.81), our single-cell data resolved this paradox: CFH expression was highest in fibroblasts and lower in malignant cells; thus, the overall downregulation likely reflects reduced fibroblast abundance, whereas high CFH expression marks activated CAFs that suppress complement-mediated antitumor immunity, fostering an immunosuppressive microenvironment (21,22). This interpretation is supported by our immune infiltration analysis, in which the high-risk group exhibited increased macrophage and MDSC infiltration, and CFH correlated with macrophage abundance (ρ=0.560). Fibroblast subclustering further showed that CFH was enriched in apCAF and iCAF rather than myCAF, supporting its association with immunoregulatory rather than purely matrix-remodeling CAFs and providing a cellular basis for the immunosuppressive TME observed in the high-risk group.

SPINT2 served as a protective gene in our prognostic model, with lower expression significantly associated with poor prognosis. Single-cell analysis revealed that SPINT2 was primarily expressed in malignant cells and cholangiocytes. SPINT2 inhibits HGFA and thereby restrains HGF/c-MET signaling (23); its downregulation may promote an invasive phenotype, consistent with EMT enrichment in the high-risk group and its tumor-suppressor roles in other cancers (24-26). Furthermore, SPINT2 expression was negatively associated with infiltration of NKT cells, monocytes, and NK cells, suggesting that loss of SPINT2 may be linked to an inflammatory microenvironment associated with tumor progression.

Beyond their individual roles, CFH and SPINT2 may cooperatively shape an immunosuppressive TME through a dual-compartment stroma-tumor axis. In the high-risk phenotype (CFH-high/SPINT2-low), fibroblast-derived CFH (enriched in immunoregulatory CAF subsets) and malignant-cell SPINT2 loss likely act in parallel, and these two programs may converge on a myeloid-enriched, inflammatory yet immunosuppressive TME, although direct causal interactions require experimental validation.

Consistent with this cooperative model, the high-risk group (CFH-high/SPINT2-low) exhibited increased macrophage and MDSC infiltration, together with an inflammatory component associated with SPINT2 loss. In terms of therapeutic significance, the high-risk group showed enhanced sensitivity to gemcitabine, cisplatin, PARP inhibitors, and mTOR inhibitors, whereas the low-risk group was more sensitive to lapatinib. These findings support a risk-stratified treatment strategy and provide a rational basis for clinical application of the CFH + SPINT2 signature.

The model showed lower discriminative ability in the E-MTAB-6389 cohort, with a C-index =0.592 (vs. 0.721 in TCGA-CHOL), likely due to its higher event rate (66.7%), indicating more advanced-stage patients. Despite this, the model demonstrated robust performance in both the training and external validation cohorts, and the risk score was an independent prognostic factor after adjustments for sex and vascular invasion, supporting its utility for risk stratification in patients with CCA. Additionally, the model includes only two genes and thus is low cost and feasible for clinical implementation. Finally, single-cell RNA sequencing revealed a cell-type-specific expression network between fibroblast-derived CFH and malignant cell-derived SPINT2, suggesting that therapeutic strategies targeting CAFs or the HGF/c-MET pathway may have better efficacy in high-risk patients.

Several limitations to this study should be acknowledged. First, all data were derived from public databases, and independent validation with prospective clinical cohorts was not conducted. Second, the study was primarily based on transcriptomic analysis, and validation at the protein level and functional experiments are needed to further elucidate the biological functions of CFH and SPINT2. Third, the specific regulatory mechanisms underlying CFH secretion in CAFs and the compensatory upregulation of SPINT2 in CCA remain to be further clarified. Fourth, the single-cell dataset originated from a single dataset, and there was no external single-cell validation; therefore, the inference of reduced fibroblast abundance should be confirmed with quantitative evidence in future studies incorporating spatial transcriptomics or immunohistochemistry. Fifth, although we supplemented bulk drug predictions with GDSC2 biliary cell line data and pan-cancer correlations, only five biliary-specific lines were available, and experimental validation of drug response in CCA remains necessary.

Conclusions

We constructed and validated a CFH- and SPINT2-based dual-gene prognostic signature for postoperative risk stratification in CCA. Single-cell analyses revealed a stroma–tumor expression pattern involving fibroblast-derived CFH and malignant cell-derived SPINT2, which together relate to an immunosuppressive TME. This signature may facilitate prognostic assessment after resection and help inform adjuvant treatment selection.

Acknowledgments

None.

Footnote

Reporting Checklist: The authors have completed the TRIPOD reporting checklist. Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0798/rc

Peer Review File: Available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0798/prf

Funding: This work was supported by the Project of Fujian University of Traditional Chinese Medicine (Clinical Special Project; grant No. XB2024091).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://jgo.amegroups.com/article/view/10.21037/jgo-2026-0798/coif). All authors report that this work was supported by the Project of Fujian University of Traditional Chinese Medicine (Clinical Special Project; grant No. XB2024091). The authors have no other conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. The study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Liu J, Fang W, He X, et al. Global burden of early-onset gallbladder and biliary tract cancer from 1990 to 2021. BMC Gastroenterol 2025;25:461. [Crossref] [PubMed]
2. Li X, Guan R, Zhang S. Factors Contributing to the High Malignancy Level of Cholangiocarcinoma and Its Epidemiology: Literature Review and Data. Biology (Basel) 2025;14:351. [Crossref] [PubMed]
3. Hussain MM, Wang JM, Zhai AQ, et al. Comparison of prognostic factors and their differences in intrahepatic, hilar, and distal cholangiocarcinoma: A systematic review and meta-analysis. World J Gastrointest Oncol 2025;17:107995. [Crossref] [PubMed]
4. Liu CX, Wong CC. Stratifying cholangiocarcinoma: tumor microenvironment, molecular drivers, and novel immunotherapeutic approaches. Clin Mol Hepatol 2026;32:127-55. [Crossref] [PubMed]
5. Beyersdorf J, Zhang ML. Extrahepatic cholangiocarcinoma: Current concepts in histopathology, immunohistochemistry, and molecular diagnostics. Semin Diagn Pathol 2025;42:150949. [Crossref] [PubMed]
6. Qurashi M, Vithayathil M, Khan SA. Epidemiology of cholangiocarcinoma. Eur J Surg Oncol 2025;51:107064. [Crossref] [PubMed]
7. Laurenzi A, Brandi G, Greco F, et al. Can repeated surgical resection offer a chance of cure for recurrent cholangiocarcinoma? Langenbecks Arch Surg 2023;408:102. [Crossref] [PubMed]
8. Macarulla T, Neuzillet C, Prager GW, et al. Opportunities and Approaches to Optimising Advanced Cholangiocarcinoma Outcomes in the Era of Targeted Therapies: A Narrative Review. Oncol Ther 2025;13:939-62. [Crossref] [PubMed]
9. Khan SA, Rushbrook SM, Kendall TJ, et al. Guidelines Development Group for the British Society of Gastroenterology guidelines for the diagnosis and management of cholangiocarcinoma. Gut 2025;74:504-5. [Crossref] [PubMed]
10. Czekay RP, Higgins CE, Aydin HB, et al. SERPINE1: Role in Cholangiocarcinoma Progression and a Therapeutic Target in the Desmoplastic Microenvironment. Cells 2024;13:796. [Crossref] [PubMed]
11. Minini M, Fouassier L. Cancer-Associated Fibroblasts and Extracellular Matrix: Therapeutical Strategies for Modulating the Cholangiocarcinoma Microenvironment. Curr Oncol 2023;30:4185-96. [Crossref] [PubMed]
12. Huang F, Liu Z, Song Y, et al. Bile acids activate cancer-associated fibroblasts and induce an immunosuppressive microenvironment in cholangiocarcinoma. Cancer Cell 2025;43:1460-1475.e10. [Crossref] [PubMed]
13. Lin X, Zhou J, Yonemori K, et al. Recent advances in targeted and immunotherapeutic strategies for cholangiocarcinoma: a comprehensive review. Genes & Diseases 2026; In Press. [Crossref]
14. Xue J, Zhang L, Zhang K, et al. Immunotherapy in biliary tract cancer: reshaping the tumour microenvironment and advancing precision combination strategies. Front Immunol 2025;16:1651769. [Crossref] [PubMed]
15. Nishida N. Biomarkers and Management of Cholangiocarcinoma: Unveiling New Horizons for Precision Therapy. Cancers (Basel) 2025;17:1243. [Crossref] [PubMed]
16. Ding Z, Zeng Y. Artificial intelligence and multi-omics driven models would be the future of intrahepatic cholangiocarcinoma prediction research. Hepatobiliary Surg Nutr 2024;13:560-1. [Crossref] [PubMed]
17. Choi WJ, Walker R, Rajendran L, et al. Call to Improve the Quality of Prediction Tools for Intrahepatic Cholangiocarcinoma Resection: A Critical Appraisal, Systematic Review, and External Validation Study. Ann Surg Open 2023;4:e328. [Crossref] [PubMed]
18. Zhang Z, Wu G. Considerations for improving generalizability and robustness of predictive models in perihilar cholangiocarcinoma. Hepatobiliary Surg Nutr 2026;15:19. [Crossref] [PubMed]
19. Charoentong P, Finotello F, Angelova M, et al. Pan-cancer Immunogenomic Analyses Reveal Genotype-Immunophenotype Relationships and Predictors of Response to Checkpoint Blockade. Cell Rep 2017;18:248-62. [Crossref] [PubMed]
20. Chu D, Huang R, Shi J, et al. NETs-related genes predict prognosis and are correlated with the immune microenvironment in osteosarcoma. Front Oncol 2025;15:1551074. [Crossref] [PubMed]
21. Saxena R, Gottlin EB, Campa MJ, et al. Complement factor H: a novel innate immune checkpoint in cancer immunotherapy. Front Cell Dev Biol 2024;12:1302490. [Crossref] [PubMed]
22. Daugan MV, Revel M, Thouenon R, et al. Intracellular Factor H Drives Tumor Progression Independently of the Complement Cascade. Cancer Immunol Res 2021;9:909-25. [Crossref] [PubMed]
23. Liu F, Cox CD, Chowdhury R, et al. SPINT2 is hypermethylated in both IDH1 mutated and wild-type glioblastomas, and exerts tumor suppression via reduction of c-Met activation. J Neurooncol 2019;142:423-34. [Crossref] [PubMed]
24. Pereira MS, Celeiro SP, Costa ÂM, et al. Loss of SPINT2 expression frequently occurs in glioma, leading to increased growth and invasion via MMP2. Cell Oncol (Dordr) 2020;43:107-21. [Crossref] [PubMed]
25. Ma Z, Liu D, Li W, et al. STYK1 promotes tumor growth and metastasis by reducing SPINT2/HAI-2 expression in non-small cell lung cancer. Cell Death Dis 2019;10:435. [Crossref] [PubMed]
26. Ye P, Yang X, Wang H, et al. SPINT2 inhibits NEDD4L-mediated ACSL4 ubiquitination to promote ferroptosis and suppress gallbladder cancer progression. Int J Biol Macromol 2025;330:148141. [Crossref] [PubMed]

(English Language Editor: J. Gray)