MDS represents a heterogenous cluster of myeloid neoplasms, featured by dysplasia, ineffective hematopoiesis, with or without excessive blasts. Leukemic transformation is one of the main causes of death in MDS patients (Platzbecker et al., 2021). However, evaluation of transformation risk for individual patients, remains an essential but difficult aspect in clinical investigations. By far several prognostic systems have been established, such as IPSS (International Prognostic Scoring System (Greenberg et al., 1997)), IPSS-R (Revised International Scoring System (Greenberg et al., 2012)) and WPSS (WHO Prognostic Scoring System (Malcovati et al., 2007)), which are based on clinical (cytopenia, blast percentage, disease subtypes, etc.) or genetic variables (cytogenetics). The risk categories correlated with duration of leukemic transformation without treatment, according to which the risk-adapted treatment strategy was used in clinical practice. Nevertheless, it is recognized that patients with the same genetic signature (mutation/cytogenetic variation) frequently have distinct clinical and prognostic features. Genes involving in RNA splicing, epigenetic modification and signaling transduction, were most frequently mutated in MDS (Lindsley and Ebert, 2013), which lead to dysregulation of gene expression. Beyond DNA variation, the emerging transcriptomic investigations unraveled the gene expression signature for MDS (Pellagatti et al., 2006; Pellagatti et al., 2010; Gerstung et al., 2015), and a few prognostic models were built. The work of Moritz Gerstung et al. even proved that the predictive power of expression signature was superior to that of traditional markers (IPSS, genetics, cytogenetic, etc.) on AML-free survival (Gerstung et al., 2015). While these studies were lack of external validation by independent cohorts, and the robustness of models was questioned due to the heterogeneity between experimental platforms and cell origins. Additionally, these models were either too elaborate (included too many variables), or derived from CD34⁺ sorted cells, which limited the feasibility in practice (unsorted samples instead of CD34⁺ sorted were used in most clinics).

New methods had emerged to analyze the transcriptomic datasets, which made it possible to improve the accuracy of the prediction model. An updated transcriptomic analysis method, WGCNA, can recognize co-expressed gene clusters according to scale-free network theory and correlated clinical or genetic traits with MEs (the first principal component representing gene co-expression modules) (Zhang and Horvath, 2005; Langfelder and Horvath, 2008). By WGCNA, the co-expressed gene clusters, significantly associating with AML transformation or high-risk factors, were unraveled using GSE58831 dataset in the present study. MDS15 model was established using multiple-iteration LASSO, consisting of 15 gene expression variables. The MDS15 model comprises expression markers from 15 genes: NEAT1, LYSMD2, SLC4A1AP, KMT2A, PHC1, ADHFE1, TFAP2E, TPBG, TRIP11, GAS6-AS1, KCNMB4, ZNF225, LOC100506730, WT1, and STARD9. The NEAT1 gene yields a lncRNA known to influence AML progression (Zhao et al., 2019; Feng et al., 2020; Yan et al., 2021; Rostami et al., 2022). KMT2A also stands as a prognostic determinant in MDS/AML (Li et al., 2013; Kotani et al., 2019), with its partial tandem duplication or rearrangement traditionally recognized as harbingers of diminished survival. GAS6-AS1, another long noncoding RNA, is notably overexpressed in AML and correlates with adverse survival outcomes; curbing its expression has been shown to decelerate AML progression (Zhou et al., 2021). Several studies have illuminated that elevated WT1 expression detrimentally affects the survival rates in MDS, and its expression in peripheral blood astutely forecasts progression-free survival (Rautenberg et al., 2019). However, the remaining genes delineated in the MDS15 model have not been definitively linked to either MDS prognosis or the biology of myeloid malignancies.

Then to address the versatility of MDS15, predictive power was validated in both GSE58831 and GSE15061 datasets. Furthermore, GSEA was implemented to uncover the possible related biological process to MDS15 risk scores. The flowchart of this study was shown in Figure 1.

The significantly enriched pathways relating with MDS15 scores were illustrated in Figure 13. The activated pathways consisted of MYC (Myc proto-oncogene protein) targets, E2F (Transcription factor E2F1) targets, Oxidative Phosphorylation and DNA repair pathways, etc., significantly suppressed pathways constituted apoptosis, IL6-JAK-STAT3 pathway, IL2-STAT5 pathway, Interferon alpha response pathway, apoptosis and p53 pathway, etc.

The transcriptomic features have been widely investigated for bone marrow samples of MDS patients, to link the mRNA expression signature with clinical outcomes or eligibility of target therapy (Pellagatti et al., 2006; Pellagatti et al., 2010; Kondo et al., 2011; Heinrichs et al., 2013; Xu et al., 2016; Huang et al., 2019). The updated analytic methods have provided new clues on the ‘old data’. WGCNA has been used in digging from omics data of hematological malignancies (Fu et al., 2020; Liu et al., 2020; Adebayo et al., 2022; Chen et al., 2022). In comparison with traditional differential expression analysis, WGCNA focused on the interplay of co-expressed modules based on scale-free network theory, so that the gene clusters instead of individual genes were investigated to establish complicated, but robust correlation with clinical variables.

The current MDS risk stratification systems profiling survival and leukemic transformation, were challenge by emerging omics data, including genomics, epigenetics (methylation or chromosomal accessibility), transcriptomic, and proteomics. The omics data profiled the actual MDS cell status and further predict prognosis or discern therapeutic targets. In the present study, the co-expressed gene modules were recognized, among which the brown module was identified as the key gene cluster associating with AML transformation and high-risk clinical factors (RAEB phenotype, etc.. 1,301 genes constituting brown module (Supplementary Table S1), was demonstrated to be predominantly enriched in epigenetic regulating pathways, by biological process (BP) analysis based on GO database. One node of this pathway in brown module, KMT2A gene encodes histone methyltransferase, which plays an essential role in hematopoiesis (Nakamura et al., 2002; Dou et al., 2005; Patel et al., 2009; Avdic et al., 2011), rearrangement or abnormal expression of which had been linked to poor prognosis (Tsai et al., 2022). Other histone methylating genes, such as KMT5B/EZH1/KMT2D/KMT2E/SETD5, were also included in brown module, and related to pathological process or prognosis of myeloid neoplasms (MDS/AML) (Mochizuki-Kashio et al., 2015; Jin et al., 2020; Liquori et al., 2020; Janusz et al., 2021; Lin et al., 2021; Zhong et al., 2021).

MDS15 model was derived from LASSO analysis based on expression signature of brown module, and greatly improved the prediction of leukemic transformation (Figure 3). For GSE15061, only 3.09% (3/97) MDS15 low-risk patients transformed eventually, while 9.09% (6/66) of IPSS low-risk group, 24.39% (10/41) of IPSS int-1 group, progressed into AML eventually. Moreover, 57.14% (20/25) patients in MDS15 high-risk group finally transformed, superior to the result of IPSS (26.32% and 33.33% transformation rate in IPSS int-2 and high-risk group, respectively). Identical results were found in GSE58831, irrespective of MDS15 vs. IPSS (Supplementary Figure S6A) or MDS15 vs. WPSS (Supplementary Figure S6B). The demonstratable advantage in distinguishing actually risky patients in Kaplan-Meier analysis (Figure 5A, B) confirmed prognostic value of MDS15 model.

Time-dependent ROC analysis revealed 12/24/36-month AUC were 0.759/0.792/0.792 for training set, and 0.838/0.835/0.819 for testing set, respectively (Figure 6A, B). The riskscore-survival curves were steady, continuous, and well-fitted in both sets. Notably, GSE58831 and GSE15061 were not only sampled from independent cohorts, but also from different cell origins (GSE58831 derived from CD34⁺ bone marrow cells, GSE15061 from unsorted bone marrow cells). Robust and comparable predictive power was exhibited by both datasets, suggesting MDS15 model bridged between CD34⁺ sorted and unsorted samples. Since only unsorted samples were available in most hospitals or clinics, MDS15 is a promising risk system in clinical practice. In comparison with traditional prognostic scoring systems, superior predictive power (AUC and Harrison’s C index) of this model was discerned by both training and testing set (Figures 9A, B). Transformed patients had significantly higher baseline MDS15 scores (Figure 10A, B). Moreover, the independent prognostic variables were screened in covariate Cox analysis, in which diagnostic subtype, IPSS, WPSS and MDS15 scores were identified to be significant prognostic factors. While only MDS15 scores were demonstrated to be an independent prognostic factor by multi-Cox analysis. The abovementioned results elucidated the prognostic power and future potential in clinical practice.

In spite of the rapid advancement of omics techniques in the realm of myeloid neoplasms, including RNAseq and microarray, few findings have been integrated into clinical practice owing to the prohibitive costs and laboratory requirements. The MDS15 model offers precise prognostic predictions and reduces costs by incorporating fewer variables, rendering it practical for clinical applications.

Moreover, the whole-genome expression analysis profiled cell signaling pathways in relation with MDS15 scores, combining with GSEA. The transcript factors, such as MYC and E2F signaling was activated along with ascending MDS15 scores (core enrichment of which included MYC, TRAF6 and EZH2) (Figure 14). Overexpression of MYC was associated with early AML progression in MDS patients (Gajzer et al., 2021). One of the MYC regulating genes in both core enrichment and brown module, TRAF6 (TNF receptor-associated factor 6 protein) is generally overexpressed in MDS and was attributed to hematopoietic progenitor cell defects, but its expression declined and dysfunction in a subset of AML patients (Muto et al., 2022). Further investigation indicated TRAF6 mediates ubiquitination of MYC protein and plays a role of tumor suppressor in myeloid neoplasms (Muto et al., 2022), suggesting the decreasing gradient of TRAF6 expression was related to leukemic transformation risk. Consistent with this result, expression of TRAF6 was significantly negatively correlated with MDS15 risk scores (p = 0.02, Pearson’s coefficient = −0.21, Supplementary Figure S8).

EZH2 is a downstream target of E2F signaling (Bracken et al., 2003; He et al., 2019), the loss-of-function mutation of which is an independent prognostic factor for MDS. Despite of EZH2 expression declines in loss of 7q karyotype, EZH2 or SFSR2 mutant patients (Nagata and Maciejewski, 2019; Sakhdari et al., 2022), the expression profile in general MDS patients had not been investigated. After we dichotomized GSE58831 cohort into EZH2-high and EZH2-low expression subgroups by the median expression level, a trend toward LFS difference was found in Kaplan-Meier plot (log-rank p = 0089, Supplementary Figure S9A). Considering that harboring 7q abnormality (−7 or del7q) is an independent prognostic factor and lead to haplo-insufficiency of EZH2 (located in 7q36), 5 such cases were then excluded in LFS analysis, resulting a significant LFS difference between EZH2-high and EZH2-low subgroups (log-rank p = 0.0017, Supplementary Figure S6B). The expression of EZH2 was also correlated with MDS15 risk scores (p = 0.003). Abovementioned results indicated the transcriptomic signature of oncogenic signaling (e.g., MYC/E2F), associating with MDS15 scores, contributed to tumor cell proliferation and disease progression.

Notably, inflammatory signaling pathways (interferon alpha/gamma signaling, TNF alpha signaling via NF-kB, inflammatory response, etc.) were suppressed when the risk score increased (Figures 13, 14). The interplay between inflammatory cytokines and MDS stem cells was complicated and paradoxical (Hemmati et al., 2017). Both interferon alpha and gamma induce apoptosis and cell growth inhibition in hematopoietic cells (Snoeck et al., 1994; Mayer et al., 2001). TNF alpha was reported to be upregulated in bone marrow plasma and positively correlated with cell apoptosis in low-risk MDS (Kerbauy and Deeg, 2007). While inflammation signaling (e.g., TLR signaling) is specific in low-risk MDS, which promotes tumor cell death (Paracatu and Schuettpelz, 2020). These results suggested that suppression of inflammatory signaling attenuated cell death in MDS, which was in favor of leukemic transformation. MYD88 and TLR2, in the core enrichment gene list of inflammatory response signaling, was previously reported to be associated with low-risk MDS, while downregulated in high-risk MDS (Dimicoli et al., 2013; Wei et al., 2013). MYD88 is a downstream mediator following Toll-like receptor activation, the mRNA of which was reported to be over-expressed in CD34⁺ cells of low-risk MDS compared to healthy donors in previous investigation, inhibition of which restored erythroid colony formatting capacity of HSPC. Nevertheless, higher expression level of MYD88 was displayed in high-risk MDS than that in low-risk MDS (Dimicoli et al., 2013). Consistently, expression of MYD88 was significantly negatively correlated with MDS15 risk scores (p = 2.15e-3, Pearson’s coefficient = −0.28, Supplementary Figure S10), suggesting the association of MDS15 model with dysregulated immunome. The signaling analysis shed light on potential therapeutic targets for MDS transformation, paving the way for future drug development.

The primary innovation of the MDS15 lies in its enhanced prognostic accuracy compared to traditional systems like IPSS and WPSS. This enhanced precision can assist clinicians in identifying truly high-risk patients, enabling them to implement more proactive measures such as increased monitoring or interventions that can alter the disease’s natural progression, such as demethylating agents. Furthermore, by incorporating fewer variables in the MDS15, costs can be significantly reduced in contrast to comprehensive genome transcriptome analyses. Despite our study’s insights, there are limitations. The patient populations from GSE15061 and GSE58831 displayed heterogeneity in terms of age, diagnostic subtype, and various clinical and genetic variables (Mills et al., 2009; Gerstung et al., 2015). And the treatment profile was not described in both sets, which may produce potential bias. Moving forward, it will be imperative to develop a refined model that incorporates individualized patient information, and this should be validated using larger cohort studies.

By the updated transcriptomic analysis method, WGCNA, leukemic transformation correlated gene cluster was identified, and a novel prediction model (MDS15) was established by LASSO. Then, the predictive power was shown to be superior to that of traditional prognostic systems (IPSS, WPSS, etc.), which was validated by datasets derived from different cell origins and independent cohorts. Moreover, disrupted MYC/E2F signaling, and inflammatory pathway were demonstrated to be associating with MDS15 risk scores.