We obtained Transcripts Per Million (TPM)-normalized RNA sequencing (RNA-seq) data for Colon Adenocarcinoma (COAD) via The Cancer Genome Atlas (TCGA) data portal (https://portal.gdc.cancer.gov/). A total of 524 samples were downloaded, comprising 455 tumor samples and 69 normal samples (23). Corresponding clinical data were retrieved from the UCSC Xena database (https://xenabrowser.net/), which initially included 530 samples. Duplicate samples were removed, and only those overlapping with the RNA-seq dataset were retained. Samples lacking clinical information were excluded, resulting in a final cohort of 480 patients, including 431 tumor samples and 49 normal samples. For validation, we obtained the GSE39582 dataset from the Gene Expression Omnibus (GEO) database (www.ncbi.nlm.nih.gov/geo/). After removing duplicates and samples with missing values, the dataset comprised 566 tumor cases and 19 normal cases. All data processing was performed using R software (v4.5.1), without applying batch correction. TCGA data were extracted and processed using the R package tidyverse (2.0.0); gene expression levels were in TPM format and were finally processed by taking the log2 + 1 value. The GSE39582 dataset was extracted and processed using the R package GEOquery (v2.76.0), and probe identifiers were converted using the corresponding GPL570–55999 platform annotation files obtained from the official GEO website.

We secured the gene set “HALLMARK BILE ACID METABOLISM” from the Molecular Signatures Database (MSigDB; https://www.gsea-msigdb.org/gsea/msigdb) (24), which is a curated compilation of genes specifically associated with bile acid metabolic mechanisms. Using this gene set, we calculated the bile acid metabolism gene-set score for each specimen. To identify molecular subtypes of COAD associated with bile acid metabolism, we performed unsupervised consensus clustering on filtered COAD samples via the ConsensusClusterPlus (k = 6, 1,000 resampling iterations) package (v 1.72.0) in R (25). Taking into account the above analyses and the statistical power required for subsequent survival analyses, immune infiltration comparison, and other statistical tests for each subgroup, we ultimately classified COAD patients into distinct molecular subtypes.

To assess the robustness of the molecular typing results, we used two independent Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway gene sets for validation: the bile secretion pathway and the primary bile acid biosynthesis pathway (data obtained using the R package KEGGREST v1.48.1). Pathway activity scores were calculated as the average expression levels of genes within each pathway, after which correlations with the original Hallmark pathway scores were evaluated. Survival analysis was then conducted based on stratification using these gene sets. All three enrichment scores were calculated using GSVA implemented in the R package GSVA (v2.2.0).

Survival and immune infiltration analyses of molecular subtypes based on bile acid metabolism were performed. Overall survival (OS) was defined as the study endpoint, and Kaplan–Meier survival curves were generated for each molecular subtype.

The log-rank (Mantel–Cox) test was used to compare survival differences among subtypes, and a two-sided p < 0.05 was considered statistically significant. Moreover, the connection between the continuous bile acid metabolism score and patient survival was evaluated using multivariate Cox proportional hazards regression, with results reported as hazard ratios. All statistical analyses and visualizations were conducted in R using the packages ggplot2 (v3.5.2), survminer (v0.5.0), and survival (v3.8-3).

We employed the Cell-type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORT) algorithm to estimate the relative percentages of 22 immune cell types (using the LM22 signature matrix; downloaded from https://cibersort.stanford.edu/). Only samples with a deconvolution p-value < 0.05 from the CIBERSORT output were retained (perm = 100). The main differences in immune cell proportions among subgroups were compared using a two-sided Wilcoxon test, with p < 0.05 indicating statistical significance. All computations were implemented in R using the packages e1071 (v1.7-16), preprocessCore (v1.70.0), ggplot2 (v3.5.2), pheatmap (v1.0.13), and parallel (base package). To verify the stability of the above results, we further applied the xCell algorithm to 431 samples using the xCell package (v1.1.0) and then extracted the immune cell-related data for downstream analysis.

We identified differential expression genes (DEGs) by comparing gene expression differences across subgroups using the limma package (v3.64.3) in R. The selection criteria were |log2(fold change)| > 1 and an adjusted FDR < 0.05. KEGG enrichment analysis and Gene Ontology (GO) pathway enrichment analysis of DEGs were performed using the clusterProfiler (v4.16.0) and enrichplot (v1.28.4) R packages (26, 27).

To prioritize biologically relevant and prognostically significant genes, we employed a stringent two-step filtering strategy. First, 70 DEGs were submitted to the STRING online database (28) (https://string-db.org), and a PPI network was constructed using Cytoscape software (v3.9.1). The top 20 genes with the highest degree values were selected using the “Degree” algorithm within the PPI network. Moreover, Cox analysis (R package survival, v3.8-3) was performed on the 70 identified DEGs to screen for genes significantly associated with OS, using a two-sided Wald test with p < 0.05 as the significance threshold (without adjustment for multiple testing). The final hub genes for validation were defined as the intersection of two gene sets.

Based on the screened key genes, all patients were stratified into high- and low-expression groups. Survival analysis was performed between the two subgroups. The extent of immune cell infiltration in each subgroup was analyzed using a two-sided Wilcoxon test (p < 0.05 was considered statistically significant). In addition, the bile-acid metabolism score (single-sample Gene Set Enrichment Analysis (ssGSEA); same method as described in Section 2.1) was compared between the two groups. Spearman correlation analysis (implemented using the stats package in R) was employed to examine correlations between gene expression levels and key T-cell chemotactic factors (C-X-C motif chemokine ligand 9 (CXCL9), C-X-C motif chemokine ligand 9 (CXCL10)). The R packages used for survival analysis and the Wilcoxon test were the same as those described previously.

All samples’ standardized expression matrices (gene expression z-scores) were uploaded to the online Tumor Immune Dysfunction and Exclusion (TIDE) platform (link: http://tide.dfci.harvard.edu/), and relevant metrics, including the TIDE score, dysfunction score, and exclusion score, were obtained for each sample. According to the developer’s recommendations, a TIDE score > 0 indicated a predicted nonresponder, whereas a score ≤ 0 indicated a predicted responder. Differences in TIDE scores between expression-defined groups were assessed using the Wilcoxon rank-sum test, and correlations between TIDE scores and gene expression levels were evaluated using Spearman correlation analysis (29). In addition, the associations between the expression levels of CLCA1, UGT2A3, and ZG16 and bile acid metabolism pathway activity (ssGSEA score) were examined using a two-sided Wilcoxon rank-sum test and a two-sided Spearman correlation analysis.

Continuous variables were analyzed according to the number of groups. For comparisons between two groups, the Wilcoxon rank-sum test was used; for comparisons involving more than two groups, the Kruskal–Wallis test was applied. Categorical variables were compared using the Chi-square test (χ² test). Correlations between variables were assessed using Spearman’s rank correlation analysis. All statistical tests run were two-sided, and a p-value < 0.05 was considered statistically significant.

To establish a molecular subtyping system based on bile acid metabolism, we applied unsupervised consensus clustering to evaluate clustering schemes with k ranging from 2 to 6. As k increased, the diagonal of the consensus matrix became progressively clearer (Figure 2a). The CDF curve exhibited a distinct separation and clear upward trend at k = 2 (Figure 2b), while both cluster-consensus and item-consensus analyses indicated that k = 2 and k = 4 had the highest stability (Figures 2c, d). Although k = 4 demonstrated comparable stability, k = 2 provided the advantage of more balanced subgroup sizes (> 100 patients per group), ensuring sufficient statistical power for downstream survival and immune analyses and avoiding potential overstratification. When k = 2, the patients were divided into two subgroups with balanced sample sizes (bile-high group, n = 215; bile-low group, n = 216), which enhanced the statistical power of subsequent analyses. The baseline characteristics of the two subgroups are presented in Table 1. No significant differences were observed in age (p = 0.774, Wilcoxon test), gender (p = 0.598, Chi-square test), and tumor stage (p = 0.292, Chi-square test) between the groups. However, tumor location differed significantly between the subgroups (p = 0.001, Chi-square test). Therefore, this grouping scheme was adopted for further analyses. Survival analysis indicated that overall survival was significantly shorter in the bile-low group compared with the bile-high group (p = 0.0049, Figure 3a). Multivariate Cox regression analysis, adjusting for key clinical confounding factors, confirmed that this subtyping model was an independent prognostic factor (HR = 0.13; 95% CI: 0.02–0.73; p = 0.020; Figure 3b). Independent validation analysis demonstrated that the activity scores of the KEGG bile secretion pathway and the primary bile acid biosynthesis pathway were significantly positively correlated with the original Hallmark bile acid metabolism score (Figures 3c, d; Spearman’s ρ = 0.337 and 0.411, respectively; p < 0.001 for both).

Based on the results of CIBERSORT, a total of 319 samples (74% for total samples) had p < 0.05 and were included in the subsequent immune infiltration analysis. We examined the immune cell infiltration landscape of these samples (Figure 4a). Correlation analysis revealed that cytotoxic lymphocytes (CD8⁺ T cells and NK cells) were significantly positively correlated (Figure 4b), which may play an important role in antitumor immunity. Interestingly, the abundance of CD8⁺ T cells was negatively correlated with the abundance of immunosuppressive cells, suggesting a potential antagonistic relationship between effector and suppressive immune components. Figure 4c depicts the detailed distribution of specific immune cell types between the two groups (CIBERSORT). Compared with the bile-high group, the bile-low group exhibited significantly higher levels of CD8⁺ T cells and M1 macrophages (p < 0.05, Figure 4c). In addition, the xCell analysis confirmed that, in addition to the consistent infiltration trends of CD8⁺ T cells (p < 0.001) and M1 macrophages (p < 0.01) observed in the CIBERSORT results, B cells (p < 0.05) and CD4⁺ naive T cells (p < 0.05) also showed increased infiltration in the bile-low group compared with the bile-high group (Figure 4d). These findings provide important information for further exploring the mechanisms by which bile acid metabolism regulates the function of specific immune cell types.

Intergroup comparison identified 70 significantly DEGs, of which 59 were upregulated, and 11 were downregulated. The distribution of these genes is illustrated in a volcano plot (Figure 5a), while a heatmap in Figure 5b shows the top 50 genes. GO analysis indicated that the DEGs were primarily associated with the biological processes “killing of cells of other organisms,” “antimicrobial peptide-mediated humoral immune response,” and “response to nutrients” (Figure 5c). KEGG pathway analysis revealed that the DEGs were notably enriched in “retinol metabolism,” “drug metabolism–cytochrome P450,” and “fat digestion and absorption” (Figure 5d).

Subsequently, the top 20 hub genes with the highest degree values were selected from the PPI network (Figures 6a, b). Univariate Cox regression analysis revealed 10 DEGs significantly associated with OS (p < 0.05; Figure 6c). Given that the Cox regression analysis was not adjusted for multiple testing, these genes were considered exploratory candidate prognostic genes. The intersection of hub genes from the PPI network and DEGs significantly correlated with OS yielded three genes: ZG16, CLCA1, and UGT2A3 (Figure 6d).

Regarding immunocyte infiltration, CLCA1 expression was associated with distinct infiltration patterns of various immune cell types compared with low expression (Figure 7a). Similarly, high expression of UGT2A3 and ZG16 showed significantly different infiltration scores across multiple immune cell types compared with low expression (Figures 7b, c). Survival analysis indicated that patients with high CLCA1 expression had notably better overall survival than those with low expression (p = 9e−04; Figure 7d). Although the survival curves for UGT2A3 (p = 0.23; Figure 7e) and ZG16 (p = 0.17; Figure 7f) appeared visually separated, the differences did not reach statistical significance; these genes should therefore be considered hypothesis-generating candidates requiring further validation. Notably, patients with high expression of CLCA1, UGT2A3, or ZG16 also exhibited significantly higher bile acid metabolism scores (all p < 0.001; Figures 7g–i), confirming that these hub genes are functionally linked to bile acid metabolic activity in CRC. Spearman correlation analysis revealed that the expression of these three genes was not significantly associated with CXCL9 or CXCL10 expression (all Spearman’s correlation coefficients |R| < 0.19; Figure 7j).

TIDE analysis indicated that the low-expression groups of CLCA1, UGT2A3, and ZG16 had significantly higher TIDE scores compared with the high-expression groups (p < 0.0001, p = 0.0029, p = 0.029; Figures 8a–c). Remarkably, Spearman correlation analysis revealed a significant negative correlation between the expression of these three genes and the TIDE score (CLCA1: R = − 0.24, p < 0.001; UGT2A3: R = − 0.15, p = 0.0022; ZG16: R = − 0.14, p = 0.0039; Figures 8d–f).

Our study reveals a paradoxical immune landscape in the bile-low subgroup: despite elevated infiltration of CD8⁺ T cells and M1 macrophages, patients exhibit significantly worse survival. This suggests that immune cell quantity alone is insufficient for effective antitumor immunity, and functional impairment may predominate. Notably, we found no significant correlation between CLCA1/UGT2A3/ZG16 expression and key T-cell chemokines CXCL9/CXCL10, implying that the observed T-cell accumulation is unlikely to be driven by canonical chemotactic signaling. Instead, it may reflect T-cell exhaustion or spatial exclusion—a phenotype captured by the TIDE score but not by bulk infiltration metrics. The strong association of low bile acid metabolism with right-sided tumors further suggests site-specific metabolic-immune crosstalk, possibly mediated by microbiota or bile acid gradients (16, 30–33).

A notable finding of this study is that CLCA1, UGT2A3, and ZG16—three genes highly expressed in normal colorectal tissues—are consistently downregulated in tumor tissues, suggesting that their loss may represent a coordinated event in CRC progression. Low expression of these genes is significantly associated with higher TIDE scores, indicative of immunotherapy resistance. Among the hub genes, only CLCA1 exhibited a robust association with favorable overall survival, whereas UGT2A3 and ZG16 serve as hypothesis-generating candidates that warrant further validation. Importantly, these hub genes also display strong positive correlations with bile acid metabolism pathway activity, highlighting their potential role as functional links between metabolic dysregulation and immune dysfunction in colorectal cancer.

CLCA1 is an accessory protein for chloride ion channels and plays a well-established role in mucus secretion and the maintenance of intestinal epithelial barrier integrity. Reduced expression of CLCA1 may cause “damage to the intestinal mucosal barrier,” increasing local inflammation and attracting immunosuppressive cells, which can facilitate tumor development (34). Based on our findings, we hypothesize that CLCA1 loss in CRC not only impairs the physical barrier but also generates proinflammatory milieu that attracts immunosuppressive cells. Consistent with this, lower CLCA1 expression was associated with higher TIDE scores, suggesting that CLCA1 deficiency contributes to a TIME characterized by T-cell dysfunction, potentially due to the failure to maintain a mucosal microenvironment supportive of effective immunity.

UGT2A3 functions as an essential second-stage enzyme in bile acid metabolism, catalyzing the glucuronidation of bile acids to facilitate their excretion. Reduced expression of UGT2A3 may lead to the “accumulation of toxic bile acids within cells, activation of oncogenic signaling pathways, and induction of endoplasmic reticulum stress. This, in turn, could recruit immunosuppressive cell populations through the release of specific chemokines, further remodeling the TIME (35). In this study, UGT2A3 expression was significantly downregulated in tumors. We hypothesize that decreased UGT2A3 expression may result in the accumulation of specific bile acid species, triggering stress-response or protumorigenic signaling pathways. The precise mechanisms underlying these effects require further investigation in future studies.

ZG16 is a lectin, and its reduced expression may impair host–microbiota interactions, increase the risk of pathogenic bacteria translocation, trigger chronic inflammation, and create an immunosuppressive environment conducive to tumor growth (36). We speculate that decreased ZG16 expression in CRC could destabilize the selective interface with the gut microbiota, facilitating the translocation of pathogenic bacteria or altering microbially derived signals, thereby promoting a state of chronic, low-grade inflammation.

All three genes illustrate how the interruption of the “metabolism–barrier–immunity” axis contributes to colorectal cancer progression and T-cell dysfunction, offering a new perspective for exploring the mechanisms of immune escape in CRC. Naturally, this hypothesis requires further validation.

The observed association between low CLCA1, UGT2A3, and ZG16 expression and elevated TIDE scores highlights a potential link between bile acid metabolism and T-cell dysfunction in CRC. The consistency with TIDE (a tool specifically focused on T-cell exhaustion and exclusion) suggests that these genes may play a role in functional immune impairment. From another perspective, the TIDE score is designed to model T-cell dysfunction and immune rejection and is particularly effective for predicting responses to ICIs. This indicates that these genes may influence the functional state of T cells rather than overall immune infiltration or cytotoxic potential. This interpretation aligns with the phenotype observed in our study, characterized by “high CD8⁺ T-cell infiltration but poor prognosis”. Thus, rather than functioning as standalone predictive biomarkers, CLCA1, UGT2A3, and ZG16 are more likely mechanistic mediators linking bile acid metabolism to T-cell dysfunction in CRC. Animal experiments have confirmed that a circular single-stranded DNA (CSSD) can upregulate CLCA1 expression (37). Similarly, the UGT2A3 expression can be modulated using specific antibodies (38). Additionally, interventions targeting the intestinal microbiota can indirectly upregulate the expression of genes like ZG16, thereby reshaping the TIME (39). However, studies investigating the regulation of these three key genes remain at the preclinical stage or have been conducted in other tumor types.