Integrative Analysis of Gene Expression Through One-Class Logistic Regression Machine Learning Identifies Stemness Features in Multiple Myeloma

Tumor progression includes the obtainment of progenitor and stem cell-like features and the gradual loss of a differentiated phenotype. Stemness was defined as the potential for differentiation and self-renewal from the cell of origin. Previous studies have confirmed the effective application of stemness in a number of malignancies. However, the mechanisms underlying the growth and maintenance of multiple myeloma (MM) stem cells remain unclear. We calculated the stemness index for samples of MM by utilizing a novel one-class logistic regression (OCLR) machine learning algorithm and found that mRNA expression-based stemness index (mRNAsi) was an independent prognostic factor of MM. Based on the same cutoff value, mRNAsi could stratify MM patients into low and high groups with different outcomes. We identified 127 stemness-related signatures using weighted gene co-expression network analysis (WGCNA) and differential expression analysis. Functional annotation and pathway enrichment analysis indicated that these genes were mainly involved in the cell cycle, cell differentiation, and DNA replication and repair. Using the molecular complex detection (MCODE) algorithm, we identified 34 pivotal signatures. Meanwhile, we conducted unsupervised clustering and classified the MM cohorts into three MM stemness (MMS) clusters with distinct prognoses. Samples in MMS-cluster3 possessed the highest stemness fractions and the worst prognosis. Additionally, we applied the ESTIMATE algorithm to infer differential immune infiltration among the three MMS clusters. The immune core and stromal score were significantly lower in MMS-cluster3 than in the other clusters, supporting the negative relation between stemness and anticancer immunity. Finally, we proposed a prognostic nomogram that allows for individualized assessment of the 3- and 5-year overall survival (OS) probabilities among patients with MM. Our study comprehensively assessed the MM stemness index based on large cohorts and built a 34-gene based classifier for predicting prognosis and potential strategies for stemness treatment.


INTRODUCTION
Multiple myeloma (MM), a clonal B-cell malignancy (Bataille and Harousseau, 1997), is characterized by the aberrant proliferation of bone marrow plasma cells and the overproduction of lightchain or monoclonal immunoglobulin (Röllig et al., 2015;Dhakal et al., 2016). Among the diseases with abnormal plasma cells, MM is the second most common hematologic malignancy (Röllig et al., 2015). Owing to the development of immunomodulators and proteasome inhibitors, patients with MM have significantly improved survival over the past decade. However, progressive disease remains a common cause of fatal outcomes (Martinez-Lopez et al., 2011). As a result, it is essential to understand the underlying mechanisms that contribute to disease relapse and progression in MM, as well as novel targets for therapeutic improvements or prognostic prediction.
Stemness was defined as the potential for differentiation and self-renewal from the cell of origin (Seguin et al., 2015;Prasetyanti and Medema, 2017). It has been reported that the gradual loss of differentiation capacity and acquisition of stem cell-like characteristics are important factors that promote tumor progression (Seguin et al., 2015). A growing body of research has confirmed the existence of cancer stem cells in multiple tumors, including hematological malignancies (Lapidot et al., 1994;Singh et al., 2004;O'Brien et al., 2007). Cancer stem cells play a key role in tumorigenesis, progression, metastasis, and drug resistance. Therapies targeting cancer stem cells are of great value for tumor prevention and treatment. For example, the use of retinoic acid to induce the maturation and differentiation of malignant proliferating cells has achieved great success in the clinical treatment of promyelocytic leukemia (Nowak et al., 2009). With the application of serial transplantation models and clonogenic in vitro assays, MM stem cells have been suggested to be part of a subset of CD38 − CD19 + CD27 + B-cell precursors that do not express the classic MM markers CD38 or CD138 (Matsui et al., 2008). These cells with clonogenic potential could mediate tumor regrowth and chemoresistance. The one-class logistic regression (OCLR) machine learning algorithm is an effective method of quantifying the cancer stemness index using two independent indices (Malta et al., 2018;Lian et al., 2019). One is the stemness index [mRNA expression-based stemness index (mRNAsi)] based on gene expression that reflects mRNA expression, and the other is mDNAsi, which reflects epigenetic characteristics (Sokolov et al., 2016). Previous studies have proven the effective application of stemness indices calculated by the OCLR algorithm in a variety of malignant tumors. However, there are few studies on the role and prognostic value of MM stemness (MMS); hence,it is an urgent need to develop prognostic or predictive biomarkers associated with stemness index.
In this study, we collected a total of 1,095 newly diagnosed or pre-treatment MM patients with expression data and clinical information and systematically evaluated the MM stemness index (mRNAsi) using the OCLR algorithm. By combining weighted gene co-expression network analysis (WGCNA) with MM mRNAsi, we searched for key genes related to stemness in 1,095 MM patients. The analysis of gene and module functions showed significance in MM. We classified MM patients into different subgroups based on the expression of stemness-related genes. Finally, mRNAsi was integrated with other clinicopathological characteristics to construct a nomogram to predict the prognosis of patients with MM.

Data Collection and Pre-processing
The workflow of this study is shown in Figure 1. In this study, two MM cohorts were downloaded from the Gene Expression Omnibus (GEO) 1 database. The related information contained clinical, molecular, and microarray datasets. Samples in the GSE4204 cohort (N = 538) were pre-treated with bone marrow aspirates from MM patients (Driscoll et al., 2010). The GSE24080 cohort consisted of 559 newly diagnosed MM patients (Shi et al., 2010).
The gene expression values in the two cohorts were analyzed using the robust multi-array average algorithm for background correction, quantile normalization, and final summarization. After excluding samples without survival information, a total of 1,095 samples were selected for subsequent analysis (558 cases from the GSE24080 cohort and 537 cases from the GSE4204 cohort). Malta et al. (2018) provided a novel OCLR machine learning algorithm to assess oncogenic dedifferentiation that considered mRNAsi. We utilized the method to calculate the mRNAsi for each MM patient. The gene expression-based stemness index was mapped from zero to one. We assessed the prognostic value of mRNAsi and its relationship with other clinical features in different cohorts. First, we utilized univariate and multivariable Cox proportional hazard regression to calculate the hazard ratios (HRs) for overall survival (OS). Second, each cohort of MM patients was divided into low and high mRNAsi groups based on the same cutoff value of mRNAsi (0.4124166, identified by the Survminer package), and the Kaplan-Meier method and log-rank tests were used to determine the significance of survival differences.

Identification of Key Modules and Genes Associated With mRNAsi by WGCNA
The WGCNA was developed to discover correlations among genes by constructing significant modules (Langfelder and Horvath, 2008). The mRNAsi and expression modules were calculated using WGCNA to identify key modules and genes associated with mRNAsi. Module eigengenes (MEs) were defined as the first principal component of each gene module and were adopted as the representative of all genes in each module. Gene significance (GS), as the mediator p-value (GS = lgP) for each gene, represented the degree of linear correlation between gene expression of module and mRNAsi or other clinical features. MRNAsi-related modules were defined according to P < 0.01, and a higher GS value was used for further analysis. Finally, we selected the intersection of the key genes identified in different data sets.

Screened Differentially Expressed Genes
The limma R package was used to screen differentially expressed genes (DEGs) between the high and low mRNAsi groups (Ritchie et al., 2015). The selection criteria were as follows: adjusted P-value < 0.0001 and | log2 fold change| > 0.5.

Functional Annotation and Pathway Enrichment Analysis
We first selected the intersection between the modular genes identified by WGCNA and the DEGs identified by differential expression analysis and then adopted a univariate Cox regression model to assess the prognostic significance of the overlapping genes. Only genes with a P-value < 0.05 were considered as prognosis-related hub genes. Functional annotation and pathway enrichment analysis of these hub genes was performed using Metascape (Zhou et al., 2019). 2 2 https://metascape.org/gp/index.html#/main/step1

Consensus Clustering and Prognostic Analysis Based on mRNAsi-Related Hub Genes
Protein-protein interaction enrichment analysis was conducted, and the molecular complex detection (MCODE) algorithm (Bader and Hogue, 2003) was used to identify densely connected network components. Based on the 34 mRNAsi-related hub genes identified by the MCODE algorithm, we subsequently performed unsupervised clustering on MM patients. The R package ConsensusClusterPlus was used to conduct unsupervised clustering (1,000 iterations, 80% resampling), and k-means and Euclidean distances were used as the clustering algorithm and distance metric, respectively (Wilkerson and Hayes, 2010). We then assessed the prognosis in each MMS-cluster via Kaplan-Meier analysis. Meanwhile, mRNAsi was also compared in distinct MMS clusters.

Development and Validation of a Nomogram for Prognosis Prediction of MM Patients
To develop a clinically applicable method of predicting MM prognosis, we integrated mRNAsi and other clinicopathological covariates to build a nomogram. Predictive factors included mRNAsi, cytogenetic abnormalities, albumin (ALB), beta-2 microglobulin (B2M), and lactate dehydrogenase (LDH). The nomogram was verified using receiver operating characteristic (ROC) curves and calibration curves.

Statistical Analyses
The Wilcoxon rank-sum test was used to compare two groups with non-normally distributed variables, while Student's t-test was conducted to compare groups with normally distributed variables. The one-way analysis of variance or Kruskal-Wallis tests were used to compare differences between three or more groups. Based on the correlation between mRNAsi and patient survival, the cutoff point of mRNAsi was determined using the Survminer R package. The "surv-cutpoint" function, which repeatedly tested all potential cut points to find the maximum rank statistic, was applied to dichotomize mRNAsi, and the patients were then divided into high and low mRNAsi groups based on the maximally selected log-rank statistics to decrease the batch effect of calculation. The survival curves for the prognostic analysis were generated using the Kaplan-Meier method, and log-rank tests were used to determine the significance of differences. Independent prognostic factors were ascertained using the multivariable Cox regression model.

The Association Between mRNAsi and Patient Survival
We collected 1,095 MM samples with clinical information and corresponding expression data to systematically characterize the stemness features of MM. The overall characteristics of the MM cohort are listed in Table 1. MM samples were sorted according to their mRNAsi values (from low to high stemness index) and examined whether any clinical characteristic/molecular/demographic was associated with either a low or high stemness index (Figures 2A, 3A). We did not find significant differences in mRNAsi in terms of gender and race, but we found that patients who died and patients who experienced adverse events occurred often showed a higher mRNAsi.

Correlations of mRNAsi With Clinical Prognosis in MM
We further explored the prognostic value of mRNAsi in patients with MM. First, the multivariate Cox regression model analysis in different cohorts confirmed mRNAsi as an independent and robust prognostic biomarker for evaluating MM patient outcomes ( Table 2 and Supplementary Tables 1, 2). Second, based on the cutoff point (0.4124166) identified by the Survminer R package, each MM cohort was separated into low and high mRNAsi groups. The Kaplan-Meier analysis confirmed that the high-mRNAsi group had shorter OS and event-free survival (EFS) than the low mRNAsi group in the GSE24080 cohort ( Figures 2B,C). Consistent with this, in the GSE4204 cohort, the Kaplan-Meier analysis indicated that the high-mRNAsi group had shorter OS than the low mRNAsi group ( Figure 3B).
Surprisingly, we found that in both the GSE24080 and GSE4204 cohorts, the treatment regimen (TT2 or TT3) did not improve the survival of MM patients (Figures 2D, 3C). Regardless of the treatment plan (TT2 or TT3), it failed to improve the survival disadvantages of patients in the high-mRNAsi group (Figures 2E,F, 3D,E). Consistently with previous reports (Zhan et al., 2006), our analysis found that among the seven subgroups of MM, the PR subgroup had the worst prognosis, followed by the MS and MF subgroups ( Figure 3F). Finally, we investigated the association between mRNAsi and other clinical characteristics. We found that the difference in mRNAsi among the groups was statistically significant when comparing groups of patients with B2M ≥ 3.5 and those with B2M < 3.5 ( Figure 4A). There was no significant difference in mRNAsi between the different isotype groups ( Figure 4B). However, for the seven subgroups of MM, consistent with poor survival, patients in the PR subgroups possessed the highest mRNAsi ( Figure 4C). In terms of cytogenetic abnormalities, patients with cytogenetic abnormalities tended to have higher mRNA levels ( Figure 4D). Since different therapies (TT2 or TT3) failed to improve the survival of MM patients, we further investigated whether there were any differences in mRNAsi among MM patients who received different therapies. Consistent with our conjecture, different treatment regimens failed to effectively affect mRNAsi (Figures 4E,F). Based on the ESTIMATE algorithm, we calculated the individual stromal and immune scores to evaluate the level of infiltrating stromal and immune cells in any given MM sample (Yoshihara et al., 2013). We found that  the high-mRNAsi group had lower immune and stromal scores than the low mRNAsi group, which suggested that MM patients with high mRNAsi had a lower level of infiltration of tumor microenvironment (TME) cells (Figures 4G,H).

Weighted Gene Co-expression Network Analysis: Identification of the Most Significant Modules and Genes
Weighted gene co-expression network analysis was used to build a gene co-expression network to classify all genes into biological gene modules based on average linkage hierarchical clustering and further identify genes strongly associated with MM stemness (Figures 5A,D, 6A,D). The soft-thresholding powers in the WGCNA were determined based on scale-free R2 (Figures 5B,C, 6B,C). The pink module had the highest correlation with mRNAsi in the GSE24080 cohort (Figures 5E,F). Additionally, the tan and green modules were highly associated with mRNAsi in the GSE4204 cohort (Figures 6E-G). The pink module contained 643 genes (Figure 5F), and tan and green contained 313 and 623 genes, respectively (Figures 6F,G). There were 379 overlapping genes in the three modules ( Figure 7A). A total of 379 overlapping genes were retained for further analysis.

Differentially Expressed Genes Between Low and High mRNAsi Groups of MM Patients
Since the prognosis of MM patients in the high and low mRNAsi groups differed significantly, we conducted a differential expression analysis between high and low mRNAsi groups to identify the differentially expressed key genes that regulate stemness. We identified 370 DEGs, of which 326 were upregulated and 44 were downregulated in the high-mRNAsi group. After merging the 370 DEGs with the 379 co-expressed genes identified by WGCNA, 128 hub genes were retained for further analysis (Figure 7B).

Pathway and Process Enrichment Analysis of Hub Prognostic Genes
The prognostic evaluation of the 128 hub genes was performed using the Cox proportional hazards regression analysis in the GSE24080 cohort dataset. In total, 127 hub genes were identified as significantly prognostic of the outcome (Supplementary Table 3, P < 0.05). Functional enrichment analysis was performed using Metascape to elucidate the biological functions of the 127 hub prognostic genes ( Figure 7C).  The enrichment analysis suggested that these genes were significantly enriched in pathways and processes related to cell cycle, cell differentiation, and DNA replication and repair ( Figure 7D, hypergeometric test P < 0.01). Finally, based on the MCODE algorithm, the 127 hub prognostic genes were divided into four densely connected network components. The MCODE networks identified for the individual gene lists are shown in Figure 7E. A total of four major gene modules (MCODE 1-4) were identified, including 34 hub prognostic genes (Supplementary Table 4).

Consensus Clustering to Distinguish Different Stemness Prognostic Subtypes
First, we identified different stemness prognostic subgroups in the GSE24080 cohort. We utilized the 34 stemness-related signatures of MCODE 1-4 to conduct unsupervised clustering and identified the stemness molecular subtypes of MM for prognostic analysis. The R package of ConsensusClusterPlus was used to iterate 1,000 times for the stabilization of classification categories (parameters: pItem = 0.8, reps = 1,000, pFeature = 1), and three distinct stemness molecular subgroups were eventually identified using unsupervised clustering (Figures 8A,B). Unsupervised clustering is a useful technique in tumor research, where intrinsic groups sharing biological characteristics may exist but are unknown. A Kaplan-Meier analysis was conducted across the three clusters, where MMS-cluster3 had the worst OS prognosis and MMS-cluster2 had the most favorable prognosis (Figure 8C, log-rank test P < 0.0001). In agreement with these results, MMS-cluster3 had the highest mRNAsi, and MMS-cluster2 had the lowest mRNAsi ( Figure 8D). In addition, to verify our stemness prognostic subgroups, based on the 34 stemness-related signatures, we also performed unsupervised clustering and prognostic analysis on the GSE4204 cohort. Consistent with the above results, the MM samples could be classified into three clusters with different prognoses ( Figure 8F); MMS-cluster3 had the worst prognosis with the highest mRNAsi, while MMS-cluster2 had the most favorable prognosis with the lowest mRNAsi ( Figure 8G). We also investigated whether there were any differences in immune and stromal scores among the three clusters. The results showed that MMS-cluster2 had higher stromal and immune scores, while MMS-cluster3 had lower stromal and immune scores (Figures 8E,H). In total, these analyses indicated that the 34 MMS-related signature could guide molecular classifications, by which the MMS clusters possessed different immune infiltration and had different prognostic outcomes.

Development and Validation of a Nomogram
Based on the multivariate Cox regression analysis, a nomogram was built to determine MM prognosis (Figure 9A), and a timedependent ROC curve was used to evaluate the effect of the nomogram. As the data with complete clinical information was limited, we used the GSE24080 cohort as the whole dataset (training set) to construct the nomogram. To verify our nomogram, we set the random seed to 123 and used the caret R package to randomly cut the whole set and divided the samples into two different data sets (50%/50%) to verify the nomogram.
In the whole dataset, the area under the curve was 0.759 and 0.740 for the 3-and 5-year survival, respectively (Figures 9B,C). The performance of the nomogram was far superior to that of mRNAsi, cytogenetic abnormalities, ALB, B2M, and LDH alone for assessing patient prognosis (Figures 9B,C). In addition, the nomogram showed good prediction performance in the verification cohorts, and the 3-and 5-year ROC curves of test sets 1 and 2, respectively, are shown in Figures 9D-G (Figures 9D,E,  test set 1; Figures 9F,G, test set 2). A calibration curve was constructed to assess the accuracy of the nomogram. As shown in Figures 9H,I (training set) and Figures 9J-M (Figures 9J,K, test  set 1; Figures 9L,M, test set 2), the combined nomogram showed good performance in predicting the 3-and 5-year survival rates of patients, and the prediction probability was close to the actual observed situation.

DISCUSSION
Cancer stem cells are commonly considered to be responsible for tumor persistence and progression, tumor recurrence, and resistance to traditional therapies (Seguin et al., 2015;Prasetyanti and Medema, 2017;Saygin et al., 2019). We comprehensively analyzed data from large cohorts to determine the cancer stemness of MM and found that mRNAsi was tightly associated with OS and EFS, possessing good fitness and quantification of stemness in MM. In addition, low-and high-mRNAsi groups had different features, including B2M, cytogenetic abnormalities, ALB, subgroups, and immune infiltration. This evidence indicates that stemness, as a yardstick to assess the degree of gradual loss of a differentiated phenotype and gain of progenitor and stem cell-like characteristics, strongly affects the prognosis of MM. Our researches showed that mRNAsi could effectively classify patients with MM into groups with low and high risks of poor prognosis. Additionally, the proposed mRNAsi provided additional prognostic value to existing clinicopathological prognosticators for MM. Of particular importance, this is the first study to demonstrate the clinical utility of the stemness signature as a prognostic tool in patients with MM. We observed that the prognoses of patients with MM who received different therapies (TT2 or TT3) were not significantly different. Similarly, mRNAsi did not differ significantly; however, in both the TT2 and TT3 groups, patients with high mRNAsi had shorter OS than patients with low mRNAsi. This suggests that the TT3 regimen did not significantly improve the prognosis of patients compared with the TT2 regimen. This might be because the TT3 regimen failed to remove the stemness, thus failing to influence the degree of oncogenic dedifferentiation, which might provide valuable insights for targeted therapies aimed at tumor differentiation of MM. Moreover, the stemness index-related modules and genes were identified using WGCNA and differential expression analysis. A total of 127 hub prognostic genes were found to be most significantly associated with stemness index. Pathway and process enrichment analysis revealed that these hub prognostic genes were involved in various biological functions related to the initiation and progression of MM, including the cell cycle, cell differentiation, and DNA replication and repair. Interestingly, these hub prognostic genes allowed for the discovery of innovative targets and possible targeted therapies aimed at tumor differentiation. Based on the MCODE algorithm, the most important connected network components were identified from these hub prognostic genes. We considered tumor heterogeneity and classified the MM cohorts into three MMS clusters based on the 34 stemness-related signatures. We observed distinct survival outcomes across the three MMS clusters and compared  the differential mRNAsi among the three MMS clusters. In particular, MMS-cluster3 samples possessed the highest mRNAsi and had the worst OS outcomes compared with other clusters. We also applied the ESTIMATE method to evaluate the level of immune infiltrating cells in the TME of MM and found different infiltration patterns across the three clusters. These findings revealed the patterns of intra-tumor molecular heterogeneity and the different patterns of TME infiltration within MM and supported the negative regulation of the stemness index and anticancer immunity. Lastly, by integrating mRNAsi and other clinical characteristics, we proposed a prognostic nomogram that allows for individualized estimations of the 3-and 5-year OS probabilities among MM patients. Stemness signatures have been identified in many malignancies and have different prognostic values in gastric cancer, acute myeloid leukemia, prostate cancer, and breast cancer (Ng et al., 2016;Malta et al., 2018;Wang et al., 2018;Miranda et al., 2019;Chang et al., 2020;Zhang et al., 2020). However, there are few studies on the stemness of MM and its prognostic value in MM patients. In our study, mRNAsi could stratify MM patients into two groups with notably different prognoses, which was convenient for risk stratification of MM in clinical practice. Integrated analyses also confirmed that mRNAsi is an independent prognostic biomarker in MM. In our screening in multiple cohorts, some of the stemness-related genes have never been reported in other previous researches that could predict the outcomes of MM, which was convenient to implement in clinical practice. BUB1B encodes a kinase involved in spindle checkpoint function. The protein is localized to the kinetochore and plays a role in the inhibition of the anaphase-promoting complex/cyclosome. An increasing body of literature has verified that aberrant expression of BUB1B is highly involved in the tumorigenesis and the development of various tumors (Wan et al., 2012;Qiu et al., 2020). A previous study revealed that BUB1B could promote MM cell proliferation through the CDC20/CCNB axis (Yang et al., 2015). MCM2, MCM3, MCM5, and MCM6, members of the minichromosome maintenance (MCM) family, are highly involved in DNA replication and are vital in limiting replication in the cell cycle (Freeman et al., 1999;Forsburg, 2004). Some investigations have shown that the expression of MCM family plays an important role in the prognosis of MM, and MCM2 is an independent risk factor for MM (Quan et al., 2020). Additionally, MCM2 is associated with many types of cancer, including acute lymphocytic leukemia, gallbladder cancer, and glioma (Hua et al., 2014;Liu et al., 2016;Li et al., 2018). Trichostatin A, a classical histone deacetylase inhibitor, could downregulate the expression of MCM2, and the silencing of MCM2 in colon cancer cells could induce cell cycle arrest and apoptosis as reported in a previous study (Liu et al., 2013). Therefore, MCM2 could be a potential therapeutic target for the treatment of MM. ZWILCH kinetochore protein is an essential part of the Rod-Zw10-Zwilch complex and is important in maintaining the normal function of mitotic checkpoints (Karess, 2005;Gama et al., 2017). The abnormal function of mitotic checkpoints is related to the appearance of chromosomal instability, a consensus sign of many human cancers. In MM, chromosomal instability contributes to the acquisition of tumor heterogeneity and thereby to drug resistance, disease progression, and eventual treatment failure (Fujibayashi et al., 2020;Neuse et al., 2020). Aberrant BUB1 overexpression promotes mitotic segregation errors and chromosomal instability in MM (Fujibayashi et al., 2020), and the synergistic mechanism of BUB1B and ZWILCH in the occurrence and development of MM requires further investigation.
Moreover, we proposed prognostic nomogram that contribute to individualized evaluations of the 3-and 5-year OS probabilities among patients with MM. Taken together, mRNAsi and the associated nomogram might serve as a clinically helpful tool to improve surveillance and guide decision-making regarding the administration of adjuvant chemotherapy.
Collectively, mRNAsi could effectively classify patients with MM into groups with different risks of outcomes, thereby raising the possibility that stemness might be supplementary to the conventional clinicopathological risk factors as a prognostic scheme. The 34-gene based MMS-related signature could be a good molecular classifier for uncovering distinct stemness clusters. Additionally, the proposed nomogram incorporating mRNAsi and existing clinical prognosticators might facilitate personalized surveillance and management of patients with MM.

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/s.

AUTHOR CONTRIBUTIONS
JH designed the study, interpreted the data, and revised the manuscript. CB and FY performed the statistical analysis, interpreted the data, and wrote the manuscript. MW, QL, and JW provided the bioinformatic data and contributed to the manuscript revision. LC, LTL, DX, and LL provided the key scientific insights and contributed to the manuscript revision. All authors have read and approved the final version of the manuscript.

FUNDING
This work was financially supported by grants from the Selffunded Scientific Research Project of the Guangxi Zhuang Autonomous Region Health Commission (contract number: Z20200114) and the "The role and mechanism of quercetin in acute myeloid leukemia" project (contract number: lryjj201909).