Comparative Genomic Landscape of AML

We performed exome sequencing on 622 of the specimens from the cohort representing 531 unique patients. The final, high-confidence variant list (Supplementary Information, Table S7) revealed a range of 1–80 somatic variants per patient (cohort median of 13 somatic variants) (Extended Data Fig. 1). Comparison of the top 33 most commonly mutated genes across Beat AML and TCGA⁶ showed generally similar frequencies. Higher frequency of mutations in serine and arginine rich splicing factor 2 (SRSF2) were seen in Beat AML than in TCGA, and this difference was conserved comparing only the de novo cases in Beat AML with TCGA (Fig. 1A). In contrast, mutational events that were seen at less than 2% frequency across Beat AML and TCGA were much more divergent with variants in 221 mutant genes called in common, 939 mutant genes called only in TCGA, and ~1,500 mutant genes called only in Beat AML, irrespective of whether we compared only de novo or non-de novo cases (Fig. 1B,​,C).C). Most of these divergent mutational events were observed in only single patients, however, there were mutations in 11 genes that were called in 1% or more of patients in Beat AML, but were not observed in prior AML sequencing studies (Extended Data Fig. 1). Finally, co-occurrence and exclusivity of the most frequent variants were computed revealing significant patterns of mutational co-segregation suggesting biological cooperativity of certain mutational events (Fig. 1D).

Open in a separate window

Fig. 1.

Comparative Genomic Landscape of AML.

(A) Frequency of the 33 mutational events that were cumulatively most frequent in Beat AML (n=531 patients) and TCGA (n=200 patients). Top row represents the full Beat AML cohort and the middle bar represents only the de novo Beat AML cases. Mutations were summarized by gene as was done by TCGA whereas the FLT3-ITD mutations were kept separate in the rest of this manuscript.

(B) Mutational events at 2% frequency or less in the de novo cases of Beat AML and TCGA were compared for overlap. Venn diagram displays the overlap with the small circles within each compartment representing a size-scaled frequency of each mutational event.

(C) Analysis as in panel B with only the non-de novo Beat AML cases versus TCGA.

(D) Co-occurrence or exclusivity of the most recurrent mutational events in the Beat AML cohort (n=531 patients) were assessed using the DISCOVER⁴⁰ method with a dot plot indicating the odds ratio of co-occurrence (blue) or exclusivity (red) using color-coding and circle size as well as asterisks that indicate FDR-corrected statistical significance.

Drug Response and Gene Expression

RNA sequencing was performed on 451 specimens from 411 patients. Clustering of the 2,000 most variably expressed genes across the cohort revealed gene expression signatures that associated with many of the prominent genetic and cytogenetic disease groups (Extended Data Fig. 2). To understand the profile of sensitivity and resistance to a variety of small-molecule inhibitors, we profiled primary tumor cells from 409 specimens derived from 363 patients against a panel of 122 small-molecule inhibitors using an ex vivo drug sensitivity assay³⁷. Drug sensitivity patterns were analyzed with respect to clinical and genetic features of tumors (Extended Data Fig. 3). We compared the average area-under-the curve (AUC) values for each drug between the de novo and transformed samples via a series of single factor analysis of variance (ANOVA) tests. Generally, transformed cases showed less sensitivity to most drugs than de novo cases. Of the 122 drugs tested, 64 were significantly (false discovery rate (FDR) < .1) more sensitive in the de novo samples, while only 1 drug Panobinostat (HDAC inhibitor) was significantly more sensitive in the transformed cases (Extended Data Fig. 4). In addition, we analyzed the concordance of drug sensitivity patterns with respect to predicted target gene, gene family, or pathway for each drug (drug assignments to target families found in Supplementary Information, Table S11). This analysis revealed drug targets/families that were highly concordant amongst constituent members as well as drug families that were quite discordant (Extended Data Fig. 5). To create a global view of overall sensitivity or resistance of each case, we generated a heatmap of binary sensitive/resistant calls for each sample to each drug. We then annotated the sensitivity or resistance fraction of each case against European Leukemia Net (ELN) 2017 (Extended Data Fig. 6A) and World Health Organization (WHO) 2016 (Extended Data Fig. 6B) classifications.

Whole Exome and Custom Capture Validation Sequencing

For exome sequencing, Illumina Nextera RapidCapture Exome capture probes and protocol were utilized, which provided coverage of 37 Mb of genomic DNA coding regions. Briefly, following initial QC on a TapeStation (Agilent), 50 ng of intact genomic DNA was fragmented and tagged (tagmentation) in a single step. Following clean-up, the tagmented DNA was amplified by 10 cycles of PCR, which added the indexed adaptors for clustering and sequencing. Libraries were hybridized to capture pools in 12 sample sets with two rounds of hybridization performed to increase specificity. Libraries recovered with streptavidin magnetic beads were amplified by 10 cycles of PCR, unincorporated reagents were removed with AMPure beads (Agencourt), and validated on the TapeStation. Quantification of capture pools was done using real time PCR (Kapa). Libraries were denatured, flow cells set up using the cBot (Illumina), and run on a HiSeq 2500 using paired end 100 cycle protocols. For the AML tumor samples, 5 or 6 lanes were run per capture group. For the matched skin biopsy samples, 3 lanes (or equivalent) were run per capture group. The instrument and chemistry for all capture groups are provided in the Supplementary dashboard tables (Supplementary Information, Table S12,13).

For validation of sequencing results, an 11.9 megabase custom capture library was developed that provided coverage of all variants previously reported in AML as well as all new variants detected from exome sequencing in this study (Genes, Variants, and Capture Regions for this custom library are found in Supplementary Information, Table S14). This capture library was then applied to sequence 96 specimens that had previously been subjected to whole exome sequencing for validation of variant calls.

Whole Exome Sequencing Data Processing

We developed customized analytical pipelines that combined published algorithms with novel filtering, curation, and quality control steps. Detailed depictions of analytical workflows are found in Supplementary Information (Table S6) and on our online browser (http://vizome.org/additional_figures_BeatAML.html). Initial data processing and alignments were performed with commonly used analytical tools. For each flowcell and each sample, the FASTQ files were aggregated into single files for read 1 and 2. During this process these reads were trimmed by 3 on the 5’ end and 5 on the 3’ end. BWA MEM version 0.7.10-r789⁴¹ was used to align the read pairs for each sample-lane FASTQ file. As part of this process, the flowcell and lane information was kept as part of the read group of the resulting SAM file. The Genome Analysis Toolkit (v3.3) and the bundled Picard (v1.120.1579) were used⁴² for alignment post-processing. The files contained within the Broad’s bundle 2.8 were used including their version of the build 37 human genome (These files were downloaded from: ftp://ftp.broadinstitute.org/bundle/2.8/b37/). The following steps were performed per sample-lane SAM file generated for each CaptureGroup:

• The SAM files were sorted and converted to BAM via SortSam

• MarkDuplicates was run, marking both lane level standard and optical duplicates

• The reads were realigned around indels from the reads--RealignerTargetCreator/IndelRealigner.

• Base Quality Score Recalibration

The resulting BAM files were then aggregated by sample and an additional round of MarkDuplicates was carried out at the sample level. Quality control reports were generated using the ReportingTools⁴³ and qrqc⁴⁴ Bioconductor R packages along with sequencing core and alignment output files. Each AML sample BAM was paired with its skin biopsy pair and an additional round of indel realignment was carried out to ensure consistency of genotypes between the two samples. If an AML sample did not have a pair, the indel realignment was instead done at the sample level.

Whole Exome Sequencing Variant Detection

For genotyping, each AML/skin biopsy pair was realigned at the sample level and then genotyped for single nucleotide variations using Mutect v1.1.7⁴⁵ and Varscan2 v2.4.1⁴⁶. Indels were produced using Varscan2. Each variant call format (VCF) file was annotated using the Variant Effect Predictor v83 against GRCh37⁴⁷. The resulting VCF files were filtered to include only those annotated to a gene and were converted to mutation annotation format (MAF) format using the vcf2maf v1.6.6 tool⁴⁸.

Mutect v1.1.7⁴⁵ was run using default parameters except that no limit was placed on the number or frequency of the alternative allele frequency in the normal to help address normal contamination.

Varscan2 v2.4.1⁴⁶ was run in somatic mode with the recommended filtering scheme⁴⁹ except as shown in Supplementary Information, Table S23.

Indels and single nucleotide variants (SNVs) were produced for the tumor-only samples again using Mutect without a specified normal for consistency and VarScan2 in mpileup2indel or mpileup2snp mode respectively.

These variants were assigned to their most deleterious effect on Ensembl transcripts using Ensembl VEP v83 on GRCh37. This assignment was done using the same VEP parameters as the vcf2maf (v1.5.0) program.

The TCGA AML variants⁶ in MAF form were downloaded from the GDC archive site: https://portal.gdc.cancer.gov/legacy-archive/files/c410d927-d49c-4d4f-8356–601bee563ebe. The MAF was converted to VCF files using the vcf2maf suite⁵⁰. The resulting VCF files were lifted over from NCBI36 to GRCh37 of the human genome using CrossMap⁵¹. Only those variants that successfully lifted over were kept. Variants from Table S2 in Jaiswal et al 2014¹⁴ were extracted and further processed removing variants that were ambiguous in terms of external sources and could not be found in their whole exome sequencing (WES) variants. The unique set of BeatAML variants was annotated relative to RefSeq transcripts using Ensembl VEP similar to above and all consequences were kept. This set of variants and consequences was searched against the set of processed Jaiswal variants.

Using the runs from MuTect and VarScan, these data were next filtered to keep only the protein impacting SNVs and indels from Mutect and VarScan2 and filtered requiring that the variants had at least 5 reads and either not be seen in the exome aggregation consortium (ExAC)⁵² or be seen at a frequency < .01. These data present several additional challenges. First somatic calls cannot be obtained directly from the tumor-only samples, second there is always a possibility of tumor contamination of the skin samples for those samples that were paired. To address these issues and maximize comparability we used an iterative approach. The following was done separately for the two genotypers:

• An initial set of higher confidence somatic mutations were retrieved from the paired tumor/normal samples requiring tumor variant allele frequency (VAF) >= 8% and normal VAF <= 5% in addition to the significance tests already performed by the programs.

• A list of all candidate mutations was collated requiring that a mutation was either seen in the high confidence somatic set, the set of variants from Jaiswal et al or from the lifted over set of variants from the TCGA AML paper.

• Mutations from the overall set were kept if:

  • ○ The overall number of calls in the paired samples was not more than twice the number of high confidence somatic calls

  • ○ The tumor-only frequency for the calls was less than 50% greater than the number of calls in the paired samples

  • ○ The mutation was seen in Jaiswal/TCGA list

• High Confidence somatic mutations were kept regardless

The data from the two genotypers were combined along with FLT3-ITD calls from Pindel⁵³. Comparing our variant lists from whole exome sequencing versus our custom capture validation sequencing, we noticed, similar to others⁵⁴, that low allele frequency C->A variants (< 15%) tended to have poor concordance (7.7%; data not shown) between the initial run and the technical validation run. These variants were removed in these data and along with a curated ‘blacklist’ (Supplementary Information, Table S15) of known problematic variants/genes including mitochondrial DNA variants. In addition, all variants that were seen in a cumulative list of BeatAML normal samples at a frequency greater than 1% were removed. Cumulatively, of this set, 94% of covered single-nucleotide variants were validated with 82% of insertion/deletion calls also being confirmed with validation sequencing.

Manual review was then carried out in the following steps:

• a.)

  The addition back of all Jaiswal flagged rows.

• b.)

  Reviewed all TCGA flagged rows for VAF pattern that matched or did not match with known drivers in same specimen. Some TCGA variants were added back based on convincing VAF pattern and known pathogenic role, other TCGA variants were kept excluded based on VAF pattern unlike known drivers in same specimens.

• c.)

  Other variants were added back based on other specimens that had the same variant that was still on the include list and VAF pattern looked convincing for inclusion.

• d.)

  All Jaiswal genes with only frameshift/nonsense variants were manually reviewed and missense mutations were manually removed.

• e.)

  Genes/variants that were on both the include and exclude lists were manually reviewed and were removed if c to a with over 15% VAF, did not validate, and/or VAF pattern unlike known drivers in same specimen

• f.)

  Further review of all genes in summary sheet with cohort frequency of 8 or more (1% of more). Any that were not familiar from knowledge of AML literature were manually reviewed for VAF patterns that did or did not match known drivers within same specimens. Those that did not match were manually removed.

After this manual review, additional curated mutations from the UnifiedGenotyper run were added back in along with a curated set of variants from tumor-only patients in Jaiswal et al genes. The tumor-only AML samples utilized in Zhang et al (submitted) were removed and utilized for that study.

Internal FLT3-ITD and NPM1 mutation detection

A subset of samples was tested for FLT3-ITD and NPM1 mutation status using an internally run PCR assay and capillary electrophoresis. Genomic DNA was extracted from fresh blood or bone marrow aspirates of AML patients was used to detect the presence or absence of FLT3-ITD and NPM1 4 bp insertion mutations^55,56. Primers for FLT3 spanned approximately 330 bp to include the common internal duplication site⁵⁵. Primers for NPM1 spanned approximately 170 bp to cover the clustered multiple insertional or insertion/deletion sites^56,57. Primers were HPLC-purified by the manufacturer. The multiplex PCR reaction solution consisted of 100 ng gDNA, 10 pmol of the respective forward and reverse primers for FLT3 and NPM1, 25 mmol/l MgCl₂, 2·5 mmol/l dNTP, 5 μl 10 × PCR buffer, 0·2 μl AccuTaq DNA polymerase and water in a total volume of 50ul⁵⁸. The PCR conditions were: initial denaturing for 5 min at 94°C, followed by 30 cycles at 94°C for 30 s, 56°C for 45 s and 72°C for 30 s with a final cycle of 10 min at 72°C. The PCR products were diluted 1:10 and analyzed by capillary electrophoresis on a QIAxcel high resolution DNA cartridge according to the manufacturer’s protocol.

Forward Primer FLT3: 5′ - AGCA ATT TAG GTA TGA AAG CCA GCTA - 3′

Reverse Primer FLT3: 5′ - CTT TCA GCA TTT TGA CGG CAA CC - 3′

Forward Primer NPM1: 5′ - GTT TCT TTT TTT TTT TTT CCA GGC TAT TCA AG - 3′

Reverse Primer NPM1: 5′ - CAC GGT AGG GAA AGT TCT CAC TCT GC - 3′

Derivation of FLT3-ITD and NPM1 consensus calls

Consensus FLT3-ITD and NPM1 mutation calls found in the clinical summary table (Supplementary Information, Table S5) were determined by comparing the internal capillary PCR test (internal; method described above), the Clinical Laboratory Improvement Amendments/College of American Pathology (CLIA/CAP) laboratory run test (Sequenome, GeneTrails, Foundation Medicine, Genoptix, Illumina). The internal test was used for the sample consensus call when available as it was performed on the exact sample that was used for further ex vivo drug sensitivity assays. Where discordance existed between the internal test versus the CLIA lab test results, the sample was flagged for manual review. The trace file for the internal test was visually inspected and if discordance with the CLIA/CAP test results persisted, the whole exome sequencing data was then used to help determine the consensus call.

Derivation of CCAAT enhancer binding protein alpha (CEBPA) biallelic consensus calls

N-terminal and C-terminal mutations have been described to occur on opposing alleles and patients harboring CEBPA biallelic mutations have been shown to fall into a favorable risk category⁵⁹. Patients were scored positive for biallelic CEBPA mutation if described in the clinical notes as biallelic or double positive. Patients were also scored as CEBPA biallelic if both N-terminal and C-terminal mutations were identified in the whole exome sequencing data.

RNA-Sequencing and Data Processing

All samples were sequenced using the Agilent SureSelect Strand-Specific RNA Library Preparation Kit on the Bravo robot (Agilent). Briefly, poly(A)+ RNA was chemically fragmented. Double stranded cDNAs were synthesized using random hexamer priming with 3’ ends of the cDNA adenylated then indexed adaptors were ligated. Library amplification was performed using three-primer PCR using a uracil DNA glycosylase addition for strandedness. Libraries were validated with the Bioanalyzer (Agilent) and combined to run 4 samples per lane, with a targeted yield of 200 million clusters. Combined libraries were denatured, clustered with the cBot (Illumina), and sequenced on the HiSeq 2500 using a 100-cycle paired end protocol. In addition to the AML samples, there was also a sample of purified CD34 molecule (CD34)+ cells from healthy control bone marrow, which was included in each sample group (for a total of 12 times sequencing this control RNA). This control served as both a healthy comparator and a quality check on inter-group batch effects. 21 additional, individual healthy bone marrow samples were also included, two of which were CD34-selected (17–00053 and 17–00056) with the other 19 being whole mononuclear bone marrow cells from healthy donors.

Workflows for processing and analysis of RNA sequencing data to generate gene counts and gene fusions for each sample are shown on our online browser (http://vizome.org/additional_figures_BeatAML.html) with processed gene expression values for each specimen listed in Supplementary Information (Table S8,9). For each flowcell and each sample, the FASTQ files were aggregated into single files for read 1 and read 2 (if not already done by the sequencing core). During this process these reads were trimmed by 3 on the 5’ end and 5 on the 3’ end. Alignments of reads was performed using the subjunc aligner (1.5.0-p2)⁶⁰. BAM files obtained from subjunc were used as inputs into featureCounts (1.5.0-p2)⁶¹ and gene-level read counts were produced. For a reference genome, the GRCh37 build provided by the Broad as part of the GATK bundle was used. Gene assignments were based on the Ensembl build 75 gene models on GRCh37. See the parameters below for software usage.

subjunc -i /path/to/reference/ -u -r fastq1 -R fastq2 -o outputBAMFilename -I 5 -T 7 -d 50 -D 600 -S fr

featureCounts -a Homo_sapiens.GRCh37.75.gtf -o output -F GTF -t exon -g gene_id -s 2 -C -T 10 -p -P -d 50 -D 600 -B BAM_files

The data was collated from featureCounts matrices and all genes with no counts across the samples were excluded. Genes with duplicate gene symbols and those where the counts were < 10 for 90% or more of the samples were additionally removed prior to normalization similar to the approach suggested for weighted gene correlation network analysis (WGNCA)⁶². Samples for which their median expression was less than 2 standard deviations below the average were removed from the dataset (N=10). Normalization was performed using the conditional quantile normalization procedure⁶³, which produced GC-content corrected log2 reads per kilobase per million mapped reads (RPKM) values. This procedure produces both offsets to be used in conjunction with edgeR as well as a matrix of log2 normalized RPKM values for clustering.

In addition, the subjunc BAM files were processed using the RNA-sequencing genotyping protocol (as of GATK v3.3) which was similar to the WES protocol described in the ‘Whole Exome Sequencing’ section including for each sample:

• MarkDuplicates

• SplitNCigarReads

• RealignerTargetCreator/IndelRealigner.

• Base Quality Score Recalibration

The resulting BAM files were used to produce RNA genotypes using the UnifiedGenotyper for the purposes of QC and ethnicity estimation.

Gene fusion data was additionally generated using the TopHat-Fusion (v2.0.14) program using default parameters⁶⁴.

Coexpression network formation

We formed coexpression modules using the WGCNA procedure on the RNA-sequencing data from the ‘RNA-Sequencing and Data Processing’ section. All RNA-sequencing samples were used to form the set of modules. Due to the heterogeneity of clinical expression data we generated ‘signed hybrid’ networks using ‘bicor’ correlation setting the max proportion of outliers to .1⁶⁵. We ran the procedure multiple times, varying several parameters to choose the most relevant set for further analysis. The WGCNA procedure was run on datasets formed from the top 2,000 and 5,000 most variable genes. For each dataset we set the ‘power’ variable to either 2 or 3. For each of these runs we varied the module detection parameters of dynamicTreeCut⁶², namely the deepSplit parameter was set to 0 or 2 and the pamStage parameter was set to TRUE or FALSE. For each of these sets of modules we computed a series of module quality statistics⁶⁶, mean correlation, mean adjacency, mean MAR, mean KME, proportion of variability explained and the mean cluster coefficient. Significance of modules was determined by computing a Zscore of each of these values relative to the mean and standard deviation of those from 100 random assignment of modules. We chose the set of modules to use in our analyses as those that were most correlated with the ‘specimentSpecificDx’ using the module quality as a tie-breaker. The analysis set of modules was chosen to be the version using the 5,000 most variable genes, power set to 2 and modules formed using deepSplit=2 and pamStage=F. Of this set of modules only the grey module didn’t have a summary Z statistic (median across the 4 density measures) of at least 2. Additionally after correcting the data using the estimated principal components⁶⁷, the module structure didn’t change appreciably (data not shown).

Quality Control

The UnifiedGenotyper runs for both the WES and RNA-sequencing were combined into a single VCF file using the GATK CombineVCFs functionality. This combined VCF file was converted to a GDS file using SNPRelate (1.12.2)⁶⁸. Note the version is an upper bound as several versions were used across the entire project). The overall similarity of the genotypes of each pair of samples were computed, termed identity by state (IBS) and a hierarchical clustering was performed using one minus this similarity. From this clustering and visualization we had devised hard cutoffs for further inspection based on the types of data being compared. For instance samples not meeting the specified IBS thresholds (DNA-DNA=.9; RNA-RNA=.83; DNA-RNA=.89) were subject to manual review. Based on the dendrogram structure as well as the clinical/lab information, samples were either excluded, assigned to a different patient ID or in rare cases assigned to a different sample. It was observed that bone marrow transplants between sample collections produced a noticeable but milder effect in these dendrograms and such samples were flagged for removal in RNA-sequencing analysis and for treatment as tumor-only samples in the WES analysis as is described in the ‘WES Variant Detection’ section.

Fusion annotation for analysis

Fusions calls were determined from a consensus of three datatypes, a specific diagnosis categorization at the time of sample acquisition, current set of clinical karyotypes and fusions detected in RNA-sequencing data by TopHat-Fusion. All sources were limited to the same set of known fusions: RUNX1-RUNX1T1, CBFB-MYH11, MLLT3-KMT2A, DEK-NUP214, GATA2-MECOM and PML-RARA. It was determined that the RNA-sequencing calls did not provide additional resolution in detecting these known fusions and was not performed on all the samples so the consensus was limited to the clinical karyotype calls as well as the specific diagnosis categorization (which was determined based on karyotype and other cytogenetic clinical tests). Overall the calls were based on the karyotype data except in 10 cases, 3 where the karyotype and diagnosis was sufficiently complex to warrant a separate ‘Complex’ categorization. The remaining 7 of these cases were set to the specific diagnosis classification. It should be noted there was additional support from the RNA-sequencing data for several of these cases.

Ethnicity

The combined RNA and WES VCF from the ‘Quality Control’ step was merged with a set of Hapmap genotypes⁶⁹ lifted over to build 37. The SNPRelate package was used to convert the VCF to GDS, perform LD pruning using an LD threshold of .2, MAF cutoff of .05 and allowing a missing rate of .3 and calculation of the principal components. The methodology of Zheng and Weir 2015⁷⁰ was used to assign admixture proportions relative to the HapMap samples using the principal components. Each sample was assigned to an ethnicity group based on the group with the maximum admixture proportion. If the maximum was 50% or less, we labeled it Admixed. As we have observed previously that the clustering of ethnicities for RNA-sequencing samples are more diffuse than exome sequencing, we assigned the final inferred ethnicity to each patient based on the distinct WES calls if available, deferring to RNA-sequencing only if not available. If multiple exome sequencing samples were present with discrepant calls, the patient was subject to manual review. The only patient where this occurred was patient 4043, a self-identified Hispanic who had 2 RNA samples and an exome sequencing sample inferred as White and one exome sample inferred to be HispNative. Only the White call for the exome sequencing had an admixture proportion over .5. The patient was kept consistent with the self-identification and labeled as HispNative.

Sex

For DNA, coverage is first computed over the Y chromosome and the counts for each sample are added up and log10 transformed (after adding 1 to all the counts). K-means clustering is used to assign samples to 2 clusters with the cluster with the lower mean labeled as the Female cluster.

For RNA-sequencing, counts were converted to counts per million (CPM) after applying the Trimmed Mean of M scaling normalization⁷¹. A set of 28 genes chosen to successfully discriminate the genders via DE analysis over multiple studies (data not shown) were used in conjunction with Kmeans clustering to form two clusters. The Female cluster was labeled based on high XIST expression.

ELN 2017 Classification

This procedure is based off of the categorization in Table 5 of the 2017 ELN update paper⁵

Karyotypes in the clinical file were first cleaned and parsed into clones/subclones and distinct abnormalities using standard conventions⁷². The current representation was corrected for nomenclature type (e.g. idem vs sl) in a basic manner. For instance, ambiguous events such as chromosomal loss (e.g. −15) was not corrected for whether the preceding clone had a counteracting gain. Also, additional +/− symbols in conjunction with valid karyotype operators in a separate clone (e.g. +del(12)(q?15) or -del(12)(q?15)) ) were treated separately with gains (+) being kept in the unique count of events and losses (−) being removed.

Abnormalities were first checked for the following categories:

1. RUNX1-RUNX1T1

2. CBFB-MYH11

3. MLLT3-KMT2A

4. DEK-NUP214

5. KMT2A-*

6. BCR-ABL1

7. GATA2-MECOM

8. −5/del(5q) (termed minus_5)

9. −7 (termed minus_7)

10. −17 (termed minus_17)

Note the bold were further considered to be WHO recurrent fusions.

The number of unique abnormalities (across clones) was then computed.

Whether or not a karyotype was considered to be polyploid was also recorded (at least 60 chrs or ‘(>=3)n’ or label).

NPM1, FLT3-ITD, and biallelic CEBPA were derived from consensus calls. FLT3-ITD allelic ratios were determined solely for the samples with an internal assay. The MAF values of the internal assay were converted to a ratio using the formula MAF/(1-MAF).

RUNX1, ASXL1 and TP53 were derived from the clinical genotypes.

Abn_17 calls were manually curated from the karyotype data and clinical genotype calls respectively.

The determination of ELN 2017 categories proceeds by assigning TRUE/FALSE/NA values to one or more of the 5 columns (3 ELN and 2 ambiguous) in the following way:

1. isFavorable is considered TRUE if a sample has at least one of the following:

   1. RUNX1-RUNX1T1

   2. CBFB-MYH11

   3. Positive NPM1 and negative FLT3-ITD

   4. Positive NPM1 and positive FLT3-ITD with allelic ratio < .5

   5. Biallelic CEBPA

2. isFavorableOrIntermediate

   1. NPM1 is positive and FLT3-ITD is positive but the allelic ratio is not available

3. isAdverse is considered TRUE if a sample has at least one of the following:

   1. DEK-NUP214

   2. KMT2A-*

   3. BCR-ABL1

   4. GATA2-MECOM

   5. minus_5

   6. minus_7

   7. minus_17

   8. abn_17

   9. 3 or more abnormalities and not a WHO recurrent fusion

   10. One monosomy (autosomal) and at least one additional abnormality except for CBFB-MYH11

   11. Positive RUNX1 or ASXL1 and not considered to be isFavorable or isFavorableOrIntermediate

   12. Positive TP53

4. isIntermediate is considered TRUE if a sample has at least one of the following:

   1. MLLT3-KMT2A

   2. NPM1 positive and positive FLT3-ITD with allelic ratio >= .5

   3. NPM1 negative and negative or low allelic ratio (<.5) FLT3-ITD

   4. At least one abnormality and is not considered isFavorable or isAdverse

5. isIntermediateOrAdverse

   1. NPM1 negative and FLT3-ITD positive without an allelic ratio

NAs occur in the absence of FLT3-ITD or NPM1 calls.

Samples where the specific diagnosis at inclusion indicated ‘Acute promyelocytic leukaemia with t(15;17)(q22;q12)’ were automatically set to ‘Favorable’.

Any overlaps of the categories were resolved based on manual expert review.

Ex vivo Functional Drug Screens

Ex vivo functional drug screens were performed on freshly isolated mononuclear cells from AML samples. Briefly, 10,000 cells per well were arrayed into three, 384-well plates containing 122 small-molecule inhibitors. This panel contained graded concentrations of a drugs with activity against two-thirds of the tyrosine kinome as well as other non-tyrosine kinase pathways, including mitogen activated protein kinases (MAPKs), phosphatidylinositol-4,5-bisphosphate 3-kinase/AKT serine/threonine kinase 1/mechanistic target of rapamycin kinase (PIK3C/AKT/MTOR), protein kinase AMP-activated (AMPK/PRKAA), ATM serine/threonine kinase (ATM), Aurora kinases, calcium/calmodulin dependent protein kinases (CAMKs), cyclin-dependent kinases (CDKs), serine/threonine protein kinase 3 (GSK3), I-kappaB kinase (IKK), cAMP dependent protein kinase (PKA), protein kinase C (PKC), polo-like kinase 1 (PLK1), and RAF proto-oncogene serine/threonine kinase (RAF). In addition, the library contained small molecule inhibitors with activity against the BCL2 family, bromodomain containing 4 (BRD4), Hedgehog, heat shock protein 90 (HSP90), NOTCH/gamma-secretase, proteasome, survivin, signal transducer and activator of transcription 3 (STAT3), histone deacetylase (HDAC), and WNT/beta-catenin. Drug plates were created using inhibitors purchased from LC Laboratories and Selleck Chemicals and master stocks were reconstituted in dimethyl sulfoxide (DMSO) and stored at −80 °C. Master plates were created by distributing a single agent per well in a seven-point concentration series, created from three-fold dilutions of the most concentrated stock resulting in a range pf 10 μM to 0.0137 μM for each drug (except dasatinib, ponatinib, sunitinib, and YM-155 which were plated at a concentration range of 1 μM to 0.00137 μM). DMSO control wells and positive control wells containing a drug combination of Flavopiridol, Staurosporine and Velcade were placed on each plate, with the final concentration of DMSO ≤0.1% in all wells. Daughter plates were created using a V&P Scientific 384-well pin tool head operated by the Caliper Sciclone ALH 3000 and equipped with 0.457mm diameter, 30 nanoliter, slotted stainless steel pins (cat num: FP1NS30). Daughter and destination plates were sealed with pealable thermal seals using a PlateLoc thermal sealer. Destination plates were stored at −20 °C for no more than three months and thawed immediately before use. Primary mononuclear cells were plated across single-agent inhibitor panels within 24 h of collection. Cells were seeded into 384-well assay plates at 10,000 cells per well in Roswell Park Memorial Institute (RPMI) 1640 media supplemented with fetal bovine serum (FBS) (10%), l-glutamine, penicillin/streptomycin, and β-mercaptoethanol (10−4 M). After 3 d of culture at 37 °C in 5% CO2, MTS reagent (CellTiter96 AQueous One; Promega) was added, optical density was measured at 490 nm, and raw absorbance values were adjusted to a reference blank value and then used to determine cell viability (normalized to untreated control wells).

Ex vivo Functional Drug Screen Data Processing

A workflow summarizing data normalization, curve fit parameters, and quality assurance/quality control (QA/QC) steps is found on our online browser (http://vizome.org/additional_figures_BeatAML.html) with processed drug response data for each specimen listed in Supplementary Information (Table S10). A given sample was run on one or more panels and within each panel, the majority of drugs were run without within-panel replicates. Two steps were performed to harmonize these data prior to model fitting:

1. A ‘curve-free’ AUC (integration based on fine linear interpolation between the 7 data points themselves) was calculated for those runs with within-panel replicates after applying a ceiling of 100 and a floor of 0 for the normalized viability. The maximum change in AUC amongst the replicates was noted and those runs with differences > 100 were removed.

2. Remaining within-plate replicates had their normalized viability averaged and subject to a ceiling of 100 and floor of 0. An additional set of ‘curve-free’ AUCs was computed for sample-inhibitor pairs run on multiple panels. The maximum change in AUC amongst the across-panel replicates was noted and those runs with differences > 75 were removed.

At this point, the within and across plate replicates for the normalized viability were averaged together and a ceiling of 100 was applied. From the steps above, the floor was already at 0.

Based on the methodology used in our prior drug combination study⁷³, a probit regression was fit to all possible run groups using the model:

Where for all groups there were N=7 dose-response measurements.

The summary measures of curve fit were inspected and cutoffs were devised removing all runs with an AIC > 12 and deviance > 2. For inhibitors that were run using multiple concentration ranges, only the latest concentration range was kept. Finally, these data were compared to the AUC values from third order polynomial fits. Those runs that were discrepant in terms of sensitive/resistant calls were manually reviewed as subject to removal.

Ex vivo Functional Drug Screen Analysis

For all drug analyses requiring a call of sensitivity or resistance (e.g. the gene expression signatures), sensitivity/resistance was determined by the lowest and highest 20% of the AUC values for each drug.

Expression Analysis and Integration with Ex Vivo Functional Drug Screen

For all the below analyses the earliest sample was chosen for each patient.

Mutation Analysis and Integration with Ex vivo Functional Drug Screen

For all the below analyses where groups of samples were compared, the earliest sample was chosen for each patient.

We thank all of our patients at all sites for donating precious time and tissue. DNA and RNA quality assessments, library creation, and short read sequencing assays were performed by the OHSU Massively Parallel Sequencing Shared Resource. Sheenu Sheela, Catherine Lai, Katherine Lindblad and Kary Oetjen assisted in study coordination at NIH. Barry Sawicki and Christina Cline assisted in study coordination at the University of Florida. S. Ravencroft assisted with patient sample shipping and data entry and K. Schorno provided project management and support of activities at the University of Kansas Cancer Center. Jack Taw assisted with patient sample shipping and Shyam Patel assisted with data entry at Stanford University.

Funding Funding for this project was provided in part by a Leukemia Lymphoma Therapy Acceleration Grant to B.J.D. and J.W.T. and by support provided by the Knight Cancer Research Institute (Oregon Health & Science University, OHSU). Supported by grants from the National Cancer Institute (1U01CA217862, 1U54CA224019, 3P30CA069533–18S5) and NIH/NCATS CTSA UL1TR002369 (SKM, BW). ASB was supported by the National Library of Medicine Informatics Training Grant (T15LM007088). JWT received grants from the V Foundation for Cancer Research, the Gabrielle’s Angel Foundation for Cancer Research, and the National Cancer Institute (1R01CA183947). CRC received a Scholar in Clinical Research Award from The Leukemia & Lymphoma Society (2400–13), was distinguished with a Pierre Chagnon Professorship in Stem Cell Biology and Blood & Marrow Transplant and a UF Research Foundation Professorship. This work was supported in part by the Intramural Research Program of the National Heart, Lung, and Blood Institute of the National Institutes of Health.

The authors declare competing financial interests detailed below.

Competing Interests JWT receives research support from Agios, Aptose, Array, AstraZeneca, Constellation, Genentech, Gilead, Incyte, Janssen, Seattle Genetics, Syros, Takeda; JWT is a co-founder of Vivid Biosciences. MWD serves on the advisory boards and/or as a consultant for Novartis, Incyte, and BMS and receives research funding from BMS and Gilead. CSH receives research funding from Merck. TLL consults for Jazz Pharmaceuticals and receives research funding from Tolero, Gilead, Precient, Ono, Bio-Path, Mateon, Genentech/Roche, Trovagene, Abbvie, Pfizer, Celgene, Imago, Astellas, Karyopharm, Seattle Genetics, and Incyte. DAP receives research funding from Pfizer and Agios and served on advisory boards for Pfizer, Celyad, Agios, Celgene, AbbVie, Argenx, Takeda and Servier. BJD serves on the advisory boards for Gilead, Aptose, and Blueprint Medicines. BJD is principal investigator or coinvestigator on Novartis and BMS clinical trials. His institution, Oregon Health & Science University, has contracts with these companies to pay for patient costs, nurse and data manager salaries, and institutional overhead. He does not derive salary, nor does his laboratory receive funds from these contracts. The authors certify that all compounds tested in this study were chosen without input from any of our industry partners.