CCA sequencing data were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/). All the datasets used in this study are listed in Supplementary Table S1.

The single-cell datasets were downloaded from the GEO database (Supplementary Table S1). The raw gene expression matrix was imported and processed using the Seurat R package (version 3.1.5) (Butler et al., 2018). Cells with unique molecular identifier counts below 200 or mitochondrial content above 35% were removed. Normalization and dimension-reduction processes were performed using Seurat R package (3.1.5). Clusters were computed using the FindClusters function (resolution = 0.8) and visualized using uniform manifold approximation and projection (UMAP), as implemented in Seurat. Differential expression between clusters was calculated using a likelihood-ratio test for single-cell gene expression implemented in Seurat, with a family-wise error rate of 5%. Cell types were defined according to lineage-specific marker genes.

Interactions are frequent among TME cells depending on the ligand and receptor pairs. We evaluated the contribution of cell-interaction genes to cell communication by measuring the importance of their participant-enriched pathways in each cell type. Meanwhile, transcription factors can enhance cell signaling by promoting the expression of interacting genes. The importance and critical TF regulons were calculated using empirical Bayes loopy belief propagation (eLBP). The cell clustering results and differential genes in each cell type were used as the input data for eLBP. The algorithm consists of the following five steps based on cell interaction genes (CIGs).

Human leukemia monocytic THP-1 cells were cultured in Gibco™ RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS), 100 units/ml penicillin, and 100 g/ml streptomycin. THP-1 cells were seeded at a density of 6 × 10⁵ cells/well into 24-well plates and treated with 200 ng/ml phorbol12-myristate13-acetate (PMA; Sigma-Aldrich, St. Louis, MO, United States) for 48 h. The cells were cultured at 37°C in a humidified atmosphere with 5% CO₂, where THP-1 monocytes were differentiated into macrophages.

Small interfering RNA (siRNAs) against IRF1 and SPI1 were used to validate the regulation of the expression of LAGLS9 in THP-1 cells. IRF1 and SPI1 siRNAs were designed and synthesized by GenePharma (Shanghai, China). THP-1 cells were transfected with GP-transfect-Mate, according to the manufacturer’s instructions. Briefly, THP-1 cells were seeded at a density of 6 × 10⁵ cells/well in 24-well plates and treated with 200 ng/ml for 48 h. GP-transfect-Mate was diluted in Opti-MEM (50 ml; Invitrogen, Waltham, MA, United States) for 5 min before mixing with an equal volume of Opti-MEM containing siRNA (40 pmol). After 20 min of incubation, 100 μL of the resulting GP-transfect-Mate/siRNA mixture was added directly to the cells. After 24 h of incubation at 37°C in a 5% CO₂ atmosphere, cells were harvested for real-time PCR. The sequences of the primers used were as follows: SPI1-siRNA: CAG​GCA​GCA​AGA​AGA​AGA​UTT and AUC​UUC​UUC​UUG​CUG​CCU​T. IRF1-siRNA: GGG​CUC​AUC​UGG​AUU​AAU​ATT, UAU​UAA​UCC​AGA​UGA​GCC​CTT.

Total RNA was isolated using TRIzol reagent (Sigma-Aldrich), according to the manufacturer’s instructions. cDNA was synthesized using 5× All-in-One RT-Master Mix (Abm, Canada). Real-time PCR was performed in a LightCycler 96 using 2× chemQSYBR QPCR Master Mix (Vezyme, China), according to the manufacturer’s instructions, in a total volume of 10 μL. The sequences of the primers used were SPI1-Fw and GCGTGCAAAATGGAAG GGTTT. SPI1-Rev: GGT​ATC​GAG​GAC​GTG​CAT​CT. IRF1-Fw, GCT​GGG​ACA​TCA​ACA​AGG​AT. IRF1-Rev: CCT​GCT​CTG​GTC​TTT​CAC​CT; LGALS9-Fw: AAG​GTG​ATG​GTG​AAC​GGG​AT. LGALS9-Rev: ACT​GTC​TGG​GTA​ATG​GGA​GC; CD68-Fw: TCC​AGG​GAA​GCT​GTG​AGG​GT. CD68-Rev: AGC​CGA​GAA​TGT​CCA​CTG​TGC. CD163-Fw: TTT​GTC​AAC​TTG​AGT​CCC​TTC​AC. CD163-Rev: TCC​CGC​TAC​ACT​TGT​TTT​CAC. VEGFA-Fw: TCC​TCA​CAC​CAT​TGA​AAC​CA. VEGFA-Rev: TTT​TCT​CTG​CCT​CCA​CAA​TG.

Genes from eLBP were selected for consistent clustering using the R software package ConsensusClusterPlus to sort the immune molecular CCA subtypes (Wilkerson & Hayes, 2010). Correlations between subtypes, clinical features, immunity, and prognosis were analyzed.

The somatic mutation data of all patients with CCA was obtained from the International Cancer Genome Consortium (ICGC, https://dcc.icgc.org/). Mutations were analyzed and visualized using maftools (Mayakonda et al., 2018). The enrichment scores of the hallmark genes were evaluated through single-sample GSEA (ssGSEA) using the “GSVA” R package (Hänzelmann et al., 2013). Hallmark gene sets were obtained from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/).

The “limma” package was used to perform differentially expressed gene (DEG) analysis. An empirical Bayesian method was applied to estimate the DEGs between two clusters, which were identified using a consistent clustering method based on moderated t-tests (Ritchie et al., 2015). The adjusted p-value for multiple testing was calculated using the Benjamini-Hochberg correction. Genes with an absolute log₂ (fold change) greater than one and false discovery rate (FDR) < 0.05 were identified as DEGs between the two subtypes.

We downloaded the “CIBERSORT” scripts (https://cibersort.stanford.edu/) to estimate the immune composition of patients with CCA using a normalized express matrix (Chen et al., 2018). Immune, stromal, and tumor purity scores were calculated using the “estimate” R package (Yoshihara et al., 2013).

Differential DNA methylation analysis was conducted on both normal and tumor tissues, as well as on the two tumor subtypes. Data were pre-processed using the “minfi” package (Aryee et al., 2014). Hypomethylation probes with β < 0.5 hypo-methylation probes were used to conduct the analysis. We considered CpG probes to be hypo-methylated if they met the following criteria (Song You et al., 2018): β < 0.5 (Tirosh et al., 2016); M-value was significantly different between the two subtypes (q-value < 0.05); and (Stuart et al., 2019) mean β-value difference (Δβ) of >0.2 between the two groups. Copy number profiles were derived from the signal intensity values of methylated and unmethylated probes using the “conumee” package in R (http://bioconductor.org/packages/conumee). Copy number variations (CNVs) were derived from the log2-ratios of the tumor samples to the average value of matched normal tissues.

All computational and statistical analyses were performed using R software (https://www.r-project.org/). The unpaired Student’s t-test was used to compare two groups with normally distributed variables, while the Mann-Whitney U-test was used to compare two groups with non-normally distributed variables. Survival analysis was performed using the Kaplan–Meier “survival” R package. The log-rank test was used to determine whether the survival curves were significantly different. A p-value < 0.05 was considered statistically significant.

We used the least absolute shrinkage and selection operator (LASSO) method with 10-fold cross-validation and the Cox proportional hazards model with Akaike information criterion (AIC) selection criteria to build a nine-gene prognostic risk model based on the 54 TF-CIGs using the “glmnet” R package. The tumor mutation-burden-related signature (TMBRS) was calculated using Eq. 11: T M B R S = ∑ ( β i × E X P i ) (11) where β_i can be derived from multivariate Cox analysis. The prognostic accuracy of the classifiers in the training and testing sets was evaluated using the Kaplan–Meier curve and log-rank test.

We predicted the IC₅₀ values for the CCAs in GSE89749 via OncoPredict software using the 54 genes obtained from the eLBP algorithm (Maeser et al., 2021). In this step, we adopted the gene expression matrix from Genomics of Drug Sensitivity in Cancer (https://www.cancerrxgene.org/) and acquired drug sensitivity scores for the two subtypes. Next, we conducted a differential analysis of drugs between the two subtypes using the Wilcoxon rank-sum test. Drugs with p-values < 0.001 and drug sensitivity scores <1 were considered potential candidates.

A schematic of the study design is shown in Figure 1A. To investigate the composition of the TME in CCA, we reanalyzed the single-cell data from GSE138709 (Zhang et al., 2020). We deconstructed the cell types in tumors and adjacent tissues. We observed evident differences between the two tissues. We obtained 16 cell types in the adjacent tissues, including CD8_GZMA, CD8_MKI67, and CD69-positive T cells, M1-like macrophage subtypes (CD68_C1QC, CD68_S100A9), and two DC subtypes (Figure 1B; Supplementary Figure S1A). GSEA showed that cytokine production and inflammatory responses were enriched in M1-like macrophage subtypes (CD68_C1QC and CD68_S100A9). Additionally, T cell activation and T cell-mediated immunity pathways were enriched in CD8_GZMA cells (Figure 1B). In tumor tissues, we obtained three T cell subtypes, including memory-like T cells (CCR7+ and CCR7_T) and exhausted CD8 (TIGIT+ and TIGIT_CD8). We also obtained M2-like macrophages (VEGFA+ and VEGFA_MACRO; Figure 1C; Supplementary Figure S1A). GSEA results showed that pathways including cell communication, cell motility, and response to stimulus were enriched in VEGFA_MACRO, while leukocyte activation and T cell differentiation were enriched in the TIGIT_CD8 subtype (Figure 1C). Interestingly, the sample distribution results showed that T cells and macrophages were mostly from samples 23T and 24T. In contrast, the cells from 18T to 20T samples were almost all epithelial cells (Methods; Supplementary Table S2). On this basis, the tumor samples in this dataset were divided into high- and low-TIL tumor subtypes. These results indicate that T cells and macrophages perform different functions in tumor and para-cancerous tissues.

To explore the mechanisms that induce the establishment of different immune states between normal and tumor tissues, we developed the eLBP algorithm to evaluate the correlation among immune-infiltrating, storm, and epithelial cells in normal and tumor single-cell data by calculating the intensity scores among the cell types (Figure 2A; Methods). First, we analyzed 7,071 activated cell type-specific pathways in normal tissues and 8,096 pathways in tumor tissues, in which 813 ligand or receptor genes in tumor tissues and 559 in normal tissues were computed based on the CellChat dataset (Supplementary Table S3). Second, we constructed protein-protein interaction (PPI) networks and calculated the cosine similarity between the PPI genes in each cell type. The degree of ligand/receptor (L/R) genes in the network is shown in Supplementary Figure S2A. The expression of ligand and receptor genes in CAF, CD68_C1QC, DC_RAMP3, and END was higher in normal tissues. T, 24T, TIGIT-CD8, and VEGFA_MACRO had higher degrees in the tumor tissues, indicating these cell types were in hyperactive status in cell communication among the TME (Supplementary Figure S2B). The levels of ligand and receptor genes indicate that macrophages with high levels of VEGFA play an important role in tumor tissues. Next, we used the eLBP algorithm to calculate the interaction probability of each ligand and receptor (Supplementary Table S3). Using this probability, we obtained the interaction intensity scores for each L/R pair among the cell types in the TME (Supplementary Tables S4, S5). The results demonstrated significantly different communities in tumors and normal tissues. In normal tissue, cells from DC_RAM3 and CD68_C1QC interacted with activated CD8 cells from CD8_GZMA according to NAMPT and chemokines, respectively (Figure 2B). In contrast, there were high levels of TIL infiltration in the tumor tissues. In addition, CAF could interact with memory-like T cells in CCR7_T, tumor cells in T_24T, and macrophages from VEGFA_MACRO. Moreover, VEGFA_MACRO interacted with exhausted CD8 + cells from TIGIT_CD8 (Figure 2C). The interaction genes between CAF and others are mainly related to cell adhesion, such as glycoprotein CD44 and collagen genes (Figure 2C). Meanwhile, VEGFA_MACRO, with a high level of LGALS9, interacted with TIGIT and HAVCR2 (Figure 2C). These results illustrate that VEGFA_MACRO plays an important role in CD8 exhaustion.

Subsequently, we analyzed the key TFs that promote the expression of ligand genes in normal and cancerous tissues using the eLBP algorithm (Methods). In normal tissues, CD68_C1QC cells interacted with CD8_GZMA via CXCL12/CXCR4 (Figure 2B). We next obtained seven TFs (STAT1, EGR1, KLF4, and JUND; Figure 2D; Supplementary Table S6) from ENCODE and one TF (KLF2; Supplementary Table S6) from ChEA3. Moreover, KLF4 and STAT1 regulated the ligand gene CXCL12 (Figure 2D), which may promote the expression of ligand genes and influence cell signaling. In addition, cells in DC_RAMP3 interact with those in CD8_GZMA via CLEC1B/KLRB1 and TF STAT3 in DC_RAMP3 could regulate CLEC1B (Supplementary Table S6). Notably, in the tumor tissues, we found that only the cells in VEGFA_MACRO interacted with those in TIGIT_CD8 via LGALS9/HAVCR2 (Figure 2B). We then obtained 18 TF regulons from the three datasets, in which IRF1 was enriched in both ENCODE and TRANSFAC, and SPI1 was enriched in both ENCODE and ChEA3 (Supplementary Table S6). Moreover, these two TFs also regulated the ligand gene LGALS9 (Figures 2E,F). In addition, the profiles of LGALS9, HAVCR2, IRF1, and SPI1 in CCA patients from another single-cell dataset GSE151530 also verified their functions related to TAM and exhausted T cells (Supplementary Figure S3) (Ma et al., 2021). Moreover, immunohistochemical analysis further confirmed the high expression of HAVCR2 in human CCA patients and the co-expression of LGALS9 and these two genes (Figure 3A). The significantly low expression of LGALS9 in siRNA-IRF1/SPI1 VEGFA-positive macrophages further validated the transcript regulon of SPI1/IRF1-LGALS9 (Figures 3B,C). Finally, we derived the gene trimmed-PPI network from the mutually connected cell types in which L/R genes participate in normal and tumor tissues via eLBP (Supplementary Figure S4).

Concerning the potential function of VEGFA_MACRO to induce the expression of exhausted genes such as HAVCR2, we further overlapped the genes in TF regulons with those in ligand or receptor-interaction trimmed-networks from eLBP from VEGFA_MACRO and obtained 54 genes (Supplementary Figure S4C; Supplementary Table S7). We then investigated the clinical features of patients with CCA using these genes. We grouped patients into different subtypes in two bulk RNA datasets using consistent clustering (Methods; GSE89749; GSE76297; Supplementary Table S1) (Chaisaingmongkol et al., 2017; Jusakul et al., 2017). Patients in GSE89749 were grouped into two subtypes (S1 and S2; Figure 4A; Supplementary Figure S5). Moreover, we found that the samples from the S2 group had a significantly better prognosis (p = 0.027; Figure 4A). Notably, the tumor purity results showed that patients in S2 had higher infiltrative immune and stromal scores, but lower tumor purity (Figure 4B). The immune infiltration results showed that patients in S2 had a higher ratio of CD8, activated memory CD4, S1, and M2-like macrophage subtypes (Figure 4C). In particular, CD8A, CD8B, CD68, HAVCR2, and LGALS9 were highly expressed in S2 (Figure 4D). On this basis, S2 was classified as the inflammatory “hot tumor” subtype and S1 as the “cold tumor”. We further identified the clustering results from other RNA sequencing datasets (GSE76297). The patients in S2 from GSE76297 had higher expression levels of LGALS9, IRF1, SPI1, and immune genes and similar TIL profiles with the GSE89749 data (Supplementary Figure S6). These results further confirmed the existence of high-TIL tumor subtypes in CCA and indicated the exhaustion of CD8 cells.

Next, we analyzed the heterogeneity of somatic frequencies across the two CCA subtypes. The results showed that the cell replication genes TP53 and genes in the PI3K/Akt pathway, such as BAP1, KRAS, EPHA2, ARID1A, and SMAD4, exhibited high mutation ratios in the two subtypes (Figure 4A). S1 had high ratios of ADAMTS20, KMT2C, and APC (Figure 5A, 13%, 13%, and 11%, respectively). S2 had a higher mutation frequency of IDH1 and MUC16 (Figure 5A, 14% and 12%, respectively). Tumor cells with TP53 mutations are generally identified as being more immunogenic (Chasov et al., 2020), and we found that immune genes, such as CD3D and HAVCR2, were highly expressed in the MUC16 mutated group. Subsequently, we investigated the heterogeneity of methylation profiles between normal tissues and tumors. We obtained 15,520 differential probes between normal and tumor tissues, including 5,219 downregulated and 10, 321 upregulated probes (Methods, Figure 5B). In normal tissues, hypomethylated genes were mainly enriched in pathways such as T cell differentiation, cell adhesion, and inflammatory response, which were also enriched in normal single-cell datasets (Figure 5B). This indicated a high immune infiltration and activated immune status of the normal. Among the two CCA subtypes, we obtained 697 downregulated probes and 1,115 upregulated probes in S2. In addition, there were 684 downregulated hypomethylated probes enriched in 271 genes in S2. These genes were enriched in MAPK, CAMP, cGMP-PKG, and PI3K-AKT pathways (Figure 5C). We then conducted CNV analysis using the methylation data. The results showed that significant CNV signaling occurred across chromosomes 1, 6, and 8 in both subtypes (Figure 5D). Interestingly, we found that S2 had more CNV regions on chromosomes 6 and 8 (Figure 5D). We then analyzed the hypomethylated genes on chromosomes 6 and 8 in S2. These genes were enriched in antigen presentation and immune-toxic suppression, as well as T cell differentiation and G protein-coupled receptor-associated signaling pathways (Figure 5E).

LASSO-COX regression was used to construct the eLBP-COX TME risk score (TMRS) using the following formula in dataset GSE89749: TMRS = −0.03085 × IRF1 − 1.03491 × RHOA + 0.5868 × PLAUR − 0.2188 × NCF2 − 0.2273×CXCR4 − 0.1201 × HCK − 0.1773 × LYZ − 0.1244 × RGS1 + 0.2259 × VCAN (Methods; Supplementary Figure S7). We found that patients with low TMRS had a better prognosis with higher immune scores and presented with IRF1, LGALS9, and other immune genes in GSE89749 (Figures 6A,B), indicating an effective indicator of RS score in prognostic risk prediction. Furthermore, we verified the RS model in two other datasets (EMTAB-6389; GSE107943 Supplementary Table S1; Methods; Figure 6A) (Ahn et al., 2019). Meanwhile, a nomogram was drawn to visualize the results of the eLBP-COX regression analysis (Figure 6C). Subsequently, we conducted drug susceptibility analysis and explored seven drugs (AZD7762, BI 2536, CDK9 5576, gemcitabine, mitoxantrone, teniposide, and topotecan) with IC50 values that were significantly lower in S2 than in S1 (Methods; Figure 6D). Specifically, BI 2536, teniposide, and mitoxantrone were used as potential therapeutic drugs for HCC that could be used in the treatment of CCA.

To take part in the Resource Identification Initiative, the corresponding catalogue numbers were listed in Table 1.