We utilized the R package TCGAbiolinks (15) to download the 36 CHOL RNA-seq data [log2(FPKM+1)] and relevant clinical information from the TCGA database (https://portal.gdc.cancer.gov/). The 36 samples with both available RNA-seq data and clinical information were used as the training set for hypoxia-related markers. The GSE107943 dataset was obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which were used as the external validation set, containing available microarray expression data and clinical data of bile duct cancer samples. The microarray data preprocessing process was as follows: remove null probes and probes corresponding to multiple genes at the same time, and use the median expression value of these probes as the expression value of the gene when multiple probes correspond to the same gene.

The 200 potential hypoxia-related genes were retrieved from MSigDB’s HALLMARK HYPOXIA pathway (https://www.gsea-msigdb.org/gsea/index.jsp) (v7.4 version). All prognosis-associated hypoxia genes were screened using univariate Cox hazard analysis, with a p-value of 0.1 chosen as the cutoff. Utilizing the R package glmnet, the Lasso Cox regression model was built using all prognosis-related hypoxia genes, and candidate genes were penalized to exclude redundant components and reduce the risk of overfitting (16). A 10-fold cross-validation under the constraint of low bias is used to calculate the model’s penalty parameter (λ).

To determine the best prognosis-associated hypoxia signature, Lasso regression analysis was used, and the regression coefficients were calculated for each gene separately, and then the sample risk scores for multivariate Cox regression were calculated using the R package survminer (17), with the formula as follows:

βi denotes the weight of each gene in Lasso Cox regression; χi denotes the expression level of each gene.

The CHOL patients were divided into high-risk and low-risk groups based on Risk score. Kaplan-Meier curves of survival differences between samples in the high-risk and low-risk groups were plotted, and significant differences in OS between groups were assessed based on log-rank tests.

The univariate and multivariate Cox hazard regression analysis were performed separately for risk scores and patient clinical characteristics (including age, gender, and stage) to determine whether the predictive power of the prognostic model was independent of clinical characteristics. Risk ratios (HR) and 95% confidence intervals were calculated for each candidate prognostic factor, with a p-value < 0.05 as the threshold of significance. In addition, differences in risk scores across age, sex, and stage subgroups were assessed based on the Wilcoxon rank-sum test.

The expression correlation between genes in the prognostic model was assessed based on Pearson correlation analysis. In addition, the interactions of proteins encoding hypoxia-related genes were obtained based on the STRING database (http://www.string-db.org/), and protein-protein interactions (PPI) network of hypoxia-related genes were constructed, and the network graphs were visualized using Cytoscape (v3.9.1) to assess the interactions of the genes in the prognostic model.

Gene set enrichment analysis (GSEA) was used to identify biological processes that differed significantly between high- and low-risk group samples (18). We examined biological processes from the GO database and KEGG signaling pathways for significant variations between high- and low-risk group samples using the R tool ClusterProfiler (19). To screen the collection of substantially linked genes, a corrected p-value of less than 0.05 was employed as a criterion, and the p-value correction technique Benjamini Hochberg was applied.

TCGA provided the tumor mutational data. The mutation genes were calculated using tumor mutational burden (TMB). The mutational data in both the high-risk and low-risk groups were analyzed using the R maftools package (20). P-values less than 0.05 were deemed statistically significant.

The cancer immune cycle has emerged as the primary research framework for cancer immunotherapy. It describes a cyclical process in which the immune system eradicates cancer: Antigen Processing and Presentation (Category 1), Atimicrobials (Category 2), Chemokines (Category 3), Cytokine Recepors (Category 4), Cytokines (Category 5), Matrix remodeling (Category 6), and Nature Killing Cell Cytoxicity (Category 7) (21). In the low and high hypoxia risk groups, we looked at the expression of genes that negatively regulate these processes. Tracking Tumor Immunophenotype (http://biocc.hrbmu.edu.cn/TIP/index.jsp) was used to identify immunonegative regulation-related genes (22).

From September 2021 to March 2022, 10 paired fresh cholangiocarcinoma tissues and neighboring non-tumor tissues were taken from patients who underwent radical surgery at Xi’an Jiaotong University’s First Affiliated Hospital. In addition, between December 2013 and April 2016, 62 patients with cholangiocarcinoma who had a curative procedure were enrolled in a cohort. The protocol was authorized by the First Affiliated Hospital of Xi’an Jiaotong University’s Institutional Medical Ethics Committee (No. XJTU1AF2021LSK-261).

The RNeasy Mini Kit (Qiagen) was used to extract total RNA, and the reverse transcription was done with the High Capacity cDNA Reverse Transcription Kit (Thermo Fisher Scientific). A LightCycler480 PCR equipment was used to perform quantitative PCR using SYBR Green dye (Roche). GAPDH was used to standardize the relative fold change in expression.

Following deparaffin deparaffinization, rehydration and antigen retrieval, the paraffin-embedded slides were incubated with primary anti-PPFIA4 antibody (1:100, Sigma) at 4°C overnight. HRP-conjugated secondary antibody and DAB peroxidase substrate were utilized for immunostaining.

All statistical analyses were performed using R software (version 4.1.2; https://www.R-project.org). Wilcoxon rank-sum test was used to compare significant differences in continuous variables between the two groups. Pearson correlation analysis was used to analyze the expression correlation between genes. Differentially expressed genes between high and low risk groups were identified based on the R package limma (23) with a differentially expressed gene threshold of |log2FC| > 0.58 and p-values < 0.05. The univariate Cox regression analysis was performed using the R package survival, and multi-factor Cox regression analysis was performed on prognostic model scores to explore the the independent prognostic value of the score. Kaplan-Meier curves were plotted using the R package survminer, and the significance of differences between groups was assessed based on the log-rank test. And p-values < 0.05 were used as the threshold of significance for statistical analysis.

We acquired 200 hypoxia-related genes from the MSigDB database to study the predictive efficiency of hypoxia-related genes in CHOL. EGFR, PFKL, CHST2, ADORA2B, EDN2, HMOX1, CP, PPFIA4, PLAC8, IDS, ALDOC, STBD1, BCL2 were among 13 hypoxia-related genes found to be substantially linked with overall survival (OS) in the TCGA CHOL cohort using univariate Cox hazard analysis (p-value < 0.1, Figure 1A ). Based on the median value of gene expression (EGFR, PFKL, CHST2, ADORA2B, EDN2, HMOX1), we separated the samples into two groups, and Kaplan-Meier curve analysis revealed that gene expression affected patient survival (log-rank test p-value < 0.05, Figure 1B ).

Although we initially observed genes with prognostic efficacy in CHOL patients by univariate Cox regression analysis and log-rank test, we further performed Lasso regression analysis based on these 13 genes in order to integrate these informative genes to obtain a more optimal prognostic model, while removing redundant factors and controlling the risk of overfitting. Ten-fold cross-validation was performed under optimal conditions to determine the penalty parameter (λ) of the model. We screened out 6 most predictive factors for OS ( Figures 1C, D ). PPFIA4 and CP were retained as a valid prognostic risk factor, while BCL2, PLAC8, PFKL and EGFR were retained as a prognostic protective factor, Together, they constituted a prognostic risk model associated with hypoxia in CHOL. Based on the linear combination of the expression levels of these genes and the corresponding weights, we could assess the hypoxia-related prognostic risk score for each patient ( Figure 1E ). The corresponding association between the expression of 6 hypoxia-related signature in the prognostic model and the vital status and survival of patients were visualized in Figure 1F .

The risk score for each patient with CHOL in the TCGA was determined using the hypoxia-related prognostic model for CHOL that was constructed. After obtaining each patient’s risk value, the samples were sorted from low to high in terms of risk score, and scatter plots were used to display the results. The median of all sample risk scores was used as the cutoff, and the samples were divided into two groups: high risk and low risk, the low-risk group’s OS was likewise considerably lower than the high-risk group’s (log-rank test p-value < 0.001, Figure 2A ). We also calculated heat maps for display to investigate the expression of different genes in the risk model across all genes in different samples, and we discovered that the expression of six hypoxia-related genes in the model differed considerably in the high- and low-risk groups ( Figure 2B ).

To assess the robustness and generalizability of the hypoxia-related prognostic model, we collected CHOL samples from GSE107943 as an independent validation cohort. Patients were divided into high-risk and low-risk groups according to the risk score of the prognostic model. We obtained similar results to the training cohort through an external validation cohort ( Figures 2C, D ).

We then explored the independence of the prognostic model in patients with CHOL from TCGA. The univariate Cox regression analysis revealed that the hypoxia-related risk score was the only prognostic factor of OS (HR = 5.94, p-value < 0.001, Figure 3A ). Further, the multivariate Cox regression analysis confirmed the independent prognostic value of the hypoxia-related risk score (HR = 74.71, p-value = 0.004, Figure 3B ).

In the GEO cohort, univariate Cox analysis showed that hypoxia-related risk score (HR = 2.65, p-value = 0.046), stage (HR = 3.48, p-value = 0.015) were prognostic factors in patients with CHOL ( Figure 3C ). Hypoxia-related risk score (HR = 3.06, p-value = 0.032) and stage (HR = 3.88, p-value = 0.017) were found to be independent predictive markers in patients with CHOL in a multifactorial Cox analysis ( Figure 3D ), which is consistent with the TCGA cohort’s findings. The prognostic risk scores established in our study might be used as independent risk factors affecting patient prognosis based on the above investigations.

To further corroborate the independence of prognostic risk scores from clinical characteristics, we compared the differences in the distribution of risk scores by grouping patients based on their age, gender and stage, respectively. We found no significant differences in risk scores by age and tumor stage grouping (p-value > 0.05, Figures 3E, G ), but slightly higher in male patients than in female patients (p-value < 0.05, Figure 3F ).

Although the 6 signature genes we obtained from the prognostic model were all involved in the tumor hypoxic state hallmark, but their interactions with other hypoxic genes were not clear, so we examined the co-expression of these genes. Overall, there was no strong co-expression among the hypoxia-related genes involved in the prognostic model (Pearson |r| < 0.34, Figure 4A ), indicating that the redundant genes with strong associations had been filtered out through LASSO Cox analysis. In addition, we looked at how the hypoxia-related genes in the model interacted with other hypoxia-related genes. We built a PPI network connecting the genes in the prognostic model and other hypoxia-related genes using the interactions of gene-encoded proteins retrieved from the String database online tool ( Figure 4B ).

We found extensive interactions between EGFR, CP, PFKL genes and other hypoxia-associated genes in the model. The model includes hypoxia-related genes, which might play a key role in tumor hypoxia regulation.

To explore the underlying molecular mechanisms of hypoxia-related prognostic models, we first examined somatic mutations in different risk groups in the TCGA-CHOL cohort. Due to the small sample size of the TCGA-CHOL dataset, the comparison of differences in mutation frequencies failed to present significant results ( Figure 5A , p-value < 0.05). However, we found four patients in the high-risk group enriched for the IDH1 mutation, while only one patient in the low-risk group was enriched for this mutation ( Figure 5A ).

Using GSEA enrichment analysis, we looked into the biological pathways altered by the hypoxia-related predictive risk score in subgroups. Metabolic processes such as “Carbon metabolism”, “Glycolysis/Gluconeogenesis”, “Biosynthesis of amino acids”, “Cholesterol metabolism”, “PPAR signaling pathway” and “Fatty acid degradation” were dramatically increased in the high-risk group (NES > 0, FDR < 0.0001), whereas T cell receptor signaling route, natural killer cell mediated cytotoxicity, antigen processing and presentation, chemokine signaling pathway, and cytokine-cytokine receptor interaction were all significantly downregulated in the high-risk group (NES < 0, FDR < 0.0001, Figure 5B ). In addition, GSEA enrichment analysis of GO-BP showed that metabolic processes such as “protein-lipid complex subunit organization”, “protein-lipid complex subunit organization”, “protein-lipid complex remodeling”, “protein-containing complex remodeling”, “organic acid catabolic process”, “requlation of lipid catabolic process”, “carboxylic acid catabolic process”, “neutral lipid metabolic process”, “sterol metabolic process”, “small molecule catabolic process”, “cellular amino acid metabolic process”, “secondary alcohol metabolic process”, “steroid metabolic process”, “organic acid biosynthetic process”, “lipid catabolic process” were increased dramatically in the high-risk group (NES > 0, FDR < 0.0001). We observed that “leukocyte migration”, “leukocyte mediated immunity”, “regulation of cell-cell adhesion”, “mononuclear cell migration”, “I-kappaB kinase/NF-kappaB signaling”, “adaptive immune response”, “leukocyte cell-cell adhesion”, “interleukin-12/10 production “, “cell activation involved in immune response”, “antigen processing and presentation”, “positive regulation of cell activation” were significantly downregulated in the high-risk group (NES < 0, FDR < 0.0001; Figure 5C ). Taken together, we deciphered the tumor molecular signature of the differences between hypoxia-related risk subgroups, suggesting potential molecular mechanisms for the different risk subgroups.

Tumor-immune interactions are crucial in controlling the immune milieu because they incorporate the effects of several regulators. We looked examined how immune regulatory genes were expressed differently in high-risk and low-risk people. We discovered that genes involved in immune modulation were expressed differently in high and low risk groups ( Figure 6A ). We looked examined how divergent genes were expressed in distinct immune regulatory groups.

The genes HLA-DPB1, HLA-DQA1, HLA-DRA, HLA-DRB1 and HLA-DRB5 in the “Antigen Processing and Presentation” category were considerably downregulated in the high-risk group. PLA2G2A, KNG1, HRG, and other genes associated to “antimicrobials” were considerably downregulated in the high-risk group. CCL17, a gene associated to “chemokines,” was considerably downregulated in the high-risk group. The gene EGFR, which belongs to the “Cytokine Receptors” category, was considerably downregulated in the high-risk group. The genes CTF1, GDNF, HAMP, IGF2, SAA1, SAA2, TGFA, which belong to the “cytokines” category, were considerably downregulated in the high-risk group. CA9 and PLOD2, which are connected to “matrix remodeling,” were considerably downregulated in the high-risk group. The gene PLCG2, which is connected to “Natural Killer Cell Cytotoxicity,” was considerably downregulated in the high-risk group. The above results are shown in figures ( Figures 6A, B ). These findings show that immune checkpoint expression differs under the effect of CHOL’s hypoxic condition, affecting the immunological microenvironment’s control.

To further corroborate these findings, PPFIA4 mRNA levels were measured in 10 matched fresh CHOL tissues with surrounding non-tumor tissues. When compared to non-tumor tissues, the amount of PPFIA4 mRNA in cholangiocarcinoma tissues became much higher ( Figure 7A ). In addition, we investigated PPFIA4 expression in a cohort of patient with cholangiocarcinoma samples from our hospital by IHC assay, which demonstrated predominant expression of PPFIA4 within tumor cells ( Figure 7C, D ). The expression of PPFIA4 was found in 26 of the 62 carcinoma samples (41.9%), whereas the remaining 36 samples (58.1%) were negative. Table 1 outlined the link between PPFIA4 expression and clinicopathological characteristics. The presence of PPFIA4 was significantly linked to large tumor size, advanced TNM stage, and lymph node metastasis, while there was no correlation of PPFIA4 expression with gender, age, differentiation and tumor stage. Additionally, individuals with high PPFIA4 expression had a poorer 5-year OS prognosis ( Figure 7B , p-value = 0.03).