Clonal architecture predicts clinical outcomes and drug sensitivity in acute myeloid leukemia

Abstract

The impact of clonal heterogeneity on disease behavior or drug response in acute myeloid leukemia remains poorly understood. Using a cohort of 2,829 patients, we identify features of clonality associated with clinical features and drug sensitivities. High variant allele frequency for 7 mutations (including NRAS and TET2) associate with dismal prognosis; elevated GATA2 variant allele frequency correlates with better outcomes. Clinical features such as white blood cell count and blast percentage correlate with the subclonal abundance of mutations such as TP53 and IDH1. Furthermore, patients with cohesin mutations occurring before NPM1, or transcription factor mutations occurring before splicing factor mutations, show shorter survival. Surprisingly, a branched pattern of clonal evolution is associated with superior clinical outcomes. Finally, several mutations (including NRAS and IDH1) predict drug sensitivity based on their subclonal abundance. Together, these results demonstrate the importance of assessing clonal heterogeneity with implications for prognosis and actionable biomarkers for therapy.

Fig. 1. Cohort curation and summary.

a A systematic literature review was performed to identify studies reporting clinically annotated samples and/or molecularly profiled cohorts of adult AML. 12 cohorts met inclusion criteria and were curated for mutations, drug screening results, clinical features, and survival outcomes. These studies were then aggregated into a uniformly annotated database with an admixture of overlapping data types available for analysis. b Oncoprint for the most frequently mutated genes in our cohort. Each column is an individual sample (n = 2914) and the color of the vertical line represents the type of mutation reported. Sex, cohort, ELN 2017 risk group, subset, and survival information is indicated on the bottom of the plot. Source data are provided as a Source Data file.

Fig. 2. Distinct patterns of mutation co-occurrence associate with overall survival.

a Co-occurrence and mutual exclusivity of the most frequent mutations present in de novo AML was performed using a two-sided Fisher’s Exact test. Significant associations (FDR < 0.05) are colored according to the odds ratio of co-occurrence (red) or mutual exclusivity (blue). Points are sized based on the number of patients with co-occurring mutations for each genotype. b Summary of pairwise mutations and their association with prognosis based on Cox proportional-hazards regression modeling. Significant genotypes (FDR < 0.05) are colored according to the log-transformed hazard ratio compared to wild-type patients, with green depicting better prognosis (HR ≤ 1) and purple representing worse prognosis (HR ≥ 1). Points are sized based on the number of patients with co-occurring mutations for each genotype. c Scatterplot of the correlation between the odds ratio and hazard ratio of co-occurring mutations from a, b. Percentages indicate the fraction of genotypes per quadrant which associate with significant (Bonferroni FDR < 0.05) survival associations. For each error band, the measure of center is the line of best fit as derived from linear regression between the odds ratio and hazard ratio for each group. Shaded bands represent 95% confidence intervals for each linear regression. Points are sized based on the number of patients with co-occurring mutations for each genotype and colored according to the log-transformed hazard ratio compared to wild-type patients, with green depicting better prognosis (HR ≤ 1) and purple representing worse prognosis (HR ≥ 1). d Frequency distribution of the number of de novo patients with the most frequent 3-way mutation combinations. Bars are colored based on the association with a significant survival correlation (p ≤ 0.05) compared to patients with only two genes mutated: red = a significant survival association, gray = no significant association. e Forrest plot (left) and Kaplan–Meier plots (right; Bonferroni FDR ≤ 0.15) depicting survival analysis between triple-mutated and double-mutated genotypes. For the forest plot (left), points represent the hazard ratios calculated between triple vs. double-mutated patients using a Cox proportional-hazards model. Significant genotypes (two-sided log-rank p ≤ 0.05) are colored: green represents cases where the presence of all three mutations correlated with improved survival, while purple hits represent genotypes where all three mutations correlated with worse survival. q-values were calculated in terms of the false discovery rate using Bonferroni correction. Points are sized relative to the number of patients with all three mutations and bars represent the 95% confidence intervals of the hazard ratios. Source data are provided as a Source Data file.

Fig. 3. Variant allele frequency associates with features of disease presentation and overall survival.

a VAF distribution for frequent mutations in our de novo cohort (n = 1636 patients). Mutations were assigned as high or low VAF based on the median copy number-corrected VAF for each gene. Red points represent cases where the VAF is above the median while blue points represent those below the median VAF. For each distribution, the boxplot represents the boundaries for the first and third quartiles with a line at each median; whiskers delimit the highest data point below the third quartile +1.5× the interquartile distance and the lowest data point above the first quartile −1.5× the interquartile distance. b Distribution of effect sizes for differences in clinical features of AML based on high or low VAF for each mutation. Points represent the effect size (Cohen’s d) between high and low VAF for each genotype across the different clinical variables (n patients: WBC = 1507; Hemoglobin = 1280; Platelet = 1292; LDH = 1228; BM blasts = 1472; PB blasts = 1385; Age = 1643). Significant associations are colored based on the level of significance (Bonferroni FDR < 0.2); error bars represent the 95% confidence intervals of the effect sizes. c Scatterplots of effect sizes for WBC levels and peripheral blood blast percentages between mutated and wild-type patients versus effect sizes calculated between high and low VAF for each mutation. Points are sized based on the number of patients analyzed and colored based on VAF effect size significance (FDR < 0.1). d Forest plot summarizing univariate Cox proportional-hazards regression modeling of common mutations based on VAF thresholds in the de novo cohort. Points represent the hazard ratio for overall survival between high and low VAF groups based on VAF thresholds calculated using maximally selected rank statistics. Points are sized based on the number of patients above the VAF threshold. Error bars represent the 95% confidence intervals of the hazard ratios. Green hits represent cases where higher VAF correlated with improved survival, while purple hits represent genotypes where increased VAF correlated with worse survival. e Scatterplot of hazard ratios calculated between mutated and wild-type patients versus hazard ratios calculated between high and low VAF for each mutation. Hazard ratios are calculated using a standard Cox proportional-hazards model. Points above the dotted line indicate mutations where greater VAF associates with worse outcomes compared to patients with lower VAF for that mutation. Points are colored by significance of VAF hazard ratio calculations (red points = Bonferroni FDR < 0.1) and sized relative to the VAF threshold for each genotype. Source data are provided as a Source Data file.

Fig. 4. Clonal dominance of co-occurring mutations stratifies survival.

a Summary of pairwise clonal relationships for recurrent mutations present in de novo AML. Color represents the fraction of patients where mutations in gene 1 occurred before gene 2 (i.e. are clonally dominant; green = gene 1 before gene 2; brown = gene 1 after gene 2). Size is scaled to reflect the total number of patients with co-occurring mutations. b Schematic depicting how pairwise VAF relationships were used to bin patients into groups (left) and forest plot of hazard ratios calculated from univariate Cox proportional-hazards modeling based on the order of pairwise mutations (right). Points represent the hazard ratio, as calculated using standard Cox proportional-hazards regression, between patients with different mutation order. Bars represent the 95% confidence intervals of the hazard ratios and points are sized based on the number of patients with defined mutation ordering. c Scatterplot (left) and Kaplan–Meier plot (right) showing how the order of mutation acquisition in patients with co-occurring mutations in NRAS and GATA2 robustly improved patient stratification. Green points/line represent cases where NRAS mutations occur before those in GATA2. Brown points/line represent cases where GATA2 mutations occur before those in NRAS. The reported p-value was calculated using a two-sided log-rank test. d A Bradley–Terry model was used to assign the relative order of global mutation acquisition from pairwise relationships as determined in a. Only patients with at least two mutations were considered in this model. Density plots (left) represent the VAF distribution for each mutation in the analysis (corrected for copy number and X-linkage in males) and are ordered on the y-axis based on their relative order of acquisition compared to all other genes in the analysis. Points and error bars (right) represent the Bradley–Terry model results for the point estimate and 95% confidence interval, respectively, for relative gene ordering in temporal acquisition. e Schematic depicting how pairwise VAF relationships were used to bin patients into groups (left) and forest plot of hazard ratios calculated from univariate Cox proportional-hazards regression modeling based on the order of mutation category acquisition (right). Bars represent the 95% confidence intervals of the hazard ratios. f, g Kaplan–Meier plots for significant pairs from e. p-values in c, f, g were calculated between nonambiguous groups using a two-sided log-rank test. Source data are provided as a Source Data file.

Fig. 5. Architecture of clonal evolution associates with survival outcomes.

a Schematized workflow for modeling clonal architecture in a cohort of WES patients. Briefly, a deep-sequenced cohort was assembled and analyzed using PyClone to generate robust clonal populations and cellular prevalences. These cellular prevalence estimates were then leveraged to model the temporal acquisition of mutations and clonal architecture using ClonEvol. b Heatmap summarizing PyClone results of the per-patient cellular prevalence for the most common clonal genotypes. Each column is an individual patient sample grouped by hierarchical clustering based on similarity in clonal patterns. For patients with multiple mutations in the same gene, only the mutation with the largest cancer cell fraction (CCF) is shown. c Correlation between mutation burden and the number of unique clones derived from PyClone in the de novo cohort (Kruskal–Wallis p = 1.6e−33; n = 409 patients). Points are colored by broad clonal evolution architecture as determined by ClonEvol (blue = branched evolution; gray = linear evolution). For each distribution, the boxplot represents the boundaries for the first and third quartiles with a line at each median; whiskers delimit the highest data point below the third quartile +1.5× the interquartile distance and the lowest data point above the first quartile −1.5× the interquartile distance. d Kaplan–Meier plot showing the association between higher median mutational burden per clone (red curve) and worse outcomes in de novo patients (two-sided log-rank test). e Kaplan–Meier plot showing the association of improved outcomes in patients exhibiting branched evolutionary architecture (blue curve = branched evolution; gray curve = linear evolution; two-sided log-rank test). f Forrest plot depicting univariate Cox proportional-hazards ratios for various aspects of the clonal architecture analyses. g Forrest plot depicting univariate Cox proportional-hazards ratios for clonal and mutational burden risk stratification based on linear or branched architecture. h Schematic depicting the different genetic and clinical features associated with evolutionary architecture. Source data are provided as a Source Data file.

Fig. 6. Clonal abundance predicts drug sensitivity in primary AML samples.

a Schematic depicting how linear regression of drug response (AUC) against VAF can identify correlations between drug sensitivity and the clonal prevalence of mutations. Red points represent sensitivity trends while blue points represent resistance trends. b Volcano plot of drug response between mutated and wild-type samples for de novo samples from the Beat AML study. Points are sized based on the number of samples analyzed and colored by significance (Bonferroni FDR < 0.1; red = sensitive, blue = resistant). c Copy number-corrected VAF distribution for mutations with paired drug data in the de novo cohort of the Beat AML study (n^{mut + drug} ≥ 5; n = 202 biologically independent patient samples). For each distribution, the boxplot represents the boundaries for the first and third quartiles with a line at each median; whiskers delimit the highest data point below the third quartile +1.5× the interquartile distance and the lowest data point above the first quartile −1.5× the interquartile distance. d Dotplot of the most significant (p < 0.05) drug-gene correlations identified through linear regression of drug AUC against mutation VAF in de novo AML samples. Points are sized based on the range of VAFs for each mutation and are colored based on the type of drug sensitivity trend (red—sensitive; blue—resistant). Asterisks represent drug-gene associations with a Bonferroni FDR < 0.1. e, f Representative binary distributions (left) and AUC-VAF scatterplots (right) for clinically relevant sensitivity and resistance VAF correlations for IDH1 (e; n^WT = 195 samples; n^Mut = 10 samples) and NRAS (f; n^WT = 177 samples; n^Mut = 14 samples)^, respectively. For each distribution (left), the boxplot represents the boundaries for the first and third quartiles with a line at each median; whiskers delimit the highest data point below the third quartile +1.5× the interquartile distance and the lowest data point above the first quartile −1.5× the interquartile distance; p-values are calculated using a two-sided Wilcoxon rank-sum test. For each scatterplot (right), shaded bands represent 95% confidence intervals for each linear regression. For each error band, the measure of center is the line of best fit as derived from linear regression between the drug AUC and VAF for each mutation-drug pair. g Schematic (top) depicting the potential correlation between the subclonal prevalence of a secondary mutation (e.g. FLT3-ITD) and sensitivity to inhibitors. AUC-VAF scatterplots (bottom) for pairwise genotypes with enough samples (n ≥ 5; DNMT3A:FLT3) where linear regression of drug AUC against VAF revealed strong resistance trends. Shaded bands represent 95% confidence intervals for each linear regression. For each error band, the measure of center is the line of best fit as derived from linear regression between the drug AUC and VAF for each mutation-drug pair. h Schematics representing possible relationships between VAF and drug response. Source data are provided as a Source Data file.