This research was carried out strictly in accordance with the predetermined process. For details, please refer to Supplementary Figure 1. The patient cohorts and datasets used for the analysis include the following parts: (1) The TCGA-CHOL (The Cancer Genome Atlas Cholangiocarcinoma) dataset, derived from the UCSC Xena platform (12). A total of 45 patients with cholangiocarcinoma were included in their transcriptome sequencing data, combined with their publicly available clinical information, including age (range 32–85 years, median 62 years), gender (23 males and 22 females), survival time and survival status. The inclusion criteria are as follows: a clear pathological diagnosis of CHOL, complete gene expression and clinical follow-up information; The exclusion criteria were: cases with missing key clinical data or follow-up data. Ultimately, 36 patients met the analysis criteria. The patients in this cohort were mainly from North America and Europe. (2) Datasets such as TCGA-LIHC (hepatocellular carcinoma), GDC TARGET-OS (osteosarcoma), and IMvigor210 (immunotherapy for bladder cancer) were respectively used for the analysis of pan-cancer HRD. They were all downloaded from public databases. For specific inclusion and exclusion criteria, please refer to Supplementary Table 1. (3) The external validation cohort included E-MTAB-6389 (from the Bio-Studies database) (13) and GSE107943 (from the GEO database) (14), with a total of 75 patient samples of cholangiocarcinoma included. All cases were patients with a clear diagnosis and complete clinical information, and those with other malignant tumors or lacking follow-up data were excluded. The specific sources, enrollment times, geographical distributions, and demographic characteristics of each dataset (for details, see Supplementary Table 1) are provided to enhance the reproducibility and data transparency of the study.

The HRD score was defined as the unweighted sum of the loss of heterogeneity (LOH) (15), telomere allele imbalance (TAI) (16), and large-scale transition (LST) scores (17, 18). Neoantigen load, predicting peptides binding to major histocompatibility complex (MHC) proteins, derived from RNA-seq HLA typing. Neoantigen burden was quantified by single nucleotide variants (SNVS) and insertion/deletion (indel) mutations. Values for HRD, neoantigen load (SNV and Indel), and mutation rate were collated from the pan-cancer Atlas study by Thorsson et al. (19). We compared genome alteration fractions, mutation counts, and MSI sensor scores between the top and bottom 20% HRD score patients in CHOL. TCGA-CHOL patients were stratified into high and low HRD subgroups by median HRD score (Supplementary Table 2).

We used the Limma (V3.56.2) R package (20) to analyze the differences between HRD high-risk and low-risk patients. We set cutoffs at log2fold change > 1 and P < 0.05. Genes with log2 fold change > 1 were considered highly expressed in the HRD group, while those with log2 fold change < -1 were considered low-expressed. We visualized these differences using volcano plots, ranking maps, and a heat map. Spearman’s correlation analysis was performed to explore gene correlations.

To construct a robust prognostic model for CHOL-HRD, we followed a multi-step approach: (1) First, we integrated 10 classical algorithms: Random forest (RSF), least absolute shrinkage and selection operator (LASSO), gradient boosting machine (GBM), Survival support vector machine (survival-SVM), supervised principal component (SuperPC), ridge regression (ridge), Cox Partial least squares regression (plsRcox), CoxBoost, Stepwise Cox, and elastic network (Enet). Among them, RSF, LASSO, CoxBoost, and Stepwise Cox have dimensionality reduction and variable screening functions, leading to 37 machine-learning algorithm combinations. (2) Using TCGA-CHOL as the training cohort, we applied these algorithm combinations to screen key genes and construct a prognostic model based on previously identified feature genes. (3) We validated the model using two test cohorts (E-MTAB-6389 and GSE107943), calculating the CHOL-HRD risk score based on features obtained from the training cohort. Based on the average C-index of the two test cohorts, we selected the best prognostic model for CHOL-HRD and calculated the risk score. Patients were divided into high-risk and low-risk HRD groups based on the median score. Survival and multivariate Cox analyses were used to evaluate the independent prognostic significance of the risk score.

To further evaluate the predictive value of the HRD Score, we compared it with previously reported CHOL biomarker prognostic models (21–25). The predictive power of the other biomarkers was compared with that of the HRD Score in the two test cohorts (E-MTAB-6389 and GSE107943).

Given the importance of HRD in cancer development, we explored whether the HRD Score has a stable prognostic effect on other types of tumors. Based on the overall survival of patients in the TCGA-LIHC, TCGA-OV, and TARGET-OS datasets, we used the survival R package (https://CRAN.R-project.org/package=survival) for survival analysis to verify the prognostic efficiency of the HRD Score.

We performed Gene Set Enrichment Analysis (GSEA) using the clusterProfiler 4.8.2 R package (26), calculating normalized enrichment scores for each gene set. The gene set “c2.cp.kegg.symbols.gmt” from the MSigDB database was used for GSEA to evaluate the impact of high- and low-risk groups on tumor KEGG pathways, with FDR < 0.25 considered significant.

To analyze the Single Nucleotide Polymorphisms (SNPs) in different risk score groups of TCGA-CHOL patients, maftools 2.16.0 package was used to analyze the frequently mutated genes in the high and low risk groups. To evaluate the interactions between drugs and gene mutations, the drugInteractions function was used for analysis. The goal was to identify genetic mutations associated with susceptibility or resistance to specific drugs. In addition, biological oncogenic pathway analysis was performed on the mutation data to understand which biological oncogenic pathways are affected by gene mutations. We used the OncogenicPathways and PlotOncogenicPathways functions to achieve this goal. Meanwhile, we calculated the TMB data of CHOL patients through the maftools 2.16.0 R package. CHOL, patients with MSI - Sensor data is from cBioportal database (https://www.cbioportal.org).

The CIBERSORT (https://cibersort.stanford.edu/) algorithm evaluated immune cell infiltration in different CHOL sample datasets. Differences in immune cell infiltration among disease subgroups were tested using the Wilcoxon test, with p < 0.05 considered statistically significant. We also assessed the differential expression of Immune Checkpoint Proteins (ICP) and Immunogenic Cell Death (ICD) modulators between high- and low-risk groups. The “estimate” R package was used to analyze differences in tumor immune score, stromal score, and tumor purity.

We used the Cancer Treatment Response Portal (CTRP, https://portals.broadinstitute.org/ctrp/) and the PRISM dataset to identify potential drugs for CHOL patients in the high-risk HRD group. Following the protocol by Yang et al. (27) 2021: (1) We obtained drug sensitivity data for CCLs from the CTRP and PRISM datasets and expression data from the Cancer Cell Line Encyclopedia (CCLE) database. (2) The CTRP and PRISM datasets provided AUC values, where lower AUC values indicate increased sensitivity to the compound. (3) Using the Wilcoxon rank-sum test, we performed a differential analysis of drug response between high- and low-risk groups based on the HRD Score. We set a threshold of log2FC > 0.03 to identify compounds with low AUC values in the high-risk HRD group. (4) Next, we used Spearman’s correlation to further screen compounds with a negative correlation coefficient between the Area Under the Curve (AUC) value and the HRD Score (threshold R < -0.1). (5) Finally, we identified potential drugs for high-risk HRD patients by intersecting the compounds obtained from steps (3) and (4). Furthermore, we used CMap to validate the results obtained from the CTRP and PRISM databases based on differential expression profiles.

To identify potential drugs with significant associations with CHOL, we performed a comprehensive analysis using CTRP, PRISM, and Cmap. We downloaded the molecular structures of the drugs from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and the target proteins from the Protein Data Bank (PDB; http://www.rcsb.org/). The connection between potential drugs and targets was visualized using CB-Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php).

This study also included clinical samples from 6 patients with primary intrahepatic cholangiocarcinoma (ICC) diagnosed by pathology. All of them underwent surgical resection in the Department of Interventional Oncology of Putuo Hospital affiliated with Shanghai University of Traditional Chinese Medicine from January 2023 to June 2023. The patients’ ages ranged from 63 to 78 years old (median 70.5 years old), including 4 males and 2 females. The inclusion criteria were: (1) ICC was clearly diagnosed by imaging and/or histology before the operation; (2) No history of neoadjuvant therapy, initial surgical resection; (3) No history of other malignant tumors or serious underlying diseases; (4) Postoperative pathology confirmed ICC, and both the tumor and the paired adjacent tissues could be obtained. The exclusion criteria are: (1) Previous or combined with other malignant tumors; (2) Receive radiotherapy, chemotherapy or immunotherapy/targeted therapy before the operation; (3) There are serious complications that affect the quality of the samples. All patients signed the informed consent form, and the sample collection was approved by the hospital ethics committee (Approval Number: PTEC-A-2024-7-1). Tumor tissues and paired adjacent non-tumor tissues were collected from each patient respectively. They were immediately frozen in liquid nitrogen after the operation and stored at -80°C for subsequent molecular detection and analysis. Detailed demographic and clinical information can be found in Supplementary Table 2.

Total RNA was extracted from cancerous and adjacent tissues using TRIzol reagent (Invitrogen) according to the manufacturer’s protocol. cDNA was synthesized with the RevertAid First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). Real-time PCR (RT-PCR) was performed using the Maxima SYBR Green/ROX qPCR Master Mix under the following thermal cycling conditions: 1) 95°C for 10 min; 2) 95°C for 15 s, then 60°C for 45 s, repeated for 40 cycles; 3) 95°C for 15 s; 4) 60°C for 1 min; 5) 95°C for 15 s; 6) 60°C for 15 s. Relative miRNA levels were analyzed using ABI Prism7300 SDS Software (Foster City, CA, USA). The primer sequences used for the analysis are shown in Supplementary Table 3.

All statistical analyses were performed using R software (version 4.3.1 and 3.6.0). The Wilcoxon rank-sum test was used to compare differences between two groups, while the Kruskal–Wallis test was applied to evaluate differences among more than two groups. Spearman’s correlation test was employed for correlation analysis. p < 0.05 (two sides) were considered statistically significant.

The high-risk HRD group of CHOL is characterized by genomic instability, with significant differences in genomic alterations and MSI observed in TCGA-CHOL. A longer OS indicates that the HRD score can serve as a reliable prognostic biomarker for CHOL. We performed 240 identified candidate genes through differential gene screening. The CoxBoost and Lasso combination, with the highest average C-index among 37 machine-learning algorithms, was selected as the best model and validated using two external datasets. Finally, multi-omics analysis, including genome, proteome, metabolome, and single-cell analysis, was performed for cancer screening, diagnosis, drug response, and prognosis prediction (Supplementary Figure 1).

To explore the correlation between HRD score and genomic instability, patients with CHOL were divided into high- and low-scoring groups based on the median HRD score. Significant differences were observed in genomic alterations (Figure 1A), somatic mutation count (Figure 1B), and MSI (Figure 1C) between the groups. Additionally, TAI (Supplementary Figure 2A), LST (Supplementary Figure 2B), and LOH (Supplementary Figure 2C) were higher in the high HRD score group. Survival analysis showed significant differences in OS (log-rank p = 0.043, HR (95%Cl): 2.0108 (1.415-2.858), high: 181, low: 182, Figure 1D) and PFS (log-rank p = 0.043, HR (95%Cl): 3.5163 (1.547-7.995), high: 42, low: 43, Figure 1E) between the groups, with worse prognosis in the high HRD score group, suggesting HRD score as a good prognostic indicator.

Identifying key HRD-related genes is crucial for developing CHOL predictive biomarkers and understanding the source of HRD. To identify these key genes, we analyzed the differences between the high and low-expression groups. We identified 240 candidate genes (Figure 2A). The distribution of the logFC values for these genes is shown in Figure 2B. We then focused on constructing the best prognostic model. The combination of CoxBoost and Lasso, with the highest average C index (0.724) among the 37 machine-learning algorithms, was selected as the final model (Figure 2C). In this combined model, 17 significant genes were screened, and 16 were selected for model construction after multivariate analysis. Correlation analysis of these 16 genes showed multiple strong correlations (Figures 2D, E).

Based on the C-index comparison of 37 machine-learning model combinations, the CoxBoost and Lasso algorithms were identified as the optimal combination for constructing an HRD-related prognostic prediction model for CHOL patients. From 240 key genes, 16 critical genes were selected to construct the LASSO prognostic model (Figure 3A). These genes’ median expression values were used to divide patients into subgroups for separate prognosis analysis. Seven genes (ANXA2P1, ATPB, BBOX1, KLHL33, MN1, OR51A4, and TRDN) exhibited significant prognostic differences (Supplementary Figures 3A–G). Using the penalty coefficients calculated by multivariate Cox analysis (Supplementary Table 4), we established a risk score by multiplying gene expression values by their corresponding coefficients. This allowed us to calculate the final risk score for each sample, dividing patients into high- and low-risk groups. A risk factor heat map was created based on patients’ risk scores and the expression values of the 16 key genes. In addition, a higher number of patients with high HRD scores died (Figures 3B, C).

The Kaplan–Meier curve showed a significantly worse prognosis for high-risk patients, with substantial survival differences between the high- and low-risk groups (Figure 3D). This finding was validated in two other datasets: TCGA-CHOL (Supplementary Figure 2D) and GSE107943 (Supplementary Figure 2E). Time ROC analysis showed that the area under the curve (AUC) for 1-, 3-, and 5-year survival predictions were 0.774, 0.976, and 0.922, respectively, demonstrating the HRD risk score model’s strong predictive performance (Figure 3E). Multivariate Cox Regression analysis showed that the HRD risk score serves as an independent prognostic factor for CHOL patients, similar to bile duct necrosis and infection (Figures 3F, G). According to these results, bile duct necrosis, infection, and the HRD risk score were incorporated into a nomogram for further prognosis prediction (Figure 3H). The diagnostic calibration curve indicated a good fit (Figure 3I).

To validate the prognostic ability of the HRD score model, we compared its performance with that of published CHOL prognostic models. The results showed that the HRD score had high and stable predictive performance in the validation sets E-MTAB-6389 and GSE107943 (Figure 4A). The average C-index of the HRD score model was significantly better than that of the other models (Figure 4B). HRD significantly affects tumor progression not only in CHOL but also in other types of tumors. To explore the generalization of the prognostic value of the HRD score, we evaluated its performance in the TCGA hepatocellular carcinoma and TARGET-OS (osteosarcoma) datasets. The Kaplan–Meier survival curves showed that the HRD score had significant prognostic differences not only in CHOL patients but also in hepatocellular carcinoma (Figure 4C) and osteosarcoma (Figure 4D). There were significant prognostic differences between high and low HRD scores.

To investigate differences in signal pathway enrichment between high- and low-risk HRD groups, we performed GSEA on both groups. The top five gene pathways significantly enriched in the high-risk HRD group were cardiac muscle contraction, drug metabolism, long-term potentiation, proximal tubule bicarbonate induction, and retinol metabolism (Figure 5A). In the low-risk HRD group, the top five enriched pathways were the chemokine signaling pathway, cytokine receptor interaction, graft versus host disease, hematopoietic cell lineage, and Toll-like receptor signaling pathway (Figure 5B).

We further evaluated the effect of HRD score on genetic variants, including SNPs and copy number variations (CNVs), in patients with CHOL. Analysis of single-nucleotide mutations in common driver genes during tumorigenesis showed a higher mutation rate in the high-risk group compared to the low-risk group (Figures 5C, D). From mutation data in CHOL patients, we identified mutations associated with susceptibility or resistance to specific drugs, revealing strong links between mutated genes and multiple proprietary pathways (Figure 5E). Furthermore, we identified cancer pathways enriched with mutated genes in CHOL patients, with the RTK-RAS signaling pathway showing the highest enrichment (Figure 5F).

A differential analysis of the TMB (Figure 5G) and MSI scores (Figure 5H) demonstrated significantly higher genomic instability in the high-risk group. Although statistically significant differences were observed in the mutation counts (Supplementary Figure 3H), no significant differences were found in the genomic alteration scores (Supplementary Figure 3I). Patients in the high-risk HRD score group showed higher mutation counts, indicating greater genomic instability.

Immune cells are crucial in the development and progression of CHOL. We evaluated the distribution of 28 immune cell types using ssGSEA and compared their abundance to determine differences in immune cell infiltration between high- and low-risk HRD groups. Immune cell infiltration was more pronounced in the low-risk group (Figures 6A, B). Notably, active B cells and regulatory T cells, which are vital in tumor immunity, showed higher infiltration levels in the low-risk group.

We also evaluated the differential expression of ICP and ICD modulators between the two groups. Most ICD modulators were significantly up-regulated in the low-risk group (Figure 6C). Among ICP modulators, some genes remained significantly up-regulated in the low-risk group (Figure 6D). This suggests that low-risk patients may be more responsive to immunotherapy, while high-risk patients may develop immune-tolerant subtypes. The ESTIMATE analysis was performed to examine the relationship between tumor immune score, stromal score, and tumor purity in both groups. The ESTIMATE score was significantly higher in the low-risk group (Figure 6E), indicating a lower tumor purity (Figure 6F). Both the immune score (Figure 6G) and stromal score (Figure 6H) were also higher in the low-risk group, consistent with ssGSEA results. Overall, the low-risk HRD group exhibited higher immune infiltration, suggesting that these patients are more suitable for immunotherapy.

To assess differences in immunotherapy responses between high- and low-risk HRD groups, we used the TIDE approach to predict immunotherapy responses, favoring the low-risk group (Figure 7A). To validate our results, we conducted differential analyses using the IMvigor210 bladder cancer immunotherapy dataset. The results indicated that the risk score was significantly lower in the CR/PR group, suggesting that low-risk patients are more likely to benefit from immunotherapy (Figure 7B).

To identify potential therapeutic drugs for CHOL patients, we used a comprehensive analysis combining CTRP and PRISM with Cmap. This analysis identified six CTRP-derived drugs (Figure 7C) and three PRISM-derived drugs (Figure 7E). The estimated AUC value of these drugs was negatively correlated with the HRD rating (Figure 7D) and was significantly lower in the high-risk groups (Figure 7F). Cross-referencing the CTRP and PRISM results, we identified two candidate compounds for the high-risk CHOL group: gemcitabine and panobinostat. The corresponding targets and mechanisms of action (MOA) are shown in Supplementary Table 5. The high sensitivity of CHOL patients to these compounds suggests their potential as therapeutic agents for the high-risk group.

We performed molecular docking analysis using CB - Dock2 (https://cadd.labshare.cn/cb-dock2/php/index.php) to examine the interactions between potential drugs and receptor molecules. This analysis included target identification and visualization of the interactions (Supplementary Figure 4).

In the prognostic model combining CoxBoost and LASSO machine-learning, ANXA2P1, BBOX1, KLHL33, MN1, OR51A4, and TRDN showed significant prognostic differences. We assessed the mRNA expression levels of these six genes in cancerous and adjacent tissues from six patients diagnosed with CHOL. RT-PCR results showed significantly higher mRNA expression of ANXA2P1, BBOX1, KLHL33, MN1, OR51A4, and TRDN in CHOL tissues compared to adjacent tissues (Supplementary Figure 5).