Machine learning-based identification of telomere-related gene signatures for prognosis and immunotherapy response in hepatocellular carcinoma

Abstract

Telomere in cancers shows a main impact on maintaining chromosomal stability and unlimited proliferative capacity of tumor cells to promote cancer development and progression. So, we targeted to detect telomere-related genes(TRGs) in hepatocellular carcinoma (HCC) to develop a novel predictive maker and response to immunotherapy. We sourced clinical data and gene expression datasets of HCC patients from databases including TCGA and GEO database. The TelNet database was utilized to identify genes associated with telomeres. Genes with altered expression from TCGA and GSE14520 were intersected with TRGs, and Cox regression analysis was conducted to pinpoint genes strongly linked to survival prognosis. The risk model was developed using the Least Absolute Shrinkage and Selection Operator (LASSO) regression technique. Subsequently, evaluation of the risk model focused on immune cell infiltration, checkpoint genes, drug responsiveness, and immunotherapy outcomes across both high- and low-risk patient groups. We obtained 25 TRGs from the overlapping set of 34 genes using Cox regression analysis. Finally, six TRGs (CDC20, TRIP13, EZH2, AKR1B10, ESR1, and DNAJC6) were identified to formulate the risk score (RS) model, which independently predicted prognosis for HCC. The high-risk group demonstrated worse survival outcomes and showed elevated levels of infiltration by Macrophages M0 and Tregs. Furthermore, a notable correlation was observed between the genes in the risk model and immune checkpoint genes. The RS model, derived from TRGs, has been validated for its predictive value in immunotherapy outcomes. In conclusion, this model not only predicted the prognosis of HCC patients but also their immune responses, providing innovative strategies for cancer therapy.

Introduction

Liver cancer, which is globally the third leading cause of mortality from oncology poses a formidable challenge to efforts to extend life expectancies in all countries [1]. The year 2022 saw nearly 865,269 new liver cancer diagnoses, which is 4.3% of all novel cancer cases globally [2]. Additionally, liver cancer led to roughly 757,948 deaths, making up 7.8% of the total cancer mortality rate [2]. Within this context, hepatocellular carcinoma (HCC) stands out as the primary form, constituting 75–85% of all primary liver cancer cases according to GLOBOCAN 2020 data [3]. Currently, the main strategies for treating cancer encompass surgical resection, ablation therapy, interventional therapy, systemic anti-tumor immunotherapy, and other comprehensive treatments that are focused on optimizing therapeutic benefits [4]. Despite the considerable advancements in clinical management, HCC remains highly heterogeneous, insensitive to chemotherapy, and slow to respond to targeted therapy [5]. Consequently, HCC patients continue to face high mortality and recurrence rates, resulting in suboptimal five-year survival rates and median survival times [6]. Hence, the importance of in-depth investigation to uncover fresh predictive “signatures” for the prognosis and the efficacy of immunotherapy in HCC sufferers cannot be overstated.

Telomere shortening in somatic cells normally leads to senescence or apoptosis, acting as a barrier to unlimited proliferation and tumor formation [7]. Nevertheless, malignant cells circumvent this restriction by securing a mechanism for telomere preservation, known as TMM [8]. Telomere maintenance can be maintained by telomerase, alternative length replication, telomere-binding proteins, and DNA damage repair pathways, resulting in unlimited proliferation of tumor cells [8,9,10]. Telomerase activity is observed in the majority of HCC cases, nearly 90%, and this enzyme’s resurgence is associated with alterations to the promoter of the TERT gene [11]. Analyzing telomere lengths in hepatic tissues from 978 individuals with HCC, scientists identified that extended telomeres were a hallmark of a highly malignant subtype of the disease. These longer telomeres correlate with the G3 subclass of the transcriptome, TP53 gene mutations, and an unfavorable outcome [12]. The construction of a prognostic profile of telomere-related genes (TRGs) in lung cancer and renal cancer has been recently recorded in the literature, both with favorable predictive results [13, 14]. However, prognostic prediction models for TRGs have not been well discussed in HCC. Previous studies have also utilized similar bioinformatics approaches to develop potential biomarkers [15, 16].

Within this investigation, we harnessed the power of machine learning algorithms to detect TRGs that correlate with the progression of HCC, and from these findings, formulated a predictive genetic signature. Our findings not only provided potential targeted therapeutic strategies but and also provided predictive indicators for the evaluation of HCC survival outcome.

Materials and methods

Data and resources

The overall workflow is summarized in Fig. 1. Utilizing TCGA, we accessed and downloaded RNA-sequencing results and clinical data pertaining to HCC patients. From this database, we gathered RNA-sequencing data along with clinical profiles from 374 patients with HCC and also collected samples from 50 adjacent non-tumorous tissues. GSE14520 was secured from the GEO database through a web-based retrieval process. These included 225 cases of HCC and 220 cases of paraneoplastic liver tissue. We also downloaded ICGC database(https://dcc.icgc.org). Before its June 2024 closure, we accessed its data retrieval page, filtered for hepatocellular carcinoma - related content, and downloaded gene expression, donor, sample, and specimen data. Moreover, 2093 TRGs were initially retrieved from the TelNet database (http://www.cancertelsys.org/telnet), and a final set of 2086 genes was obtained after removal of duplicates.

Fig. 1

Flowchart of the study design illustrating data acquisition, TRGs screening, and model development

Identification of TRGs associated with HCC and functional analyses

Utilizing the “ DESeq2” (v 1.44.0) and “limma”(v 3.60.2) packages, we aimed to discern gene expression disparities within tumor tissues as compared to their surrounding normal counterparts [17, 18]. For the TCGA dataset, a higher threshold was set with|log2FC| ≥ 2 and p < 0.05 to minimize false positives. For the GSE14520 dataset, considering the smaller sample size and to ensure that biologically relevant genes were not overlooked due to a stringent threshold, we applied a lower threshold with|log2FC| ≥ 1 and p < 0.05. Subsequently, we retrieved TRGs from the Telnet database and performed an intersection analysis with genes that showed varied expression levels in both the TCGA-LIHC and GSE14520 datasets, thereby generating a pool of prospective candidate genes. This process resulted in a final list of 34 differentially expressed TRGs. To explore the functional characteristics and involvement in biological processes and signaling pathways of these 34 TRGs, we used the clusterProfiler R package (v 3.14.3) via the SangerBox platform [19], we executed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses for the differentially expressed TRGs.

Filtering the optimal models through machine learning

We crossed variably expressed genes in TCGA and GSE14520 with TRGs and then identified genes significantly associated with survival prognosis (p < 0.01). To split the TCGA dataset into training and testing sets, we used the caret package in R. Specifically, we employed a 50:50 split ratio, which is a common practice in machine learning to balance the trade-off between model training and testing. The follow-up analysis leveraged several machine learning approaches, including the Random Survival Forest using the rfsrc function from the “randomForestSRC” package (v 3.2.3) [20], CoxBoost using the CoxBoost function from the “CoxBoost” package (v 1.5) [21], Elastic Net using the glmnet function from the “glmnet” package(v 4.1.8) [22], and LASSO-Cox using the glmnet function from the “glmnet” package with alpha set to 1 [22]. The LASSO-Cox regression model (implemented via the glmnet package, v4.1.8) was trained with 10-fold cross-validation. The optimal regularization parameter (λ) was selected based on the minimum cross-validation error. To ensure reproducibility, a random seed (set.seed(10210)) was applied during dataset partitioning.

Development and confirmation of a risk assessment model utilizing TRGs

To identify TRGs suitable for constructing a prognostic signature, we conducted a LASSO regression analysis with the aid of the “glmnet” package (v 4.1.8). Afterwards, the risk model was built using the gene expression levels and their associated regression weights. The TRGs risk score (RS) was computed based on the weighted LASSO coefficients derived from individual gene expression levels in the following manner: TRGs risk score= ∑ (expression of gene_i × Coefficient of gene_i). Patients in the TCGA discovery set and validation cohort, including GSE14520 and ICGC, were stratified into high and low-risk groups utilizing the median RS as the cutoff, and survival outcomes was assessed using Kaplan-Meier. Detailed patient information can be found in Supplementary Tables 1–3. Distinctions in survival curves were evaluated with the log-rank test among the various cohorts. According to “survivalROC” R package, specificity and sensitivity of our risk model for prognostic prediction were assessed through Receiver operating characteristic (ROC) analysis.

Exploration of the immune cell infiltration features within the model

CIBERSORT—a computational approach for appraising the relative frequencies of different cell populations—is shorthand for Cell-type Identification by Estimating Relative Subsets of RNA Transcripts [23]. TCGA-LIHC dataset was employed to determine the infiltration degrees of immune cells.“Hmisc” package (v 5.1.3) was implemented for correlation of model genes with immune checkpoint genes [24].

Variations in pathway enrichment among patients belonging to distinct risk categories

Hallmark and KEGG pathway enrichment analyses were performed using the msigdbr package (v 7.5.1) in R [25]. Gene sets for Homo sapiens from MSigDB were organized into lists, and pathway activity was quantified with GSVA (parallel.sz = 10). Enrichment scores were analyzed using the limma package to identify significant pathways between risk groups. Heatmaps visualized the top 20 enriched pathways.

Drug sensitivity analysis and analysis of patients’ response to immunotherapy

“oncoPredict” is an R language package that predicts drug responsiveness based on publicly available drug susceptibility data and patient gene expression data [26]. Creating a heat map visualizes and box plots the relationship between the transcript levels of model genes and drug response as indicated by IC50. Generate histograms to illustrate the distribution of patients’ immune therapy responses across high and low-risk categories.

Statistical analyses

With R software (v 4.2.3) employed for statistical assessments, the Wilcoxon signed-rank test was used for individual sample comparisons, and the Kruskal-Wallis test was applied for comparisons among multiple samples. A p-value under 0.05 indicated a statistically significant result.

Results

Extraction and biological function analysis of differentially expressed TRGs in HCC

Extraction of transcriptomic data and clinical records of patients from the TCGA and GEO repositories. The TCGA-LIHC database comprises 374 tumor samples and 50 normal samples, while GSE14520 consists of 225 cancer samples and 220 non-cancer samples. Following the differential analysis, 2113 and 769 DEGs were identified, respectively. Subsequently, these genes were intersected with 2086 TRGs, leading to the discovery of 34 TRGs with varying expression levels (Fig. 2.A-C). To elucidate the biological roles of these 34 TRGs, we conducted GO and KEGG, which revealed that these genes have significant functions during the cell cycle and mitosis, involving key events in chromosome organization, cytoplasm and nucleus, such as DNA binding, chromatin binding, and protein kinase binding. Key pathways implicated were the Cell cycle, Oocyte meiosis, Cellular senescence, and the p53 signaling pathway, further underscoring their importance in cell division and the transmission of genetic material (Fig. 2.D-G).

Fig. 2

Assessment of the expression profiles of TRGs in HCC. (A-B) Volcano plots depicting differentially expressed genes from the TCGA-LIHC and GSE14520 databases. (C) Venn diagram illustrating the intersection between differentially expressed genes(DEG) in TCGA-LIHC and GSE14520 datasets with TRGs. (D-G) GO and KEGG analyses of the variably expressed TRGs

Filtering the optimal models through machine learning

Among the 34 intersecting TRGs, we further extracted 25 TRGs significantly associated with survival prognosis (p < 0.01). In creating prognostic models for the TCGA-LIHC cohort, a range of machine learning techniques, including Random Survival Forest, CoxBoost, Elastic Net, and LASSO-Cox, were employed. Upon model development, comparison analysis demonstrated that the predictive model based on LASSO-Cox exhibited the highest predictive power, achieving the highest C-index score among the different models evaluated (Table 1).

Table 1 C-index of machine learning methods

Developing and confirming a TRG-Based risk model

The LASSO-Cox regression analysis identified six genes, namely CDC20, TRIP13, EZH2, AKR1B10, ESR1, and DNAJC6, to form the optimal prognostic signature related to telomere (Fig. 3.A-B). To calculate the RS, the expression and coefficiention of these six genes were extracted. The equation is as follows: TRGs RS=(0.1391* CDC20 expression)+ (0.1574*TRIP13 expression) + (0.0109*EZH2 expression) +(0.0739*DNAJC6 expression) + (0.0203*AKR1B10 expression) + (-0.0134*ESR1 expression). The cohort was stratified into low and high-risk groups using the median RS. Subsequent, KM survival curve analysis of six genes within the TCGA-LIHC cohort disclosed pronounced variations in survival times contingent on gene expression levels. High expression of CDC20, TRIP13, EZH2, AKR1B10, and DNAJC6 was associated with decreased overall survival (OS) when contrasted with patients showing low expression. Conversely, high expression of ESR1 was associated with prolonged OS time (Fig. 3.C-H). Additionally, we conducted a visual exploration of the expression profiles of six TRGs, comparing low-risk and high-risk groups within both validation cohorts, GSE14520 and ICGC, while also analyzing the relationship between risk factors and survival (Fig. 3.I-K).

Fig. 3

Building and confirming a risk prediction model utilizing TRGs.(A) LASSO regression plots independent variable coefficients against the log-transformed regularization parameter λ on the x-axis.(B) For each λ, LASSO regression provides confidence intervals for the coefficients, depicted on the y-axis.(C-H) Analysis of KM survival curves for CDC20, TRIP13, EZH2, AKR1B10, and DNAJC6 in the TCGA-LIHC dataset suggested significant differences in overall survival (OS) time.(I-K) Survival profiles and the distribution of TRGs risk scores across the three cohorts are examined. Patients are ordered by their TRGs risk scores, with the vital status of each individual patient categorized by their risk score as illustrated in the central figure(I: TCGA-LIHC cohort; J:GSE14520;K: ICGC cohort). P< 0.05 is statistically different

Evaluating of model robustness and prognostic efficacy

To evaluate the model’s generalization capability and robustness, we performed validation using TCGA-LIHC, GSE14520, and ICGC cohorts. Following a similar methodology as described before, the samples were divided into two categories: high and low-risk groups. Patients in the high-risk group had a markedly unfavorable outcome relative to the low-risk group within the cohorts (Fig. 4.A-C). Subsequently, to appraise the model’s predictive performance, ROC curves were developed for patient prognosis prediction. As shown in Fig. 4. D-F, For the TCGA-LIHC cohort, the AUC values at one, three, and five years were 0.85, 0.78, and 0.73, respectively. Corresponding values in the GSE14520 cohort were 0.70, 0.69, and 0.71. The AUC of ICGC cohort were 0.86, 0.69 and 0.71. In conclusion, these results demonstrated that the prognostic model we developed possesses considerable robustness and prognostic efficacy in HCC.

Fig. 4

Assessing the durability of the model and its predictive effectiveness in various sets. (A-C) Survival curves of sufferers in the different risk groups from the TCGA cohort, GSE14520 and ICGC cohorts. E-F) AUC curves at one, three, and five years for models trained on TCGA data, as well as cohorts from GSE14520 and ICGC. P< 0.05 is statistically different

Association between RS model and clinical features of HCC

To determine the independent prognostic significance of RS, we conducted both univariable and multivariable Cox regression analyses of patients’ clinical profiles, encompassing age, stage, and gender. Combining the results, age and gender showed no significant impact on survival outcomes. In contrast, disease stage emerged as a significant predictor of survival, with higher stages associated with increased mortality risk. Furthermore, RS demonstrated independent prognostic significance, with higher scores indicating elevated mortality risk (Fig. 5.A-B). These findings underscored the importance of disease staging and RS assessment in predicting survival outcomes in the studied population.

Fig. 5

The prognostic value of RS is distinct and significant for understanding clinical attributes. (A-B) The forest plot presents the findings from Cox regression analysis conducted on clinical features from the TCGA dataset

Analysis of immune cell infiltration associated with TRGs RS

The tumor microenvironment (TME) is a pivotal part in either controlling or promoting tumor growth [27]. Hence, utilizing CIBERSORT analysis, we compared the TME variations between patients categorized as high risk versus low risk, focusing on the comparative richness of 22 immune cell types in the tumor-infiltrating population in HCC cases. Histogram demonstrated the profile of the 22 immune cell types in patients belonging to the different risk groups (Fig. 6.A). We observed variations in immune cell infiltration, including Macrophages M0, resting NK cells, regulatory T cells (Tregs), Monocytes, activated CD4 memory T cells, resting Mast cells, T follicular helper cells, and resting CD4 memory T cells. Particularly, macrophages M0 and Tregs were found at elevated levels in the high-risk group when contrasted with the low-risk group (Fig. 6.B). These findings validate the elevated presence of Macrophages M0 and Tregs, suggesting their potential involvement as infiltrating immunosuppressive cells in HCC characterized by high TRGs RS. In parallel, Immune checkpoints play a crucial role in HCC, enabling malignant cells to exploit the signaling pathway to evade immune system attacks, thus promoting tumor growth and spread [28]. To determine the correlation between the six TRGs in the model and the immune checkpoint genes, correlation heat maps were drawn. The results suggested a significant correlation between CDC20, TRIP13, EZH2, and DNAJC6 with immune checkpoints(Fig. 6.C).

Fig. 6

Inspecting the model for its immune cell infiltration characteristics. (A) The histogram demonstrated the ratio of 22 immune cell infiltrations in the different risk groups.(B) Box line graph showing the results of different subgroups of each type of immune cell to make a comparison. (C) Correlation heatmap between risk model signature and immune checkpoint genes. ^*p < 0.05, ^**p < 0.01, ^***p < 0.001

Variations in enriched pathways among patients of distinct risk groups

To continue exploring the discrepancy in enrichment pathways between cases at high- and low-risk, we computed the scores for HALLMARK and KEGG pathway enrichment analysis. These analyses will help unravel the fundamental mechanism of HCC development. Heat mapping of different risk groups was drawn by utilizing the “limma” R package. In high-risk patients, we have observed that TGF-Beta signaling, cell cycle pathway, and DNA repair signaling pathway were crucial to HCC, and significantly impacted the development and progression of HCC (Fig. 7.A-B).

Fig. 7

Comparative analysis reveals differences in pathway enrichment within the both groups. (A-B) Heatmap showing HALLMARK and KEGG pathway enrichment analyses in distinct risk groups

Predicting sensitivity to drugs and immunotherapy response

Using the TCGA-LIHC database, we forecasted the IC50 values of agents, which were sourced from GDSC. The boxplot illustrated significant disparities in the response to chemotherapeutic agents between high- and low-risk groups using the variances in IC50 value distribution. Among high-risk HCC patients, the prognosis was notably worse, so we clearly found that MK-1775_1179, AZD7762_1022, Telomerase Inhibitor IX_1930, Pevonedistat_1529, Buparlisib_1873, BI-2536_1086 were the drugs with higher sensitivities for the high-risk group. (Fig. 8.A). Next, we visualized the association between the gene transcription of six TRGs and the reaction of drugs through a correlation heatmap (Fig. 8.B). Meanwhile, to investigate the potential of model for immunotherapy, we employed the IMvigor210CoreBiologies package to procure a dataset comprising of a group of immunotherapy responders and categorized them into high- and low-risk groups based on their RS. As depicted in Fig. 8.C, immunotherapy was more available in high-risk group, with a complete response or partial response in 58% of cases and stable disease or progressive disease in 42%, compared to the low-risk group where CR/PR was 31% and SD/PD was 69%.

Fig. 8

Prediction of sensitivity and analysis of immunotherapeutic response of drugs across various risk groups. (A) Box plots represented sensitivity analyses for drugs in both risk groups. (B) Relationship heatmap between model genes and the reaction of drugs. (C) Histograms illustrated the distribution patterns of immunotherapy outcomes across various patient risk groups

Discussion

Recent studies have revealed that targeting telomeres in therapy effectively eliminates cancer-triggering cells and boosts the host immune response against tumors [29]. Hence, telomeres are crucial in the progression of malignant neoplasm. Telomere shortening contributes to genomic instability through multiple mechanisms, including the exposure of chromosomal ends and the triggering of DNA damage responses (DDR), such as the formation of γ-H2AX foci and the activation of ATM kinase. These DDR activations lead to chromosomal end-to-end fusions, breakage-fusion-bridge cycles, and large-scale genomic rearrangements, which are hallmarks of cancer initiation and progression [7, 8]. These mechanisms play a crucial role in the development and progression of HCC and are closely related to TERT promoter mutations and telomere lengthening mechanisms, such as alternative lengthening of telomeres(ALT) [11, 12].Previous literature has reported the favorable prognostic value of TRGs-based prognostic models in renal cell carcinoma [14], lymphoma [30], and lung adenocarcinoma [13], but there is limited evidence regarding HCC. The objective of this research was to assess the prognostic significance of TRGs in HCC. We developed a predictive model utilizing TRGs and determined that it could be instrumental in guiding the selection of immunotherapy for HCC patients.

In this study, using differentially expressed TRGs, including CDC20, TRIP13, EZH2, DNAJC6, AKR1B10, and ESR1, were developed for construction of a prognostic signature. These genes were documented to exert a substantial impact on the progression of HCC. CDC20 was a pro-oncogenic factor that was remarkably up-regulated in HCC and contributed to radioresistance in P53-mutant HCC cells by regulating the Bcl-2/Bax signalling pathway, thereby affecting tumor outcomes [31, 32]. TRIP13 interacted with and amplified ACTN4 expression in HCC, thereby initiating the AKT/mTOR signaling cascade, which in turn stimulated cell proliferation, migration, and invasion, as well as facilitated epithelial-mesenchymal transition (EMT), collectively contributing to tumor progression [33]. EZH2, functioning as an epigenetic modulator, upregulated CXCR4 expression by inhibiting miR-622, which subsequently triggered the activation of the CXCR4/CXCL12 and PI3K/AKT signaling axes, thereby enhancing the proliferation and metastatic potential of HCC cells [34]. DNAJC6 triggered the TGF-β signaling pathway, contributing to the EMT and furthering the advancement of HCC [35]. AKR1B10 promoted HCC proliferation through glycolytic pathway activation [36]. Furthermore, its serum proteins served as superior markers for HCC compared to AFP, offering higher sensitivity and specificity with a shorter half-life [37]. Overexpression of ESR1 promoted hepatocarcinogenesis and cell proliferation, while also influencing the behavior of tumor cells through feedback regulation of the CCL23 and AKT signaling pathways [38]. In summary, these signature genes are crucial in HCC development and could be valuable markers for prognosis and treatment strategies.

Patients within the validation set were stratified into high- and low-risk categories based on the median RS threshold. The low-risk group exhibited notable higher survival prognosis than the high-risk group, suggesting that our model had superior prognostic predictive efficacy, as demonstrated by a higher area under the curve (AUC). Moreover, RS was recognized as a standalone prognostic determinant for HCC outcome. To investigate the TME in HCC sufferers from various risk groups, we performed immune infiltration analyses. We noted elevated macrophage M0 and Treg levels in the high-risk cohort relative to the low-risk cohort, aligning with existing studies that have documented their heightened infiltration within HCC tissues, which is linked to poorer prognostic outcomes [39,40,41]. These cells actively promoted HCC growth, metastasis, and immune evasion by suppressing other immune cell activity and attenuating the anti-tumor response [39, 41]. Our findings revealed a positive association between macrophage M0 and Tregs levels and the risk scores.

Molecules that inhibited immune responses, including PD-1 and CTLA-4, facilitated tumor cell escape from immunosurveillance, thereby impeding their efficient eradication [42, 43]. Hence, we conducted an in-depth analysis exploring the association between immune checkpoints and the six TRGs of the model. Intriguing, we found a obvious association between CDC20, TRIP13, EZH2, DNAJC6 and immune checkpoint. It has been reported in the literature that inhibition of CDC20 in tumor models increases CD8⁺T lymphocyte infiltration and enhances anti-PD-1 immunotherapy [44]. EZH2 inhibited CD4⁺T, CD8⁺T, and NKT cell differentiation, migration, and enhanced Treg inhibitory activity to promote immune escape, thereby affecting immune checkpoint expression and the efficacy of immunotherapy [45, 46]. Our research findings also demonstrated a correlation between these genes and immunotherapy, indicating a greater effectiveness of immunotherapy in the high-risk group. As a result, immunotherapy emerges as a beneficial choice for high-risk patients.

Of course, our study has still some limitations. Firstly, this study is currently limited to bioinformatics analyses based on public datasets. Future collaborations with clinical institutions will be prioritized to validate the RS model in prospective cohorts and real-world samples. All data was acquired from retrospective studies in public database, which may result in some bias. Hence, our findings need to be validated through prospective studies. Secondly, further investigation is required to clarify the intricate regulatory pathways and functional roles of the TRGs in HCC.

Conclusion The RS model, concentrating on TRGs, successfully forecasted patient outcomes and immune reactions in HCC, offering novel therapeutic approaches for oncology.