A Phenogenetic Axis that Modulates Clinical Manifestation and Predicts Treatment Outcome in Primary Myeloid Neoplasms

Although the concept of “myeloid neoplasm continuum” has long been proposed, few comparative genomics studies directly tested this hypothesis. Here we report a multi-modal data analysis of 730 consecutive newly diagnosed patients with primary myeloid neoplasm, along with 462 lymphoid neoplasm cases serving as the outgroup. Our study identified a “Pan-Myeloid Axis” along which patients, genes, and phenotypic features were all aligned in sequential order. Utilizing relational information of gene mutations along the Pan-Myeloid Axis improved prognostic accuracy for complete remission and overall survival in adult patients of de novo acute myeloid leukemia and for complete remission in adult patients of myelodysplastic syndromes with excess blasts. We submit that better understanding of the myeloid neoplasm continuum might shed light on how treatment should be tailored to individual diseases. Significance: The current criteria for disease diagnosis treat myeloid neoplasms as a group of distinct, separate diseases. This work provides genomics evidence for a “myeloid neoplasm continuum” and suggests that boundaries between myeloid neoplastic diseases are much more blurred than previously thought.


Introduction
Myeloid neoplasms comprise a compendium of rare diseases that ranges from neoplasms that are more common [acute myeloid leukemia (AML) and chronic myeloid leukemia (CML)] to neoplasms that are much rarer [myelodysplas-into the 2017 European LeukemiaNet (ELN) scheme for risk stratification by genetics (17). Note, however, that-despite supporting evidence from various reports (9,(11)(12)(13)(14)-the 2017 ELN risk stratification did not consider the mutational status of BCOR, DNMT3A, EZH2, IDH2, SF3B1, SRSF2, STAG2, U2AF1, WT1, or ZRSR2; nor did the 2017 ELN risk stratification consider partial tandem duplication of the KMT2A gene. The Revised Acute Myeloid Leukemia Composite Model (AML-CM) proposed by Sorror and colleagues in 2019 followed the 2017 ELN scheme faithfully when counting cytogenetic/molecular aberrations and did not incorporate additional gene mutations (18). Some AML cases are the results of transformation from previously diagnosed MDS, MPN, or MDS/MPN, and such patients with "secondary AML" have poorer prognosis than de novo AML (19) and often carry so-called "secondarytype" gene mutations in SRSF2, SF3B1, U2AF1, ZRSR2, ASXL1, EZH2, BCOR, or STAG2 (13). In addition, some of the de novo AML cases carry multi-lineage dysplasia that are suggestive of MDS; such "AML with myelodysplasia-related changes" (AML-MRC) cases have lower probability of complete remission and shorter survival (20,21), and furthermore their mutational spectra are statistically indistinguishable from "MDS with excess blasts" (MDS-EB) (22). It is also possible for de novo AML cases without complex karyotypes or myelodysplasiarelated changes to carry "secondary-type" gene mutations, and it has been reported that such cases have poorer prognosis (13,23). Although there is significant correlation between AML-MRC and "secondary-type" gene mutations, a high percentage of AML cases that carry "secondary-type" mutations are neither secondary AML (13) nor AML-MRC (23). In other words, either some of the patients with "de novo" AML are transiting through a prodrome that does not meet the current World Health Organization (WHO) criteria for MDS and AML-MRC, or they may be true de novo cases that represent hitherto unrecognized biological overlap between de novo AML and MDS; better understanding of clonal hematopoiesis of indeterminate potential (CHIP) might shed light on this issue in the future (24,25).
The blurred distinction between AML and MDS has received ever more attention nowadays. Traditionally, MDS cases (even high-grade MDS) were less likely to be treated with AML-type induction therapy (2,26), so it was often difficult to conduct "apple-to-apple" comparison of the clinical courses of AML and MDS. Nonetheless, a Seattle group has been treating patients with MDS-EB2 with AML-type therapy in some of their locally initiated studies (27). Recent analysis of the Seattle dataset showed that most of the prognostic distinctions between AML and MDS-EB2 disappeared when covariates such as age, performance status, 2017 ELN cytogenetic/molecular risk class, secondary-ness, induction therapy intensity, and allogeneic hematopoietic stem cell transplantation were taken into consideration (5). In other words, the 20% blast percentage threshold that is dictated by the current WHO guideline to distinguish between the two diagnoses might not be ideal, and some have even argued that reclassification of diseases along the AML-MDS continuum based on genetic abnormalities would be more clinically relevant than the current standard practice (5,13). Interestingly, one small-scale study lately reported that patients with NPM1-mutated MDS and MDS/MPN might be better treated with AML-type than MDS-type therapy (28). Recent studies that investigated fine-grained cytogenetic/molecular subtypes of MDS (29,30), MPN (31), and MDS/MPN (32) are important first steps towards precision medicine along the myeloid neoplasm continuum.
Phenotypic features have always been the bed rock for risk stratification in myeloid neoplasms, as exemplified by the AML-CM algorithm (18,33) and the Revised International Prognostic Scoring System (IPSS-R) for MDS (34).
Incorporation of cytogenetic information and gene mutations has been shown to further improve prognostic accuracy (14,30,31,35), and integration of multimodal data may enable machine learning-driven personalized prescriptions routinely someday (36). Clearly, if we were to rethink the myeloid neoplasm continuum, we would also need to characterize how genetic aberrations correlate with phenotypic features throughout the continuum, i.e., across diseases. Numerous previous single-disease studies have reported correlations between gene mutations and clinical profiles. In AML, for instance, mutations in CEBPA, CSF3R, FLT3, GATA2, IDH, KIT, NPM1, and WT1 have been shown to be associated with high white blood cell count or high blast count (6,(37)(38)(39)(40)(41)(42)(43). In MDS, previous studies have reported that SF3B1 and DNMT3A mutations are associated with high platelet count (44,45), that GATA2 mutations are associated with monocytosis (46), and that mutations in RUNX1, TP53, and NRAS are associated with severe thrombocytopenia and higher blast count (47). It has also been reported that cell morphology profiles and gene mutation patterns are correlated with each other in MDS, MDS/MPN, and secondary AML (48).
In this study, we sought to understand the distribution of genetic aberrations across primary myeloid neoplasms as well as how they correlate with clinical manifestations. The patient cohort analyzed in this study is hereafter referred to as "THAMP"' (Tianjin Hematologic Atlas of Myeloid neoplasm Phenogenetics). THAMP included 730 consecutive newly-diagnosed cases of primary  Table S2). The outgroup was included in our analysis so as to distinguish between myeloid neoplasmspecific and lymphoid neoplasm-specific genetic aberrations. All the patients in THAMP were karyotyped and sequenced at a panel of 171 genes ("THAMP-Illumina-171"; Supplementary Table S3) before receiving therapy. We integrated the multi-modal data from the patients of primary myeloid neoplasms that encompassed gene mutations, karyotypic abnormalities, and peripheral blood and bone marrow parameters to identify higher-order patterns that had clinical implications in target genes and treatment outcome.

Informed Consents and Ethics Approvals
This work was entirely retrospective. At admission to the IHCAMS, each of the patients included in this study signed an informed consent form for peripheral blood sample usage and another informed consent form for bone marrow sample usage; in both of the signed forms, the patient explicitly authorized the physicians at the IHCAMS to "utilize the test samples for reasonable medical usage (including medical research)". If the primary care oncologist of a patient had not ordered a DNA sequencing test for the patient, then this patient was eliminated from this study; that is, there was no effort made during this study to collect any additional genotypic data beyond what had been deemed necessary for routine clinical care.

Variant Calling after NGS
Sequencing artifacts or low-quality bases were discarded after quality Germline samples were not available for most of the patients in THAMP. Somatic mutations were identified through the following pipeline: step 1, reject variants whose estimated probability of the base call being wrong was >0.01 (Q score < 20); step 2, reject variants where total sequencing depth (including wildtype and other variants) was <30 at the locus; step 3, reject variants whose variant allele fraction (VAF) was <0.02 or >0. 45 For the patients who carried mutations at CEBPA, those who carried more than one variants were considered "doubly mutated" (CEBPA-dm); the rest of the cases were considered "singly mutated" (CEBPA-sm).
FLT3-ITD was detected separately using a customized version of Genomon IT-Detector (53). All positive cases underwent DNA fragment analysis, and then the allele ratio (AR; defined as the ratio of the area under the curve "FLT3-ITD" divided by area under the curve "FLT3-wild type") was calculated. When AR was ≥0.5, the variant was classified as FLT3-ITD high ; otherwise, it was classified as FLT3-ITD low .

Experimental Validation of Variant Calls
Called somatic mutations with VAF ≥ 0. 15

Dirichlet Process-Driven Clustering of Genetic Aberrations
Unlike previous single-disease studies that utilized Dirichlet process to cluster patients based on their genetic profiles (14,30), we sought to cluster genetic aberrations based on their distributions across hematologic neoplasms. Suppose there were p genetic aberrations of interest and q diseases of interest.
Suppose there were n j patients who were diagnosed with disease j, out of which m ij patients carried genetic aberration i. Suppose every gene in cluster k had true (unobserved) prevalence θ jk in disease j. Let π k be the number of genetic aberrations in cluster k. Let C i be the cluster label for genetic aberration i. We set up the following model: This model was then fitted via Markov chain Monte Carlo, i.e., reiterating between the following two steps: 1. Update C i : Suppose the latest value for C i was C (t ) i . Sample a new value for C i according to Repeat this procedure for all i.
2. Update θ jk : Sample a new value for θ jk according to Repeat this procedure for all j, k.
We set the hyperparameter α = 0.0001. Model-fitting was initiated with each genetic aberration being its own standalone cluster. The burn-in period was 5,000 cycles, followed by 500 cycles of sampling from the posterior distribution. At the end of computation, each genetic aberration was assigned a cluster label that had the highest posterior probability.

Two-Dimensional Embedding of the THAMP Patient Cohort
We trained a three-layered fully-connected neural network [input layer: 339 units; layer 1: 100 units with "rectified linear unit" (ReLU) activation; layer 2: 65 units with ReLU activation; layer 3: eight output units with softmax activation) using the "keras" application programming interface for deep learning to transform the raw phenogenetic data-including the mutational status at the 171 genes [174 variables (CEBPA mutations were classified into "CEBPA-sm" and "CEBPA-dm", and FLT3 mutations were classified into "FLT3-ITD", "FLT3-TKD (tyrosine kinase domain)", and "FLT3-Other")], presence/absence of the 72 chromosomal abnormalities, age, complete blood counts (36 variables), and peripheral blood/bone marrow aspiration smear analysis (56 variables)-into discrete diagnosis for each patient. After fitting the neural network, the output values of layer 2 were used as the input for principal component analysis (PCA) to embed each patient as a dot onto a two-dimensional phase space.

Benchmarking THAMP Against AACR GENIE 10.0-Public
Our gene panel "THAMP-Illumina-171" was a commercial panel (KingMed Diagnostics) routinely used for genotyping patients at initial diagnosis and relapse at the IHCAMS cancer clinics since February 2020. The panel originally contained 175 genes and was designed to provide 100% coverage of genes used in standard clinical practice of genetic risk stratification for hematologic neoplasms, along with an expanded set of additional genes selected based on literature search. Four of the genes in the panel (ARID1B, CRLF2, FBXO11, and KMT2C) did not pass quality control and were excluded from analysis; therefore, the final result contained mutational information regarding 171 genes (Supplementary Table S3). In our quality control experiments, 99.3% of the sequence variant calls in the THAMP dataset could be validated with Sanger sequencing (Supplementary Table S4).
We then benchmarked THAMP against AACR GENIE 10.0-public (54), an extensive database of gene mutations in cancer patients accessible at the AACR website (https://www.aacr.org/professionals/research/aacr-project-genie/aacrproject-genie-data/). Despite the overall high concordance between THAMP and AACR GENIE 10.0-public, there were notable differences: First, CML had the lowest concordance between the two datasets ( Fig. 1E). For example, while the mutation rate of TET2 (whose mutation was known to be more associated with advancedphase CML than chronic-phase CML; refs. 56, 57) was 3.7% in CML in THAMP, it was 26.3% in CML in AACR GENIE 10.0-public (Fig. 1F). This could be due to difference in the timing of genotyping; while THAMP only included newlydiagnosed cases before treatment, the AACR GENIE 10.0-public dataset did not provide information related to disease phase. Second, while the mutation rate of NRAS in AML in THAMP (22.9%) was comparable to the numbers reported in other studies (58)(59)(60), it was considerably different from the observed mutation rate in AML in AACR GENIE 10.0-public (11.2%) (Fig. 1F).

Clustering Genetic Aberrations According to their Prevalence Spectra
In THAMP, the median number of observed gene mutations in a patient was 3, We observed notable differences in gene mutation rates between the pediatric (age ≤16) and adult (age >16) cases in THAMP: For AML, FLT3-ITD and mutations in DNMT3A, NPM1, TET2, CEBPA, IDH2, FLT3-TKD, RUNX1, etc.
were more prevalent in the adult cases, while non-ITD and non-TKD FLT3 mutations ("FLT3-Other") and MPL mutations were more prevalent in the pediatric cases (Fig. 2B)
A Dirichlet process-based Markov chain Monte Carlo algorithm was applied to THAMP and classified the 171 genes into 12 clusters according to their mutational spectra (Fig. 3A). Gene mutation clusters G01, G02, G04, G08, and-to a lesser extent-G09 were relatively prevalent in both myeloid and lymphoid neoplasms. In contrast, the other clusters were skewed towards one of the two lineages: mutations in clusters G06 and G10 were in general more associated with lymphoid neoplasms than myeloid neoplasms; on the other hand, mutations in clusters G03, G05, and G07 were more associated with myeloid neoplasms than lymphoid neoplasms.
Clustering of karyotypical abnormalities was also conducted, resulting in seven clusters of chromosomal abnormalities (Fig. 3B). As expected, cluster K01 [which contained one sole abnormality t(9;22)(q34;q11), the cytogenetic signature for the BCR/ABL1 fusion gene] was primarily associated with CML and to a lesser extent ALL. Clusters K03 and K06 (which contained many of the cytogenetic abnormalities listed in the 2017 ELN risk stratification scheme for AML) was associated with AML primarily. Likewise, cluster K04 (which contained many of the cytogenetic abnormalities listed in the IPSS-R risk stratification scheme for MDS) was associated with MDS primarily. Cluster K05 [which contained the cytogenetic signature t(12;21)(p13;q22) for the TEL/AML1 fusion gene, a common genetic aberration in ALL] was associated with ALL primarily.
Finally, cluster K02 (which contained one sole abnormality +8) was associated with ALL and most of the myeloid neoplasms.
We calculated Pearson correlation coefficients between the gene mutation clusters and the karyotypical abnormality clusters across the eight diseases (Fig.  3C). The magnitude of correlations was small (ranging from −0.12 to 0.15), but there were distinct patterns. For instance: gene mutation clusters G04 and G08 (associated with both myeloid and lymphoid neoplasms) and karyotypical abnormality cluster K02 (associated with both myeloid and lymphoid neoplasms) were positively correlated; gene mutation cluster G07 (associated more with myeloid neoplasms) was positively correlated with karyotypical abnormality cluster K03 (associated more with myeloid neoplasms); finally, gene mutation clusters G01 and G09 (associated with both myeloid and lymphoid neoplasms) was positively correlated with karyotypical abnormality cluster K03 (associated more with myeloid neoplasms) and K05 (associated more with lymphoid neoplasms).

Correlations between Gene Mutations and Phenotypic Features in Myeloid Neoplasms
Since we were interested in characterizing the phenogenetic patterns along the entire myeloid neoplasm continuum in this study, we also calculated Pearson correlation coefficients between gene mutations and a wide range of phenotypic features across AML, MDS, MDS/MPN, and MPN (Fig. 4A).
Interestingly, we found that there was a prominent gradient of phenogenetic correlations between two extremes of clinical manifestations. Reading Fig. 4A from top to bottom: blast counts (in both bone marrow and peripheral blood) and white blood cell count (peripheral blood) were prominently positively correlated with the total number of mutations in cluster G07 and (to a lesser extent) clusters G09 and G01-consistent with previously reported associations in CEBPA, CSF3R, FLT3, GATA2, IDH, KIT, NPM1, and WT1 (6,(37)(38)(39)(40)(41)(42)(43). This was followed by lymphocyte count (peripheral blood), monocyte percentage (peripheral blood), and promonocyte percentage (bone marrow), which were the most positively correlated with gene mutation clusters G07, G09, G01, and G02; previously, GATA2 mutation (classified as G07) has been reported to be correlated with monocytosis (46). Then, lymphocyte percentage (peripheral blood) and basophilic myelocyte percentage (peripheral blood) were positively correlated with clusters G06 and G08. Orthochromatic erythroblast percentage (bone marrow), erythrocytic cell percentage (bone marrow), immature platelet fraction (peripheral blood), etc. were correlated with cluster G05; previously, SF3B1 mutation (classified as G05) has been reported to be positively correlated with high platelet count (44). Finally, neutrophil percentage (peripheral blood), platelet count (peripheral blood), eosinophil percentage (peripheral blood), band neutrophil percentage (bone marrow), etc. were the most positively correlated with gene mutations in cluster G03. We also observed prominent negative correlations between platelet count (peripheral blood) and mutations in clusters G07, G09, G01, and G06; previously, mutations in TP53 (classified as G06) and NRAS (classified as G01) have been reported to be associated with severe thrombocytopenia (47).
We also investigated if the global-scale correlation between clinical manifestations and gene mutations was spurious due to Simpson Paradox, wherein an   observed global pattern might disappear or even become reversed when analysis was limited to a subgroup of patients (63). If Simpson Paradox was in effect, we anticipated that the phenogenetic correlation we observed on the global scale in Fig. 4A would cease to stand when we focused on a subset of the patients who were diagnosed with the same disease.
Keeping the same order of phenotypic features and the same order of gene mutation clusters as in Fig. 4A, we found that even when we limited our correlation analysis to AML only, phenotypic features at the top tended to be more correlated with gene mutation clusters on the right, while phenotypic features at the bottom tended to be more correlated with gene mutation clusters on the left (Fig. 4B). Similar phenomenon could also be observed in MDS (Fig. 4C). Therefore, global-scale correlations between phenotypic features and gene mutations were also reflected locally, i.e., in subgroups of patients.

Computation of the Pan-Myeloid Axis through Integration of Multidimensional Phenogenetic Data
To further understand how genetic and phenotypic aberrations jointly related to disease diagnosis, an artificial neural network (ANN; Fig. 5A Supplementary Fig. S6. multidimensional/multimodal clinical and genetic profiles) were then used to project all the patients onto a two-dimensional phase space via PCA to visualize how the 1,192 patients related to each other phenogenetically ( Fig. 5A and B).
We found that the myeloid and lymphoid cases in THAMP were distinctly grouped into two diverging axes, which we termed the "Pan-Myeloid Axis" and the "Pan-Lymphoid Axis" (Fig. 5B). Interestingly, CML formed a distinct group by itself, apart from the Pan-Myeloid and Pan-Lymphoid axes; this was consistent with the long-known fact that CML-despite being a myeloid neoplasm-in blast crisis can sometimes transform to ALL, a lymphoid cancer. When we projected onto the two-dimensional phase space the center of probabilistic mass of each gene mutation and each karyotypical abnormality (Fig. 5C), we found that-consistent with Fig. 3-the Pan-Myeloid and the Pan-Lymphoid axes each had its associated prevalent genetic aberrations, and they appeared to lie in particular orders along the two axes.

Simplifying Gene Mutation Groupings for the Pan-Myeloid Axis
The gene mutation clustering reported in Fig. 3A was calculated using all the 1,192 patients. Since we were primarily interested in myeloid neoplasms in this study, in light of the new information generated by the one-dimensional rendering of myeloid neoplasms (Fig. 5), we further simplified gene mutation clustering for myeloid neoplasms by designating cluster G03 as "myeloid.gene.L"; grouping G05, G12, G04, and G08 genes that had >1% mutation rate in the myeloid neoplasms in THAMP into "myeloid.gene.ML"; grouping G06 and G02 genes that had >1% mutation rate in the myeloid neoplasms into "myeloid.gene.MR"; and grouping G01, G09, and G07 genes that had >1% mutation rate in the myeloid neoplasms into "myeloid.gene.R" (Fig.   6A). All the other genes were designated as "other genes".
Using bootstrapping experiments, we confirmed that the orders of genes, phenotypic features, and patients along the Pan-Myeloid Axis were stable with statistical significance (Fig. 6B). In other words, the patterns we reported in Fig. 5D and E, and Fig. 6A were not dependent on one particular composition of patients and therefore less likely to be spurious discoveries.
Note that the ranking order of cancer genes along the Pan-Myeloid Axis was not a naïve sorting of genes based on their mutation rates in myeloid neoplasms (Fig. 6C). While the average mutation rate of genes in AML indeed formed a decrescendo from the right end of the Pan-Myeloid Axis to its left end ("myeloid.gene.R", 10.1% in AML; "myeloid.gene.MR", 5.8% in AML; "myeloid.gene.ML", 2.6% in AML; "myeloid.gene.L", 1.0% in AML; "other genes", 0.3% in AML), mutation rates of individual genes were a lot more scrambled. For instance, CSF3R, IKZF1, and RAD21 were in group "myeloid.gene.R", and yet their mutations rates in AML were lower than the average of "myeloid.gene.MR". Similarly, PTPN11 and RUNX1-both in "myeloid.gene.ML"-had higher mutation rates in AML than the average of "myeloid.gene.MR".

Clinical Relevance of the Pan-Myeloid Axis in Adult De Novo AML
We subsequently asked if the Pan-Myeloid Axis we uncovered in Figs. 5 and 6 had any relevance to treatment response, even though no clinical outcome data had been used in any computation so far. The ideal experiment to truly answer how clinical course varies along the Pan-Myeloid Axis would be administering the same array of treatment regimens to more than one disease along the axis-similar to what was recently proposed by Estey and colleagues (5). Only through such clinical trials can one determine which patients should be treated with "AML-type" regimens, which patients should be treated with "MDS-type" regimens, etc. In the real world, however, most patients are treated according to their diagnoses based on the contemporary diagnostic criteria; thus, patients with "AML" are usually treated with "AML-type" regimens, and patients with "MDS" are usually treated with "MDS-type" regimens. In other words, the clinical course of a patient cannot be evaluated independently of diagnostic criteria. Still, it would have been a missed opportunity to not investigate if the Pan-Myeloid Axis could also help with risk stratification in a single disease. Our reasoning was that some of the patients might not respond well to standard treatment regimens because they deviated from the canonical cases, and one evidence for this could be that they were situated at some distance from the canonical cases on the Pan-Myeloid Axis.
We first examined the AML cases in THAMP (Supplementary Table S9). To eliminate confounding influence from age and the PML-RARA fusion gene (which has favorable prognosis and is usually treated with a chemotherapy regimen that is entirely different from other AML subtypes; refs. 64, 65), we limited clinical outcome analysis to the adult (age >16) de novo cases that were not M3.
Since the patients with AML in THAMP were all relatively recent, our analysis focused on MRD after induction therapy.
We recognized that even these relatively homogeneous adult de novo non-M3 AML cases could yet still be more heterogeneous than ideal. For instance, we were aware of at least two additional confounding factors that could interfere with our analysis: First, 9 (4%) of the remaining 204 de novo AML cases carried myelodysplastic-related changes ("AML-MRC"); they could be secondary AML cases in reality, but somehow their MDS prodromes had escaped diagnosis (21). Second, 35 (17%) additional cases carried gene mutations (SRSF2, SF3B1, U2AF1, ZRSR2, ASXL1, EZH2, BCOR, or STAG2) that were commonly associated with secondary AML ("secondary-type AML"; ref. 13), even though they did not carry myelodysplastic-related changes according to the current WHO criteria. Close examination of these two subgroups of the patients showed that their gene mutational spectra and phenotypic features were indeed substantially different from the other adult de novo non-M3 AML ("bona fide") cases (Fig. 7A).
To validate the proposed "2017 ELN enhanced" scoring system and also to further investigate the prognostic significance of the Pan-Myeloid Axis, we queried cBioPortal in search of de novo AML patient sequencing data that included survival information. There was one hit-referred to as the "Beat-AML-2018" dataset (55) in this study-that contained 140 adult (age >16) patients with non-M3 de novo AML who had clearly defined 2017 ELN genetic risk scores, had disease samples sequenced at initial diagnosis, were not recorded as AML-MRC, did not carry "secondary-type" gene mutations, were treated with standard chemotherapy, and had known overall survival (OS) status.
Applying the same "2017 ELN enhanced" risk scoring method described in Fig.   7C, we found that the three subgroups of patients with de novo AML had mutually distinct survival curves that passed statistical significance, with much better intersubgroup separations than 2017 ELN risk stratification (Fig. 7F). In addition, the concordance score for predicting OS improved from 0.62 (2017 ELN) to 0.65 (2017 ELN enhanced).
Therefore, although the Pan-Myeloid Axis was discovered through a crossdisease analysis that was meant to quantitatively characterize the transitions between the different diseases, this axis nonetheless carried statistically significant information related to how the patients with "de novo AML" who were diagnosed with the current diagnostic criteria responded to standard "AML-type" treatment.

Clinical Relevance of the Pan-Myeloid Axis in Adult MDS-EB
Finally, we examined the MDS cases in THAMP (Supplementary Table S10) to investigate if there was also evidence for correlation between the Pan-Myeloid Axis and treatment response in MDS. 42 adult patients (age >16) of MDS with excess blasts (MDS-EB) in THAMP underwent chemotherapy that included both hypomethylating agents and cytotoxic drugs. CR (defined as ≤5% myeloblasts in bone marrow nucleated cells and decrease by ≥50% over pretreatment) after six cycles of chemotherapy was used as the criterion for treatment response. We found that high mutation load in "L + ML" (defined as carrying >2 mutations in "myeloid.gene.L" and "myeloid.gene.ML" combined) was significantly correlated with CR in MDS-EB ( Fig. 8A and B). In particular, AUROC increased from 0.42 (based on IPSS-R) to 0.67 (based on the total mutation load in "L + ML") (Fig. 8C). Thus, the Pan-Myeloid Axis carried sta-tistically significant information related to how the patients with "MDS-EB" who were diagnosed with the current diagnostic criteria responded to standard "MDS-EB-type" treatment.

Discussion
We acknowledge several limitations in our study: First, the THAMP cohort only contained patients that were treated at a single medical center; therefore, conclusions of this study might not extrapolate to other patient groups. Second, our study cohort was composed of primarily patients of Han ethnicity (selfidentified ethnicities: Han, n = 1,105; Manchu, n = 37; Mongol, n = 26; Hui, n = 16; Korean, n = 3; Kazakh, n = 2; Uyghur, n = 1; Tibetan, n = 1; unidentified, n = 1), and our observed spectra of genetic aberrations might be specific to this particular composition of racial backgrounds, as interracial differences have been reported previously (66)(67)(68). Third, our study of mutational spectra was based on targeted sequencing and might have missed crucial driver mutations or concomitant mutations. Fourth, the patients with AML and MDS in THAMP were primarily treated with "AML-type" and "MDS-type" regimens, respectively, and therefore we were limited to examining the clinical impact of the computed Pan-Myeloid Axis "locally" (i.e., examining how the axis modulates AML and MDS patients' response to standard AML and MDS therapies, respectively), rather than "globally" (i.e., examining how the axis modulates the responses of a continuum of patients with myeloid neoplasm to an array of commonlyadministered therapy regimens; refs. 5,22,28). Fifth, our data did not allow us to examine the temporal order of mutations, which had been reported to be influential on phenotype, disease course, and treatment outcome in myeloid neoplasms (69)(70)(71)(72). Sixth, we did not have preclinical molecular data for the THAMP patients; as such, we did not investigate how CHIP affected these patients' trajectories of clinical manifestations and disease progress, a subject that has received more attention lately (24,25,73). Seventh, the THAMP cohort was very recent, and our study was limited by its relatively short follow-up period.
Nonetheless, our study had the advantage that all the patients were subjected to the same protocols for genotypic and phenotypic data collection. Our comparative genomics analysis provided strong evidence for the existence of a phenogenetically-defined "Pan-Myeloid Axis", by which commonly mutated genes in myeloid neoplasms were classified into right ("AML-ish") versus left values of the total number of mutations in "R + MR", "L + ML", and "others" in the adult MDS-EB patients were 1, 2, and 0, respectively. B, CR rates for the patients with MDS-EB stratified according to the total mutation load in "L + ML" ["High" (>2 mutations), n = 13; "Low" (≤2 mutations), n = 29]. P values were calculated using the Fisher exact test (two-sided). C, AUROCs for predicting CR in THAMP. Only one trial was conducted for each stratification model; therefore, there was no SE or CI. *, P < 0.05; **, P < 0.01.
We observed some correlations between the Pan-Myeloid Axis and gene function. For instance, numerous leukemic genes involved in receptor tyrosine kinase pathways such as FLT3, KIT, KRAS, and NRAS (all classified as "myeloid.gene.R") were situated at the far right of the Pan-Myeloid Axis, while preleukemic genes such as TET2 (classified as "myeloid.gene.MR"), DNMT3A (classified as "myeloid.gene.MR"), and ASXL1 (classified as "myeloid.gene.ML") that were associated with CHIP (79-82) were all placed more to the left. In addition, all the gene mutations that have been previously reported to be associated with secondary AML (SRSF2, SF3B1, U2AF1, ZRSR2, ASXL1, EZH2, BCOR, or STAG2; ref. 13) were all classified as "myeloid.gene.ML" and situated at the left side of the Pan-Myeloid Axis. Note, however, that the "2017 ELN enhanced" risk stratification scheme continued to perform substantially better than the 2017 ELN scheme even when all the AML-MRC and "secondary-type" AML cases were excluded from our clinical outcome analysis (to be more precise, "2017 ELN enhanced" exhibited even stronger advantage over "2017 ELN" after we eliminated the AML-MRC and "secondary-type" AML cases). Therefore, the Pan-Myeloid Axis was not a mere recapitulation of the molecular distinction between de novo AML and secondary AML.
There were apparent discrepancies between the THAMP-defined Pan-Myeloid Axis and the 2017 ELN genetic risk stratification scheme. Most notably, while NPM1 mutation was labeled as favorable in 2017 ELN (especially when FLT3-ITD allelic ratio was low), NPM1 mutation was classified as part of "myeloid.gene.R"-a poor-prognosis cluster for AML-according to the Pan-Myeloid Axis. We, however, were not alone in finding this inconsistency regarding NPM1; previous researchers have published clinical outcome analyses that contradicted 2017 ELN with respect to the role of NPM1 mutation in prognosis (83)(84)(85). Similarly, CEBPA biallelic mutation was designated as favorable in 2017 ELN, but CEBPA mutation was classified as "myeloid.gene.R" in our study. The "2017 ELN enhanced" scheme reclassified some of the 2017 ELN-designated "favorable" patients in THAMP as intermediate-risk ("group 2") because of their high total mutation load (>2) in "R + MR", i.e., they had concomitant mutations in the other "high-risk" genes. It has been reported that not all CEBPA mutations carry favorable prognosis for AML (42); moreover, concomitant mutation in WT1 (classified as "myeloid.gene.R") has been reported to predict poor prognosis in AML even when CEBPA biallelic mutations were present (86). Taken together, our pan-myeloid analysis provided comparative genomics evidence that primary myeloid neoplasms indeed form a continuum. In addition, we showed that variation of treatment response in individual patients is influenced by how the diseases are organized along the phenogenetically-defined Pan-Myeloid Axis. With the ever growing trend of redefining hematologic neoplasm entities in molecular terms (21,96), diagnostic boundaries between the diseases and decision trees for treatment planning will continue to evolve and change shape in light of new genomics evidence. The ultimate test for the concept of "myeloid neoplasm continuum" would be systematic testing of successively graded and measured customization of treatment regimens along the Pan-Myeloid Axis.