Mutation Patterns Predict Drug Sensitivity in Acute Myeloid Leukemia

Acute myeloid leukemia (AML) is a hematologic malignancy characterized by maturation arrest and uncontrolled proliferation of myeloid progenitor cells (1). Data from several large cohorts of patients with AML have been analyzed to understand the mutational landscape and how the genetic diversity defines the pathophysiology of AML (2–5). According to bulk molecular profiling, the acquisition of mutations in leukemogenesis follows a stepwise pattern, where mutants with high variant allele frequencies emerge early, while mutations with lower variant allele frequencies are thought to occur later (6–8). Single-cell mutational profiling showed AML is dominated by a small number of clones, which frequently harbor co-occurring mutations in epigenetic regulators (9). Mutations in signaling genes often occur in distinct subclones from the same patient (9). The complexity of malignant cell evolution and presence of heterogeneous subclones make it challenging to stratify patients and optimize treatment.

There are only a limited number of FDA-approved targeted inhibitors for AML, such as the FLT3 inhibitors midostaurin, gilteritinib, and quizartinib, the IDH1 inhibitors ivosidenib and olutasidenib, the IDH2 inhibitor enasidenib, tyrosine kinase inhibitors (TKI) for ABL or KIT mutations, and the BCL2 inhibitor venetoclax. There are several newer agents under study, but the overall number of such inhibitors remains relatively small compared with the known number of mutated genes. Novel methods for inhibiting mutated genes, such as targeted protein degradation (10) and antibodies to drug–peptide complexes (11), appear promising, but remain under development. Thus, the path to optimize individual treatment remains complex with inadequate options.

Different factors may affect drug sensitivity, such as the genetic/epigenetic alterations, subclonal evolution, and the bone marrow stromal and immune microenvironment (12–15). Nevertheless, the majority of patients lack the genetic mutations targeted by the approved drugs. Co-occurrence of gene mutations has been widely reported in many cancer types, including AML (3). It has been reported that leukemia stemness and co-occurring mutations drive resistance to isocitrate dehydrogenase (IDH) inhibitors in AML (16). The co-occurrence of gene mutations can be due to clonal evolution with acquisition of additional mutations in the same cell or to distinct new clones with other mutations (9), which may require different treatments or combinations to optimize response.

Several groups have performed ex vivo drug screening for AML patient-derived samples, including the Beat AML project (5) and the functional precision systems medicine study (17). These datasets allow us to explore the genetic heterogeneity of AML cells and predict drug sensitivity for patients with given genetic alterations. However, more targeted inhibitors are needed to treat patients with AML with clear genetic alteration patterns for which there are still no specific agents.

In this study, we sought to characterize the heterogeneous genotypes of patients with AML with different co-occurring or mutually exclusive mutation patterns and to associate such patterns with drug response to develop clinical prediction models. Using the abundant public data resources for AML, we have built co-mutation graphs and identified subgraphs that allow us to understand the potential functional impact of different mutation patterns and predict drug sensitivity. We also characterized the mutational landscape of patients with AML using targeted sequencing of 194 AML-related genes and carried out ex vivo drug screening using a library of 207 drugs for patients with AML. We identified potential drugs for patients with various genetic alterations. We also performed single-cell genomic profiling to confirm the co-occurrence of mutations at the single-cell level and constructed machine learning models for the prediction of drug sensitivity using the genetic variants. Our newly generated targeted sequencing, ex vivo drug screening, and single-cell genomic profiling datasets, together with the accompanying analysis, provide options for more rational and potentially more effective selection of therapies.

Patients were enrolled with written informed consent on a protocol approved by the University of Washington/Fred Hutchinson Cancer Center Cancer Consortium Institutional Review Board. The studies were conducted in accordance with the U.S. Common Rule and the Belmont Report. Samples from blood and/or bone marrow were collected from 99 patients with AML. Deidentified clinical data were provided for data analysis, and deidentified samples underwent laboratory testing. Targeted gene sequencing and ex vivo drug screening were performed to investigate the drug sensitivity or resistance-associated biomarkers and to build the prediction models. Integrative analysis from large datasets was carried out to characterize the mutation co-occurrence patterns, which we further utilized as features for the drug sensitivity prediction models.

Mononuclear cells were isolated from AML blood and marrow samples by density gradient centrifugation using lymphocyte separation medium. Blasts were enriched to >80% if needed by either positive selection for CD34⁺ AML, or depletion of non-myeloid populations (e.g., CD3⁺ lymphocytes and CD235a⁺ erythroid progenitors) when CD34⁻, using magnetic beads (Miltenyi Biotec). The final blast percentage was determined by multicolor flow cytometry.

The cell populations were analyzed after a 72-hour exposure to 8–12 customized drug concentrations (within the range of 5 pmol/L to 100 μmol/L) of each drug spanning 4–5 logs. After exposure viability was determined using CellTiter-Glo luminescent reagent (Promega) per manufacturer's protocol, then the plates were analyzed with the EnVision Multilabel plate reader (Perkin Elmer). XLFit (IDBS), a Microsoft Excel Add-in, was used to analyze the data and generate dose–response curves based on standard 4-parameter logistic fit [i.e., fit = (A + (B/(1 + ((x/C)^D)))) where A and B equal minimum and maximum asymptotes, C equals IC₅₀ and D equals slope]. The AUC values were calculated using the XLFit software utilizing minimum and maximum concentrations of drugs/compounds within the panel as the limits of the AUC calculation. For each plate, data were normalized to DMSO 100% viability and blank controls. The samples were cryopreserved for the first 30 samples used to establish the assay, then all the rest of the samples were fresh. Clinical Laboratory Improvement Amendments standards require that thawed cell line data be repeated for standardization of the assay every 6 months. The assay results were highly reproducible.

An Illumina DNA prep kit was used to isolate 0.5–1 μg of genomic DNA from the patient samples enriched for leukemic blasts and the purified DNA was submitted to Invivoscribe, Inc. for sequencing. The method described on their website is briefly summarized as follows: The MyAML platform was used to perform next-generation sequencing (NGS) to analyze the 3′ and 5′ untranslated region and exonic regions of 194 genes (at the time of these analyses) and somatic gene fusion breakpoints known to be associated with AML. Fragmented genomic DNA (∼3.4 Mb) is captured with a customized probe design, and sequenced with 300 bp paired end reads on an Illumina MiSeq instrument to an average depth of coverage >1,000×. Using a custom bioinformatics pipeline, MyInformatics, single-nucleotide variants (SNV), insertion/deletions (indel), inversions, and translocations are identified, annotated, characterized, and allelic frequencies calculated. The commonly associated variants in dbSNP (18) and 1000 Genomes (19) were eliminated.

The most commonly mutated genes in AML are targeted with an average depth of coverage of 975× (range = 417× to 1,370×). The overall analytic sensitivity is >96% and specificity >99.99% for SNVs, and >95% and >99.98% for indels. There is 100% sensitivity and 100% specificity for known pathogenic mutations, including missense and nonsense mutations in FLT3, DNMT3A, IDH1, IDH2, KIT, NRAS, KRAS, and TP53. The method detects SNVs and indels with allelic frequencies as low as 2.5% with >95% reproducibility.

We utilized the MissionBio Tapestri Platform and manufacturer provided protocols (Tapestri Single-Cell DNA + Protein Sequencing User Guide) with the MissionBio 45-gene myeloid panel for targeted sequencing, with 312 amplicons, covering genes in myeloid disorders including AML, myelodysplastic syndrome, myeloproliferative neoplasms, and chronic myelomonocytic leukemia. Cryopreserved mononuclear cell samples from primary patient blood or bone marrow were thawed, incubated in Iscove's modified Dulbecco's medium with 15% FCS and 15% horse serum to enhance viability, then subjected to mononuclear cell isolation using lymphocyte separation media. Samples with less than 95% viable cells also underwent enrichment using the dead cell removal kit (Miltenyi Biotec). The viable cell fraction was verified by flow cytometry using viability dyes. The cells were incubated with Trustain FcX, then the TotalSeq-D Human Heme Oncology Cocktail, V1.0 (BioLegend) panel of 42 antibodies to surface markers. A total volume of 35 mcl of 3–4,000 cells/mcl was loaded onto the Tapestri microfluidics DNA cartridge. Single cells were encapsulated, lysed, and barcoded. DNA and protein products were purified with AMPure XP beads and DNA and protein libraries prepared. DNA analysis and sequencing was performed by the University of California Irvine High Throughput DNA Sequencing Core (NCI Shared Resource) on an Agilent Bioanalyzer and Illumina NovaSeq 6000. A total of 294 million to 893 million reads were obtained for each sample. The sequences were analyzed using Tapestri Insights software (MissionBio) to provide the characterization of the clonal mutations and distribution. The phenotype was analyzed using Mosaic software (MissionBio).

Initial steps for filtering low-quality cells or genotypes were carried out in Tapestri Insights software with default parameters, which included removing genotype in cells with quality smaller than 30 or read depth smaller than 10, or alternative allele frequency smaller than 20, cells with smaller than 50% of genotypes present, and variants mutated in <1% of cells. We then selected known clinical variants that are annotated as “Pathogenic” or “Likely Pathogenic” in the ClinVar database (20) for subclone annotation or variants that have been reported to be in AML. Subclones were identified using Tapestri Insights 2.2 using the selected variants. Major clones are defined by the subclones that show genotype in at least 1% of the cells. Allele drop-outs were identified where there was a matching pair of a homozygous clone and wildtype clone.

We filtered out all synonymous mutations, SNVs with minor allele frequency smaller than 2.5%, SNVs with population base minor allele frequency greater than 0.1% from the Genome Aggregation Database (gnomAD v2.1.1, reference genome GRCh37/hg19, both exon and whole genome-based datasets; ref. 21) and ExAC database (ref. 22; GRCh37/hg19 reference genome). GenVisR package was used to visualize the variants. We further extracted mutations in genes that are highly frequently mutated in the The Cancer Genome Atlas (TCGA)-LAML (2), Beat AML (5), or German-Austrian AML Study Group (AMLSG) (3, 4) projects.

Co-mutation and mutual exclusivity of gene mutations were estimated using three different methods, including Discover (23), Fisher exact test, and permutation test for four datasets separately, namely (i) TCGA-LAML (2), (ii) Beat AML (5), (iii) AMLSG (3, 4), and a combined cohort from both TCGA-LAML and Beat AML. Genes that have been mutated in at least five samples were taken into account. All hypothesis tests were followed by Benjamini–Hochberg (BH) multiple testing correction. Co-mutation pairs were selected using the following criteria: (i) gene pairs with FDR less than 0.05; and (ii) the number of samples showing co-mutation greater than 3. Mutually exclusive gene pairs were selected on the basis of FDR values (<0.05).

Ranked co-mutation or mutually exclusive gene pairs (FDR < 0.05) resulted from different resources (TCGA-LAML, Beat AML, AMLSG, and combined cohort of TCGA and Beat AML) and methods were used to compute the aggregated score using the aggregateRanks function in the R package of RobustRankAggreg (24). This algorithm allows one to merge lists of different lengths. In the final list, it detects the items that are ranked consistently better than expected under the null hypothesis of uncorrelated inputs and assigns a significance score for each item (aggregated score).

Igraph was used to compute communities of co-mutated genes in the co-mutation graph using the random walk method (cluster_walktrap function in the igraph R package; ref. 25). The co-mutation network is treated as a weighted undirected graph with weights for each node as 1 − log₁₀ (aggregated score) resulting from the rank aggregation method.

BeWITH is also used to identify modules with different combinations of mutation and interaction patterns (26). We used the BeME-WithCo algorithm to identify modules which have co-occurring mutations within the modules while having mutual exclusivity between the modules. In its optimization function, the weight for the co-mutation and mutual exclusivity is defined weight (w) = max(1 − log P, 7) and the sensitivity is 10⁻⁶ (any P value smaller than 10⁻⁶ considered as 10⁻⁶) where P is the aggregated score from the rank aggregation method.

Random Forest classifier in Sklearn python library was used to predict drug sensitivity using the mutation features, co-occurrence/mutual exclusivity features and variant allele frequency. We classified the samples into drug sensitive and resistant groups using either the 50% of the maximum plasma drug concentration or the median of the IC₅₀, whichever is smaller. Plasma concentration of a drug indicates the pharmacokinetic aspect of a drug in the in vivo environment which results from the drug absorption, distribution, metabolism, and excretion. Comparing the IC₅₀ with half of the Cmax (maximal plasma concentration) provides a reasonable way to simulate drug response effect in the in vivo environment. The Cmax was obtained from the literature. The features include: (i) the mutation status of genes which show mutation frequency of at least 3% of samples (value range: 0 or 1), (ii) the mutually exclusive pattern for genes in category (i) (0 or 1), and (iii) the variant allele frequency for variants of genes in category (i). For each drug, a stratified sampling approach is used to select a training set and an independent test set from both the resistant and sensitive groups, with a ratio of 4:1. The training process was repeated 20 times with randomly selected training sets after which the resulting models were used to predict the balanced accuracy scores for the remaining test sets. The median values of the balanced accuracy scores were used to rank the drugs. Models with balanced accuracy greater than 0.6 were selected to analyze the most important features that contribute to the prediction of drug sensitivity/resistance. We averaged the feature importance and selected the top important features.

We performed Spearman correlation analysis for all the feature and drug IC₅₀ values. Multiple testing correction was performed using the BH method. Features with negative correlation suggest the observation of the event or high allele frequency of the variants is associated with sensitivity to the drug while features with a positive correlation suggest they are associated with resistance to the drug.

For the combinatory features that were predictive to drug sensitivity, a Kruskal–Wallis test was performed to test whether the drug response (AUC) showed differences among the four groups of patients defined by the combinatory features in an independent cohort (27).

Drugs with IC₅₀ values less than 200 nmol/L in at least 10% of samples were selected to perform correlation analysis. Spearman correlation was performed to measure the correlation coefficient between the IC₅₀ values for each drug. Hierarchical clustering was used to cluster the correlation coefficient with “complete” linkage and the clusters are visualized using pheatmap R package.

We used Wilcoxon rank-sum test to compare the drug IC₅₀ value distribution between the samples with or without a gene mutation. Multiple testing correction was performed using the BH adjustment. Genes with a P value smaller than 0.05 in the rank sum test were designated as signature genes for drug sensitivity or resistance.

We used the clinical feature of de novo or relapsed AML. We compared the distribution of drug IC₅₀ values between the two groups using the Wilcoxon rank-sum test followed by BH adjustment. The difference between the median log-transformed IC₅₀ value of the relapsed group and the de novo group were determined.

All data are accessible from the Supplementary Tables. Additional data, including deidentified targeted sequencing data, is available from the corresponding author upon reasonable request.

This study includes a cohort of 99 patients with AML of which >70% had relapsed, with a median age of 58 (range: 19–83). The distribution of clinical features such as antecedent hematologic disorder (AHD), European LeukemiaNet (ELN2022) classification, complete remission (CR), duration of overall survival, and the status of patients (de novo or relapse) are shown in Table 1. More comprehensive clinical data are shown in Supplementary Table S1, including duration of the first CR, cytogenetics, blast percentage, and others. Blast populations from the AML patient blood or marrow samples (enriched to >80% using magnetic beads when needed) were subjected to ex vivo drug screening and targeted sequencing. We collected drug screening data for all 99 samples and full targeted sequencing data for 75 out of the 99 samples (Fig. 1A) owing to the limitation in cell numbers. With the targeted sequencing and ex vivo drug screening datasets, we then analyzed the mutation patterns for patients with AML in our study as well as public datasets, built machine learning models to predict drug sensitivity for AML patient-derived samples, and identified the features that are predictive of sensitivity and resistance (Fig. 1A). Single-cell genomic sequencing for eight samples was performed to confirm the mutation pattern at the single-cell level and reveal mutational heterogeneity.

Genetic variants for common AML-related genes were sequenced using the MyAML gene panel assay, comprised of 194 genes, including 36 translocations with an average sequencing depth greater than 1,000×, allowing the detection of variants with low allele frequency. By excluding variants with population-based allele frequency >0.01 observed in ExAC (28) and gnomAD (21), we found AML-related genes such as DNMT3A (with mutation frequency: 28%), WT1 (25%), NRAS (25%), NOTCH2 (24%), TET2 (23%), RUNX1 (21%), FLT3 (21%), STAG2 (19%), TP53 (17%), NPM1 (17%), PTPN11 (13%), NF1 (12%), and CEBPA (12%), as the most frequently mutated genes (Fig. 1B; Supplementary Table S2). The mutation frequency distribution is similar to previous studies including TCGA (2), Beat AML (5), and AMLSG (3, 4). The alterations of NRAS and KRAS are shown with amino acid substitutions at positions 12 and 61 (Supplementary Fig. S1). The alterations in WT1, NPM1, and FLT3 were typically due to deletions or insertions. The distribution of VAF for specific sites in several other highly frequently mutated genes is shown in Supplementary Fig. S1. We found the mutation frequency is associated with AHD or disease transformation from both our dataset and the Beat AML dataset. For example, NPM1 shows lower mutation frequency in patients with AHD or with disease transformation. TP53 shows higher mutation frequency in patients with ADH or with disease transformation (Supplementary Fig. S2).

The ex vivo drug screening library includes 207 drugs or inhibitors that target a broad spectrum of pathways and molecular aberrations (Supplementary Table S3, the list of abbreviations shown in Supplementary Table S4). Among them, 103 drugs are FDA-approved drugs. The current FDA-approved drugs used in AML include idarubicin, venetoclax, daunorubicin, azacitidine, cytarabine, decitabine, fludarabine, cladribine, enasidenib, ivosidenib, midostaurin, and mitoxantrone. Log-transformed IC₅₀ (μmol/L, inhibitory concentration at half maximal effect; Supplementary Table S5) and AUC values (Supplementary Table S6) are measured for each drug in different patients.

The large drug screening library allows us to select potential drugs for patients with AML. With the criteria of drugs having an IC₅₀ less than 200 nmol/L in at least 10% of the tested samples, 91 drugs were revealed that show potential utility in either most or some of the samples (Fig. 2A; Supplementary Table S7). Among them, we observed the FDA-approved AML drugs idarubicin (with IC₅₀ <200 nmol/L in 96% samples), venetoclax (83%), mitoxantrone (83%), daunorubicin (81%), clofarabine (76%), ponatinib (45%), and volasertib (27%). The drug concentration threshold for the selection of sensitive or resistant samples needs to be compared with the concentration achieved in patients to be clinically relevant. For example, high-dose cytarabine has peak plasma concentration of approximately 115 μMol/L and 94% of patient samples are sensitive to cytarabine if the threshold is raised to 115 μMol/L. Azacitidine and decitabine are hypomethylating agents (HMA), known to affect the methylation of the CpG islands that affect gene expression (29). In addition, HMAs do not exhibit much cytotoxicity as single agents in the 72-hour assay and often take 1–6 months in patients to achieve remission as their action is related to changes in gene expression. There are alternative approaches to determine drug sensitivity to azacitidine (30). The IDH inhibitors ivosidenib and enasidenib exhibit low drug sensitivity in the assay as they also do not have in vitro direct cytotoxicity as single agents due to their mechanism of action of suppressing the high levels of 2-hydroxyglutarate that occur in IDH1-mutated AML that block differentiation (31).

Beyond the observation of different drug responses of FDA-approved AML drugs, we also found other drugs or chemicals that show potential sensitivity in AML samples, including drugs that have been approved for the treatment of other diseases or undergoing clinical investigation (Supplementary Table S7). Those potential drugs include the histone deacetylase inhibitors romidepsin, panobinostat, mocetinostat, and belinostat, the 20S proteasome inhibitors approved for multiple myeloma, bortezomib, and carfilzomib, the HSP90 inhibitors luminespib and ganetespib, the CDK inhibitors dinaciclib and alvocidib, and the XPO1/CRM1 inhibitors KPT185 and selinexor. More than half of the samples show sensitivity to the BCL-2 inhibitors venetoclax, navitoclax, and ABT737. We also observed that 43% to 37% of the samples show sensitivity to the apoptosis inhibitor YM155 and the NFκB inhibitor BAY 11-7085. Drugs such as MEK inhibitors, PI3k/mTOR inhibitors, VEGFR/FGFR inhibitors, CHK1 inhibitors, BMP inhibitors, JAK inhibitors, and BET inhibitors show high activity in a smaller percentage of samples.

Hierarchical clustering of the Spearman correlation coefficient matrix for the drug sensitivity data reveals that drugs with similar actions or targets are often highly correlated, indicating a similar drug response. For example, the MEK inhibitors selumetinib, trametinib, pimasertib, AZD8330, and mirdametinib show similar response patterns (Supplementary Fig. S3; Supplementary Table S8). The BCL-2 family inhibitors venetoclax, navitoclax, and ABT737 are also clustered together (Supplementary Fig. S3).

As many of the drugs show sensitivity in a subset of the AML samples, we investigated the genes or mutations that could serve as biomarkers for patient stratification for each of the drugs tested. We applied rank-sum tests to test the significance of each gene that contains genetic variants in at least five samples by measuring the difference of median IC₅₀ values in the samples with and without mutations in that gene. Among the top significant gene-drug sensitivity/resistance associations, we find the TP53 mutation is associated with resistance to AML drugs cytarabine and mitoxantrone (P value <0.01; Fig. 2B). The mutations in CEBPA, ASXL2, and ELL are associated with sensitivity to venetoclax (P value <0.01). The association between the CEBPA mutation and sensitivity to venetoclax can also be supported from a clinical trial observation that suggested the association of CEBPA biallelic mutation with a favorable response to venetoclax (21). Our result also indicates that mutations in PRDM9, KMT2A, SRRM2, and STAG2 are associated with drug sensitivity to daunorubicin (P value <0.01) and an SRRM2 mutation is associated with sensitivity to azacitidine (P value <0.01).

The most frequently mutated genes in AML as observed in our datasets as well as other publications are associated with sensitivity or resistance to different types of drugs (Fig. 2C). The mutations of NRAS or KRAS are associated with a wide range of MEK inhibitors, including trametinib, pimasertib, selumetinib, mirdametinib, and AZD8330. Similar results are observed from the Beat AML dataset (5). We also found that NRAS mutations are associated with sensitivity to PI3K/mTOR inhibitors AZD-8055 and PF04691502. TET2 mutations may be associated with resistance to trametinib (P value <0.01). The genetic alterations of FLT3 are associated with sensitivity to receptor tyrosine kinase (RTK) inhibitors including sorafenib, PIK75, TKI258 and lenvatinib, and PIM kinase inhibitor SGI1776 (Fig. 2C). Similarly, NPM1 and DNMT3A mutations are associated with sensitivity to certain RTK inhibitors, such as TKI258 (dovitinib) and dasatinib. The associations between the mutated gene and drug sensitivity for all significant results are shown in Supplementary Table S9. The power of single-gene alterations for the association with drug sensitivity is low. When applying multiple testing corrections, none of these associations show significance with the threshold of FDR smaller than 0.05. This may be due to the limitation of sample size, but it also suggests that comprehensive gene mutation patterns in which mutations are aggregated on the basis of their co-occurrence or mutual exclusivity may work better for the prediction of drug response than single gene mutations.

We also found that certain gene fusions were associated with drug sensitivity. RUNX1-RUNXT1(t(8;21)) and CBFB-MYH11(inv16) are the core binding factor (CBF) mutations which play critical roles in most hematopoietic lineages (32). We found several drugs that show significantly different activities for the CBF class (Supplementary Fig. S4). With the MyAML panel, we detected gene fusions/translocations for some of the samples (Supplementary Table S10), but given the small numbers of each type, we did not find significant associations with drug sensitivity.

To attain a more systematic understanding of the mutation landscape of AML and how the genetic patterns are associated with drug sensitivity, we analyzed public datasets (2–5). A total of 15 genes in common show mutation frequency in over 5% of patients with AML, including FLT3, NPM1, DNMT3A, IDH2, RUNX1, IDH1, TET2, NRAS, TP53, CEBPA, WT1, PTPN11, KRAS, ASXL1, and STAG2 (Supplementary Fig. S5A and S5B). The mutations of the 15 commonly mutated genes dominate the AML samples (Supplementary Fig. S5C). To investigate the co-occurring and mutually exclusive patterns in patients with AML, we used multiple statistical methods (Fisher exact test, permutation test, and DISCOVER; ref. 23) on multiple datasets (TCGA-LAML, Beat AML, AMLSG, or a combined dataset of TCGA-LAML and Beat AML), and integrated the statistical results using the rank aggregation method (24). We further transformed the co-mutation pairs into a co-mutation graph by assigning a weight for each co-mutation pair or mutually exclusive pair as a negative log₁₀-transformed ranking aggregation score (adding 1 to the score to make sure the weight for the graph will be a positive number). The co-mutation graph is shown in Fig. 3A and supporting data are provided in Supplementary Table S11.

We detected the community of the mutated genes using a graph-based approach (ref. 25; Fig. 3A). The NRAS and KRAS mutations are not covered by the co-mutation graph, but they show high mutation frequency (16%) in AML (33). We then added NRAS and KRAS to the graph as an independent group. As shown in Fig. 3B, five co-occurring mutation groups were defined and named as five groups: the RUNX1 group, the CEBPA group, the NPM1 group, and the TP53 group as these key genes have the highest degrees in the co-mutation subgraph, and the remaining RAS group includes only RAS (NRAS, KRAS) mutations. To confirm the co-occurring mutation clusters, we also applied another algorithm, Be-With (26), which considers both the co-mutation and the mutual exclusivity of mutations to obtain the clusters that are enriched in co-mutated genes in each cluster, but have mutually exclusive mutations in different clusters. This method gives similar results, but does not include all of the genes. The results from the Be-With algorithm show NPM1, FLT3, DNMT3A, and IDH1 are in one cluster, and RUNX1, SRSF2, STAG2, ASXL1, IDH2, BCOR, and SFRS2 are in another cluster. TET2 is clustered with CEBPA and WT1 using the Be-With algorithm as it has a mutual exclusivity relationship with the other genes in the RUNX1 cluster.

By grouping the highly frequently mutated genes by function, we found that different co-mutation subgraphs show different features (Fig. 3C). The RUNX1 co-mutation group is associated with the alteration of epigenetic modifiers and splicing factors, as well as the chromatin remodeling factors ASXL1, EZH2, the splicing factors U2AF1, SFRS2, and SRSF2, the epigenetic modulators TET2, IDH2, KDM5A, MLL, and the protein ubiquitination and degradation factors BCOR and BCORL1, and PHF6 (Fig. 3C). The TP53 group and CEBPA groups are dominated by dysregulation of different transcription factors. The CEBPA group includes mutations of CEBPA, GATA2, and WT1 and is associated with abnormal transcriptional regulation and epigenetic dysregulation. The NPM1 co-mutation group features alterations in different factors, including the epigenetic modifier DNMT3A, the genomic stability modulator NPM1, kinase FLT3, and the DNA damage-related gene RAD21. FLT3 mutations alter the PI3K-RAS-MAPK signaling pathway. The KRAS and NRAS mutations comprising the RAS group also affect the RAS/MAPK signaling pathway, downstream of FLT3.

We observed that the genes that play similar roles in a signaling pathway show mutual exclusivity, such as KIT and FLT3; both of them are receptor tyrosine kinases, which initiate the signaling in MAPK pathways. Our results show that the FLT3 and PTPN11 belong to the NPM1 co-mutation group. However, each of them shows co-occurrence with NPM1, yet mutual exclusivity is observed between FLT3 and PTPN11 (Fig. 3C). IDH1 and IDH2 show mutual exclusivity, but show co-occurrences with the RUNX1 group and NPM1 group, respectively. An independent cohort of 60 patients with AML also shows the most frequently found mutations that co-occurred with IDH1/2 mutations were DNMT3A, SRSF2, ASXL1, and RUNX1 (12).

At the single-cell level, the co-occurrence of gene mutations arises from the clonal hierarchy during the evolutionary history of the tumor's development, while at the patient level, the co-occurrence of gene mutations results from the coexistence of two subclones (34). We hypothesize that the single-cell level co-occurrence may have a higher probability of co-occurrence, as the double mutations may contribute to cancer cell proliferation; on the other hand, the patient level coexistence of subclones may show a lower probability of co-occurrence as it would be more random. From the aggregated co-mutation graph (Supplementary Table S11), we found the co-mutation of FLT3-NPM1, DNMT3A-NPM1, and DNMT3A-FLT3 are among the top of the co-mutated gene pairs. Single-cell genotyping studies have proved that the three genes often co-occur in the same AML cell (34, 35). Previous single-cell genotyping studies also reported the coexistence of FLT3, NRAS, KRAS, PTPN11, and KIT in the same patient while they were often mutually exclusive at the cellular level. It is suggested that the sequence of mutation acquisition of distinct patterns of clonal evolution in AML follows linear and branching models of evolution (27).

To understand the clonal heterogeneity of the tested patient samples, we performed single-cell genomic sequencing and simultaneous phenotyping using cell surface markers for selected patients using the MissionBio platform. We focused on eight representative samples (corresponding to the reagent kits that provide materials for eight samples) with multiple AML-related genetic alterations. We detected specific clones with distinct mutations using the MissionBio 45-gene myeloid panel. The analysis of the phenotype for these samples with the BioLegend TotalseqTM D Heme Oncology cocktail for 42 surface antigens demonstrated heterogeneity at the single-cell level, but no differences in surface markers by subclone as detected in the eight samples. The variant allele frequency as measured in the bulk genome sequencing (NGS by MyAML) is consistent with the variant allele frequency estimated using the single-cell DNA sequencing (Fig. 4A). The number of cells detected ranges from 4,167 to 28,723, with total reads ranging from 294 to 893 million. Many subclones can be detected using the MissionBio gene panel. However, we detected some major clones that constitute a large proportion of cell populations (Fig. 4B). By defining subclones using selected variants (Supplementary Table S12), different subclones can be observed for each patient-derived sample (Fig. 4C). Although multiple subclones exist in the same patient, we observe dominant clones exhibiting co-occurring variants in the same cells. In the bulk targeted sequencing, we observed mutual exclusivity of FLT3 and NRAS; however, in some patients, such as AML186, we can still observe that both variants exist in the same patient. By single-cell genome sequencing, we observed cells with an FLT3 mutation or an NRAS mutation are commonly found in different cells, or different clones, which still show mutual exclusivity at the single-cell level. As different clones may respond differently to drugs, it will be important to consider the co-occurrence of mutations within clones and the relative proportion of clones to provide better predictions for drug response.

We investigated whether the gene mutation patterns can predict drug sensitivity for different drugs. We selected drugs that have greater than 40 measurements for predictive model construction using Random Forest classifiers (see Materials and Methods). Samples were labeled to sensitive or resistant groups using either the 50% of the maximum plasma drug concentration from literature (Supplementary Table S13) or the median of the IC₅₀, whichever is smaller. The results suggest that gene mutation patterns are more predictive for some of the drugs, such as the MEK inhibitors trametinib and selumetinib, chemotherapy drugs daunorubicin and mitoxantrone, CDK inhibitor alvocidib, and FLT3 inhibitors quizartinib and gilteritinib (Fig. 5A).

We then selected models with balanced accuracy greater than 0.6 in the test sets to evaluate the most important features for the drug sensitivity prediction. The most important features for several selected drugs are shown in Fig. 5B–D. The combinations of features such as TET2-wildtype-and-NRAS-mutant or the allele frequency of NRASp.G12D are predictive of sensitivity to the MEK inhibitor trametinib (Fig. 5B). This is confirmed by the Spearman correlation between the genetic alteration events or their combinations and the drug response data (Supplementary Fig. S6A and S6B). It is also confirmed through a statistical analysis using Beat AML data (Supplementary Fig. S7), which indicates samples with both NRAS and TET2 mutation show highest sensitivity to trametinib, followed by samples with NRAS mutation and TET2 wildtype. The association between all the drugs in Fig. 5A and genetic features can be found in Supplementary Table S14. We found TP53 mutation and U2AF1 mutation are associated with resistance to mitoxantrone (Supplementary Fig. S6C), and NPM1-wildtype-and-TP53 mutation is predictive for the resistance to mitoxantrone (Fig. 5D; Supplementary Fig. S6C). The prediction accuracy of the venetoclax model just using the gene mutation is still low, and we did not find variants that significantly correlated with the sensitivity of venetoclax using the threshold we chose. It may suggest other features still need to be taken into consideration when predicting the drug sensitivity of venetoclax, such as data at the expression level.

Our ex vivo drug screening included both patients with relapsed AML (N = 72) and de novo AML (N = 27). We then identified the drugs that show higher sensitivity in patients with relapsed AML than in patients with de novo AML. We found patients with relapsed AML show higher sensitivity to several drugs, including arsenic trioxide (ATO), melphalan, rigosertib, afatinib, navitoclax, and azacitidine (Supplementary Fig. S8A). The samples with top ATO sensitivity exhibit either FLT3 or TP53 mutations. We also observed that navitoclax shows significantly higher sensitivity in patients with relapsed AML (Supplementary Fig. S8A). We further compared the drug response for de novo versus relapsed patients with specific mutations. In our analysis of patients with FLT3 mutations, we observed that drugs like afatinib demonstrate greater sensitivity in samples obtained from relapsed patients, compared to those from de novo patients (Supplementary Fig. S8B). On the other hand, dactolisib shows higher sensitivity in the de novo patient samples. A phase II clinical trial of pazopanib has been conducted in patients with AML (36). Our results suggest it shows higher sensitivity in patients with relapsed AML, which may guide future clinical trials.

The prediction of drug responses based on the molecular profiling of patients with AML now opens a rational path for choosing the best therapeutic strategy. Our study provides a rich resource including the genetic variants and correlative ex vivo drug screening data on individual patient samples for more than 200 drugs for community use and provides clues for the selection of potential drugs or drug targets for further drug development for the application of AML. The statistical associations between ex vivo drug screening and genetic information reveal biomarkers for the prediction of drug sensitivity or resistance for a variety of drugs. For example, our analysis reveals that AML samples with RAS mutations show higher sensitivity to several MEK inhibitors, as has been reported clinically for the individual drug trametinib (37). We expect that the systematic analysis of the genetic patterns and definitions of co-mutation groups will facilitate more robust prediction of drug sensitivity or resistance.

Leveraging large datasets in AML, we observed co-occurring gene mutation patterns in patients with AML, which may have developed during the clonal evolution characteristic of pathogenesis of the disease. The co-occurring gene mutations can originate from coexisting distinct clones which are defined by the occurrence of individual mutations in different cells, or the accumulation of mutations in the same cell. The single-cell analysis performed for this study confirms the coexistence of different clones in the same patient, as well as accumulated mutations in one single clone. Because cells with different mutations respond differently to drugs, co-occurrence data may provide clues to predict the drug response given the heterogeneity of cell populations in each individual patient. With the co-occurrence of different clones, we may look into drug combinations that can be effective against multiple clones. With the mutation pattern exhibited in a single clone, we may expect to identify one drug or combination that may be effective against that clone.

The cohort in our study has a higher fraction of relapsed patients than other large datasets, such as TCGA-LAML and Beat AML. Our results suggested several drugs that could elicit better responses in the relapsed patients than in the patients with de novo AML, providing clues to find more effective drugs for relapsed AML.

It is important to emphasize that all the data for drug sensitivity and resistance reported here as well as in the other databases were based on ex vivo studies. However, the current ex vivo drug sensitivity testing methods are now showing predictive value or correlations with response in clinical trials (17, 38, 39). Moreover, an ex vivo assay for venetoclax drug sensitivity exhibited very high positive and negative predictive capability in a multicenter trial (40). While our current study utilized ex vivo drug screening of AML patient samples enriched for blasts using magnetic beads for CD34 or lineage depletion, we acknowledge the critical role of leukemia stem cells (LSC) in AML propagation and relapse, as noted in our prior study (41). We have also separately described differential sensitivity of the LSC and blast populations derived from the same patients demonstrating resistance to the common drugs used in AML in the LSCs as compared with the blasts (42). Future clinical trials that incorporate multicolor flow cytometry analysis (43) of not only viability but also apoptosis for LSCs and differentiated cell populations may improve the predictive models. This method that enables gating on the populations of interest allows them to be assayed in the presence of all bone marrow cells, including the components of the immune and stromal microenvironment that has been shown to confer drug resistance. A multi-omic approach considering the data presented here correlating drug sensitivity associated with specific mutations or networks of mutations in combination with functional drug screening, incorporating other parameters such as gene expression (44) proteomics and metabolomics, may further optimize outcomes of precision medicine clinical trials.

R.J. Monnat reports a patent for U.S. Patent Application 17/257,945 pending. I. Shmulevich reports being a cofounder of Twinome Health. P.S. Becker reports grants from NCI and University of Washington Foundation, and other support from Invivoscribe, Inc., during the conduct of the study as well as grants from Notable Labs, Glycomimetics, Pfizer, and GPCR Therapeutics and personal fees from Accordant Health Services outside the submitted work; in addition, P.S. Becker has a patent for High Throughput Drug Screening of Cancer Stem Cells, PCT/US2019/040617, submitted July 3, 2019, by University of Washington pending. No disclosures were reported by the other authors.

G. Qin: Conceptualization, resources, data curation, software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. J. Dai: Resources, data curation, software, methodology. S. Chien: Resources, data curation, methodology. T.J. Martins: Resources, data curation, formal analysis, methodology, writing–original draft. B. Loera: Formal analysis, investigation, writing–original draft, writing–review and editing. Q.H. Nguyen: Data curation, formal analysis, writing–original draft, writing–review and editing. M.L. Oakes: Resources. B. Tercan: Formal analysis, investigation, methodology, writing–original draft, writing–review and editing. B. Aguilar: Formal analysis. L. Hagen: Resources, data curation, writing–review and editing. J. McCune: Conceptualization, resources, data curation. R. Gelinas: Writing–review and editing. R.J. Monnat: Conceptualization, supervision, funding acquisition, writing–review and editing. I. Shmulevich: Conceptualization, resources, data curation, supervision, funding acquisition, investigation, methodology, writing–original draft, writing–review and editing. P.S. Becker: Conceptualization, resources, data curation, formal analysis, supervision, funding acquisition, investigation, methodology, writing–original draft, project administration, writing–review and editing.

This project was supported by NCI, NIH grants P01CA077852 (R.J. Monnat), and R01CA270210 (I. Shmulevich). G. Qin, R.J. Monnat, B. Tercan, I. Shmulevich, and P.S. Becker were supported by P01CA077852. G. Qin, B. Tercan, and I. Shmulevich were supported by NCI CTD² project U01CA282109 and U01CA217883. G. Qin, B. Aguilar, B. Tercan, and I. Shmulevich are supported by NCI grant R01CA270210. This work was made possible, in part, through access to the Genomics Research and Technology Hub Shared Resource of the NIH Cancer Center Support Grant (P30CA-062203) at the University of California, Irvine, and NIH shared instrumentation grants 1S10RR025496-01, 1S10OD010794-01, and 1S10OD021718-01, and the Clinical Trials Support office of the University of Washington/Fred Hutchinson Cancer Center Cancer Consortium Cancer Center Support Grant (P30 CA015704). P.S. Becker was supported by NIH Cancer Center Support Grant (P30CA-062203). We acknowledge information stored in a RedCap database supported by an ITHS grant to University of Washington UL1 TR002319. We are grateful for support from private donors to the University of Washington Foundation. We appreciate early collaboration with Invivoscribe to optimize mutation testing in AML. The content is solely the responsibility of the authors and does not necessarily represent the official views of the NIH.