The samples source with single-cell RNA-sequencing (scRNA-seq, GSE116256) and expression profiles (GSE63270 and GSE12417) were downloaded from Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) of The National Center for Biotechnology Information (NCBI) (22). We enrolled three normal samples and ten patients with acute myeloid leukemia from GSE116256 for analysis of scRNA-seq (23–25). There were 104 normal and acute myeloid leukemia (42 populations and 62 leukemic populations) included in GSE63270 dataset (26). GSE12417 dataset contained the analysis of 79 samples of bone marrow or peripheral blood mononuclear cells from adult patients with untreated acute myeloid leukemia (27, 28). In addition, the expression profiles and clinical information were acquired from TCGA-LAML cohort, including 151 patients with acute myeloid leukemia (29, 30).

First, the data of single cell was filtered by setting each gene to be expressed in at least three cells, and each cell to express at least 500 genes, resulting in 9891 cells. We calculated the percentage of mitochondria and Ribosomal RNA (rRNA) through the PercentageFeatureSet function of Seurat package (31). The number of genes expressed in each single cell was greater than 100 and less than 5000, and we ensured the percentage of mitochondria was less than 20%. Furthermore, the Unique Molecular Identifier (UMI) of the single cell was at least greater than 100, resulting in 9886 cells. Subsequently, we used the method of log-normalization to standardize the single-cell data from each of the 13 samples. The highly variable features were identified by FindVariableFeatures function (32) based on variance stabilization transformation (VST). The genes were then scaled by using the ScaleData function for all genes. We utilized RunPCA function for PCA dimension reduction to find anchors. The FindNeighbors function with dim=15 and FindClusters function with Resolution=0.1 was used to luster cells. Ulteriorly, the RunTSNE function was used to conduct t-SNE (T-Distribution Stochastic Neighbour Embedding) dimension reduction analysis and the RunUMAP function was used to conduct UMAP (Uniform Manifold Approximation and Projection) reduction analysis. The marker genes for single cell were supplied by SingleR package (33) and the classical marker from the published literature (25).

Monocle R package was applied for the data of single cell to explore the correlation of aggrephagy-related genes and pseudotime trajectories (34). The graphs for the pseudotime trajectories of specific cell with aggrephagy-related genes were plotted using the function from Monocle R package, such as plot_pseudotime_heatmap and so on.

Based on the expression matrix of the scRNA-seq, dimension reduction analysis of aggrephagy-related genes in each cell clusters were conducted employing NMF (Nonnegative Matrix Factorization) R package (35, 36), thus displaying the effect of aggrephagy-related genes in single cell for acute myeloid leukemia.

FindAllMarkers function was applied to identify he marker genes of single cell for acute myeloid leukemia (31). The aggrephagy-mediated cell clusters were identified based on differentially expressed genes (DEGs) with log Fold Change (logFC) and aggrephagy-related genes. The NK cell subtypes were summarized from the published literature of Huan Liu et al (37).

SCENIC was a tool for simultaneously reconstructing gene regulatory networks and identifying stable cell states from single-cell RNA-seq data (38). The gene regulatory network was inferred based on co-expression and DNA motif analysis, and then network activity was analyzed in each cell to identify cell status (38). We carried out analysis of transcription factors (TFs) for aggrephagy-mediated cell clusters for acute myeloid leukemia through pySCENIC package (39–41). RcisTarget R package and two gene-motif rankings (hg19-tss-centered-10 kb and hg19-500 bp-upstream) was used to identify binding motifs of TFs in the gene list for acute myeloid leukemia (42, 43). The threshold value for the TFs was set as Benjamini–Hochberg false discovery rate (BH-FDR) <0.05.

The signaling inputs and outputs among the cell types and aggrephagy-mediated cell clusters were assessed by applying CellChat package (44) and CellChatDB database (45). The netVisual_circle function was utilized for evaluating the strength of cell–cell communication networks among cell subsets (44, 45). In addition, the ligand-receptor interactions among the specific cell subsets were estimated via the netVisual_bubble function (44, 45).

Based on the data of scRNA, Gene Set Variation Analysis (GSVA) was applied to compute the signature scores involved in aggrephagy for public database (46). We carried out Cox proportional hazard regression to evaluate the prognosis for the aggrephagy-related cell clusters (47). The Kaplan–Meier curves was plotted through the survminer R package.

We used TIDE (Tumor Immune Dysfunction and Exclusion) algorithm to analyze the immune checkpoint blockade immunotherapeutic for the aggrephagy-related cell clusters (48). We also reviewed the published literature to validate the prognostic and therapeutic effects of each cell subtype using real-world immunotherapy cohorts (49–60).

The continuous or category variables were compared using Student’s t-test, Wilcoxon rank sum test, Kruskal–Wallis’s test, or Chi-square test. The log-rank test was used for survival analyses.

We carried out dimensionality reduction and annotation of single cell for acute myeloid leukemia as described in the materials and methods section. We ensured that the number of genes expressed in each single cell was greater than 100 and less than 5000, the percentage of mitochondria was less than 20%, and the Unique Molecular Identifier (UMI) of the single cell was at least greater than 100, resulting in 9886 cells. Supplementary Figures S1A–B was the statistical diagram of cell filtration, which could be seen to meet all thresholds set above (Supplementary Figures S1A, B). The highly variable features were identified via FindVariableFeatures function based on VST, and the top ten highly variable genes among the single cell were marked out in the volcano plot, including IGLL5, HBB, JCHAIN, HBG1, HBG2, HBD, HBA2, CLC, CA1, and HBA1 (Supplementary Figure S1C). Ulteriorly, PCA analysis was carried out on the highly variable genes, we used the Elbow algorithm to carry out the Standard Deviation based on the highly variable genes (Supplementary Figure S1D). The RunTSNE function was used to conduct t-SNE dimension reduction analysis and the RunUMAP function was used to conduct UMAP reduction analysis, thus identifying a total of 18 cell subsets (Supplementary Figure S1E). Afterwards, the marker genes were used to annotate the specific cell types, thus identifying 11 specific cell types, including T cell, Natural Killer (NK) cell, Common Myeloid Progenitor (CMP) cell, Myeloid cell, Granulocyte Monocyte Progenitor (GMP) cell, Megakaryocyte Erythroid Progenitor (MEP) MEP, Promonocyte (Promono) cell, Plasma cell, Hematopoietic Stem Cell (HSC) cell, B cell, and Erythroid cell (Figure 1A). We plotted the correlation network for the number of interactions among the 11 specific cell types (Figure 1B). Figure 1C visually showed the proportion of different cell types in each sample. Finally, we created the heat map to show the expression of the aggrephagy-related genes in different cell types (Figure 1D). We could see that there were several aggrephagy-related genes that were highly expressed in the 11 specific cell types, including HSP90AA1, RPS27A, UBA52, UBB, UBC, and VIM (Figure 1D). We displayed the global view of the expression pattern for marker genes gained as described in the methods section, reflecting the dynamic features of each cell subsets (Supplementary Figure S2).

There was a total of 2397 cells in the T cell type. Using UMAP reduction analysis, the 2397 cells in the T cell type could be clustered into eight cell clusters (Supplementary Figure S3), and the global view of the expression pattern for marker genes of the eight cell clusters was displayed in Supplementary Figure S3. Further, the eight cell clusters could be re-clustered into nine cell subsets (Supplementary Figure S3). Based on the aggrephagy-related genes, six clusters were identified using Monocle analysis (Figure 2A), including Cluster 1 (DYNC1LI2 and TUBB1), Cluster 2 (PRKN, TUBB2A, DYNLL2, RPS27A and UBA52), Cluster 3 (TUBA1A, UBE2V1, IFT88 and VCP), Cluster 4 (TUBA1C, TUBA3C, TUBA3D, HSP90AA1, VIM, DYNC1I2, PARK7, UBE2N, DYNLL1 and HSF1), Cluster 5 (TUBA4A, TUBB4B, TUBA1B, ARL13B and PCNT), and Cluster 6 (DYNC1LI1, UBB, HDAC6, DYNC1H1 and UBC). From the heatmap generated by Pseudotime analysis, the critical role of the aggrephagy-related genes in the trajectory process of T cells was observed (Figure 2A). Four subgroups of T cells were obtained by re-clustering annotation using t-SNE dimension reduction analysis, including CD8⁺ T cell, CD4⁺ T cell, natural killer (NK) cells, and Regulatory T (Treg) cells (Figure 2B). Among the four cell subgroups, we found that CD8⁺ T cell and CD4⁺ T cell had a higher percentage both in tumor samples and normal samples than the other two cell subgroups (Figure 2C), and CD8⁺ T cell occupied the highest proportion among the four cell subgroups (Figure 2C). The result of CD8⁺ T cell revealed that the whole trajectory could be divided into three segments (State 1, State 2, and State 3) on the basis of the developmental order (Figure 2D). Ulteriorly, NMF clustering analysis of the aggrephagy-related gene set for the trajectories showed that these cells aggregated into nine clusters (Figure 2E). In addition, the results of UMAP reduction analysis indicated that the NMF cell types were clustered into six aggrephagy-related CD8⁺ T clusters, including TUBA1B+CD8⁺ T−C1, DYNC1H1+CD8⁺ T−C2, UBE2V1+CD8⁺ T−C3, UBE2N+CD8⁺ T−C4, Unc−CD8⁺ T−C5, and Non−Aggre−CD8⁺ T−C6 (Figure 2F). The number of ligand-receptor links among the six aggrephagy-related CD8+ T clusters was computed by Cell-Chat analysis (Figure 2G). The weights and strength of ligand-receptor links among the six aggrephagy-related CD8⁺ T clusters was computed by Cell-Chat analysis (Supplementary Figure S3). Lastly, the discrepancies in the exhausted CD8⁺ T (CD8⁺_exhau), cytotoxic CD8⁺ T (CD8⁺_cyoto), and TFs (BTN3A1, BTN3A2, BTN2A2, LGALS9, TIGIT, CD274, BTLA, CTLA4, IL10, LAIR1, CD247, TGFB1, SLAMF7, CD160, CD244, HAVCR2, LAG3, CD96, ADORA2A, PDCD1, and CD48) among the six aggrephagy-related CD8⁺ T clusters were visually displayed in the pathway heatmap (Figure 2H). TUBA1B+CD8⁺ T−C1 tended to be exhausted CD8⁺ T (CD8⁺_exhau), while UBE2N+CD8⁺ T−C4 tended to be cytotoxic CD8⁺ T (CD8⁺_cyoto) as shown in Figure 2H. It is noteworthy that TFs of LGALS9, TIGIT, BTLA, and CTLA4 were upregulated in the TUBA1B+CD8⁺ T−C1, TFs of IL10 and CD160 were upregulated in the DYNC1H1+CD8⁺ T−C2, ADORA2A was upregulated in the UBE2V1+CD8⁺ T−C3, BTN3A1, BTN3A2, CD274, CD247, SLAMF7, LAG3, and PDCD1 were upregulated in the UBE2N+CD8⁺ T−C4, BTN2A2 and LAIR1 were upregulated in the Non− Aggre−CD8⁺ T−C6 (Figure 2H).

There was a total of 1067 cells in the NK cell type. Using UMAP reduction analysis, the 1067 cells in the NK cell type could be clustered into eleven cell clusters (Supplementary Figure S4). Based on the aggrephagy-related genes, six clusters were identified using Monocle analysis (Figure 3A), including Cluster 1 (TUBA1B, TUBA4A, UBA52, UBB, TUBB4B, and UBE2V1), Cluster 2 (TUBA1C, HDAC6, VCP, DYNC1LI1, UBC, PRKN, TUBB1, HSP90AA1, RPS27A, PCNT, DYNC1H1, and DYNLL2), Cluster 3 (TUBA1A, VIM, HSF1, and UBE2N), Cluster 4 (TUBB6, ARL13B, DYNC1I2, DYNC1LI2, and TUBB2A), Cluster 5 (IFT88 and TUBA8), Cluster 6 (TUBA3C, DYNLL1, and PARK7). The results of UMAP reduction analysis indicated that the NMF cell types were clustered into six aggrephagy-related NK clusters, including UBE2N+NK−C1, UBE2V1+NK−C2, DYNC1H1+NK−C3, PARK7+NK−C4, Unc+NK−C5, and Non−Aggre−NK−C6 (Figure 3B). The number of ligand-receptor links among the six aggrephagy-related NK clusters was computed by Cell-Chat analysis (Figure 3C). The number, weights and strength of ligand-receptor links among the aggrephagy-related NK clusters was computed by Cell-Chat analysis (Supplementary Figure S4). Lastly, the discrepancies in the NK−CD56bright, NK−CD56dim, NK−HIA, LrNK−FCGR3A, LrNK−XCL1, KIR2DS1, NCR1, NCR2, NCR3, TLR3, TLR9, KIR3DL1, KIR2DL3, KLRB1, LILRB1, LILRB2, KLRG1, CEACAM1, CD244, LAIR1, CD96, TIGIT, and LAG3 among the six aggrephagy-related NK clusters were visually displayed in the pathway heatmap (Figures 3D, E). It is noteworthy that NK−CD56bright was upregulated in PARK7+NK−C4 (Figure 3D), NK−HIA was upregulated in UBE2N+NK−C1 (Figure 3D), LrNK−FCGR3A was upregulated in UBE2N+NK−C1 (Figure 3D), LrNK−XCL1 was upregulated in PARK7+NK−C4 (Figure 3D), NCR1 was upregulated in UBE2N+NK−C1 (Figure 3E), NCR3 was upregulated in PARK7+NK−C4 (Figure 3E), TLR9 was upregulated in PARK7+NK−C4 (Figure 3E), KIR3DL1 was upregulated in UBE2V1+NK−C2 (Figure 3E), LILRB1 was upregulated in UBE2V1+NK−C2 (Figure 3E), LILRB2 was upregulated in DYNC1H1+NK−C3 (Figure 3E), KLRG1 was upregulated in Non−Aggre−NK−C6 (Figure 3E), CEACAM1 was upregulated in PARK7+NK−C4 (Figure 3E), CD96 was upregulated in DYNC1H1+NK−C3 (Figure 3E), TIGIT was upregulated in PARK7+NK−C4 (Figure 3E), LAG3 was upregulated in PARK7+NK−C4 (Figure 3E).

There was a total of 3167 cells in the Myeloid cell type. Based on Myeloid cell type, the PCA analysis was carried out on the highly variable genes, the Elbow algorithm to carry out the Standard Deviation based on the highly variable genes (Supplementary Figure S5). Using UMAP reduction analysis, the 3167 cells in the Myeloid cell type could be clustered into eleven and twelve cell clusters (Supplementary Figure S5). Based on the aggrephagy-related genes, six clusters were identified using Monocle analysis (Figure 4A), including Cluster 1 (VCP, HDAC6, TUBA1A, DYNC1H1, and TUBB1), Cluster 2 (UBC, TUBA4A, DYNC1LI1, VIM, TUBB2A, UBB, HSP90AA1, and TUBB6), Cluster 3 (TUBA3C and TUBA3E), Cluster 4 (PCNT, TUBA1B, UBE2V1, PARK7, TUBB4B, TUBAL3, RPS27A, UBA52, and UBE2N), Cluster 5 (TUBA1C, ARL13B, TUBA8, DYNLL2, DYNC1LI2, and HSF1), Cluster 6 (TUBA4B, IFT88, DYNC1I2, and DYNLL1). Three subgroups of Myeloid cells were obtained by re-clustering annotation using UMAP reduction analysis, including Mono (monocytes) cell, Macrophages (MAC) cell, and Dendritic cell (DC) cell (Figure 4B). Further, we displayed the global view of the expression pattern for marker genes of Mono (monocytes) cell and Macrophages (MAC) cell in Supplementary Figure S5. In addition, the results of UMAP reduction analysis indicated that the NMF cell types were clustered into six aggrephagy-related Mac clusters, including DYNLL1+Mac-C1, UBE2V1+Mac-C2, TUBA1A+Mac-C3, PAPK7+Mac-C4, Unc-Aggre-Mac-C5, and Non-Aggre-Mac-C6 (Figure 4C). The number of ligand-receptor links among the aggrephagy-related Mac clusters was computed by Cell-Chat analysis (Supplementary Figure S5). Ulteriorly, we used scMetabolism package to assess the correlation between the aggrephagy-related Mac clusters and metabolic pathways, and we could intuitively see the differences in metabolic pathways of each aggrephagy-related Mac cluster from the bubble map (Figure 4D). To identify M1/M2 type cells, we scored related genes, suggesting that M1 type macrophages were more active in AML (Figures 4E–H).

There were 104 normal and acute myeloid leukemia (42 populations and 62 leukemic populations) samples included in GSE63270 dataset, we compared the abundance of UBE2N+NK-C1, UBE2V1+NK-C2, DYNC1H1+NK-C3, PARK7+NK-C4, TUBA1B+CD8⁺ T−C1, DYNC1H1+CD8⁺ T−C2, UBE2V1+CD8⁺ T−C3, UBE2N+CD8⁺ T−C4, DYNLL1+Mac-C1, UBE2V1+Mac-C2, TUBA1A+Mac-C3, and PAPK7+Mac-C4 between normal and AML samples (Figure 5A). The results indicated that the higher abundance of PARK7+NK-C4, DYNC1H1+CD8⁺ T−C2, DYNLL1+Mac-C1, and TUBA1A+Mac-C3 were observed in the AML samples, while the higher abundance of TUBA1B+CD8⁺ T−C1, UBE2V1+CD8⁺ T−C3, and UBE2V1+Mac-C2 were observed in the normal samples (Figure 5A). Based on the differentially expressed genes (DEGs) generated by the DYNC1H1+CD8+ T−C2, TUBA1A+Mac-C3, and UBE2V1+CD8+ T−C3, the prognostic models were established using TCGA-LAML cohort, the poor prognosis was observed in patients with higher level of DYNC1H1+CD8⁺ T−C2, lower level of TUBA1A+Mac-C3, and higher level of UBE2V1+CD8⁺ T−C3 (Figure 5B). GSVA was used for calculating the aggrephagy-related score, the prognosis of the AML patients in the GSE12417 and TCGA-LAML cohorts were further evaluated as displayed in Figure 5C. We found that the survival rates of AML patients in GSE12417 and TCGA-LAML cohorts were significantly different among DYNC1H1+CD8⁺ T−C2, DYNC1H1+NK-C3, DYNLL1+Mac-C1, PAPK7+Mac-C4, PARK7+NK-C4, TUBA1A+Mac-C3, TUBA1B+CD8⁺ T−C1, UBE2N+CD8⁺ T−C4, UBE2N+NK-C1, UBE2V1+CD8⁺ T−C3, UBE2V1+Mac-C2, and UBE2V1+NK-C2 (Figure 5C).

We compared the response status (False or True) of immune checkpoint blockade therapy for patients with AML among the aggrephagy-related cell clusters (UBE2N+NK-C1, UBE2V1+NK-C2, DYNC1H1+NK-C3, PARK7+NK-C4, TUBA1B+CD8⁺ T−C1, DYNC1H1+CD8⁺ T−C2, UBE2V1+CD8⁺ T−C3, UBE2N+CD8⁺ T−C4, DYNLL1+Mac-C1, UBE2V1+Mac-C2, TUBA1A+Mac-C3, and PAPK7+Mac-C4) using TIDE algorithm (Figure 6A and Supplementary Figure S6). For the AML patients with True response status, the abundance of UBE2N+NK-C1, PARK7+NK-C4, and TUBA1A+Mac-C3 was higher, while the abundance of DYNC1H1+CD8⁺ T−C2 was lower (Figure 6A). We found that the OR rates of AML patients in GSE12417 and TCGA-LAML cohorts were significantly different among DYNC1H1+CD8⁺ T−C2, DYNC1H1+NK-C3, DYNLL1+Mac-C1, PAPK7+Mac-C4, PARK7+NK-C4, TUBA1A+Mac-C3, TUBA1B+CD8⁺ T−C1, UBE2N+CD8⁺ T−C4, UBE2N+NK-C1, UBE2V1+CD8⁺ T−C3, UBE2V1+Mac-C2, and UBE2V1+NK-C2 (Figure 6B).

In addition, we compared the response status (SD/PD or CR/PR) of immunotherapy for patients with AML among the aggrephagy-related cell clusters (UBE2N+NK-C1, UBE2V1+NK-C2, DYNC1H1+NK-C3, PARK7+NK-C4, TUBA1B+CD8⁺ T−C1, DYNC1H1+CD8⁺ T−C2, UBE2V1+CD8⁺ T−C3, UBE2N+CD8⁺ T−C4, DYNLL1+Mac-C1, UBE2V1+Mac-C2, TUBA1A+Mac-C3, and PAPK7+Mac-C4) based on public dataset (Figure 7A). For the AML patients with CR/PR response status, the abundance of TUBA1B+CD8⁺ T−C1 and DYNLL1+Mac-C1 was higher, while the abundance of TUBA1A+Mac-C3 was lower (Figure 7A). We also found that AML patients with low abundance of TUBA1A+Mac-C3 may have the better prognosis (Figure 7B). In addition, we observed that the expression of TUBA1A was upregulated in bone marrow cells of AML patient both in mRNA (Figure 7C) and protein (Figure 7D) levels.