Full Text Article

Introduction

Acute myeloid leukemia (AML) is a malignancy of the stem cell precursors of the myeloid lineage (red blood cells, platelets, and white blood cells other than B and T cells) [1], which is characterized by the rapid growth of abnormal myeloid cells, which are immature white blood cells [2,3]. AML exhibits genetic heterogeneity, implying that there are variations in the genetic characteristics of the disease4. Acute lymphoblastic leukaemia (ALL), a malignant disorder of lymphoid progenitor cells, affects both children and adults [4]. It is a type of cancer that affects the white blood cells, particularly the lymphocytes, which are a type of immune cell [5]. It is characterized by an uncontrolled growth of lymphoid cells that are arrested at an early stage of differentiation [6,7]. It has the ability to infiltrate the bone marrow, bloodstream, and other extramedullary sites [8]. In the United States, the estimated incidence of ALL in 2014 was approximately 1.57 cases per 100,000 individuals, with approximately 5960 newly diagnosed cases and 1470 deaths reported in 2018 [9-11].

Despite substantial advancements in diagnostic methods and treatment modalities, these leukemias represent substantial challenges for the area of oncology [12]. AML and ALL are associated with complex molecular and genetic alterations that contribute to their pathogenesis and clinical heterogeneity [13]. Understanding the underlying mechanisms and molecular characteristics of these diseases is crucial for improving patient outcomes and developing targeted therapeutic approaches [14]. Transcriptomics [15], metabolomics [16], and epigenomics [17], among other high-throughput technologies, have revolutionized cancer research in recent years by offering full molecular profiles of tumors. These omics techniques have allowed for the discovery of critical molecular changes and processes involved in leukemogenesis, allowing for the creation of new diagnostic tools and treatment options.

The integration of multiple omics data sets allowed to comprehensively characterize the molecular landscape of AML and establish a predictive model for AML classification. By leveraging machine learning algorithms and utilizing the combined information from transcriptomics, metabolomics, and DNA methylation data, we aimed to compare the differences between AML and ALL, and develop a robust and accurate model for the prediction of AML.

Methods

Data Collection

Data for ALL were collected from the GSE79533 dataset obtained from the Gene Expression Omnibus (GEO) database. Data for AML were collected from two datasets, namely GSE26294 and GSE110087, and obtained from the GEO database. Leukemia methylation data were obtained from The Cancer Genome Atlas (TCGA) database. The leukemia methylation data of patients were obtained by accessing the TCGA Data Portal website (https://tcga.data.gov/). The selected subtypes of leukemia included AML. Metabolomics data for this study were obtained from the article with the PubMed identifier (PMID: 34193978). The article was accessed through the PubMed database (https://pubmed.ncbi.nlm.nih.gov/34193978/) to retrieve the relevant metabolomics dataset [18].

Differential Expression Analysis

Following data preprocessing, the gene expression profiles of the healthy population are designated as the control group, while the experimental groups consist of the gene expression profile sequences of ALL and AML patients. Differential expression analysis was to identify genes with significant expression differences between the different leukemia types. DESeq2 was employed for statistical analysis and identification of differentially expressed genes.

Functional Enrichment Analysis

To gain insights into the biological functions and pathways associated with the differentially expressed genes, functional enrichment analysis was conducted. The Gene Ontology (GO) terms and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were utilized for this analysis. R packages such as “clusterProfiler” or “Enrichr” were used to perform the enrichment analysis, applying statistical tests such as the hypergeometric test or Fisher's exact test to identify significantly enriched terms or pathways.

Gene Set Enrichment Analysis (GSEA)

In addition to the differential expression analysis, Gene Set Enrichment Analysis (GSEA) was conducted to further explore the biological pathways and functions associated with AML and ALL datasets.

Weighted Gene Co-expression Network Analysis (WGCNA)

WGCNA was employed to construct gene co-expression networks and identify modules of co-expressed genes. The analysis was conducted using the R software package “WGCNA”. A pairwise correlation matrix was calculated using the preprocessed gene expression data. The correlation matrix was transformed into an adjacency matrix using a soft thresholding power parameter (β). The power parameter was selected based on the scale-free topology criterion to ensure a scale-free network structure. The adjacency matrix was transformed into a topological overlap matrix (TOM), which measures the interconnectedness between genes. Hierarchical clustering was performed based on TOM dissimilarity to identify modules of co-expressed genes. The dynamic tree-cutting algorithm was used to define modules based on a predefined minimum module size and a specified height threshold. The module eigengene, defined as the first principal component of the module's gene expression profiles, was calculated. The module eigengene represents the overall expression pattern within a module.

Network Visualization and Analysis

To establish the transcriptional network, we utilized Cytoscape 3.5.1, a powerful tool for visualizing biological networks. In order to gain insights into the network's structure and properties, we employed a Cytoscape plug-in specifically designed for network analysis. The transcriptional network was constructed based on our dataset, and Cytoscape allowed us to visualize the network in a visually appealing and informative manner.

Analysis of methylated genes from TCGA

Differentially methylated genes were screened by comparing cancer tissues and normal tissues, using the “MethylMix” R package with a false discovery rate (FDR) threshold of less than 0.05. The “MethylMix” R package was employed in the R Project for Statistical Computing software. The identified genes with high and low methylation levels underwent bidirectional hierarchical clustering. The “pheatmap” R package was utilized to generate a differential distribution map for the genes with the most significant methylation differences. This analysis enabled the observation of the methylation degree distribution in cancer samples compared to normal tissues. Pearson's correlation test was conducted using the cor.test function in the R language to determine the correlation between gene methylation degree and corresponding gene expression. The filtering criteria for this analysis were a correlation coefficient of less than 0.3 and a p-value of less than 0.05.

Machine learning-based gene selection

Four different machine learning algorithms, namely Generalized Linear Model (GLM), Random Forest (RF), Support Vector Machine (SVM), and XGBoost (XGB), were employed for gene selection. Four R packages, namely “stats” for GLM, “randomForest” for RF, “e1071” for SVM, and “xgboost” for XGB, were utilized to develop predictive models for distinguishing between AML and ALL samples.

Results

Differentially expressed genes screening and analysis of co-expression network

A volcano plot is a type of scatter plot that integrates measures of statistical significance (such as p-values) and magnitude of change from statistical tests, thereby facilitating rapid and intuitive identification of data points (genes, etc.) that exhibit significant and substantial changes. In the volcano plot analysis, a total of 51 genes were found to be downregulated and 50 genes were found to be upregulated in the ALL samples compared to the control samples (Figure 1A left panel). The volcano plot of AML was displayed (Figure 1A right panel). The heatmap analysis revealed distinct expression patterns among the differentially expressed genes. A subset of genes exhibited higher expression levels in the AML/ALL samples compared to the common samples, indicated by the red color in the heatmap. Conversely, another subset of genes showed lower expression levels in the AML/ALL samples, represented by the blue color (Figure 1B).

This study identified 11 genes that were found to be upregulated in both AML and ALL. These genes include: DNMT3A, CCDC198, RNF125, ERG, HHAT, ABHD17B, BMI1, PRKCI, MYRIP, ARMCX4, and MEF2C. These genes exhibited increased expression levels in both AML and ALL, suggesting their potential involvement in the development and progression of these types of leukemia. Besides, it was discovered that four genes were found to be downregulated in both AML and ALL. These genes include VOPP1, MARCO, LINC00472 and WARS1. A protein-protein interaction (PPI) network was constructed using 15 common genes (Figure 1C). Among these genes, BML1 was found to have a shared interaction with DNMT3A and MEF2C (Figure 1D).

Association between genes and protein abundance

Several genes, namely RNP25, ERG, RNF125, PRKC1, WARS1, ABHD17B, and HHAT, were found to exhibit a significant correlation with protein abundance. This suggests that variations in the expression of RNP25, ERG, RNF125, PRKC1, WARS1, ABHD17B, and HHAT may contribute to changes in protein levels in the studied context (Figure 2).

Gene Ontology and Gene Set Enrichment Analysis

The GO analysis of the 15 common genes identified several significant biological processes and molecular functions associated with them. These genes were found to be involved in diverse biological processes, including positive regulation of B cell proliferation, negative regulation of gene expression through epigenetic mechanisms, DNA alkylation, DNA methylation, regulation of gene expression through epigenetic processes, and DNA modification. Additionally, they exhibited associations with cellular components such as heterochromatin, synaptic membranes, compact myelin, and postsynaptic endosomes. The molecular functions of these genes encompassed activities such as S-methyltransferase, protein kinase C, and ubiquitin-protein transferase regulation. These findings provide insights into the potential roles of these genes in B cell proliferation, epigenetic regulation, DNA modification, and various cellular activities. Further investigations are warranted to unravel the precise mechanisms and functional implications of these genes in these processes (Figure S1A). Unfortunately, the study did not find any KEGG pathway enrichment in Gene Set Enrichment. GSEA was performed to identify enriched gene sets. The GSEA results revealed significant enrichment of gene sets of the control group in the following categories: GOLDRATH _ EFF _ VS _ MEMORY _ CD8 _ TCELL _ UP, GOLDRATH _ NAÏVE _ VS _ EFF _ CDBTCELL _ ON, GSE10325 _ BCELL _ VS _ MYELOID _ DN, GSE10325 _ CD4 _ TCELL _ VS _ MYELOID _ ON, GSE10325 _ LUPUS _ BCELL _ VS _ LUPUS _ MYELOID _ DN (Figure S1B).

In addition to the previously mentioned GSEA results, the analysis further identified enriched gene sets in the “treat” group. These findings provide additional insights into the functional differences and regulatory mechanisms associated with the analyzed dataset. The enriched gene sets discovered in this study include: GSE39556 _ CD8A _ DC _ VS _ NK _ CELL _ UP, GSE39556 _ UNTREATED _ VS _ 3H _ POLYIC _ INJ _ MOUSE _ CD8A DC _ DN, GSE39556 _ UNTREATED _ VS _ 3H _ POLYIC _ INJ _ MOUSE _ NK _ CELL _ DN, GSE6269 _ E _ COLIVS _ STAPH _ AREUS _ INF _ PBMC _ UP, GSE9509 _ 10MIN _ VS _ 30MIN _ LPS _ STIM _ IL10 _ KO _ MACROPHAGE _ DN (Figure S1B).

It was conducted to assess the enrichment of gene sets in the control and treat groups, specifically focusing on the KEGG pathway database. The GSEA results revealed significant enrichment of gene sets in the following KEGG pathways such as GOLDRATHLEFF _ VS _ MEMORYCD8TCELL _ JUP and GSE10325 _ BCELLVSMYELODDN (Figure S2A). Disease Ontology (DO) results showed that a significant association with bone cancer, acute leukemia, sarcoma, lymphoblastic leukemia and musculoskeletal system cancer (Figure S2B). The GO analysis results for AML highlight the potential dysregulation of metabolic processes, cellular organization, and signaling pathways (Figure S2C).

Weighted gene co-expression network analysis

The WGCNA analysis provides a comprehensive overview of the gene co-expression patterns and their relationships with phenotypic traits. In this study, we utilized the WGCNA approach to investigate the relationships between samples and traits in the analyzed dataset. The WGCNA analysis generated a sample dendrogram and trait heatmap, providing valuable insights into the clustering patterns and associations among the control samples and treat samples of AML (Figure 3A). To ensure the creation of a scale-free network, a soft threshold power of 17 was selected (Figure 3B). The results include a network heatmap plot that displayed the interactions among selected genes (Figure 3C). Significant correlations between module membership and gene significance for IPAH in the grey and turquoise modules are presented in Figure 3D.

Machine learning-based AML prediction model

Evidence from relevant studies indicates that the prognosis of AML is poorer[19], hence our subsequent research focuses on the selection of biomarkers for AML. This study employed four machine learning algorithms, namely GLM, RF, SVM, and XGB, to construct a predictive model for AML. Using these machine learning approaches, we identified four genes, namely IGHM, FRZB, MGP, and ULK2, as the most important features for predicting AML (Figure 4A). Among the machine learning algorithms utilized in this study, the XGB model demonstrated superior fitting capabilities, as evidenced by its exceptional residual performance (Figure 4B). The analysis of area under the curve (AUC) revealed varying levels of accuracy among the machine learning models. The RF and XGB models exhibited a perfect accuracy with an AUC of 1.000, indicating flawless classification performance. On the other hand, the SVM and GLM demonstrated slightly lower accuracies with AUC values of 0.964. Although not perfect, these models still achieved a high level of accuracy in distinguishing between the classes (Figure 4C).

Gene Set Variation Analysis

The genes associated with AML can be classified into two distinct clusters (Figure 5A). In clusters 1 and 2, there were notable differences in the expression levels of TRM40, EIF5A, and TUBB2A. Specifically, these genes exhibited higher expression levels in cluster 2 compared to cluster 1 (Figure 5B). The results of GSVA revealed noteworthy enrichment of several gene sets in the analyzed dataset. These findings shed light on the molecular functions, cellular components, and biological processes that are prominently associated with the studied dataset. The analysis identified significant enrichment of the gene set GOMF_CUPROUS_ION_BINDING, indicating a heightened affinity for binding to cuprous ions. Additionally, the gene set GOMF_HISTONE_BINDING displayed notable enrichment, suggesting a prevalent interaction with histones. Furthermore, the cellular component GOCC_RNA_POLYMERASE_COMPLEX exhibited substantial enrichment, indicating its prominence in the dataset. This finding suggests the involvement of RNA polymerase complexes in the studied biological context. Regarding biological processes, the gene set GOBP _ NUCLEOTIDE _ EXCISION _ REPAIR showed significant enrichment, indicating its pronounced involvement in the repair of nucleotide excision (Figure 5C). The KEGG pathway “Aminoacyl Biosynthesis” demonstrated a significant enrichment (P value < 0.05), indicating its potential involvement in the observed dataset. Additionally, the pathways “Base Excision Repair”, “Nucleotide Excision Repair”, and “Homologous Recombination” displayed significant enrichment, suggesting their potential roles in DNA repair mechanisms (P value < 0.01) (Figure 5C).

DNA methylation analysis and metabolite expression profiles

Through GO analysis, we found that the biological function of AML enrichment is related to methylation and metabolism, so this study carried out subsequent methylation and metabolism-related analysis. The analysis of DNA methylation patterns revealed several genes that exhibited differential methylation status in tumor samples. Specifically, the genes TRM2, GALNS, ZNF321P, and FOZD7 were found to be highly methylated in tumors. On the other hand, the genes, FZD5, PROZ, and BCL6, showed low levels of methylation in tumor samples (Figure 6A). The volcano plot was displayed in figure 6B and the PCA results revealed a substantial dissimilarity between the tumor samples and the experimental samples (Figure 6C). In the context of AML , the analysis of metabolite expression profiles revealed differential abundance of Oleic acida , palmitic acid, 1-iodo-dodecanea, 2, 6,11-trimethyl-dodecanea, and 2,3- dihydroxypropyl exhibited higher levels of expression in samples compared to normal samples (Figure 6D). It was found that these metabolites were enriched in the metabolic pathways of cysteine and methionine metabolism, as well as taurine and hypotaurine metabolism (Figure 6E). The differential expression of metabolites in AML and ALL compared to healthy controls were investigated. The analysis of LogFC (fold change) values for various metabolites was visualized using a heatmap, as shown in Figure 6F.

Discussion

In this study, we initially focused on 15 genes that are common to both ALL and AML, with 11 upregulated and 4 downregulated. These genes may be involved in shared biological processes or pathways related to the pathogenesis of AML and ALL, including cell proliferation, apoptosis regulation, and cell signaling. We further elucidated gene functions and pathways through GO, KEGG, and GSEA enrichment analyses.

Additionally, due to the higher incidence, rapid progression, complex treatment, and poorer prognosis of AML, our subsequent research primarily involves data analysis focused on AML. Utilizing methods like WGCNA, a systems biology approach, we aim to understand gene relationships within the expression matrix comprehensively. This approach helps uncover potential biological mechanisms and biomarkers, providing crucial insights for further investigation [20].

Furthermore, we constructed a machine learning model for the prediction of AML by selecting four genes, namely IGHM, FRZB, MGP, and ULK2, as predictive features. Immunoglobulin Heavy Constant Mu (IGHM) is a critical component of the immune system and is a key subtype of immunoglobulin M (IgM) [21]. One study demonstrated the presence of five Ig classes, including IGHM, in AML [22]. Another study revealed the upregulation of IGHM expression in AML patients, indicating its potential role in the disease [23]. Importantly, a significant correlation between IGHM expression and patient survival, along with age as an additional contributing factor was observed [23].

Another gene with predictive value, FRZB, acts as an antagonist of the Wnt signaling pathway, which plays a critical role in regulating various cellular processes, including proliferation, differentiation, and apoptosis [24,25]. Dysregulation of the Wnt pathway has been implicated in various cancers, including hematological malignancies such as AML [26,27]. The significant difference in FRZB expression between AML patients and the control group further highlights its potential as a biomarker for AML which suggested its involvement in the disease's molecular mechanisms [28]. It may potentially provide insights into the underlying molecular pathways driving leukemogenesis.

MGP, an extracellular matrix (ECM) protein, has been demonstrated to possess inhibitory effects on the mineralization and apoptosis of chondrocytes [29]. These observations imply that the crosstalk between leukemia cells and the bone marrow (BM) stroma induces alterations in the tumor microenvironment, promoting an undifferentiated state that may potentially contribute to the progression of the disease and the development of resistance to chemotherapy [30]. Upregulation of MGP in AML has been observed in the study [31]. This finding highlights the potential involvement of MGP in the pathogenesis of AML and suggests its potential role as a biomarker or therapeutic target for the disease.

ULK1, also known as ATG1, is a serine/threonine kinase that plays a pivotal role in the initiation of autophagy, a highly conserved cellular process involved in the degradation and recycling of damaged or unnecessary cellular components [32]. Growing evidence suggests that dysregulated autophagy is a distinguishing feature of malignant disorders and plays intricate roles in the initiation and advancement of cancers, as well as their resistance to treatment [33-35]. In addition to its role in autophagy initiation, ULK1 has been implicated in various cellular processes, including cell growth, survival, and differentiation. The diverse functions of numerous autophagy-related proteins in cancer have been extensively examined [36]. Studies have revealed that heightened expression levels of ULK1 are correlated with unfavorable prognosis in various solid tumors [37,38]. Notably, it was found that the expression of ULK1 in primary AML samples was significantly increased compared to healthy control samples [39]. Besides, a growing body of literature strongly suggested that ULK1 is a viable drug target for the treatment of AML [40].

The present study focused on the investigation of stearic acid (SA) and its potential association with AML. Stearic acid, a saturated fatty acid, is known for its involvement in various cellular processes and has been implicated in the pathogenesis of several diseases. In our study, we observed that stearic acid displayed a statistically significant higher concentration in AML plasmas, consistent with previous findings reported in the literature [41]. These results are in line with the growing body of evidence suggesting a potential link between stearic acid and AML [42]. The elevated levels of stearic acid in AML plasmas may indicate its role as a biomarker or a potential contributing factor in the development and progression of AML. However, further investigations are required to elucidate the underlying mechanisms by which stearic acid influences AML pathogenesis.

Conclusion

The comprehensive analysis of AML and ALL genes in this study has yielded valuable insights into the underlying mechanisms of leukemia and the identification of potential biomarkers. Furthermore, the application of machine learning techniques has successfully identified four biomarkers that hold promise for enhancing the diagnosis and treatment of AML. These findings contribute to the understanding of leukemia pathogenesis and offer potential avenues for personalized medicine and improved patient outcomes.

However, the clinical relevance of these findings needs further exploration. Future research should focus on developing targeted therapies that address the specific gene alterations identified in this study. For instance, targeted inhibitors could be designed based on the molecular pathways disrupted by these biomarkers. Additionally, the implementation of advanced genomic sequencing techniques could enable the detection of these gene changes in patients, allowing for early prediction of disease progression and tailoring treatment strategies accordingly. Moreover, the potential correlation between these biomarkers and patient prognosis merits further investigation. A longitudinal study tracking patients with identified biomarker alterations could provide insights into their prognostic value and influence on treatment outcomes.

It is also important to acknowledge the limitations of our analysis. While our study identifies potential biomarkers, it does not establish causation. The relationship between gene expression levels and disease manifestation requires further validation. Experimental approaches such as gene knockdown or overexpression in relevant cellular models could elucidate the functional impact of these genes on leukemia development. Additionally, in vivo studies using animal models may help confirm the role of these biomarkers in disease progression and therapeutic response. In conclusion, the findings from this study pave the way for further research into the clinical applicability of identified biomarkers, enhancing our understanding of leukemia and guiding the development of personalized treatment strategies.

Funding

This study was supported by grants from Public Welfare Research Project of Zhejiang Province (grant number LGF22H080008 and LGF20H080001), Medical and Health Science and Technology Project of Zhejiang Province (grant number 2022482097), Traditional Chinese Medicine Science and Technology Project of Zhejiang Province (grant number 2021ZB245), and Scientific Research Project of Education Department of Zhejiang Province (grant number Y202249376).

Conflict of Interest

The authors declared no conflict of interest.

Author’s Contribution

L.N. Wang: Data curation, Software, Funding acquisition, Formal analysis, Validation and Writing – original draft. W.Z. Shang: Data curation, Formal analysis, Software and Funding acquisition. Y.F. Wu: Data curation, Formal analysis and Funding acquisition. G.B. Zhang: Data curation, Formal analysis and Funding acquisition. J. Xie: Conceptualization, Resources, Supervision, Funding acquisition and Writing – review and editing.

Data Statement

Data available on request from the authors.