Computational analyses to reveal the key determinants of the high malignancy level of cholangiocarcinoma

Data

We have conducted comparative analyses of transcriptomic data of 1835 cancer tissue samples of 7 cancer types vs. those of 330 control tissues, all from the cancer genome atlas (TCGA) database via the genomic data commons (GDC) data portal (https://portal.gdc.cancer.gov).^[43] The detailed information of the datasets is given in Table 2.

For each of the following factors analyzed, we have used the expression levels of the marker genes as widely in the literature (see Supplementary Tables S1 and S2): (a) the level of hypoxia;^[44,45] (b) the level of dedifferentiation;^[46] and (c) the levels of SA biosynthesis and degradation, nucleotide de novo biosynthesis, and Fenton reactions in each of the seven cancer types.^[46]

De-batch effect

We have applied an empirical Bayesian method to estimate the batch effect by modeling the relationship between gene expression levels and batch variables.^[47,48,49] The estimated batch effect is then subtracted from the raw data^[50] to ensure that transcriptomic data collected separately can be compared directly with each other.

Single sample gene set enrichment analysis (ssGSEA) for pathway enrichment of individual samples

To perform gene set enrichment analysis (GSEA) on a single sample, we used the “gsva ()” function in R,^[51] with the analysis method set to “ssgsea”, which is specifically designed for GSEA on individual samples.^[52] The method is a non-parametric one that uses the empirical cumulative distribution function of the ranks of gene expressions within and outside a specified gene set to calculate the enrichment score. This method has been shown to be highly effective in identifying gene sets or pathways enriched by the given genes.^[53,54]

Estimation of SA accumulation on cancer cell surface

It is known that surface O-glycosylated proteins are involved in cancer metastasis.^[55,56] Hence, we have estimated the level of SA accumulation onto O-glycosylated proteins via sialyltransferases—that is, ST6GALNAC1, ST6GALNAC2, ST6GALNAC3, ST6GALNAC4, and STA8SIA6.^[57] It must be noted that only ST6GALNAC1, ST6GALNAC2, and ST6GALNAC4 are upregulated in digestive cancers; hence, these are the only ones considered in our analyses. In addition, NEU3 is the main SA degrader. We have used the following equation to estimate the average level of SA accumulation in stage-i samples for a target cancer type:

Where G₁ = {ST6GALNAC1, ST6GALNAC2, ST6GALNAC4}, G₂ = NEU3, and Ave_i (G, S) represents the expression of G averaged across all samples S in stage i. Further, the estimated level of SA accumulation up to stage I is expressed in the following manner:

where I = 1 for localized tumor samples, 2 for regional tumor samples, and 3 for distant tumor samples, assuming that the duration of each of the three stages is approximately the same.

xCell for estimating the infiltration levels by immune cells

xCell is a widely-used computational tool for calculating an enrichment score by the expressions of marker genes for each of the 64 immune and stroma cell types in a given cancer tissue sample.^[58] Using the calculated enrichment scores, we have estimated the infiltration levels of B cells, neutrophils, and macrophages in each sample by assessing each enrichment score against a total enrichment score for the tissue sample.

Mechanisms framework

We have also considered bile acids (BAs) in our analysis. BAs can activate the G-protein-coupled BA receptor GPBAR1, which is known to suppress cell proliferation, particularly in cancer cells.^[59,60] This has an important implication in our analyses, which follow a specific chain of inference, described hereafter. (1) All cancer tissues in TCGA harbor the following persistent Fenton reaction:

This reaction is harbored in the cell cytosols due to chronic inflammation (chemically, with elevated levels of H₂O₂ and ⋅O2− ), coupled with local iron accumulation, which continuously produces OH^–.^[61] (2) The continuous alkalization by the persistent Fenton reaction casts a life-threatening stress on the affected cells, thus driving the cells to start multiple acidifying reprogrammed metabolisms to keep the pH stable.^[61,62] (3) The two predominant reprogrammed metabolisms are the de novo biosynthesis of nucleotides (NTs) and biosynthesis and deployment of SAs by all cancers in TCGA, and together they have the following:

Where R (OH^–, FR) represents the rate of OH^– production by Fenton reactions in cancer tissue, R (H⁺, NT) and R (H⁺, SA) are the rates of H⁺ production by nucleotide de novo biosynthesis and SA biosynthesis, respectively, and ε represents the total s produced by all other acidifying reprogrammed metabolisms, which is a relatively small quantity compared to R (H⁺, NT) + R (H⁺, SA).

A direct implication of the above equation is that if R (H⁺, NT) in cancer cannot be as high as that in other cancers for any reason, R (H⁺, SA) must be sufficiently high so their sum matches a major fraction of R (OH^–, FR) to keep the pH stable. Our previous study has provided strong evidence that the increased SA biosynthesis and deployment drives cancer metastasis.^[41,63] The above analysis serves as the basis for our study in the results section.

Principal component analysis (PCA) for geneexpression analyses

PCA is used to capture the main direction of change in a multi-marker gene set in terms of the changes in their collective expressions.^[64] Here, we use PCA analysis on the marker genes involved in the dedifferentiation process to quantify the level of dedifferentiation.

Bayesian information criterion (BIC) for contribution estimation

The BIC is used to select the best model among the given options,^[65] it is defined in the following manner:

Where L is the maximum likelihood function of the model, n is the sample size, and k is the number of free parameters in the model.

The level of contribution to a model (model_K) by the i^th free parameter is estimated using |BIC_K|-|BIC_i|, where BIC_K is the quality of the optimal model_K that is measured using the R² score and BIC_i is the best model based on all parameters except for the i^th parameter. Hence, the percentage of contribution by the i^th parameter can be calculated in the following manner:

Multiple linear regression (MLR)

MLR is used to model the linear relationship between a single dependent variable and multiple free variables. We used BIC-based MLR to perform a linear regression on the five-year survival rate against a list of possible contributors in each sample across the seven cancer types, where BIC provides information regarding the level of contribution to the regression result by each free variable.

Cell culture

RBE and HuCC-T1 cells were obtained from Procell and cultured in RPMI-1640 medium supplemented with 10% fetal bovine serum (FBS, Moregate) and 1% pen-strep antibiotic (Procell). All cells were cultured in a humidified incubator at 37°C with 5% . First, we performed STR testing on the two cell lines to prevent misidentification or cross-contamination of cell lines. STR testing uses primers to extract repeated DNA fragments, usually two to six base pairs of tandem repeats. A DNA kit is used to extract high-purity DNA; thereafter, multiple rounds of PCR amplification are performed—that is, multiple gene loci are amplified. Finally, sequencing is performed to determine the experimental conditions.

Generating persistent Fenton reactions in cytosol in cultured cells

FeSO₄ (Sigma, F8633), H₂O₂ (Aladdin, 7722-84-1), and L-ascorbic acid (Sigma, A4544) were added to the culture to induce sustained Fenton reactions. Here we used L-ascorbic acid instead of ⋅O2− as the reducing molecule for the persistent Fenton reaction.

Cell proliferation and colony formation assays

Cell Counting Kit 8 (CCK-8, Apexbio) assays were performed to assess cell proliferation in accordance with the manufacturer’s instructions; 5000 cells were seeded in 96-well plates in triplicate containing regular fullserum media and allowed to adhere overnight. Different concentrations of lithocholic acid (LCD, MCE) were added the next day, and cells were harvested at the indicated times. Further, cells were incubated with 10% CCK-8 solution at 37°C for two hours, and then the absorbance was measured at 450 nm using a microplate reader.

For the colony formation assay, 600-800 cells were seeded in six-well plates in triplicate and adhered overnight. After 15 days, cells treated with LCD were fixed with polyformaldehyde (4%), washed twice with phosphate buffer saline (PBS), and stained with crystal violet solution (1%) for 15 min at room temperature.

SA assay

The culture medium of cancer cells was collected as the acidic medium. SA was quantified using a Sialic Acid Assay kit (Abcam, ab83375) with the colorimetric assay protocol in accordance with the manufacturer’s instructions.

Statistical analysis

Experiments were performed using a minimum of three replicates for each experimental group. For all statistical analyses, GraphPad Prism software (version 8.0) was used. Unpaired or paired Student t tests were used to compare data between the two groups. Statistical significance was defined as P < 0.05 (^*P < 0.05; ^***P < 0.001; ^****P < 0.0001). All observations were expressed as mean SEM.

CHOL is among the most hypoxic cancer types vs. controls

We used the expression level of HIF1A to reflect the level of hypoxia in a tissue sample, with a higher HIF1A expression indicating a higher level of hypoxia.^[45] Figure 1 depicts the expression levels of HIF1A in samples in localized, regionally migrated, distantly metastasized, combined cancer samples, and control samples for each of the seven cancer types, respectively. It is evident from the figure that CHOL has the largest increase in the expression of HIF1A in cancer vs. control samples across all stages, thereby indicating that CHOL is the most hypoxic cancer type among the seven.

Figure 1

Hypoxic levels of cancer vs. control samples in localized (A), regionally migrated (B), distantly metastasized (C), and all cancer samples combined (D), where the expression of HIF1A is used to reflect the hypoxia level for cancer vs. control samples at each cancer stage across seven cancer types. The y-axis in each panel represents the level of log-transformed expression, measured using TPM. CHOL: cholangiocarcinoma; COAD: colon cancer; ESCA: esophagus cancer; LIHC: liver cancer; PAAD: pancreatic cancer; READ: rectum cancer; STAD: stomach cancer.

CHOL has the highest level of B-cell infiltration

We calculated the infiltration levels of B cells, macrophages, and neutrophils in each tissue sample using the xCell program (see materials and methods) in the control, localized, regionally migrated, distantly metastasized, and combined cancer samples for each of the seven cancer types. Figure 2 presents the calculation results, from which it is evident that CHOL has the highest level of B-cell infiltration vs. the controls across all seven cancer types on average.

Figure 2

The estimated infiltration levels of (A) B cells, (B) macrophages, and (C) neutrophils in the control, localized tumors, tumors with regional migration, tumors with distant metastasis, or combined cancer samples for each of the seven cancer types, where the y-axis represents the percentage of the cells under study in a tissue sample. (D) The relationship between the infiltration level of B cells (the y-axis) and the level of SA biosynthesis and deployment (the x-axis), measured using the first principal component of the expressions of the marker genes of SA biosynthesis and deployment across the relevant samples, and the two have a correlation coefficient 0.823, given in Table S1. CHOL: cholangiocarcinoma; COAD: colon cancer; ESCA: esophagus cancer; LIHC: liver cancer; PAAD: pancreatic cancer; READ: rectum cancer; STAD: stomach cancer; SA: sialic acid.

To understand why CHOL has the highest level of B-cell concentration, a pathway enrichment analysis was undertaken to identify genes (Table S1, supplementary materials), whose expressions are highly correlated with genes involved in SA biosynthesis and deployment. Interestingly, B-cell activation is among the most correlated activities (Figure 2D). Therefore, we posit that the higher levels of SAs attract more B cells into the area, which is consistent with the results of published studies (SAs are the capping molecules of cell-surface glycans, which are part of gangliosides).^{[66,67,68,69]}

CHOL is among the most dedifferentiated cancers

Using the marker genes for the level of cell dedifferentiation given in Table S1, we estimated the level of dedifferentiation via a PCA (see materials and methods). Figure 3 depicts the estimated level of cell dedifferentiation for each of the seven cancer types. It is evident from the figure that CHOL is among the cancer types that have the highest level of dedifferentiation, suggesting that the cancer is, in general, more stressed than other cancer types.

Figure 3

The level of dedifferentiation (the y-axis) is compared between cancer and control samples across the seven cancers in (A) totally localized tumor, (B) tumor with regional migration, (C) tumor with distant metastasis, and (D) combined cancer samples, where weight is the first principal component of the expressions of marker genes for dedifferentiation, where the gene list is given in Table S1. CHOL: cholangiocarcinoma; COAD: colon cancer; ESCA: esophagus cancer; LIHC: liver cancer; PAAD: pancreatic cancer; READ: rectum cancer; STAD: stomach cancer.

CHOL has the highest level of SA accumulation on cancer cell surface

As discussed earlier, cancer cells utilize de novo biosynthesis of NTs and biosynthesis of SAs as the main H⁺ producers to neutralize OH^- produced by Fenton reactions.^[61,70]

Figure 4A and 4B depicts the level of contribution of SA towards keeping the Fenton reaction-produced OH^– neutralized across the seven cancer types. We noted that CHOL stands out in terms of the level of contribution by SA biosynthesis. To understand why SA contributes particularly highly to keeping the cytosolic pH stable, we studied the roles that BAs may have played in this picture. First, we noted that (1) the level of BA, measured using the expression data of the relevant marker genes (Table S1), strongly correlates with the level of SA biosynthesis as well as the level of SA accumulation; and (2) the level of BA negatively correlates with nucleotide de novo biosynthesis, as depicted in Figure 4C, suggesting that BA might play a role in suppressing nucleotide synthesis and hence driving up the SA biosynthesis and deployment in accordance with equation (4).

Figure 4

(A) Contributions by de novo nucleotide (NT) biosynthesis and SA biosynthesis and utilization in neutralizing OH^- produced by Fenton reactions, respectively, where the x-axis represents the weight of contribution by nucleotide biosynthesis and the y-axis represents the weight of contribution by SA biosynthesis and utilization, and a dot on the diagonal indicates that the two biosynthesis processes are equally weighted in the neutralization of the OH^- produced by Fenton reactions. (B) The predicted level of SA accumulation for each cancer type. (C) Correlations between BA biosynthesis and nucleotide de novo biosynthesis as well as between BA biosynthesis and SA biosynthesis. (D) Cell-cycle progression-related pathways enriched by genes that positively correlate with GPBAR1 in CHOL. (E) The BA’s receptor GPBAR1 negatively correlates with cell-cycle marker genes, where each dot represents a CHOL sample with its normalized GPBAR1 expression level (x-axis) and the normalized first principal component of the expression levels of cell-cycle genes (y-axis). (F) Correlation between the de novo nucleotide biosynthesis and the biosynthesis and utilization of SA in different stages. The x-axis represents the (normalized) level of de novo nucleotide biosynthesis, and the y-axis represents the (normalized) level of SA biosynthesis, where all the relevant marker genes are given in Table S1. CHOL: cholangiocarcinoma; COAD: colon cancer; ESCA: esophagus cancer; LIHC: liver cancer; PAAD: pancreatic cancer; READ: rectum cancer; STAD: stomach cancer; SA: sialic acid.

To elucidate the possible connection between the two, we examined the expression level of the BA receptor GPBAR1 and found that it is upregulated, which is consistent with the increased BA level in CHOL. Our literature search revealed that GPBAR1 can suppress cell proliferation by activating the cAMP signaling pathway, which is known to be associated with cancer bone metastasis.^[71,72,73] To examine if this is the case in CHOL, we note that the expression of GPBAR1 indeed negatively correlates with cell-cycle related pathways (Figure 4D, Table S1), and with the level of cell-cycle progression (Figure 4E), thereby confirming that GPBAR1 suppresses cell proliferation in CHOL, as reported in the literature, which further limits the rate of nucleotide biosynthesis for DNA synthesis.

As noted earlier, nucleotide de novo biosynthesis and SA biosynthesis play the predominating roles in neutralizing the OH^-s produced by Fenton reactions. The data depicted in Figure 4F explain why the SA biosynthesis level in CHOL stands out among the seven cancer types—that is, the nucleotide biosynthesis is suppressed as a result of suppressed cell proliferation by BAs via their binding with GPBAR1. Hence, it is the BA that ultimately causes the levels of SA biosynthesis and accumulation to go higher, thereby leading to the high migration rate by CHOL and, hence, the lowest survival rate among the seven cancer types. This computational prediction, the core of our model, is validated experimentally in the last section of results.

Key factors contributing to the low survival rates of CHOL

We now aim to establish a quantitative relationship between the five-year survival rate and the possible contributing factors discussed earlier. To do this, we first predicted the possible causal relationships among these factors to identify and remove those factors that are largely determined by the other factors. Specifically, we sorted all the cancer samples of each cancer type in ascending order of the number of upregulated genes in a sample vs. controls, as a means to provide a pseudo-time order of the disease progression. On this ordered list of cancer samples, we used the F-test statistic in the Granger causality test^[74] to assess if some pairs are causally connected with a significant P-value. In other words, if the P-value of factor i to factor j is less than or equal to 0.05, we consider that factor i can predict factor j—that is, there is a causal relationship between the two. Table 3 lists the calculation results.

It is evident from the table that the P-values of dedifferentiation relative to ΔSA, hypoxia, and B cells are all 0.05, indicating that there is a causal relationship between dedifferentiation and each of the three factors; hence, it was excluded from our regression analyses.

We then conducted a linear regression analysis of the five-year survival rates using the five remaining factors: ΔSA; the level of hypoxia; the infiltration levels by B cell, neutrophil, and macrophage, respectively. Figure 5 presents the four regression models for localized (model₁), regionally migrated (model₂), distantly metastasized (model₃) and combined cancer (model₄) samples with the regression results, respectively.

Figure 5

The five-year survival probabilities (y-axis) of seven types of cancers against the predicted survival probabilities (x-axis) based on SA accumulation, level of hypoxia, and the infiltration levels of B cells, macrophages, and neutrophils through regression analyses, with (A), (B), (C), and (D) representing localized tumors, tumors with regional migration, distant metastasis, and combined cancer sample set, respectively. CHOL: cholangiocarcinoma; COAD: colon cancer; ESCA: esophagus cancer; LIHC: liver cancer; PAAD: pancreatic cancer; READ: rectum cancer; STAD: stomach cancer.

In the above equations, SR is the five-year survival rate; ΔSA is the levels of SA accumulation; H is the hypoxia; and B, M, and N are the infiltration level by B cells, macrophages, and neutrophils, respectively.

Our BIC-based analyses (see materials and methods) have provided the level of contribution by each of these factors, as given in Table 4. Hence, we conclude that the SA accumulation, dictated by the BA level and the level of Fenton reaction, plays the most significant role in making CHOL the deadliest cancer among the seven that were studied. Moreover, the infiltration level of B cells is the second most significant contributor to the lowest survival rate of CHOL, which is partially because of the high SA level, as discussed earlier.

Experimental validation of Fenton reactions and cell-surface SA levels

The core logic of our predictive model is that it is the Fenton reaction that drives the high levels of nucleotide de novo biosynthesis and SA biosynthesis for survival and that it is the BAs in CHOL that suppress the level of cell proliferation, hence the level of nucleotide biosynthesis, which in turn drives the level of SA biosynthesis up to keep the intracellular pH stable for cells affected by Fenton reactions.

To verify this prediction, we selected CHOL cells (RBE and HuCC-T1 cells) for study. We first induced sustained cytosolic Fenton reactions in these cells by FeSO₄, H₂O₂, and L-ascorbic acid (VC) to HuCC-T1, where L-ascorbic acid serves the same purpose of superoxide to reduce Fe³⁺ back to Fe²⁺, to induce persistent Fenton reactions. After preliminary experiments, we found that 100 μmol/L FeSO₄, 100 μmol/L H₂O₂, and 1000 μmol/L L-ascorbic acid were the most suitable concentrations for inducing persistent Fenton reactions in cytosol. When FeSO₄ was added in a concentration of 100 μmol/L, H₂O₂ in a concentration of 100 μmol/L, and L-ascorbic acid in a concentration of 1000 μmol/L, significantly increased cell proliferation was observed at 24 h and 48 h compared with the control group as result of increased nucleotide biosynthesis induced by Fenton reactions, and the difference was statistically significant (Figure 6 and Table S3).

Figure 6

The sustained Fenton reaction promotes cell proliferation (Table S3). (A and B) The cell proliferation ability of the RBE cells was detected after adding FeSO_4′, H₂ O_2′, and L-ascorbic acid. (C and D) The cell proliferation ability of the HuCC-T1 cells was detected after adding FeSO_4′, H₂ O_2′, and L-ascorbic acid. (E) The colony formation assay was performed after adding, and L-ascorbic acid. (F) Statistical diagram of the colony formation assay (^****P < 0.0001; ^***P = 0.0002).

We have then examined how the proliferating cells will be affected when treating them with the proliferation inhibitor LCD (Lithocholic acid) as in the case of BAs in CHOL. First, we note that the cell viability was the lowest when treated with the LCD with concentration at 100 μmol/L among other levels of concentrations (Figure 7A). Then, we tested the impact of different LCD concentrations on cell proliferation and found that LCD had the strongest inhibition of cell proliferation at the same concentration of 100 μmol/L for 48 h and 72 h (Figure 7B and 7C). Similarly, the colony formation assay found the same results (Figure 7D and 7E). Therefore, we chose the LCD concentration at 100 μmol/L for further experiments.

Figure 7

LCD inhibited proliferation of CHOL cells and increased cell-surface SA level (Table S4). (A) the cell viability (%) of CHOL cells was detected after adding different concentrations of LCD. (B and C) the cell proliferation ability of RBE cells and HuCC-T1 cells was detected after adding different concentrations of LCD. (D) the colony formation assay was performed after adding LCD. (E) statistical diagram of colony formation assay. (F) SA content (nmol/μL) of different types of cancer cells. (G) the SA content (nmol/μL) of CHOL cells after adding LCD was detected. (H) SA content (nmol/μL) of CHOL cells at different time points after adding LCD. (^*P < 0.05; ^***P < 0.001; ^****P < 0.0001.) SA: sialic acid.

We then measured the SA level on the cell-surface in cells of stomach, liver, colon, and pancreatic cancer, together with CHOL cells, separately, and found that consistent with our above results, the highest SA content was found in CHOL cells (Figure 7F). Subsequently, we treated the CHOL cells with LCD at 100 μmol/L and observed that the cell-surface SA content continued to increase (Figure 7G) with time (P-value = 0.0174) (Figure 7H).

The conclusions are that (1) by inducing Fenton reactions in the cytosol of cancer cells, the cells proliferate at higher rates as predicted by our computational study and validated by our cell-based experimental study, and (2) by inhibiting experimentally the proliferation rate of such cancer cells, the cell-surface SA content increases as we have predicted, which drives cells to metastasize.^[41] This explains why CHOL is one of the most malignant cancers.

Supplementary Information

Supplementary materials are only available at the official site of the journal (www.intern-med.com).