Machine learning integrates genomic signatures for subclassification beyond primary and secondary acute myeloid leukemia

Genetic mutations (somatic or germline), cytogenetic abnormalities, and their combinations contribute to the heterogeneity of acute myeloid leukemia (AML) phenotypes.^1-3 Seminal studies have described the molecular landscape of AML and its exclusive framework and have stratified patients into prognostic subgroups.^4-7 Moreover, serial sequencing studies have delineated a stepwise acquisition of subclonal mutations accompanying AML evolution.⁸ To date, prototypic founder lesions [eg, t(8;21), inv(16), and t(15;17)] define only a fraction of AML subgroups, with specific prognoses corresponding to molecular pathogenesis.^9,10 Indeed, in a larger proportion of AML patients, somatic mutations or cytogenetic abnormalities potentially serve as driver lesions in combination with numerous acquired secondary hits.³ However, their combinatorial complexity hampers the resolution of distinct genomic classifications and overlaps across classical pathomorphologic AML subtypes, including de novo/primary (pAML) and secondary AML (sAML) evolving from an antecedent myeloid neoplasm.^11,12 These AML subtypes are themselves nonspecific because of variation in the understanding of their pathogenetic links, especially in cases without overt dysplasia.^13,14 Without dysplasia, reliance is mainly on anamnestic clinical information that might be unavailable or cannot be correctly assigned because of a short prodromal history of antecedent myeloid neoplasm. Additionally, supervised analytic strategies to reproduce current pathomorphologic entities as gold standard using molecular features have been modest. Here, we explored the potential use of distinct genomic markers, uncovered by advanced machine learning methods, to subclassify AML objectively and provide personalized prognostication, irrespective of clinicopathologic information; we propose this method as a new standard in AML assessment.

We analyzed integrated genomic data from pAML and sAML patients seen in our institution and multiple other centers over 2 decades using both standard supervised approaches and unsupervised machine learning methods that better captured the complex interactions of high-dimensional genomic features underlying AML subgroups. Machine learning was instrumental for the identification of novel AML subgroups with invariant driver genomic features.

For the purpose of this study, we combined AML patient data from the Cleveland Clinic (CC; n = 855) and the Munich Leukemia Laboratory (n = 4002) with publicly available data sets (The Cancer Genome Atlas, the BEAT AML Master Trial, and the German-Austrian Study Group; n = 1931[ cases with unavailable cytogenetics were excluded])^4-6 to form a large, well-annotated cohort of 6788 patients (supplemental Table 1). Targeted next-generation sequencing results at the time of AML diagnosis were adjusted to focus on the most recurrent somatic myeloid mutations (supplemental Table 2). Patients were followed until September 2019, with a median duration of 12.4 months. Peripheral blood and/or bone marrow samples were collected after receiving written informed consent according to protocols approved by the institutional review board at CC and other institutions in accordance with the Declaration of Helsinki. Clinical parameters were obtained from medical records after securing appropriate material transfer agreements and from resources accessible online.

For the data collected at CC, whole-exome sequencing was performed on paired tumor and germline DNAs (purified CD3⁺ lymphocytes). Whole-exome capture was accomplished according to the SureSelect Human All Exon 50Mb or V4 Kit (Agilent Technologies), and captured targets were sequenced using a HiSeq 2000 (Illumina). Reads were aligned to the human genome (hg19) by a Burrows-Wheeler aligner (http://bio-bwa.sourceforge.net/). Data were validated using a TruSeq Custom Amplicon Kit (Illumina) with a panel of 44 genes (supplemental Table 2). Variants were annotated using Annovar and filtered, and a bioanalytic pipeline developed in house^15-17 identified somatic mutations as specified in the data supplement. Variants in patients from the Munich Leukemia Laboratory cohort were called as previously reported.^18-22 The gene sequencing methods of publicly shared AML patients were previously described.^4-6 

Multivariate Cox proportional hazards modeling was used to identify genomic abnormalities associated with survival in various AML cohorts. Uni- and multivariate logistic regression (MLR) analyses were performed to find distinct genomic features of pAML and sAML. We performed unsupervised analysis to cluster AML patients into genomic subgroups by latent variable modeling. More specifically, we used Bayesian latent class analysis (BLCA) coupled with resampling to generate a consensus matrix²³ that was then hierarchically clustered using Ward’s criteria to obtain final patient cluster assignments. To validate the prognostic significance of identified clusters, we used survival analysis. To determine if AML subtype distributions differed across identified clusters, we normalized pAML and sAML samples to population proportions using the bootstrap method. To identify distinct genomic features and generate a subclassification model, we used random forest (RF) classification and extracted the variables with the highest global importance measured by mean decrease in accuracy. Additionally, we performed internal and external validation of our model. The RF subclassification model and cluster-specific survival estimates are available via a Web-based open-access resource.

Using the World Health Organization 2016 diagnostic criteria,⁹ we classified 6788 AML patients as having core-binding factor AML (CBF-AML; n = 422), acute promyelocytic leukemia (APL; n = 312), KMT2A-rearranged AML (KMT2A^R-AML; n = 371), pAML (n = 4502), sAML (n = 832), or therapy-related AML (tAML; n = 349). Patients’ baseline, clinical, and cytogenetic information are presented in Table 1. Mutational profiling identified 13 879 somatic mutations of variant allele frequency (VAF) ≥1% in the selected uniformed gene panel (supplemental Tables 2 and 3). Using multivariate Cox modeling, we identified specific genomic lesion associations with survival. This approach enabled feature partitioning into favorable vs adverse risks within diverse AML groups (supplemental Table 4). Because the role of recurrent balanced translocations in AML diagnostics and the prognosis of tAML are already well established, we focused our analyses on 5334 pAML and sAML cases without these pathognomonic lesions; thus, we excluded CBF-AML, APL, KMT2A^R-AML, and tAML. Our objective was to determine if unique configurations of specific genetic lesions can produce distinguishable diagnostic patterns of pAML vs sAML or within AML subsets, including normal karyotype (NK-AML; n = 3176) and abnormal karyotype AML (AK-AML; n = 2158). This strategy was motivated by the observation of significantly different pAML vs sAML survival (Figure 1A-C). Indeed, the supervised analyses yielded distinct clinical (supplemental Table 5) and genomic features that characterized each subtype (Figure 1D). Patterns detected by univariate logistic regression and MLR included enrichment in mutations in CEBPA (both monoallelic and biallelic), DNMT3A, FLT3 internal tandem duplication (FLT3^ITD), FLT3 tyrosine kinase domain, GATA2, IDH1, IDH2^R140, NRAS, NPM1, and WT1 in pAML and a greater prevalence of mutations in ASXL1, RUNX1, SF3B1, SRSF2, and U2AF1; −5/del(5q); −7/del(7q); −17/del(17p); del(20q); +8; and complex karyotype in sAML (Figure 1D-F). Mutational burdens were similar in both AML subtypes (median, 2 mutations per individual; supplemental Figure 1). Analyses of NK-pAML vs NK-sAML (supplemental Figure 2) and AK-pAML vs AK-sAML (supplemental Figure 3) revealed significant genetic associations that characterized each subset. In addition, clonal hierarchy analyses differentiated pAML vs sAML based on founder and subclonal hits (Figure 1G). Despite these significant findings, the genomic profiles of pAML vs sAML identified by MLR resulted in only 0.74 cross-validation accuracy of predictive performance when used to reproduce pathomorphologic AML subtypes (Figure 1H).

Because the accuracy of MLR prediction was modest, we explored other machine learning approaches as an alternative analytic strategy. BLCA of AML cases with complete mutational screens (n = 2681; supplemental Table 2) uncovered 4 novel genomic subgroups (Figure 2A) based on the highest silhouette value (supplemental Figures 4 and 5). The biologic relevance of these subgroups was reflected in significantly different survival among GC-1 (median, 34.1; 95% confidence interval [CI], 25.2-50.5 months), GC-2 (median, 26.5; 95% CI, 22.9-31.0 months), GC-3 (median, 15.8; 95% CI, 13.3-18.0 months), and GC-4 (median, 9.2; 95% CI, 7.4-11.6 months; Figure 2B; supplemental Figure 6), as well as survival probability (supplemental Table 6). Of note, the implementation of survival analyses was considered only to reflect the biologic and prognostic relevance of these clusters and not to replace current prognostic schemes. Moreover, the robustness of the BLCA clustering with respect to VAF was further validated when considering a higher cutoff of 15%, which also resulted in 4 GCs with a silhouette value of 0.86 and an adjusted Rand index of 0.84 (supplemental Figure 7A-B).

The distribution of GCs within pAML and sAML was variable (Figure 2C). For instance, pAML cases showed similar percentages of GC-1 (32%), GC-2 (33%), and GC-3 (25%) but fewer cases of GC-4 (10%). In contrast, sAML cases had higher percentages of GC-4 (22%) but lower GC-1 (5%) (Figure 2C). The few GC-1 sAML cases may be suggestive of a possible subtype misclassification on presentation or impact of an important genetic alteration (supplemental Table 7). Higher percentages of patients with good molecular prognosis were found among those with pAML (GC-1/2; 65%), whereas those with sAML had more higher-risk cases (GC-3/4; 66%). Results of reverse analysis of normalized frequencies of pAML and sAML within cluster groups were consistent with the aforementioned results (Figure 2D), showing increased pAML proportion in GC-1 (89% vs 11%) and sAML in GC-4 (67% vs 33%). In addition, survival analyses within the same prognostic group showed no significant difference between pAML and sAML cases except in GC-4 (supplemental Figure 8A-D; supplemental Figure 9 shows for P values of all pairwise comparisons of survival in our 8-cluster × pAML/sAML groups).

To link each cluster to its pathogenetic features, we generated an RF model. The resulting multiclass classifier yielded a cross-validation accuracy of 0.97 (Figure 2E). Globally, the most important genomic features of the model, quantified by mean decrease in accuracy, included mutations in NPM1, RUNX1, ASXL1, SRSF2, and TP53; −5/del(5q); DNMT3A mutations; −17/del(17p); BCOR/L1 mutations; and others (supplemental Figure 10). Comprehensive group-specific observations showed that GC-1 was characterized by the highest prevalence of NK-AML (88%) and full presence of NPM1 mutations (100%; 86% with VAF >20%) that cooccurred with DNMT3A (52%), FLT3^ITD (27%; 91% with VAF <50%), and IDH2^R140 mutations (16%; whereas IDH2^R172 mutations were absent), with depletion or absence of ASXL1, EZH2, RUNX1, and TP53 mutations and complex cytogenetics (Figure 3A). GC-2 had a higher percentage of AK-AML cases than GC-1 and the highest frequency of CEBPA biallelic (9%) and IDH2^R172K (4%), FLT3^ITD (14%), and FLT3^TKD mutations (6%) occurring without NPM1 mutations, whereas mutations in ASXL1, RUNX1, and TP53 were absent (Figure 3B). GC-3 had the highest frequency of ASXL1 (39%), BCOR/L1 (16%), and DNMT3A mutations without NPM1 mutation (19%), in addition to being highly enriched with EZH2 (9%), RUNX1 (52%), SF3B1 (7%), SRSF2 (38%), and U2AF1 mutations (12%; Figure 3C). Of note, GC-3 showed a higher degree of heterogeneity. In fact, 53 cases in GC-3 had a silhouette value <0, and of these, 15 cases were misclassified by the RF model. Further investigation of these misclassified cases showed that they were wild type for RUNX1, a GC-3–defining gene (supplemental Figure 11). Finally, GC-4 had the highest prevalence of AK-AML (96%; mostly of complex karyotype [76%]) and TP53 mutations (70%) that were associated with −5/del(5q) (68%), −7/del(7q) (35%), and −17/del(17p) (31%; Figure 3D). Signature patterns, their importance, and pairwise co-occurrence with other genomic markers, in addition to the clonal hierarchy of driver mutations in each cluster, are described in Figure 3B-E, supplemental Figure 12, and supplemental Figure 13A-D, respectively.

We also analyzed the percentages of novel groups among each genomic lesion population (supplemental Figures 14-16). GC-1 represented 97% of NPM1, 50% of FLT3^ITD, 54% of DNMT3A, 43% of IDH1, and 43% of IDH2^R140 mutations as well as 43% of NK-AML cases; GC-2 accounted for 91% of CEBPA biallelic, 46% of GATA2, and 50% of WT1 mutations; and GC-3 had 90% of ASXL1, 82% of BCOR/L1, 52% of CBL, 53% of ETV6, and 46% of IDH2^R172K mutations. GC-3 also represented a majority of splicing factor mutations (48% of SF3B1, 86% of SRSF2, 70% of U2AF1, and 65% of ZRSR2 mutations), 98% of RUNX1 mutations, and the highest portion of del(20q) (65%) and trisomy 8 (49%). GC-4 represented 94% of TP53 mutations, 62% of complex cytogenetics, 92% of −5/del(5q), 100% of −6/del(6q), 88% of del(12p), 91% of del(16q), and 92% of −17/del(17p).

When the clinical and baseline characteristics of each group were studied (supplemental Table 8), GC-1/2 were found to comprise a significantly younger age population compared with GC-3/4 (median age, 61 vs 70 years; P < .0001; supplemental Figure 17). Moreover, lower numbers of white blood cells correlated with higher-risk disease (P < .0001; supplemental Figure 18), possibly because patients in GC-3/4 harbored more dysplastic features than those in GC-1/2, who had more a proliferative AML phenotype.

Finally, we revisited the previously excluded well-defined prognostic AML groups and applied BLCA, which demonstrated that APL, CBF-AML, and tAML constituted a single GC each, whereas 2 genomic groups were revealed in KMT2A^R-AML (supplemental Figure 19A-B), including GC-A (OS, 20.3 months) and GC-B (median OS, 6.9 months) with distinct survival (supplemental Figure 19C). The most important genomic markers extracted by the RF model included TP53 mutations followed by −5/del(5q), −7/del(7q), and −17/del(17p) (supplemental Figure 19D). GC-2 KMT2A^R-AML was characterized by the enriched presence of TP53 mutations (79%), −5/del(5q) (63%), −7/del(7q) (38%), +8 (38%), and −17/del(17p) (31%; supplemental Figure 19E), whereas GC-1 KMT2A^R-AML had absence/depletion of the aforementioned genomic aberrations except for +8 (18%), in addition to frequent DNMT3A (20%) and NRAS mutations (16%; supplemental Figure 19E).

We performed internal and external validation of our genomic clustering model. The internal confirmatory cohorts consisted of randomly selected training (80%; n = 2144) and test sets (20%; n = 537). BLCA and RF classification were applied on the training set only, and survival analysis of the test set was performed separately (supplemental Figure 20A-B). The test cohort did not show significant survival differences per each GC compared with the training set when Kaplan-Meier analyses were performed (supplemental Figure 21A-D). We further evaluated how the number of the identified clusters varied across subsets of the training cohort by randomly sampling 75% of observations to prevent overfitting. Silhouette values across random samples showed 4 as the optimal number of clusters (supplemental Figure 22). External validation was then conducted using an independent cohort of 203 AML patients from the MD Anderson Cancer Center (MDACC; pAML, n = 143; sAML, n = 60), with a median follow-up of 12 (range, 0.1-35.3 months) and fully annotated characteristics (supplemental Table 9). Gene sequencing of the selected gene panel (supplemental Table 2) identified a total of 723 somatic mutations in the MDACC cases (supplemental Table 10). The Kaplan-Meier survival analyses of the original data and MDACC cases showed similar survival among each GC assigned by the RF model (Figure 4A-D). Details of additional validation approaches for hyperparameter tuning and depth selection in the RF model are provided in the data supplement.

Furthermore, the 4 GCs greatly overlapped with the ELN 2017 risk groups (supplemental Figure 23). GC-1/2 mainly represented the favorable and intermediate groups, whereas GC-3/4 highly matched the adverse group. Interestingly, when OS of the GCs were compared within each ELN 2017 risk group, significant differences were observed (supplemental Figure 24).

Finally, our genomic subclassification model is available as a Web-based open-source resource that can be accessed widely by clinicians and the public to forecast the subclassification and estimated survival of AML patients without known pathognomonic lesions, balanced translocations, or tAML (https://drmz.shinyapps.io/local_app/; Figure 4E). A conceptual framework summarizing our overall approach is illustrated in supplemental Figure 25.

Apart from certain well-defined AML subtypes (eg, CBF-AML, APL, and KMT2A^R-AML), historically, AML patients have been subcategorized into subgroups defined by pathomorphologic features and broad anamnestic clinical criteria because of the inability to precisely infer the presence of an antecedent prodromal disease.^11,24 Ubiquitous application of genomic diagnostics has provided opportunities for objective subclassification of AML, which, because of its mechanistic foundation, can direct discovery and application of molecularly targeted therapeutics and allow for tailored personalized risk stratification.¹⁰ Building on this potential and the power of modern genomic and bioanalytic approaches, we investigated whether rational genomic tools would yield precise, simple, and diagnostic AML subclasses reflective of genomic pathogenesis and prognosis. This is best illustrated by the ability of GCs to redefine historical subclasses such as pAML and sAML, which have been shown to correlate only variably with genetic patterns and pathogenesis.

By applying machine learning methods to an unprecedentedly large cohort of AML patients with detailed molecular annotations, we developed a new classification model. Credence in this model is fortified by its accurate classification (∼97% correct) and the plausibility of the distinctive genomic features that contributed to the overall assignment accuracy. The utility of this model is supported by different cluster group survival outcomes. Feature association with both pathobiology (molecular parameters) and survival suggests that this model could overcome some limitations of previous pathomorphologically based classifications of pAML and sAML.

Our results reaffirm previous studies^4-6,11 and extend them by integrating more molecular features and expanding diagnostic and pathogenetic implications. In particular, we focused on the inclusion of genomic signatures, despite their variable degree of global importance, to achieve the highest possible genomic classification accuracy. For instance, and in line with previous reports,^10,13,25-28NPM1 and TP53 mutations contributed greatly to creating the lower- and higher-risk phenotypes (corresponding to the clinical/survival risk), respectively. However, the highest cumulative accuracy was achieved only by the incorporation of the status (presence or absence) of additional genomic lesions, including RUNX1, ASXL1, SRSF2, and DNMT3A mutations; −5/del(5q); −17/del(17p); and others.

Our molecular 4-tiered model is not meant to challenge or replace previously established prognostic schemes. It mainly focuses on objectively subclassifying genomically undefined pathomorphologic AML subtypes, including pAML and sAML. However, because we compare it with traditional prognostication tools that incorporate prognostic genomic features, like ELN 2017, we would like to point out certain advantages that were concluded from our GC-based model. Our model expands a larger pool of genomic signatures and quantifies their corresponding importance. The latter is crucial when determining the probability of objective subclassification in complex heterogeneous AML cases harboring combined ELN 2017–defined favorable and adverse genomic lesions. The model describes distinct clustering of a variety of previously described genomic lesions that are known to influence AML outcome and emerged as AML cluster determinants. DNMT3A, IDH2^R140, and TET2 mutations, only when occurring without NPM1 mutations, are important genomic determinants of GC-2. Splicing factor mutations (SRSF2 and U2AF1) contributed substantially to the performance of our model and were noticeably enriched in GC-3, indicating their predicted potential to be a distinct AML subgroup.^11,22 Moreover, RUNX1 mutations had the second-highest global importance and crucially contributed to the identification of all novel groups. Specifically, RUNX1 mutations were highly prevalent in GC-3. Thus, our data confirm the substantial presence of RUNX1-mutant AML in the most recent World Health Organization classification as a provisional disease category.^9,29 Interestingly, BCOR/L1 mutations emerged as a potential genomic marker of GC-3. Although complex karyotype AML was abundant in the poorest survival group as defined in ELN 2017, the concurrent presence of other important genomic markers identified by our model [−5/del(5q), −17/del(17p), and TP53 mutations] seemed to delineate its classification. When these aforementioned markers were absent, complex karyotype AML was also seen in other groups (GC-3). Therefore, our model argues that genomic subclassification of complex karyotype AML is strongly dependent on the present/absent status of other decisive correlating genomic markers. Finally, the model is dynamic and displays flexibility and personalization by accounting for accurate probabilities of assignment to each cluster per the presence/absence of each genomic feature and its interactions with other signatures, rather than predicting a single classification. It also defines the estimated survival interval of each genomic group, which can be considered when assessing prospective AML patients’ prognoses.

Because of the mechanistic focus, the deliberate exclusion of certain clinical data may appear as a limitation of our model. Although we believe that some of the phenotypic features are a result of genomic makeup and are likely codependent, we acknowledge that selected parameters may be later incorporated, similarly to the genomic features to be discovered in the future. The latter may, for example, include some of the germline alterations, clonal/subclonal burden, or configuration of hits as demonstrated for CEBPA biallelic mutations. Furthermore, some GCs showed a higher degree of heterogeneity compared with others, which can be likely improved by future incorporation of more complex models such as neural network–based clustering or the use of infinite priors in a Bayesian setting if a larger cohort is available. Also, as all of the patients receive therapy, new effective drugs could affect prognosis and thus may have a global subgroup-specific impact on survival, and the predictive value of survival curves within subgroups may have a limited shelf life. Although the predictive accuracy of our genomic model was validated and our approach accounted for a possible generalizability limitation by including multicenter cohorts, eventually prospective external validations of longer patient follow-up durations are still warranted. Additionally, we envision that molecularly based risk assessment may have a rational implication for the use of specific therapy choices, especially when targeted agents and their combinations become more widely applied and thus purely clinical classification schemes become obsolete to provide generalizable survival predictions.

In summary, our study demonstrates that despite the tremendous heterogeneity of AML genomics, nonrandom genomic relationships, captured by machine learning methods, are capable to accurately assign objective molecular classification and prognosis irrespective of the availability of clinicopathologic or anamnestic information. It clearly indicates that the classical distinction of sAML from pAML cannot be justified at molecular levels, and rather molecular signatures/patterns should provide a prevailing impetus for classification schemes. Our model provides a personalized genomic tool for AML subclassification that is publicly shared.

The authors thank The Cancer Genome Atlas, the BEAT AML Master Trial, and the German-Austrian Study Group for data accessibility.

This work was supported by National Heart, Lung, and Blood Institute, National Institutes of Health, grants R01HL118281, R01HL123904, R01HL132071, and R35HL135795 (J.P.M.), and by grants from the Edward P. Evans Foundation and the Henry and Marilyn Taub Foundation.

Contribution: H.A. collected, analyzed, and interpreted the clinical and molecular data and wrote the manuscript; A.D. studied and analyzed the data, applied biostatistical methods, and wrote the manuscript; C.G, A.K., J.D., J. Shreve, and Y.N. edited and provided important feedback on the manuscript; C.M.K. performed and analyzed DNA sequencing data; M.M., T. Kuzmanovic, A.S.A., F.R., H.E.C., A.N, C.H., Y.S., H.K., T. Kadia, T.H., and M.A.S. provided patients and edited and added important insights to the manuscript; T.R., J. Scott, and V.V. provided important editing and wrote the manuscript; J.P.M. designed the study, conceptualized and sponsored the overall project, and wrote the manuscript; and all authors read and approved the final manuscript.

Conflict-of-interest disclosure: T.H. and C.H. are part owners of the Munich Leukemia Laboratory. M.M. is employed by the Munich Leukemia Laboratory. Y.S. is a consultant for and has equity and royalty rights with EpiDestiny. The remaining authors declare no competing financial interests.

Correspondence: Jaroslaw P. Maciejewski, Department of Translational Hematology and Oncology Research, Taussig Cancer Institute, Cleveland Clinic Foundation, 9620 Carnegie Ave, N Building, NE6-314, Cleveland, OH 44106; e-mail: maciejj@ccf.org.