Gene expression profiles of 151 bone marrow samples from AML patients were downloaded from the TCGA database. Since sequencing data of control samples were not available from the TCGA database, we obtained RNA sequencing data of 337 normal peripheral blood samples from the GTEx database. The GSE37642 (Mitra et al., 2008; Liang et al., 2016; Meng et al., 2016; Nibourel et al., 2017) was obtained from the GEO database. The GSE37642 dataset (platform: GPL570; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE37642) contains transcriptional data from 136 AML patients with complete survival information, which was used for external validation analysis of the prognostic signature. Furthermore, the GSE12417 (platform: GPL570; n = 73; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE12417) and GSE71014 (platform: GPL10558; n = 104; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE71014) datasets were downloaded as complementary external validation sets to perform validation analyses on a larger sample base.

Figure 1 shows the schematic presentation of this study.

We downloaded 259 FRGs from the FerrDb online database (http://www.zhounan.org/ferrdb/). Meanwhile, a total of 60 FRGs were also obtained from the report of Yingkai Hong et al. (Herold et al., 2014). After de-duplication (retention of unique values), a total of 268 FRGs were obtained for our study (Supplementary Table S1).

Gene expression profiles of AML (n = 151) and normal (n = 337) samples from TCGA and GTEx were subjected to normalize Between Arrays function for normalization, and subsequently, differential expression analysis (AML vs. normal) was performed using the R package limma. The significance threshold was set to |log₂ fold change (FC)| ≥ 1 and p < 0.05. A volcano map showing the distribution of the identified DEGs was plotted based on the R package ggplot2. In addition, we also verified that the identified DEGs were able to distinguish between normal and AML samples by Principal Component Analysis (PCA) to imply the applicability of the samples.

To reveal the functions of target genes, R package clusterProfiler (Papaemmanuil et al., 2016) was used to conduct Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses. The GO terms were comprised of the following three divisions: biological process (BP), cellular component (CC), and molecular function (MF). p < 0.05 was regarded as statistically significant.

The CNV data of AML patients (n = 200) and normal subjects (n = 197) were downloaded from the TCGA database (sequencing platform Affymetrix SNP 6.0). Genes in CNV regions were annotated using Genome Research Consortium Human build 38 (GRCh38) as the reference genome. CNVs alteration rates between normal and tumor samples were then compared using the Chi-square test, and CNVs data with p < 0.05 were chosen for the next analysis. Then the CNVs data and DEGs data of the same sample were merged to construct a matrix. By using the Kolmogorov-Smirnov test, those genes showing the same tendency both in CNVs and differential gene expression (CNV increase-upregulated; CNV decrease-downregulated) were selected as CNV-driven DEGs. Moreover, CNV-driven DEGs and identified FRGs were analyzed for overlap, and the common genes in both gene lists were defined as CNV-driven DE-FRGs.

The TCGA-AML dataset containing 132 AML samples with complete survival information was used for prognostic gene screening and prognostic model construction and evaluation. The GSE37642 (n = 136), GSE12417 (n = 73), and GSE71014 (n = 104) datasets were used as an independent external validation set for prognostic model validation analysis. Prognostic genes were screened by Cox regression analysis. Briefly, target genes were included in univariate Cox regression analysis, and variables satisfying p < 0.1 were included in stepwise regression multivariate Cox analysis. The variables obtained from multivariate Cox regression analysis were identified as the optimal variables for the construction of the prognostic model. The risk score for each AML patient was calculated using the regression coefficient (coef) calculated from the multivariate Cox analysis and the expression of prognostic genes. The formula for calculating the risk score as shown below: r i s k   s c o r e = c o e f g e n e 1 × e x p r e s s i o n g e n e 1 + c o e f g e n e 2 × e x p r e s s i o n g e n e 2 + ⋯ + c o e f g e n e n × e x p r e s s i o n g e n e n

The samples were classified into high- and low-risk groups based on the median value of the risk score in each dataset (TCGA dataset and independent external validation set). The difference in OS between the two risk subgroups was assessed based on the R package survival using Kaplan-Meier (K-M) analysis. Receiver operating characteristic (ROC) curves were plotted by survROC package to assess the accuracy of the risk score for prognostic prediction in patients with TCGA-AML and GSE37642/GSE12417/GSE71014 datasets-AML.

Clinical characteristics of AML (age and sex) available in the TCGA database were included in the Cox regression analysis along with the risk score. Univariate Cox regression analyses with p < 0.05 output would be performed further in multivariate Cox analyses. Ultimately, variables with p < 0.05 generated by multivariate Cox regression analysis were considered as independent prognostic factors for AML.

We selected the blood of 10 AML patients and 10 healthy people to carry on the quantitative PCR test to the genes screened in this study. The experimental verification data are obtained from the sample database of the Institute of Pediatrics of Tianjin Children’s Hospital, and no samples are obtained directly from the children’s body, so they do not need to be approved by the institutional ethics committee. All experimental operations are carried out in accordance with the relevant guidelines and regulations, and have been repeatedly verified by professional laboratory researchers.

Whole-cell RNA was extracted via the RNAiso Plus (TaKaRa, Japan). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed to detect DNAJB6 and HSPB1 expression using American Bio rad Bole T100 gradient PCR instrument. The following primers were used for qRT-PCR: DNAJB6 primers: upstream primer F: 5′-TAT​GAA​GTG​CTG​TCG​GAT​GCT​AAG​A-3'; downstream primer R: 5′-GAA​GAC​ATC​ATC​TGG​GTT​ACG​GA-3'. The conditions for PCR to DNAJB6 were: the size of the amplified product was 144bp, annealing temperature was 55–60°C. The sequences of HSPB1 primers used were: upstream primer F: 5′-GCA​GTC​CAA​CGA​GAT​CAC​CA-3′ and downstream primer R: 5′-TTA​CTT​GGC​GGC​AGT​CTC​ATC-3′. The conditions for PCR to HSPB1 were: the size of the amplified product was 97bp, annealing temperature was 55–60°C. qRT-PCR was conducted in triplicate for each sample. All gene expression levels were normalized to that of GAPDH using the 2^−ΔΔCt method. Uupaired t-test (two-tailed) was used for the comparison analyses.

All analyses in this study were performed in R software. A log-rank test was used to check the significant difference in OS between groups. An area under the ROC curve (AUC) served as an indicator of prognostic accuracy. Unless otherwise specified, a p-value less than 0.05 was considered statistically significant.

The normalized expression profiles of TCGA-AML (n = 151) and GTEx-normal (n = 337) were selected as the basis for differential expression analysis. By R package limma, we identified 6,828 DEGs between AML and normal samples. Among them, a total of 3,330 met log₂ FC ≥ 1 and p < 0.05 and 3,498 matched log₂ FC ≤ -1 and p < 0.05 (Figure 2A; Supplementary Table S2). Furthermore, PCA performed based on the obtained DEGs demonstrated that samples from different groups were clustered in the same category (Figure 2B).

To explore the potential mechanisms of the above AML-related DEGs, we performed GO and KEGG function enrichment analysis for upregulated DEGs and downregulated DEGs, respectively, using the clusterProfiler package. Figures 3A–C illustrated the top5 terms that were significantly enriched in the three categories of GO, BP, CC, and MF, based on upregulated and downregulated DEGs. In the BP category, specifically, upregulated DEGs were significantly associated with “ribonucleoprotein complex biogenesis”, “RNA splicing”, “RNA splicing, via transesterification reactions” (Figure 3A); besides, these genes were found to be closely correlated with the cell cycle (“DNA replication”, “chromosome segregation”, “mitotic nuclear division”, etc.). Moreover, upregulated DEGs might be played mainly in CCs such as “chromosomal region”, “spindle”, “chromosome, centromeric region” ( Figure 3B) for MFs such as “histone binding”, “DNA-dependent ATPase activity”, “helicase activity” (Figure 3C). The results of the detailed GO analysis for the upregulated DEGs could be reviewed in Supplementary Table S3. For the down-regulated DEGs, they were significantly correlated with immune responses (“regulation of immune effector process”, “negative regulation of immune system process”, “immune response-activating signal transduction”, etc.) and biological processes of immune cells (“neutrophil degranulation”, “T cell activation”, “T-cell differentiation”, “lymphocyte proliferation”, etc.) (Figure 3A). Also, these genes were able to perform the molecular functions of “carbohydrate binding”, “organic acid binding”, and “MHC protein binding” (Figure 3C) in “specific granule”, “secretory granule lumen”, and “cytoplasmic vesicle lumen” (Figure 3B). Detailed GO analysis results for downregulated DEGs were displayed in Supplementary Table S4. The top5 pathways significantly enriched in KEGG analysis were shown in Figure 3D, with upregulated DEGs involved in “Cell cycle”, “Spliceosome”, “Nucleocytoplasmic transport”, “Protein processing in endoplasmic reticulum”, and “DNA replication” (Supplementary Table S5); downregulated DEGs were closely associated with “Osteoclast differentiation”, “Neutrophil extracellular trap formation”, “Phagosome”, “Shigellosis”, and “Malaria” (Supplementary Table S6).

By applying the Chi-square test, 4637 CNV genes associated with AML were identified (p < 0.05; Supplementary Table S7). The distribution of AML-related CNVs in chromosomes was presented in Figure 4A. Then CNV-driven DEGs were screened using the Kolmogorov-Smirnov test. We selected 337 CNV-driven DEGs, of which 127 were upregulated in AML with increased CNV (Supplementary Table S8) and the remaining 210 were downregulated with decreased CNV (Supplementary Table S9). Subsequently, by overlap analysis (Figure 4B), a total of 6 CNV-driven DE-FRGs, namely ALOX15B, MTDH, DNAJB6, HSPB1, ATF4, and PLIN2, were identified in the above list of CNV-driven DEGs and the list of 268 FRGs (Supplementary Table S1).

To illustrate the potential functional and biological effects of these CNV-driven DE-FRGs, GO (Supplementary Table S10) and KEGG (Supplementary Table S11) analyses were performed (Figure 4C). Results showed that CNV-driven DE-FRGs were significantly enriched in the BP category in terms related to oxidative stress response and its mediated cell death, such as “negative regulation of cellular response to oxidative stress”, “negative regulation of oxidative stress-induced cell death”, “negative regulation of response to oxidative stress”, and “regulation of oxidative stress-induced cell death”. KEGG analysis showed that these genes were involved in “MAPK signaling pathway”, “Cocaine addiction”, “VEGF signaling pathway”, “Arachidonic acid metabolism”, and “Cortisol synthesis and secretion”. These results suggested that CNV-driven DE-FRGs were probably implicated in the cell cycle dysregulation and the inflammatory immune response during the disease process.

We matched the expression profiles of the identified CNV-driven DE-FRGs in 132 TCGA-AML samples containing complete survival data. Univariate Cox regression analysis was applied to identify the CNV-driven DE-FRGs associated with survival time in AML patients. p < 0.1 was set as the cut-off value, and a total of 2 variables associated with survival in TCGA-AML patients were screened, namely DNAJB6 and HSPB1 (Figure 5A). HR > 1 for HSPB1 (HR = 1.2222, 95% CI. 0.992–1.506, p = 0.059) may be an oncogenic gene in AML, whereas DNAJB6 (HR = 0.514, 95% CI: 0.345–0.767, p = 0.001) with HR < 1 was expected to be a protective factor for AML. Further, the Cox model consisting of DNAJB6 and HSPB1 (Figure 5B) was identified as the optimal prognostic signature for AML by sophisticated calculations of multivariate Cox analysis with stepwise regression.

We assessed the efficacy of prognostic signature consisting of DNAJB6 and HSPB1 for predicting AML prognosis by risk scoring system. The risk score for each TCGA-AML patient was calculated according to the following formula: risk score = (−0.65 * expression of DNAJB6) + (0.19 * expression of HSPB1). All TCGA-AML samples were divided into high- (n = 66) and low- (n = 66) risk groups according to the median risk score (median = 1.029) (Supplementary Table S12). Figure 5C demonstrated the risk profile and AML survival distribution in the TCGA database, suggesting that low-risk AML patients had relatively longer OS than high-risk patients. Moreover, the heatmap showed that DNAJB6 was relatively less expressed in the high-risk group compared to the low-risk group, while HSPB1 tended to be more highly expressed in the high-risk group (Figure 5C). K-M survival analysis confirmed that the low-risk group had better OS and the high-risk score was linked to poor outcomes (Figure 5D). The ROC curve assessed the accuracy of the risk score to predict OS in TCGA-AML patients at 1, 3, and 5 years 1 year OS had an AUC of 0.726, 3 years of 0.634, and 5 years of 0.741 (Figure 5E). These results indicated that our risk score had a more reliable performance in predicting the prognosis of AML patients.

The risk scores for AML patients in the GSE37642 (n = 136), GSE12417 (n = 73), and GSE71014 (n = 104) datasets were calculated using the above equations, and the samples were divided into high-risk and low-risk groups based on the optimal threshold for risk scores (Figures 6A–C). The predictive validity of the risk scores in the independent external validation set was consistent with that in TCGA database. Patients in the low-risk group had significantly longer survival times compared with the high-risk group (Figures 6D–F). Meanwhile, ROC curve analysis showed that the AUC of risk score in predicting patients” OS at 1, 3, and 5 years was 0.631, 0.682, and 0.675 in the GSE37642 dataset (Figure 6G). 0.640, 0.621, and 0.589, respectively, in the GSE12417 dataset (Figure 6H). 0.652, 0.782, and 0.766, respectively, in the GSE71014 dataset (Figure 6I). This evidence suggested that our risk score had more satisfactory predictive validity and some general applicability.

We explored the relationship between risk score and clinical characteristics (age and gender) in the TCGA-AML dataset by the Wilcoxon test. The results revealed that the risk score was significantly different between age subgroups (p = 0.02), and the risk score was positively correlated with age, with higher risk levels at older ages (Supplementary Figure 1A). Besides, the risk scores were higher in the male subgroup than in the female group, but which was not statistically significant (p = 0.084; Supplementary Figure 1B).

Further, Cox regression analysis was utilized to evaluate whether the risk score could predict the outcome of AML patients independently of clinical characteristics (age and gender). Univariate Cox analysis noted that age and risk score were significantly associated with prognosis in AML patients (p < 0.05; Table 1). Ultimately, multivariate Cox analysis indicated that age and risk score were the independent factors of AML prognosis (Table 2).

We collected the blood of 10 AML patients and 10 healthy individuals and analyzed the mRNA expression levels of two prognostic genes by PCR. The results showed that the expression of the DNAJB6 gene and the HSPB1 gene in AML patients were both lower than that in healthy people (p < 0.0001; Figures 7A,P < 0.0001; Figure 7B). There was a significant statistical difference in the expression level after statistical analysis.