Mengyao Sha
Jun Chen
Haifeng Hou1- 1Department of Laboratory Medicine, Suzhou Yongding Hospital, Suzhou, China
- 2Department of Hematology, Suzhou Yongding Hospital, Suzhou, China
Background: T-cell suppression in patients with Acute myeloid leukemia (AML) limits tumor cell clearance. This study aimed to explore the role of T-cell senescence-related genes in AML progression using single-cell RNA sequencing (scRNA-seq), bulk RNA sequencing (RNA-seq), and survival data of patients with AML in the TCGA database.
Methods: The Uniform Manifold Approximation and Projection (UMAP) algorithm was used to identify different cell clusters in the GSE116256, and differentially expressed genes (DEGs) in T-cells were identified using the FindAllMarkers analysis. GSE114868 was used to identify DEGs in AML and control samples. Both were crossed with the CellAge database to identify aging-related genes. Univariate and multivariate regression analyses were performed to screen prognostic genes using the AML Cohort in The Cancer Genome Atlas (TCGA) Database (TCGA-LAML), and risk models were constructed to identify high-risk and low-risk patients. Line graphs showing the survival of patients with AML were created based on the independent prognostic factors, and Receiver Operating Characteristic Curve (ROC) curves were used to calculate the predictive accuracy of the line graph. GSE71014 was used to validate the prognostic ability of the risk score model. Tumor immune infiltration analysis was used to compare differences in tumor immune microenvironments between high- and low-risk AML groups. Finally, the expression levels of prognostic genes were verified using polymerase chain reaction (RT-qPCR).
Results: 31 AMLDEGs associated with aging identified 4 prognostic genes (CALR, CDK6, HOXA9, and PARP1) by univariate, multivariate, and stepwise regression analyses with risk modeling The ROC curves suggested that the line graph based on the independent prognostic factors accurately predicted the 1-, 3-, and 5-year survival of patients with AML. Tumor immune infiltration analyses suggested significant differences in the tumor immune microenvironment between low- and high-risk groups. Prognostic genes showed strong binding activity to target drugs (IGF1R and ABT737). RT-qPCR verified that prognostic gene expression was consistent with the data prediction results.
Conclusion: CALR, CDK6, HOXA9, and PARP1 predicted disease progression and prognosis in patients with AML. Based on these, we developed and validated a new AML risk model with great potential for predicting patients’ prognosis and survival.
1 Introduction
Acute myeloid leukemia (AML) is a heterogeneous hematological neoplasm characterized by the clonal proliferation of myeloid cells in peripheral blood and bone marrow (Pollyea et al., 2021), with a 5-year survival rate of less than 30% (Narayanan and Weinberg, 2020). Cellular senescence is an irreversible state of cell cycle arrest characterized by DNA damage, telomere shortening, and secretion of pro-inflammatory cytokines (Rossiello et al., 2022). Recent studies have emphasized that T-cell dysfunction, including senescence, is a key contributing factor to immune evasion and treatment resistance in AML (Chen et al., 2025; Ikeda et al., 2025; Yuan et al., 2021), as evidenced by the poor prognosis and low survival of patients with AML (Mazziotta et al., 2024).
Previous studies have identified markers to study T-cell senescence in some patients with AML, such as increased KLRG1 expression (Shive et al., 2021; Towers et al., 2024). These studies focused on identifying senescence markers; however, knowledge of the impact of T-cell senescence and subpopulation changes on the immune microenvironment of AML tumors and AML prognosis remains inadequate. In this study, we identified AML prognostic genes related to T-cell senescence using transcriptome (GSE116256) and single-cell sequencing (GSE114868) data and constructed a prognostic risk model. We also developed a predictive line graph for viability using TCGA-LAML and another transcriptome data (GSE71014). We then analyzed the differences in the immune microenvironment between the low- and high-risk AML groups. The relationship between T cell subsets and prognostic genes, and which T cell subsets were most significantly associated with AML, was evaluated, predicting target lncRNAs, miRNAs, and drugs. Finally, we concluded that prognostic gene expression influences T cell senescence, altering the immune microenvironment. This causes tumor immune escape, thereby determining the prognosis and survival duration of patients with AML. Our findings provide new insights into the role of T-cell senescence in AML and suggest potential therapeutic targets for improving immunotherapy-based treatments.
2 Materials and methods
2.1 Data collection
We downloaded gene expression data from the GSE114868, GSE116256, and GSE71014 datasets associated with AML through the Gene Expression Omnibus (GEO) repository (http://www.ncbi.nlm.nih.gov/geo/). Eight hundred and sixty-six cellular senescence-related genes (CSRGs) were extracted from the CellAge database (Avelar et al., 2020; Chatsirisupachai et al., 2019) (https://genomics.senescence.info/cells/) (Supplementary Table S1). The survival data and gene expression profiles of the TCGA-LAML cohort were obtained from the TCGA database (https://xena.ucsc.edu/) (Table 1).
Table 1. Dataset collection.
2.2 Differential expression analysis
Differentially expressed genes (DEGs) in the AML and control groups from the GSE114868 dataset were identified using the Limma software package (v 3.54.0). All hypothesis tests were corrected using the Benjamini–Hochberg method to control the false discovery rate (FDR <0.05), and the screening criteria for differentially expressed genes were |log 2 fold change (FC)|>1, adj. p < 0.05 (Ritchie et al., 2015). The ggplot2 (Gustavsson et al., 2022) and ComplexHeatmap (v 2.14.0) (Gu and Hübschmann, 2022) packages were used to visualize the top ten significantly DEGs via a volcano plot and a heatmap, respectively.
2.3 Single-cell analysis
In the GSE116256 scRNA-seq dataset, cells were curated using Seurat’s NormalizeData function (Gu and Hübschmann, 2022), selecting high-quality cells based on stringent criteria (nFeature RNA >200, genes expressed in <3 cells removed, nCount RNA between 200–3,000, and mitochondrial proportion <10%). The FindVariableFeatures function identified the top 2,000 highly variable genes (HVGs). Subsequent normalization and PCA outlier detection were conducted using Seurat’s Scale Data function, with Elbowplot visualization and Jackstraw reclustering used to determine the top 30 significant principal components (PCs) (p < 0.05). Dimensionality reduction and clustering were performed using a UMAP with 30 PCs and Seurat’s FindNeighbors and FindClusters functions (resolution = 0.4). Notable marker genes for distinct clusters were identified using the FindAllMarkers function (logfc. Threshold = 0.5, min. pct = 0.25, return. thresh = 0.01). Cell annotation was assigned based on these markers using SingleR (v 2.4.0) (Aran et al., 2019) and the CellMarker database. Finally, to investigate the differences in gene expression at the T cell level between AML and controls, we screened the DEGs of T cells in AML and controls using the FindMarkers function. We depicted them using Manhattan plots (|log2 FC| > 0.1, adj. p < 0.05).
2.4 Identification and functional analysis of candidate genes and protein-protein interaction (PPI) networks
We used the VennDiagram package (v 1.7.3) (Chen and Boutros, 2011) to intersect T cell DEGs with CSRGs and GSE114868 DEGs. The ClusterProfiler package (v 4.7.1.3) (Yu et al., 2012) was used to perform Gene Ontology (GO) and Kyoto Encyclopedia of the Genome (KEGG) enrichment analyses (p < 0.05). We visualized the data using the GO plot (v 1.0.2) (Walter et al., 2015) and enrichment (v 1.22.0) packages (Wang et al., 2022). The interacting gene search tool (STRING database [http://www.string-db.org/]) and Cytoscape software (v 3.7.2) (Liu et al., 2020) were used to build protein-protein interaction (PPI) networks with a confidence threshold >0.4.
2.5 The screening of the prognostic genes and construction of a risk model
We refined candidate genes linked to survival and prognosis using a univariate Cox regression analysis (HR ≠ 1, p < 0.05). We validated the proportional hazards (PH) assumption (p > 0.05) using the survival package (v 3.7-0) (Lei et al., 2023) and cox. zph function on TCGA-LAML data. This was followed by multivariate and stepwise Cox regression analyses (p < 0.05) to identify prognostic genes for AML risk model development.
Where coef is the gene regression coefficient and expr is the gene expression level.
Patients in the TCGA-LAML and GSE71014 databases were categorized into high- and low-risk groups based on the median risk score. We plotted risk and Kaplan-Meier Survival Curves using the survminer package (v 0.4.9) (Ramsay et al., 2018). The AUC was determined by ROC analysis of the 1-, 3-, and 5-year survival using the survivalROC package (v 1.0.3.1) (Heagerty et al., 2000).
2.6 Independent prognostic analysis and line graph
We performed a univariate Cox regression analysis, which included risk scores, age, and type. A HR≠1 and p < 0.05 represented the prognostic significance of the model, and the model met the proportional hazards assumption (p > 0.05). We used the survival package (v 3.7-0) for the multivariate analysis and the rms package (v 6.5.0) (Heagerty et al., 2000) to construct line plots to predict the 1-, 3-, and 5-year survival of patients with AML. We plotted ROC curves to verify the reliability of the line graphs using the plotROC package (v 2.3.1) (Sachs, 2017).
2.7 Tumor immune microenvironment analysis
The estimate package (v 1.0.13) and Wilcoxon test assessed the AML immune microenvironment. The ssGSEA (Charoentong et al., 2017) and GSVA (v 1.46.0) packages (Hänzelmann et al., 2013) calculated the enrichment of 28 immune cell populations, and the Wilcoxon test compared their ratios in TCGA-LAML samples (p < 0.05). We performed a Spearman correlation analysis using the psych package (v 2.4.3) to assess the association between immune cells and prognostic genes (|cor| > 0.30, p < 0.05).
2.8 Regulatory network construction
We used the MultiMiR package (v 1.24.3) (Ru et al., 2014) to predict mRNAs targeted by miRNAs based on the TargetScan and PITA databases. We used starBase (https://rnasysu.com/encori/) to identify the upstream lncRNAs of miRNAs and visualized “lncRNA-miRNA-mRNA” networks using Cytoscape (v 3.9.1) (Liu et al., 2020). We explored prognostic gene interactions and co-expression using the geneMANIA (https://genemania.org/) database.
2.9 Chromosomal and subcellular localization
The chromosomal location of prognostic genes was determined using ENSEMBL (https://asia.ensembl.org/index.html) and mapped to chromosomes using the RCircos package (v 1.2.2) (Zhang et al., 2013). Proteins were extracted from Genecards (https://www.genecards.org/) sequences to study the subcellular localization of proteins encoded by prognostic genes, and the ggplot2 package (v 3.4.1) was used to visualize the cellular distribution of the proteins.
2.10 Drug prediction and molecular docking
We assessed drug sensitivity using the GDSC database. The pRRophetic package (v 0.5) (Geeleher et al., 2014) determined the IC50 values, and a Wilcoxon test compared drug sensitivities between different risk groups (p < 0.05). The corrplot package (v 0.92) (Wang et al., 2022) performed a Spearman’s test (|cor| > 0.30, p < 0.05) to analyze the correlation between IC50 values and gene expression. Highly correlated 3D drug structures were obtained from the PubChem and RSCB PDB databases. CB-Dock2 was used to evaluate affinity, with a binding energy below −5 kcal/mol indicating affinity.
2.11 Experimental validation by RT-qPCR
Bone marrow samples from controls and patients with AML were collected at Suzhou Yongding Hospital, and informed consent and ethical approval were obtained (202450; see Supplementary Table S3 for sample details).
2.12 Statistical analysis
Data processing and comparisons were performed using R (v 4.2.3), with statistical significance determined by Wilcoxon tests at p < 0.05.
3 Results
3.1 Five hundred and seventeen T-cell DEGs identified in the scRNA-seq dataset GSE116256 and T-cell subcluster analysis
The GSE116256 dataset contained 21,253 cells and 19,003 genes for subgroup identification and annotation afterquality control (QC) (Supplementary Figure S1A). Supplementary Figure S1B shows the inclusion of 2,000 high-variance genes (red dots). Supplementary Figure S1C shows no outlier samples in GSE116256. An elbow plot was used to visualize the 2,000 high-variance genes after principal component analysis. The optimal number of dimensions for cell clustering was determined to be 30 (p < 0.05) using linear dimensionality reduction (Supplementary Figure S1D,E). After co-clustering the cells into 20 different cell clusters (Figure 1A), bubble plots showed the expression of typical marker genes in 15 cell clusters (Figure 1B, Supplementary Table S2). The annotated cell groups are shown in Figure 1C. Among them, monocytes had three clusters (7, 13, and 14), and T cells had two clusters (6 and 11). After a differential analysis of the genes in the subclustered cells, clusters 7 and 13 were shown to contain 340 and 736 DEGs (cluster 14 did not), respectively, based on the screening according to predefined criteria (|logFC| > 0.1 and adj. p < 0.05; Supplementary Figure S2A,B). A functional analysis revealed that the two clusters of monocytes were involved in significantly different biological processes and pathways (Supplementary Figure S2C–F). Given that the impact of T-cell senescence and subpopulation changes in the AML tumor immune microenvironment is still poorly understood, we clustered the T-cell subpopulations based on T-cell-typical marker gene expression (Szabo et al., 2019) (Figures 1D,E). The annotated cell clusters are shown in Figure 1F. The Manhattan plot suggests 517 T-cells DEGs between the control and AML groups (358 upregulated and 159 downregulated genes) (Figure 1G; Supplementary Table S4).
Figure 1. Five hundred and seventeen T-cell DEGs identified in the scRNA-seq dataset GSE116256 and T-cell subcluster analysis. (A) UMAP plot of 20 different cell clusters obtained by dimensionality reduction clustering before annotation. (B) Bubble plot of gene expression of different cellular markers. (C) The 15 major cell types identified after annotation. (D–F) Further dimensionality reduction and clustering of T cell cluster 6 revealed four distinct T cell types. (G) Manhattan plot of 517 DEGs in T cell clusters.
3.2 Screening of AML candidate genes based on the GSE114868, GSE116256, and CellAge databases
The GSE114868 dataset yielded 2,864 DEGs (downregulated: 1,552; upregulated: 1,312), with the top ten DEGs highlighted in the volcano and heat maps (Figures 2A,B). The intersection between the 2,864 DEGs with 866 CSRGs and 517 T-cell DEGs yielded 31 candidate genes (Figure 2C; Supplementary Table S5). The enrichment analysis and PPI networks of the 31 candidate genes is shown in Figures 2D–F.
Figure 2. Screening of AML candidate genes based on the GSE114868, GSE116256, and CellAge databases. (A) Volcano plot of DEGs in GSE114868, with genes labeled according to significance; (B)The heatmap consists of two parts: the upper part shows the density plot of DEGs expression with lines representing five quantiles and the mean; the lower part is the DEGs expression heatmap (red for AML group, green for Control group). (C) Venn diagram showing the intersection of DEGs from GSE116256, GSE114868, and CSRGs; (D–F) GO/KEGG enrichment and PPI network analysis of 31 cross-genes.
3.3 Establishment and evaluation of a prognostic model for AML based on the TCGA-LAML cohort and validation of a prognostic model for AML based on GSE71014
Using multifactorial models that directly incorporate too many variables resulted in model instability and overfitting. Therefore, we first performed a univariate Cox regression analysis on the TCGA-LAML cohort (HR = 1: no significant effect; HR < 1: protective genes; HR > 1: risk genes; p < 0.05) to initially screen for variables that were significantly associated with AML survival. The results of the univariate Cox regression analysis were then subjected to the proportional hazards (PH) hypothesis test (p > 0.05) (Supplementary Figure S3), and six genes were finally identified: CALR, CDK6, CTSD, HOXA9, PARP1, and SAMHD1 (Figure 3A). A multifactor Cox model was constructed for the six genes to obtain a more robust prediction model (Figure 3B). A stepwise regression analysis was performed based on the multifactor Cox regression analysis (Figure 3C), which resulted in the final identification of four prognostic genes: CALR, CDK6, HOXA9, and PARP1. The formula for the risk model was
Figure 3. Establishment and evaluation of a prognostic model for AML based on the TCGA-LAML cohort and validation of a prognostic model for AML based on GSE71014. (A) Six prognostic genes identified via univariate analysis; (B,C) Four prognostic genes determined by multivariate Cox and stepwise regression analyses; (D) Risk score distribution and Survival status distribution (The x-axis of both panels represents patients sorted by risk score. The y-axis of the upper panel shows risk scores, with a dashed line indicating the median risk score and its corresponding patient count. The y-axis of the lower panel shows survival time, with a dashed line indicating the median risk score and its corresponding patient count). (E) Kaplan-Meier curve; (F) 1-, 3-, 5-year ROC curves; (G–I) Validation in GSE71014 with corresponding risk score and survival analysis, and ROC curve assessment at 1, 3, and 5 years.
The risk coefficients are shown in Table 2.
Table 2. Risk model gene coefficients.
Risk scores were calculated for each patient based on the regression coefficients and expression levels of the prognostic genes, and the TCGA-LAML samples were categorized into high- and low-risk groups (median risk score of 0.1134133, high/low-risk patients = 66/66). Risk distribution and survival status maps showed that the high-risk group had higher risk scores and shorter survival durations than the low-risk group (Figure 3D). The KM curves showed that the high-risk group had a significantly lower probability of survival (Figure 3E). The AUC values of the ROC analyses of the 1-, 3-, and 5-year risk models were ≥0.6, suggesting a strong validity in predicting the survival of patients with AML (Figure 3F). The KM curves showed that the high-risk group had a higher risk score and a shorter survival duration than the low-risk group (Figure 3D).
In the independent cohort (GSE71014) split into high- and low-risk groups (52 samples each) with a median score of 3.269811, the high-risk group had more adverse outcomes than the low-risk group with an AUC of ≥0.6 (Figures 3G–I), validating the efficacy of the risk model.
3.4 Establishing a line graph for AML survival prediction based on independent risk factors
To evaluate the prognostic model and identify independent prognostic factors for AML, a one-way Cox regression analysis (HR ≠ 1, p < 0.05) was performed using the risk score and age. M-staging (M3 and others) was also performed on the TCGA-LAML cohort. We found that the risk score, age, and M-stage were significantly associated with the survival of patients with AML (Figure 4A), and the PH hypothesis (p > 0.05) was fulfilled. The risk score and age were the independent prognostic factors in the multifactorial Cox independent prognostic analysis (Figure 4B). Based on the independent prognostic factors, a column chart was constructed with survival as the outcome event (Figure 4C). Higher scores in the column chart indicated a higher risk of death and a lower survival rate. For example, if a patient was 60 years old (points = 4.5) with a risk score of 1.5 (points = 0; total points = 4.5), the probability of surviving for 1, 3, and 5 years was approximately 50%, 20%, and 10%, respectively. ROCs were evaluated on the column line graphs, and the AUCs were all greater than 0.8, indicating good model predictions (Figure 4D).
Figure 4. Establishing a line graph for AML survival prediction based on independent risk factors. (A,B) Risk score, age, and staging were included as clinical characteristics in univariate and multivariate Cox regression analyses. (C) Line graph for risk score and age and 1-, 3-, and 5-year survivability with AML. (D) ROC curves for 1-, 3- and 5-year survivability with AML.
3.5 Differences in tumor immune microenvironment between the high- and low-risk groups and prognostic gene expression in different cell subclusters
We found a significant upward trend (p < 0.05) in the immunization and ESTIMATE scores in the high-risk group of the TCGA-LAML cohort (Figure 5A). We compared the abundance of immune cell infiltration in the high- and low-risk groups. Twelve immune cell types showed significant differences, including activated dendritic cells, memory T-cells, and macrophages (p < 0.05; Figures 5B,C). Additionally, CALR, CDK6, and PARP1 were negatively correlated with most of the differential immune cells, whereas HOXA9 was positively correlated with plasmacytoid dendritic cells, macrophages, and central memory CD4+ T cells (p < 0.05; Figure 5D).
Figure 5. Differences in tumor immune microenvironment between the high- and low-risk groups and prognostic gene expression in different T cell subclusters. (A) Stromal, immune, and ESTIMATE scores were compared between low- and high-risk cohorts. (B,C) Immune cell infiltration levels were contrasted across low- and high-risk groups. (D) The correlation between prognostic genes and immune cells was assessed.
3.6 Functional analysis of molecular regulation of prognostic genes
The “lncRNA-miRNA-mRNA” regulatory network showed that CALR, CDK6, HOXA9, and PARP1 predicted six, six, two, and one miRNAs, respectively. The miRNAs predicted 0, 40, 28, and 41 lncRNAs associated with prognostic genes, including the XIST-“hsa-miR-324-5p”-HOXA9 and XIST-“hsa-miR-105-5p”-PARP1 relationships (Figure 6A). The biological functions and co-expression networks of the prognostic genes were analyzed using the GeneMANIA database, with CANX, CDKN2C, and CCND3 having the strongest association with prognostic genes (Figure 6B). Chromosomal localization revealed that CALR, CDK6, HOXA9, and PARP1 were located on chromosomes nineteen, one, seven, and seven, respectively (Figure 6C). Subcellular localization of proteins revealed that the four prognostic genes were mainly expressed in the cytoplasm and nucleus (Figure 6D).
Figure 6. Functional analysis of molecular regulation of prognostic genes. (A) lncRNA-miRNA-mRNA interaction network. (B) Functional and co-expression analysis of prognostic genes. (C,D) Genomic and subcellular mapping of prognostic genes.
3.7 IC50 difference analysis and targeted drug prediction based on AML prognostic modeling
IC50 difference analyses assess drug sensitivity. Lower IC50 values (concentration of drug required to inhibit tumor cell growth) indicate higher drug sensitivity (Garnett et al., 2012). By comparing the IC50 values of drugs in the high- and low-risk AML groups of the prediction model, we found that the high-risk group had higher sensitivity to chemotherapeutic drugs (Dactolisb, Gemcitabine, GNE.317, and 5-fluorouracil). In contrast, the low-risk group was more sensitive to entinostat (Figure 7A). Figure 7B shows the correlation between the IC50 of drugs and prognostic genes. PARP1 was negatively correlated with the IC50 of most chemotherapeutic drugs.
Figure 7. IC50 difference analysis and targeted drug prediction based on AML prognostic modeling. (A) Comparative IC50 analysis predicts differential drug sensitivity between low- and high-risk groups. (B) Assesses the correlation between prognostic gene expression and drug IC50 values. (C–F) Employs CB-Dock2 for molecular docking of prognostic genes to elucidate potential drug interactions.
The prognostic genes CALR, CDK6, HOXA9, and PARP1 were identified as four targets for CB-Dock2 molecular docking, and the prognostic genes had good binding with chemotherapeutic drugs. Specifically, CALR had a binding energy of −9.9 kcal/mol with IGF1R, with residues E60 and K151 forming hydrogen bonds with IGF1R (Figure 7C). The binding energy of CDK6 to ABT737 was −10.1 kcal/mol, with residues E127, G152, and P114 forming hydrogen bonds with ABT737 (Figure 7D). The binding energy of HOXA9 to pevonedistat was −8.9 kcal/mol, with residues T276, V277, and K265 forming hydrogen bonds with pevonedistat (Figure 7E). PARP1 binds to IBRD9 with a binding energy of −10.6 kcal/mol, with residues M890, G888, N868, and S864 forming hydrogen bonds with I. BRD9 (Figure 7F, Table 3).
Table 3. Binding energy.
3.8 Experimental validation of CALR, CDK6, HOXA9, and PARP1
We validated the mRNA expression of prognostic genes by extracting RNA from fresh bone marrow from healthy individuals and patients with AML. Compared to the control group, the AML group had significantly reduced CALR and CDK6 expression and increased HOXA9 and PARP1 expression (Figure 8), consistent with the predictions of our risk model.
Figure 8. Experimental validation of CALR, CDK6, HOXA9, and PARP1. (A–D) RT-qPCR validation of prognostic genes CALR, CDK6, HOXA9, and PARP1 in control and refractory AML cohorts.
Insert up to 5 heading levels into your manuscript as can be seen in “Styles” tab of this template. These formatting styles are meant as a guide, as long as the heading levels are clear, Frontiers style will be applied during typesetting.
4 Discussion
Prior research correlates AML risk with T-cell senescence, characterized by chronic immune activation and dysfunction (Bailur et al., 2020; Chen et al., 2024; Miao et al., 2024; Song et al., 2024). In this study, we screened candidate genes using single-cell and bulk RNA sequencing data. We constructed a prognostic model using four genes (CALR, CDK6, HOXA9, and PARP1) and a line graph for predicting survival by single-factor and multifactor Cox stepwise regression using TCGA-LAML and RNA sequencing data. We also elucidated the relationship between the tumor immune microenvironment and prognostic genes, revealing that cellular senescence was significantly associated with AML. This study aimed to provide new insights for the early diagnosis, personalized treatment, and prognostication of AML.
Thirty-one candidate genes were first identified using independent external data from GEO, TCGA, and CellAge databases. The signaling pathways significantly enriched with the candidate genes were the nucleotide-binding oligomeric structural domain (NOD)-like receptor and IL-17 signaling pathways. Studies have shown that the NOD-like receptor protein 3 (NLRP3) pathway is overexpressed and highly activated in AML cells (Chen et al., 2024; Zhong et al., 2021). IL-17 stimulates the development of granulocytes and myeloid hematopoietic cells (Wei et al., 2024), and small molecules targeting the TNF/IL-17/MAPK pathway (OUL35, KJ-Pyr-9, and CID44216842) attenuate zebrafish bone marrow proliferation (Luo et al., 2024). Univariate analysis identified six prognostic predictor genes, which multivariate Cox and stepwise regression analyses further validated. Using the above analyses, we identified the prognostic genes (CALR, CDK6, HOXA9, and PARP1). The prognostic risk model categorized patients into high-risk and low-risk groups, with high-risk patients experiencing a shorter survival period. The AUC values from the ROC analyses confirmed the accuracy of the model. The qPCR validation results were consistent with the initial predictions, emphasizing the reliability of the model predictions. IC50 difference analysis yielded significant differences in sensitivity to drugs between the high- and low-risk groups in the prognostic model. Molecular docking techniques identified targeted drugs with strong binding to the prognostic genes, including pevonedistat and azacitidine (first-line therapy for AML) (Adès et al., 2022; Murthy et al., 2024; Short et al., 2023). ABT737 neutralizes the inhibition of Bax and Bak by Bcl-2, Bcl-xl, and Bcl-w at nano-micro molar concentrations, thereby inducing apoptosis (Yalniz and Wierda, 2019).
In this study, CALR, CDK6, HOXA9, and PARP1 expression in pDCs, Mφs, and central memory CD4+ T cells significantly differed between low- and high-risk groups. CALR induces phagocytosis by sending an “eat-me signal,” and its downregulation contributes to immune evasion in AML (Bhave et al., 2023; Liu et al., 2021; Weinhäuser et al., 2023; Zhou et al., 2024). CDK4/6 inhibitors (palbociclib) slow AML progression by reducing DNA damage and telomere shortening in T cells by inhibiting CDK4 and CDK6 kinase activities (Nebenfuehr et al., 2020; Schmoellerl et al., 2020). CDK6 knockdown stalls DN T cell differentiation at the DN3/DN4 stage, resulting in insufficient thymocytes and structural atrophy (Yabas and Hoyne, 2023). HOXA9 is a transcription factor that can be used as a prognostic indicator for AML because its dysregulation is associated with AML progression (Deshpande and Zhu, 2023; Lai et al., 2020; Salik et al., 2020; Xiao et al., 2023), specifically by recruiting EZH2 to inhibit INK4A/B (cell cycle inhibitor) expression and subsequently affecting T-cell senescence (Collins and Hess, 2016a; b). PARP1 is involved in DNA repair and gene expression, and inhibiting PARP1 has anti-tumor effects in AML (Csizmar et al., 2021; Kontandreopoulou et al., 2021; Kumar et al., 2022; Paczulla et al., 2019; Wu et al., 2023).
5 Conclusion
This bioinformatics-driven study identified prognostic genes linked to T-cell senescence in AML and confirmed the findings using RT-qPCR. Using public single-cell and bulk RNA-seq data, we identified 31 candidate genes associated with inflammation and hematopoiesis. A four-gene prognostic model (CALR, CDK6, HOXA9, and PARP1) was constructed using Cox and stepwise regression to forecast AML survival. Molecular docking highlighted ABT737, IBRD9, IGF1R, and pevonedistat as potential therapeutic agents with high affinity for these genes.
Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.
Ethics statement
The studies involving humans were approved by Medical Ethics Committee of Suzhou Yongding Hospital (202450). The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
MS: Data curation, Software, Conceptualization, Writing – original draft, Methodology, Validation, Formal Analysis, Writing – review and editing. JC: Formal Analysis, Data curation, Investigation, Writing – review and editing, Validation, Supervision, Methodology. HH: Methodology, Supervision, Project administration, Investigation, Writing – review and editing, Formal Analysis. HD: Data curation, Investigation, Software, Formal Analysis, Supervision, Writing – review and editing, Project administration, Resources. YZ: Writing – review and editing, Funding acquisition, Resources, Investigation, Conceptualization, Formal Analysis, Project administration, Supervision, Validation, Methodology.
Funding
The author(s) declare that financial support was received for the research and/or publication of this article. Funding for this research was graciously provided by the Jiangsu Nutrition Society through their “2023 Nutrition Special Research Grant” (JYX2023011).
Acknowledgments
We acknowledge the valuable support provided by Suzhou Yongding Hospital.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Generative AI statement
The author(s) declare that no Generative AI was used in the creation of this manuscript.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbinf.2025.1606284/full#supplementary-material
SUPPLEMENTARY FIGURE S1 | Quality control and dimensionality reduction clustering analysis of the single-cell sequencing dataset GSE116256.
SUPPLEMENTARY FIGURE S2 | Gene difference analysis and functional enrichment of monocyte clusters (7 and 13) in the single-cell sequencing dataset GSE116256.
SUPPLEMENTARY FIGURE S3 | PH hypothesis testing for six candidate genes.
SUPPLEMENTARY TABLE S1 | List of CSRGs.
SUPPLEMENTARY TABLE S2 | The marker gene of the different cell types.
SUPPLEMENTARY TABLE S3 | Information on the control group and AML patients used for bone marrow sample collection at Suzhou Yongding Hospital.
SUPPLEMENTARY TABLE S4 | List of 517 differentially expressed genes between the control group and the AML group.
SUPPLEMENTARY TABLE S5 | 31 cross-candidate genes, the candidate genes.
References
Adès, L., Girshova, L., Doronin, V. A., Díez-Campelo, M., Valcárcel, D., Kambhampati, S., et al. (2022). Pevonedistat plus azacitidine vs azacitidine alone in higher-risk MDS/chronic myelomonocytic leukemia or low-blast-percentage AML. Blood Adv. 6 (17), 5132–5145. doi:10.1182/bloodadvances.2022007334
Aran, D., Looney, A. P., Liu, L., Wu, E., Fong, V., Hsu, A., et al. (2019). Reference-based analysis of lung single-cell sequencing reveals a transitional profibrotic macrophage. Nat. Immunol. 20 (2), 163–172. doi:10.1038/s41590-018-0276-y
Avelar, R. A., Ortega, J. G., Tacutu, R., Tyler, E. J., Bennett, D., Binetti, P., et al. (2020). A multidimensional systems biology analysis of cellular senescence in aging and disease. Genome Biol. 21 (1), 91. doi:10.1186/s13059-020-01990-9
Bailur, J. K., McCachren, S. S., Pendleton, K., Vasquez, J. C., Lim, H. S., Duffy, A., et al. (2020). Risk-associated alterations in marrow T cells in pediatric leukemia. JCI Insight 5 (16), e140179. doi:10.1172/jci.insight.140179
Bhave, R. R., Mesa, R., and Grunwald, M. R. (2023). Top advances of the year: myeloproliferative neoplasms. Cancer 129 (23), 3685–3691. doi:10.1002/cncr.35028
Charoentong, P., Finotello, F., Angelova, M., Mayer, C., Efremova, M., Rieder, D., et al. (2017). Pan-cancer immunogenomic analyses reveal genotype-immunophenotype relationships and predictors of response to checkpoint blockade. Cell Rep. 18 (1), 248–262. doi:10.1016/j.celrep.2016.12.019
Chatsirisupachai, K., Palmer, D., Ferreira, S., and de Magalhães, J. P. (2019). A human tissue-specific transcriptomic analysis reveals a complex relationship between aging, cancer, and cellular senescence. Aging Cell 18 (6), e13041. doi:10.1111/acel.13041
Chen, H., and Boutros, P. C. (2011). VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinforma. 12, 35. doi:10.1186/1471-2105-12-35
Chen, J., Hou, Q., Chang, T., Zheng, J., Yao, C., He, J., et al. (2024). Analysis of prognostic biomarker models of TXNIP/NLRP3/IL1B inflammasome pathway in patients with acute myeloid leukemia. Int. J. Med. Sci. 21 (8), 1438–1446. doi:10.7150/ijms.96627
Chen, M., Su, Z., and Xue, J. (2025). Targeting T-cell aging to remodel the aging immune system and revitalize geriatric immunotherapy. Aging Dis. doi:10.14336/ad.2025.0061
Collins, C. T., and Hess, J. L. (2016a). Deregulation of the HOXA9/MEIS1 axis in acute leukemia. Curr. Opin. Hematol. 23 (4), 354–361. doi:10.1097/moh.0000000000000245
Collins, C. T., and Hess, J. L. (2016b). Role of HOXA9 in leukemia: dysregulation, cofactors and essential targets. Oncogene 35 (9), 1090–1098. doi:10.1038/onc.2015.174
Csizmar, C. M., Saliba, A. N., Swisher, E. M., and Kaufmann, S. H. (2021). PARP inhibitors and myeloid neoplasms: a double-edged sword. Cancers (Basel) 13 (24), 6385. doi:10.3390/cancers13246385
Deshpande, A. J., and Zhu, N. (2023). Silencing with SAFB: a new role for HOXA9 in AML. Blood 141 (14), 1653–1655. doi:10.1182/blood.2022019244
Garnett, M. J., Edelman, E. J., Heidorn, S. J., Greenman, C. D., Dastur, A., Lau, K. W., et al. (2012). Systematic identification of genomic markers of drug sensitivity in cancer cells. Nature 483 (7391), 570–575. doi:10.1038/nature11005
Geeleher, P., Cox, N., and Huang, R. S. (2014). pRRophetic: an R package for prediction of clinical chemotherapeutic response from tumor gene expression levels. PLoS One 9 (9), e107468. doi:10.1371/journal.pone.0107468
Gu, Z., and Hübschmann, D. (2022). Make interactive complex heatmaps in R. Bioinformatics 38 (5), 1460–1462. doi:10.1093/bioinformatics/btab806
Gustavsson, E. K., Zhang, D., Reynolds, R. H., Garcia-Ruiz, S., and Ryten, M. (2022). ggtranscript: an R package for the visualization and interpretation of transcript isoforms using ggplot2. Bioinformatics 38 (15), 3844–3846. doi:10.1093/bioinformatics/btac409
Hänzelmann, S., Castelo, R., and Guinney, J. (2013). GSVA: gene set variation analysis for microarray and RNA-seq data. BMC Bioinforma. 14, 7. doi:10.1186/1471-2105-14-7
Heagerty, P. J., Lumley, T., and Pepe, M. S. (2000). Time-dependent ROC curves for censored survival data and a diagnostic marker. Biometrics 56 (2), 337–344. doi:10.1111/j.0006-341x.2000.00337.x
Ikeda, H., Kawase, K., Nishi, T., Watanabe, T., Takenaga, K., Inozume, T., et al. (2025). Immune evasion through mitochondrial transfer in the tumour microenvironment. Nature 638 (8049), 225–236. doi:10.1038/s41586-024-08439-0
Kontandreopoulou, C. N., Diamantopoulos, P. T., Tiblalexi, D., Giannakopoulou, N., and Viniou, N. A. (2021). PARP1 as a therapeutic target in acute myeloid leukemia and myelodysplastic syndrome. Blood Adv. 5 (22), 4794–4805. doi:10.1182/bloodadvances.2021004638
Kumar, V., Kumar, A., Mir, K. U. I., Yadav, V., and Chauhan, S. S. (2022). Pleiotropic role of PARP1: an overview. 3 Biotech. 12 (1), 3. doi:10.1007/s13205-021-03038-6
Lai, M., Liu, G., Li, R., Bai, H., Zhao, J., Xiao, P., et al. (2020). Hsa_circ_0079662 induces the resistance mechanism of the chemotherapy drug oxaliplatin through the TNF-α pathway in human colon cancer. J. Cell Mol. Med. 24 (9), 5021–5027. doi:10.1111/jcmm.15122
Lei, J., Qu, T., Cha, L., Tian, L., Qiu, F., Guo, W., et al. (2023). Clinicopathological characteristics of pheochromocytoma/paraganglioma and screening of prognostic markers. J. Surg. Oncol. 128 (4), 510–518. doi:10.1002/jso.27358
Liu, P., Xu, H., Shi, Y., Deng, L., and Chen, X. (2020). Potential molecular mechanisms of plantain in the treatment of gout and hyperuricemia based on network pharmacology. Evid. Based Complement. Altern. Med. 2020, 3023127. doi:10.1155/2020/3023127
Liu, X., Xie, P., Hao, N., Zhang, M., Liu, Y., Liu, P., et al. (2021). HIF-1-regulated expression of calreticulin promotes breast tumorigenesis and progression through Wnt/β-catenin pathway activation. Proc. Natl. Acad. Sci. U. S. A. 118 (44), e2109144118. doi:10.1073/pnas.2109144118
Luo, H., Li, Q., Hong, J., Huang, Z., Deng, W., Wei, K., et al. (2024). Targeting TNF/IL-17/MAPK pathway in hE2A-PBX1 leukemia: effects of OUL35, KJ-Pyr-9, and CID44216842. Haematologica 109 (7), 2092–2110. doi:10.3324/haematol.2023.283647
Mazziotta, F., Biavati, L., Rimando, J., Rutella, S., Borcherding, N., Parbhoo, S., et al. (2024). CD8+ T-cell differentiation and dysfunction inform treatment response in acute myeloid leukemia. Blood 144 (11), 1168–1182. doi:10.1182/blood.2023021680
Miao, P., Yu, J., Chen, Z., Qian, S., and Chen, C. (2024). Establishment and verification of a TME prognosis scoring model based on the acute myeloid leukemia single-cell transcriptome. Sci. Rep. 14 (1), 19811. doi:10.1038/s41598-024-65345-1
Murthy, G. S. G., Saliba, A. N., Szabo, A., Harrington, A., Abedin, S., Carlson, K., et al. (2024). A phase I study of pevonedistat, azacitidine, and venetoclax in patients with relapsed/refractory acute myeloid leukemia. Haematologica 109 (9), 2864–2872. doi:10.3324/haematol.2024.285014
Narayanan, D., and Weinberg, O. K. (2020). How I investigate acute myeloid leukemia. Int. J. Lab. Hematol. 42 (1), 3–15. doi:10.1111/ijlh.13135
Nebenfuehr, S., Kollmann, K., and Sexl, V. (2020). The role of CDK6 in cancer. Int. J. Cancer 147 (11), 2988–2995. doi:10.1002/ijc.33054
Paczulla, A. M., Rothfelder, K., Raffel, S., Konantz, M., Steinbacher, J., Wang, H., et al. (2019). Absence of NKG2D ligands defines leukaemia stem cells and mediates their immune evasion. Nature 572 (7768), 254–259. doi:10.1038/s41586-019-1410-1
Pollyea, D. A., Bixby, D., Perl, A., Bhatt, V. R., Altman, J. K., Appelbaum, F. R., et al. (2021). NCCN guidelines insights: acute myeloid leukemia, version 2.2021. J. Natl. Compr. Canc Netw. 19 (1), 16–27. doi:10.6004/jnccn.2021.0002
Ramsay, I. S., Ma, S., Fisher, M., Loewy, R. L., Ragland, J. D., Niendam, T., et al. (2018). Model selection and prediction of outcomes in recent onset schizophrenia patients who undergo cognitive training. Schizophr. Res. Cogn. 11, 1–5. doi:10.1016/j.scog.2017.10.001
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi:10.1093/nar/gkv007
Rossiello, F., Jurk, D., Passos, J. F., and d'Adda di Fagagna, F. (2022). Telomere dysfunction in ageing and age-related diseases. Nat. Cell Biol. 24 (2), 135–147. doi:10.1038/s41556-022-00842-x
Ru, Y., Kechris, K. J., Tabakoff, B., Hoffman, P., Radcliffe, R. A., Bowler, R., et al. (2014). The multiMiR R package and database: integration of microRNA-target interactions along with their disease and drug associations. Nucleic Acids Res. 42 (17), e133. doi:10.1093/nar/gku631
Sachs, M. C. (2017). plotROC: a tool for plotting ROC curves. J. Stat. Softw. 79, 2. doi:10.18637/jss.v079.c02
Salik, B., Yi, H., Hassan, N., Santiappillai, N., Vick, B., Connerty, P., et al. (2020). Targeting RSPO3-LGR4 signaling for leukemia stem cell eradication in acute myeloid leukemia. Cancer Cell 38 (2), 263–278.e6. doi:10.1016/j.ccell.2020.05.014
Schmoellerl, J., Barbosa, I. A. M., Eder, T., Brandstoetter, T., Schmidt, L., Maurer, B., et al. (2020). CDK6 is an essential direct target of NUP98 fusion proteins in acute myeloid leukemia. Blood 136 (4), 387–400. doi:10.1182/blood.2019003267
Shive, C. L., Freeman, M. L., Younes, S. A., Kowal, C. M., Canaday, D. H., Rodriguez, B., et al. (2021). Markers of T Cell exhaustion and senescence and their relationship to plasma TGF-β levels in treated HIV+ immune non-responders. Front. Immunol. 12, 638010. doi:10.3389/fimmu.2021.638010
Short, N. J., Muftuoglu, M., Ong, F., Nasr, L., Macaron, W., Montalban-Bravo, G., et al. (2023). A phase 1/2 study of azacitidine, venetoclax and pevonedistat in newly diagnosed secondary AML and in MDS or CMML after failure of hypomethy lating agents. J. Hematol. and Oncol. 16 (1), 73. doi:10.1186/s13045-023-01476-8
Song, Y., Yang, Z., Gao, N., and Zhang, B. (2024). MICAL1 promotes the proliferation in acute myeloid leukemia and is associated with clinical prognosis and immune infiltration. Discov. Oncol. 15 (1), 279. doi:10.1007/s12672-024-01150-6
Szabo, P. A., Levitin, H. M., Miron, M., Snyder, M. E., Senda, T., Yuan, J., et al. (2019). Single-cell transcriptomics of human T cells reveals tissue and activation signatures in health and disease. Nat. Commun. 10 (1), 4706. doi:10.1038/s41467-019-12464-3
Towers, R., Trombello, L., Fusenig, M., Tunger, A., Baumann, A. L., Savoldelli, R., et al. (2024). Bone marrow-derived mesenchymal stromal cells obstruct AML-targeting CD8(+) clonal effector and CAR T-cell function while promoting a senescence-associated phenotype. Cancer Immunol. Immunother. 73 (1), 8. doi:10.1007/s00262-023-03594-1
Walter, W., Sánchez-Cabo, F., and Ricote, M. (2015). GOplot: an R package for visually combining expression data with functional analysis. Bioinformatics 31 (17), 2912–2914. doi:10.1093/bioinformatics/btv300
Wang, L., Wang, D., Yang, L., Zeng, X., Zhang, Q., Liu, G., et al. (2022). Cuproptosis related genes associated with Jab1 shapes tumor microenvironment and pharmacological profile in nasopharyngeal carcinoma. Front. Immunol. 13, 989286. doi:10.3389/fimmu.2022.989286
Wei, X., Wang, W., Yin, Q., Li, H., Ahmed, A., Ullah, R., et al. (2024). In vivo chemical screening in zebrafish embryos identified FDA-approved drugs that induce differentiation of acute myeloid leukemia cells. Int. J. Mol. Sci. 25 (14), 7798. doi:10.3390/ijms25147798
Weinhäuser, I., Pereira-Martins, D. A., Almeida, L. Y., Hilberink, J. R., Silveira, D. R. A., Quek, L., et al. (2023). M2 macrophages drive leukemic transformation by imposing resistance to phagocytosis and improving mitochondrial metabolism. Sci. Adv. 9 (15), eadf8522. doi:10.1126/sciadv.adf8522
Wu, F., Xu, G., Li, G., Yin, Z., Shen, H., Ye, K., et al. (2023). A prognostic model based on prognosis-related ferroptosis genes for patients with acute myeloid leukemia. Front. Mol. Biosci. 10, 1281141. doi:10.3389/fmolb.2023.1281141
Xiao, G., Yu, L., Tan, W., Yang, H., Li, W., Xia, R., et al. (2023). Propofol inhibits glioma progression by regulating circMAPK4/miR-622/HOXA9 axis. Metab. Brain Dis. 38 (1), 233–244. doi:10.1007/s11011-022-01099-x
Yabas, M., and Hoyne, G. F. (2023). Immunological phenotyping of mice with a point mutation in Cdk4. Biomedicines 11 (10), 2847. doi:10.3390/biomedicines11102847
Yalniz, F. F., and Wierda, W. G. (2019). Targeting BCL2 in chronic lymphocytic leukemia and other hematologic malignancies. Drugs 79 (12), 1287–1304. doi:10.1007/s40265-019-01163-4
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16 (5), 284–287. doi:10.1089/omi.2011.0118
Yuan, J., Duan, F., Zhai, W., Song, C., Wang, L., Xia, W., et al. (2021). An aging-related gene signature-based model for risk stratification and prognosis prediction in breast cancer. Int. J. Womens Health 13, 1053–1064. doi:10.2147/ijwh.s334756
Zhang, H., Meltzer, P., and Davis, S. (2013). RCircos: an R package for Circos 2D track plots. BMC Bioinforma. 14, 244. doi:10.1186/1471-2105-14-244
Zhong, C., Wang, R., Hua, M., Zhang, C., Han, F., Xu, M., et al. (2021). NLRP3 inflammasome promotes the progression of acute myeloid leukemia via IL-1β pathway. Front. Immunol. 12, 661939. doi:10.3389/fimmu.2021.661939
Keywords: acute myeloid leukemia, T cell, cell senescence, single-cell RNA sequencing, prognostic risk model
Citation: Sha M, Chen J, Hou H, Dou H and Zhang Y (2025) Integrated single-cell and bulk RNA dequencing to identify and validate prognostic genes related to T Cell senescence in acute myeloid leukemia. Front. Bioinform. 5:1606284. doi: 10.3389/fbinf.2025.1606284
Received: 05 April 2025; Accepted: 09 June 2025;
Published: 25 June 2025.
Edited by:
Wenlin Yang, University of Florida, United StatesReviewed by:
Sheng-Yan Lin, Huazhong University of Science and Technology, ChinaEdgar Gonzalez-Kozlova, Icahn School of Medicine at Mount Sinai, United States
Copyright © 2025 Sha, Chen, Hou, Dou and Zhang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yan Zhang, MTgyNjE4NDkzNzZAMTYzLmNvbQ==