Integrative analysis of drug response and clinical outcome in acute myeloid leukemia

Abstract

Acute myeloid leukemia (AML) is a cancer of myeloid-lineage cells with limited therapeutic options. We previously combined ex vivo drug sensitivity with genomic, transcriptomic, and clinical annotations for a large cohort of AML patients, which facilitated discovery of functional genomic correlates. Here, we present a dataset that has been harmonized with our initial report to yield a cumulative cohort of 805 patients (942 specimens). We show strong cross-cohort concordance and identify features of drug response. Further, deconvoluting transcriptomic data shows that drug sensitivity is governed broadly by AML cell differentiation state, sometimes conditionally affecting other correlates of response. Finally, modeling of clinical outcome reveals a single gene, PEAR1, to be among the strongest predictors of patient survival, especially for young patients. Collectively, this report expands a large functional genomic resource, offers avenues for mechanistic exploration and drug development, and reveals tools for predicting outcome in AML.

Keywords: JEDI; LSC17; MEGF12; cell state; eigengene; hematologic malignancy; leukemia stem cell; monocyte; targeted therapy.

Figure 1.. Genomic features and outcomes for the Beat AML cohorts are concordant.

(A) The percentage of patients with a somatic mutation or gene fusion is shown for Waves 1+2 and Waves 3+4 cohorts and percentage difference is shown (average difference is 1.037%). (B) Overall survival is shown for Waves 1+2 and Waves 3+4 cohorts (p-value 0.2 log-rank test). (C). Ex vivo drug response values are shown comparing the average AUC for each drug (points) within de novo AML patients comparing Waves 1+2 to Waves 3+4 (Pearson’s correlation r=0.965). (D) Waves 3+4 can serve as validation cohort to assess prior mutation-inhibitor associations (Tyner et al., 2018). The difference in average response of specimens that are mutated versus wild type for a given gene (effect size) are plotted for Waves 1+2 on the x-axis and Waves 3+4 on the y-axis. Effect sizes were expressed as Glass’s delta with respect to the wild type group. Significance was based on the Welch’s t-test comparing mutated vs wild type and requiring a minimum of five mutations in Waves 1+2 and three mutations in Waves 3+4. Adjusted significance in Waves 1+2, Waves 3+4, or both are also annotated (qvalue < 0.05; (Storey and Tibshirani, 2003)).

Figure 2.. Differentiation scoring to characterize expression and mutational patterns.

We generated malignant cell-type gene expression scores for each of our patients relative to 6 sets of 30 genes derived from expression signatures from single-cell sequencing (Van Galen et al., 2019). (A) By examining the correlation of expression module eigengenes and the cell-type scores we can see that WGCNA expression modules (X-axis) show correlation with differentiation cell-type scores (y-axis). Pearson’s correlation r-values are annotated. (B) Similarly, we can determine the mutational events significantly associated with each of the cell-type scores. Shown are the signed –log10(Welch’s t-test p-values) for the differences in cell-type score with respect to mutational status. Up to the top 5 most significant (qvalue < 0.05; (Storey and Tibshirani, 2003)) mutational events associated with enrichment of each cell-type score are highlighted. Full data are available in Table S2.

Figure 3.. Influence of cell-type on inhibitor response.

(A) We compared the cell-type gene expression scores, developed in Figure 2, with ex vivo drug response. Pearson’s correlation of ex vivo inhibitor response (measured by AUC) with each of the six cell-type scores reveals relationship between differentiation and resistance (blue) or sensitivity (red). Note, all drugs displayed in the heatmap showed a significant correlation with at least one cell-type (BY FDR < 0.05; (Benjamini and Yekutieli, 2001)). Panobinostat and Venetoclax are shown correlated with Monocyte-like score in Figure S3A. Asterisks indicate those inhibitors in which there is an interaction between mutation and cell-type, shown below in panel B (single asterisk, qvalue < 0.1; double asterisk, qvalue < 0.05; (Storey and Tibshirani, 2003)). Full data are available in Table S3. (B) We examined whether correlations between drug response and cell-type score changed due to mutational status. Drug response modification was quantified by examining the significance (Y-axis) of the interaction between each mutational event and cell-type (x-axis) for each inhibitor, requiring a minimum of 10 mutations. Interactions with qvalue at two thresholds (q < 0.1; q < 0.05) are called out by text and shape and distinguished by a dashed and solid line, respectively. As an example, the conditional relationships between sorafenib, FLT3-ITD, and cell-type scores are shown in Figure S3B, where sorafenib activity is robust in FLT3-ITD cases with a high progenitor-like score, but activity is lost in FLT3-ITD cases with a high monocyte-like score.

Figure 4.. Genomic associations of drug family response.

As in Figure 3A, we define sensitivity by red and resistant with blue and with white indicating intermediate. (A) Drug target information, where available, was annotated for all drugs on the panel, and these targets were mapped onto pathways to create drug families that share common targets and pathways. Full mapping of drug families is available in Table S4. Family-level drug response summaries allow us to identify cohort-level responses. Some like the Aurora Kinase family (Aurk) show overall resistance and others like Type XIII RTKs (Eph) or PIKK are overall sensitive. (B). We examined correlations between mutational status and drug family response. Mutational status associates with drug family response either in terms of significant (qvalue < 0.05; (Storey and Tibshirani, 2003)) sensitivity (red) or resistance (blue). Association tests were performed using Welch’s T-test comparing mutated vs wild type and requiring at least 5 mutations for a given mutational event. Effect size was based off of Glass’s delta using wild type as the reference. Full data are available in Table S5.

Figure 5.. Drug family response is influenced by and conditional on cell-type.

(A) Similar to the analysis of single-inhibitor responses and their correlation with cell-type score, shown in Figure 3A, we can group the drug families from Figure 4 by their association with cell-type differentiation scores based on their Pearson’s correlation (BY FDR < 0.05; (Benjamini and Yekutieli, 2001)). (B) We can also determine instances where drug family correlation with mutational state is conditional on cell-type (as we did for individual drugs in Figure 3B). The mutations that play a significant drug response modification role based on the statistical interaction between cell-type score and mutation, requiring at least 10 mutations per event is shown. The Y-axis indicates the signed -log10 (Pvalue), the X-axis indicates cell-type. The interactions are listed by text (qvalue < 0.1 and qvalue < 0.05; (Storey and Tibshirani, 2003)) and distinguished by a dashed and solid line respectively. Full data are available in Table S6.

Figure 6.. Comprehensive analysis of clinical variables identifies PEAR1 expression as a prognostic factor in AML.

(A) We analyzed univariate associations between the 6 cell-type scores and the set of 9 non-redundant WGCNA eigengenes with categorical and continuous clinical variables as well as overall survival across and within age subsets. Each variable type was grouped on the Y-axis and the displayed association values were derived from the corresponding test statistics. These include the Z statistic for categorical (logistic regression) and survival outcomes (Cox Proportional Hazards) and T-statistic for the continuous outcomes (general linear model). Average linkage hierarchical clustering was applied to rows (within groups) and columns (B; top) To understand which gene(s) among the module 3 WGCNA signature are the strongest drivers of the overall signature, we calculated the correlation of each gene with the module eigengene (kME). Of all the 203 genes in module 3, PEAR1 has the highest correlation. (B; bottom) The patient expression values (dots) of PEAR1 (x-axis) correlated highly with the Mod3 module eigengene (y-axis) and both are also correlated with the ‘HSC-like’ malignant cell-type as indicated by the cell-type purple(low)-yellow(high) gradient as defined in Figure S2A. Full module gene components are available in Table S7. (C) A strategy based on a forest of conditional inference trees (cforest) was used to determine variables that differentiate overall survival. Age groups, expression modules (denoted with Mod*), gene mutation, AML gene fusions, and cell-type scores from in the harmonized dataset were all included in the model. The age groups are young (<45), middle (45–60), older (60–75), oldest (>75). Depiction of the resulting tree where oval splits indicate variables that most significantly split overall survival. Lines indicate values that were split on, for instance, “TRUE” for the ‘young’ variable, which indicates that the subgroup is <45 years old. The resulting subgroups of patients, rectangles that denote ‘terminal nodes’, are listed with the subgroup size (denoted as n) and are colored to match corresponding survival curves.

Figure 7.. Integrative modeling shows PEAR1 expression is a single, independent predictor of poor prognosis.

(A) PEAR1 expression (y-axis) is significantly higher in Adverse compared to Favorable ELN 2017 risk categorization (x-axis) for patients 60 and over (dots; older + oldest group), but the pattern is less pronounced in patients < 60. Significance determined using a using a Welch’s T-test. (B) PEAR1 expression differentiates survival for young patients equivalently to the LSC17 signature. Categorization of samples into high and low expression groups for PEAR1 and LSC17 was determined using the ctree methodology to facilitate comparison with PEAR1. (C) We compared performance of PEAR1 and LSC17 using both ctree and the median expression values for each. Shown are the hazard ratios (HRs) and 95% confidence intervals (CIs) for PEAR1 and LSC17 in the entire cohort or in the younger patients as well as in all sample types or only in bone marrow samples. The categories were high vs low expression as determined by either median threshold or significant splits using ctree (i.e., split type). In two instances, the ctree method failed to identify a split, and these are denoted with “N/A”.