Identifying Telomere-Related Genes as Prognostic Biomarkers of Cholangiocarcinoma and Exploring Their Biological Functions Based on Integrated Analysis of Single-Cell and Bulk RNA Sequencing Data

1. Introduction

CCA is a very deadly form of adenocarcinoma that develops in the epithelial cells of the bile ducts within the hepatobiliary system [1,2,3]. Worldwide, approximately 106,000 people are diagnosed with CCA each year, accounting for 0.3% of all cancers, and the incidence is higher in Asian populations compared to European and North American populations [4]. As early symptoms of bile duct cancer are not obvious, most of the patients have already progressed to advanced stages by the time of diagnosis, losing the opportunity for surgical treatment and ultimately affecting the prognosis [5]. Globally, about 99,000 people die from bile duct cancer every year, accounting for about 0.4% of all cancer mortality. Hence, there is an immediate requirement to develop prognostic models using novel biomarkers in order to enhance the unfavorable prognosis of patients with CCA.

Telomeres, consisting of repetitive DNA sequences found at the tips of chromosomes, are essential for preserving chromosome stability and integrity [6,7]. Telomere shortening and damage are closely related to the development of many diseases such as aging and cancer. Studies have shown that telomere length in malignant cells is usually longer than in normal cells, which is associated with the unlimited proliferative capacity and resistance to apoptosis of tumor cells. Regulating this process is a group of highly conserved regulatory factors, telomere-related genes. telomere-related genes are genes related to telomere structure, function, and regulation of telomerase [6]. Abnormal expression or mutation of these genes will lead to abnormal shortening or damage of telomeres, which will trigger chromosomal instability and abnormal cellular function. In tumor development, the abnormal expression of telomere-related genes is closely related to the shortening of telomeres and chromosomal rearrangements. These abnormalities lead to genomic instability, which promotes the proliferation, metastasis, and drug resistance of tumor cells. Therefore, telomere-associated genes play an important role in regulating tumor development by modulating the stability and function of telomeres [8]. In CCA, several telomere-associated genes are differentially expressed and play a wide range of biological roles, and these genes may be closely related to the treatment and prognosis of CCA. However, the complex role and prognostic significance of telomere-associated genes in CCA are still largely unknown [9,10,11].

The contribution of single-cell sequencing technology to tumor research has been immense. It has enabled researchers to gain insight into the cellular heterogeneity and evolutionary processes of tumors [12]. Single-cell sequencing data are important in identifying cellular heterogeneity, constructing cellular trajectories, and analyzing intracell communication and ligand receptor expression. By analyzing the genomes, transcriptomes, and epigenomes of individual cells, we can reveal the differences between different tumor cells, identify potential therapeutic targets, and provide a basis for personalized therapy. The rapid development of single-cell genomics in recent years has revealed the unknown biological processes and disease causes of CCA [13].

In this study, transcriptomic data, clinical data, and mutation data of CCA were downloaded through TCGA and GEO databases. Single-cell analysis was performed based on single-cell sequencing datasets to identify different cell types and telomere-associated active cells in CCA. The pseudotime analysis to construct cell trajectories and identify branching events, and cell communication analysis to reveal the interactions and signaling pathways between CCA cells. Next, we performed an enrichment analysis of telomere activity-related differential genes and explored the functions of these genes in biological processes and pathways. We constructed prognostic models and column line plots based on the acquired prognostic-related genes, predicted the overall survival of patients, and evaluated the efficacy of the models. Furthermore, we conducted gene set enrichment analysis, gene set variant analysis, immune infiltration analysis, and drug sensitivity analysis to uncover the molecular characteristics and biological pathways underlying CCA. This comprehensive approach offers valuable insights for diagnosing and treating CCA, with the ultimate goal of discovering novel biomarkers and enhancing prognostic prediction for patients with this disease.

2.1. Transcriptome Data Download and Preparation

All data used in this study are freely available to the public, mainly from the TCGA (The Cancer Genome Atlas, https://portal.gdc.cancer.gov/) database and GEO (Gene Expression Omnibus, https://www.ncbi. nlm.nih.gov/geo/) databases. CCA expression profiling data, and clinical data were downloaded from the TCGA database via the R package “TCGAbiolinks (version 2.25.0)” [14] as well as by the Single nucleotide mutation (SNV) data predicted by the “VarScan2 Variant Aggregation and Masking” tool were downloaded from the TCGA database. 36 tumor samples and 9 paraneoplastic control samples from TCGA-CHOL and 30 tumor samples from GSE107943 were included in this study. The ComBat method using the R package “sva” [15] corrected for the batch effect of non-biotechnology bias. Principal component analysis (PCA) was used to check the extent of correction.

2.2. Downloading and Processing of Single-Cell Sequencing Data

The GEO database contains a large amount of single-cell sequencing data. In this study, the single-cell sequencing dataset GSE138709 of CCA, which contains five CCA disease samples and three paracancerous controls, was obtained from the GEO database. The single-cell raw data in GSE138709 were imported using the Seurat package (version 4.2.0) [16] for the R programming language. First, low-quality cells and genes were filtered by the following criteria: 1) Cells expressing fewer than 300 genes were removed; 2) Genes expressed in fewer than 1 cell were removed. 3) Cells with fluctuating numbers of expressed genes between 200 and 6,000 were retained. 4) Cells with a percentage of mitochondrial genes of less than 10% were retained. 5) Cells with UMI readings of less than 40,000 were retained. The data were normalized using the “normalizedata” function in the Seurat R package. After data normalization, highly variable genes in single cells were identified by balancing the relationship between average expression and dispersion. Principal component analysis (PCA) was then performed and significant principal components (PCs) were used as input for graph-based clustering. A harmony approach was used to eliminate batch effects across samples. For clustering, we used the function FindClusters, which is based on the Shared Nearest Neighbor (SNN) modular optimized clustering algorithm, to generate 18 clusters on 20 PC components at a resolution of 0.3. The “Runtsne” function is then used for t-distributed random neighbor embedding (t-SNE). Cell aggregation was demonstrated using t-SNE-1 and t-SNE-2. Subsequently, cell clusters were identified by cell type-specific biomarkers, and the proportion of cell types was calculated and evaluated.

2.3. Telomere-Related Gene Scoring

The R package “AUCell” [17] scores each cell based on gene set enrichment analysis. A gene expression ranking for each cell was generated based on the area under the curve (AUC) values of 2089 human telomere-associated genes (Table S1) obtained from the TelNet (http://www.cancertels-ys.org/telnet/) database to estimate the proportion of the gene set with high expression in each cell. Cells that express more genes in the gene set have higher AUC values. The “AUCell_exploreThresholds” function was used to determine a threshold for identifying active cells in the gene set. The AUC score of each cell was then mapped to the TSNE embedding using the “ggplot2” R package (version 3.3.5) to visualize the activated cell clusters.

2.4. Constructing Cellular Trajectories by Pseudotime Analysis

The pseudotime analysis was performed by Monocle 2 [18], which performs reverse graph embedding based on a user-defined list of genes to generate a pseudo-temporal graph that can account for branching and linear differentiation processes. For the proposed temporal analysis of cells with high telomere activity, the raw count data were normalized by estimating a size factor for trajectory inference. Pseudotemporal trajectories were constructed using genes with a high degree of discretization and high expression (discretization estimate ≥ 1 and mean expression ≥ 0.1) [19]. default values were selected for the parameters of the DDRTree algorithm. To further analyze these branching events, we used branching expression analysis modeling (BEAM) implemented in Monocle 2. It helps to identify all genes with significant branching-dependent expression. Branch-dependent expression patterns were visualized as heatmaps using Monocle 2.

2.5. Cell Communication Analysis and Ligand-Receptor Expression

CellChat objects were created by the “CellChat” R package (https://www.github.com/sqjin/CellChat) [20] based on the UMI count matrix of each group (CCA and control). Cell-to-cell communication was analyzed using the “CellChatDB.human” ligand-receptor interaction database as the reference database with default parameters. The “mergeCellChat” function was used to combine the CellChat objects of each group to obtain a comparison of the total number of interactions and the intensity of interactions. The “netVisual_diffInteraction” function was used to visualize the differences in the number or strength of interactions between different cell types between groups. Finally, we used the “netVisual_bubble” and “netVisual_aggregate” functions to visualize the distribution of signaling gene expression between groups.

2.6. GO/KEGG Pathway Enrichment Analysis

Gene Ontology (GO) [21]e nrichment analysis includes biological process (BP), molecular function (MF), and cellular component (CC) analysis. GO enrichment analysis of genes was performed using the R package “clusterProfiler (version 4.2.2)” [22]. Kyoto Encyclopedia of Genes and Genomes (KEGG) [23] is a bioinformatics resource for mining gene lists enriched with significantly altered metabolic pathways. The KEGG enrichment analysis of telomere-related differentially expressed genes (DEGs) in CCA was performed using the R package “clusterProfiler” (p-value <0.05).

2.7. Construction and Validation of Prognostic Models

The genes obtained by taking the intersection with telomere-related genes as key genes were analyzed by the results of the difference analysis of telomere high-activity and low-activity cells in the single-cell dataset, the results of the difference analysis of disease and control, and the results of the difference analysis of the CCA and normal groups in the CCA dataset of TCGA. To assess the prognostic value of the key genes, we used one-way Cox regression analysis to assess the correlation between each gene and survival in the tumor cohort, and genes with P ≤ 0.05 were considered to be significantly associated with survival for further analysis. All CCA tumor samples (n=66) with clinical information were randomly divided into training and validation sets in a 7:3 ratio. The LASSO Cox regression model (R package “glmnet” [24]) was then utilized to narrow down the candidate genes and construct prognostic models. The penalty parameter (λ) was determined by the minimum criterion. Risk scores were calculated using the following formula:

$$\mathrm{riskScore}={\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\mathrm{Coef}\left({\mathrm{gene}}_{\mathrm{i}}\right)}\ast \mathrm{Expression}\left({\mathrm{gene}}_{\mathrm{i}}\right)$$

(Coef (genei): risk factor, Expression (genei): gene expression level)

The training group samples were divided into low-risk and high-risk groups based on the optimal cut-off value determined by the surv_cutpoint function. Survival curves were then generated to assess prognostic significance utilizing the Kaplan-Meier method, and statistical significance was evaluated through the log-rank test. A subject operating characteristic curve (ROC) was used to validate the efficacy of the prognostic model. AUC values typically ranged from 0.5 to 1, with closer to 1 indicating better model efficacy. The validation set was also divided into low-risk and high-risk groups for validation of the prognostic model.

2.8. Construction and Validation of Column-Line Diagrams

Clinical information of patients (age, gender and disease stage, etc.) was extracted from the TCGA cohort. A one-way Cox regression analysis was performed combining the clinical information with the risk scores obtained from the prognostic model. In addition, we combined the prognostic and clinical characteristics and constructed a column chart using the R package “RMS” to predict overall survival at 1, 3, and 5 years. The effectiveness of the column plot model was evaluated by ROC curves.

2.9. Gene Set Enrichment Analysis (GSEA)

The Gene Set Enrichment Analysis (GSEA) algorithm assesses whether a predefined set of genes exhibits statistically significant differences between two biological conditions.GSEA is performed with the R package “clusterProfiler (version 4.2.2)” on an ordered list of all genes aligned according to log2FC, with 1000 gene set alignments performed for each analysis. In this study, c2.cp.kegg.v7.5.1.symbols was used as the reference gene set, which is stored in the Molecular Signatures Database (MSigDB) [25,26,27] database. p < 0.05 was considered to be a gene set with a P < 0.05. 0.05 of the gene set was considered significantly enriched.

2.10. Gene Set Variation Analysis (GSVA)

To investigate the differences in biological function between high and low-risk groups, Gene Set Variation Analysis (GSVA) was performed using the R package “GSVA (version 1.42.0)” based on “c2.cp.kegg.v7.5.1.symbols”. GSVA). The R package “pheatmap (version 1.0.12)” was used to visualize the results.

2.11. Immune Infiltration Analysis

Single Sample Gene Set Enrichment Analysis (ssGSEA) [PMID: 36090900], an extension of Gene Set Enrichment Analysis (GSEA), calculates separate enrichment scores for each pair of samples and gene sets. Each ssGSEA enrichment score represents the extent to which genes in a particular gene set are coordinately up-or down-regulated in the sample. Based on 28 immune cells downloaded from the TISIDB (Tumor and Immune System Interactions Database) (http://cis.hku.hk/TISIDB/index.php) [28] database, including Activated CD8 T cell, Central memory CD8 T cell, Effector memory CD8 T cell, Activated CD4 T cell, Central memory CD4 T cell, Effector memory CD4 T cell, T follicular helper cell, Gamma delta cell, T cell, and T cell. helper cell, Gamma delta T cell, Type 1 T helper cell, Type 17 T helper cell, Type 2 T helper cell, Regulatory T cell, Activated B cell, Immature B cell, Memory B cell, Natural killer cell, CD56bright natural killer cell, CD56dim natural killer cell, Myeloid derived suppressor cell, Natural killer T cell, Activated dendritic cell, Type 1 T helper cell, Type 17 T helper cell, Type 2 T helper cell, Regulatory T cell Activated dendritic cell, Plasmacytoid dendritic cell, Immature dendritic cell, Macrophage, Eosinophil, Mast cell, Monocyte and Neutrophil, Quantification of gene expression profiles for each CCA sample from each CCA disease samples to quantify the relative enrichment score of each immune cell from the gene expression profiles. Differences in immune cell infiltration levels between high and low-risk groups were visualized by the R package “ggplot2 (version 3.3.6)” [29].

2.12. Somatic Mutation Analysis

A landscape of genomic mutational variation was depicted based on mutation data. the R package “maftools” was used to display somatic variants between different subgroups, including single nucleotide polymorphisms (SNPs), insertions and deletions (INDELs), tumor mutational load (TMB), and mutation frequency [30]. Generally, frequently mutated genes (FMGs) with the top 20 mutation frequencies are considered to be the main driver genes of malignant tumors [31].

2.13. Drug Sensitivity Analysis

Based on the half-maximal inhibitory concentration (IC50) data and corresponding gene expression data downloaded from the GDSC database (Genomics of Drug Sensitivity in Cancer, https://www.cancerrxgene.org/) [32], the R package was used “oncoPredict (version 0.2)” [33] to predict potential therapeutic drug sensitivities in patients in the high and low-risk groups of CCA.

2.14. Statistical Analysis

This study was statistically analyzed using R software v4.1.2. Spearman correlation test was used to infer the correlation between two parameters. The Wilcoxon test was employed to assess disparities between two groups, while the Kruskal-Wallis test was applied to evaluate differences among three or more groups. A two-tailed p-value below 0.05 was deemed to indicate statistical significance. The detailed analysis process is shown in Figure 1.

3. Results

3.1. Single-Cell Dimensionality Reduction Clustering and Annotation

To explore the origin of highly expressed genes, we analyzed cell populations from CCA samples using the single-cell sequencing dataset GSE138709. After the initial quality control assessment, we obtained a total of 29,983 cells from the single-cell transcriptome. 8 samples were included in the study, and all cells were aggregated into 18 clusters (Figure 2A), and different cell types were annotated by cell-specific biomarkers based on the gene expression profile of each cluster (Table S2). As shown in Figure 2B, 10 cell types can be found, such as NK cells, malignant cells, and B cells. specific genes for each cell type were visualized by dot plots (Figure 2C). The proportions of different cell types in the control samples (GSM4116579, GSM4116582, GSM4116586) are shown in Figure 2D; the proportions of different cell types in the disease samples (GSM4116580, GSM4116581, GSM4116583, GSM4116584, GSM4116585) are Figure 2E shown.

3.2. Characterization of Telomere-Related Activity

A subpopulation of active cells was utilized to study the expression pattern of telomere-associated genes at the single-cell level. Active cells were identified by the optimal threshold, which showed the presence of 452 telomere-associated active cells, the population of cells with AUC values greater than 0.25 (Figure 3A). Figure 3B highlights the TSNE plot of the active cells and Figure 3C shows the telomere activity of all cells, the darker the color the higher the telomere activity, from Figure 3B and 3C it can be seen that the telomere active cells are mainly malignant cells.

3.3. Pseudotime Analysis

Using telomere-associated active cells, we built a pseudo-temporal cell trajectory to find the important gene expression programs that determine CCA progression. Indeed, the transcriptional states in the trajectory revealed different processes. malignant cells were distributed in different branches, suggesting that malignant cells have different transcriptional states (Figure 4A–D). To elucidate the molecular basis of malignant cells transformation, we explored the genes that determine the branching of malignant cells in CCA. Genes highly expressed in pre-branch were mainly enriched in GO BP pathways such as wound healing and skin development. Genes enriched in response to toxic substance-related pathways were highly expressed in cell branch 2, and genes enriched in T cell activation, antigen processing, and presentation of exogenous peptide antigen via MHC genes in class II were highly expressed in cell branch 1 (Figure 4E).

3.4. Cell Communication Analysis

To further investigate the cellular network of interactions in CCA, we used the R package “Cellchat” to reveal changes in crosstalk between CCA cells. The network plot showed an overall increase in the number of interactions between CCA cell types, an increase in the strength of interactions associated with B cells, fibroblasts, and hepatocytes, and a general decrease in the strength of interactions associated with NK cells (Figure 5A,B). We then further compared the signaling patterns between cells from the CCA group (Figure 5C,D). NK cells were the major signal receivers. Subsequently, from the probability point plots of communication signaling between malignant cells and the ligand-receptor pairs between NK cells and B cells, we could see that the stronger changing communication signaling pairs were MIF - (CD74+CXCR4), HLA-A - CD8A, etc. (Figure 6A). We also explored the interaction of malignant cells with other cell types through the MIF, MHC-I signaling pathway where the above signaling pairs are located (Figure 6B,C). These results preliminarily elucidated the potential interactions between these cell types and helped us to further investigate the integrative role of malignant cells in CCA. As can be seen from Figure 7A–D expression of specific pathway ligand receptors, ligand HLA-A was highly expressed in NK cells.

3.5. Enrichment Analysis of Telomere Activity-Related Differentially Expressed Genes

A total of 805 differentially expressed genes (DEGs 1) were identified by taking the intersection of the results of the differential analysis of telomere high and low activity cells in the single-cell dataset (n=1674, Table S3), and the results of the differential analysis of disease vs. control (n=2225, Table S4), and the differences of these genes between the two groups were statistically significant (p-value < 0.05, | Log2 fold change| > 0.25). We then compared the CCA and normal groups in the TCGA CCA dataset and identified a total of 3658 differentially expressed genes (Table S5), and intersected with telomere-associated genes (n=2089) to obtain 472 telomere-associated differentially expressed genes (DEGs 2), which were statistically different between the two groups (p-value < 0.05, | Log2 fold change| > 2). The distribution of telomere-associated differential genes was demonstrated by volcano plots, and the top 5 up-regulated genes (RALY, LCMT1, CD2BP2, FUS, CCDC115) and top 5 down-regulated genes (LINC01831, C3orf85, LINC02754, AC016682.1, SLC6A13) were illustrated by heat maps (Figure 8A,B). The two sets of differential genes were taken to intersect to obtain 47 key genes (hub genes, Table S6) as shown by the Wayne diagram (Figure 8C).

In order to investigate the biological functions associated with the differential genes, we analyzed the enrichment of telomere-associated differential genes for GO entries (Table S7) and KEGG pathway (Table S8), which are shown by lollipop diagrams. The GO results showed that these genes were mainly involved in “DNA replication”, “positive regulation of DNA metabolic process”, “regulation of DNA metabolic process”, “regulation of DNA metabolic process”, and “DNA metabolism”. “regulation of DNA metabolic process” and “chromosomal region” (BP). Enriched in biological processes such as “chromosomal region”, “chromosome, telomeric region”, “nuclear chromosome”, etc. (BP). Enriched in “chromosomal region”, “chromosome, telomeric region”, “nuclear chromosome” and other cellular components (CC). Enriched in “catalytic activity, acting on DNA”, “single-stranded DNA binding”, “single-stranded DNA helicase”. “single-stranded DNA helicase activity” were enriched in molecular functions (MF) (Figure 8D). KEGG pathway includes “Cell cycle”, “DNA replication”, “Cellular senescence”, etc. (Figure 8E).

3.6. Construction and Validation of Prognostic Models

A univariate Cox analysis was used to identify key genes associated with prognosis (p < 0.05), and disease samples in TCGA and GEO were integrated to finally identify five genes associated with CCA prognosis (DDX39A, SNRPE, TMEM109, ANXA4, ATP1A1). 7/10 CCA samples (66 in total) were selected by random sampling as the training set (n=44) and 3/10 as the validation set (n=22). LASSO regression analysis was performed on the training set to remove redundant genes, and the random number seed was set to 17. A total of 3 genes associated with the prognosis of CCA patients were identified (Table S9), and the results are shown in Figure 9A,B. To determine the robustness of the model constructed using the 3 gene features, the samples were categorized into high-risk and low-risk groups based on the median risk value as a threshold. KM survival curves were constructed for different groups of patients in the training cohort (Figure 9C) and validation cohort (Figure 9D). The results showed that patients in the high-risk group had a significantly worse prognosis compared to patients in the low-risk group in all cohorts. ROC curves were used to determine the efficacy of the model in predicting patient prognosis. In the training cohort, the AUC values for 1, 3, and 5year survival were 0.701, 0.717, and 0.769, respectively (Figure 9E). In the validation cohort, the AUC values for 1, 3, and 5 year survival were 0.779, 0.720, and 0.759, respectively (Figure 9F).

3.7. GSEA and GSVA

To further explore the potential mechanisms of differentially expressed genes, we performed a GSEA enrichment analysis. Using the pathway information in the MsigDB database as a reference, we selected the most significant pathways based on the normalized enrichment score (NES) (Table S10). The results of the GSEA analysis indicated that DRUG METABOLISM CYTOCHROME P450, RETINOL METABOLISM, COMPLEMENT AND COAGULATION CASCADES, TIGHT JUNCTION, ECM RECEPTOR INTERACTION, and FOCAL ADHESION were the most significant pathways. FOCAL ADHESION and other pathways and biological processes were more active in CCA (Figure 10A–F, P< 0.05). In addition, we performed GSVA enrichment analysis using the pathway information in the MsigDB database as a reference and selected the pathways with the most significant differences between the high- and low-risk groups to draw pathway activity maps (Figure 10G, Table S11).

3.8. Immune Infiltration Analysis

We further investigated the level of infiltration of 28 immune cell types in the high- and low-risk groups using the ssGSEA method. The heat map of correlation between immune cells showed that most immune cells were positively correlated with each other (Figure 11A). Immune cells central memory CD4 T cell, Immature B cell, etc. showed significant differences between high and low-risk groups (P value <0.05, Figure 11B). In addition, we examined the significant correlation between prognostic genes and corresponding immune cells and found that ATP1A1 was positively correlated with Central memory CD4 T cell (R=0.327, p < 0.05) and negatively correlated with Immature B cell (R=-0.25, p < 0.05) (Figure 11C,D).

3.9. Construction of Nomogram Model

In order to verify whether risk score could be used as a prognostic factor, a one-way Cox regression analysis was performed on the clinical characteristics of patients (e.g., age, disease stage, and gender). The results showed that risk score and disease stage were independent prognostic risk factors for patients (Figure 12A). The results of constructing a nomogram showed that risk score significantly predicted clinical outcomes (Figure 12B). ROC curves were used to determine the efficacy of the nomogram in predicting patient prognosis. The results showed that the nomogram model had good stability and accuracy at 1, 3, and 5 years (Figure 12C).

3.10. TMB and Drug Sensitivity Analysis

We also assessed specific gene mutations in CCA and visualized the top 20 genes with the highest mutation frequency. PBRM1 had the highest mutation frequency in both high and low-risk groups and was much more frequent in the high-risk group than in the low-risk group (Figure 13A,B). Tumor mutational load (TMB) is a key criterion for the success of immunotherapy. Therefore, somatic mutations associated with CCA were analyzed, and TMB was significantly higher in the high-risk group than in the low-risk group (Figure 13C, P<0.05). We analyzed whether the risk score accurately predicted chemotherapy sensitivity in patients with CCA. The results showed that high-scoring patients were more sensitive to Vinblastine_1004 (Figure 13D), Vinorelbine_2048 (Figure 13E), Docetaxel_1007 (Figure 13F) and other drugs.