A clinical transcriptome approach to patient stratification and therapy selection in acute myeloid leukemia

Abstract

As more clinically-relevant genomic features of myeloid malignancies are revealed, it has become clear that targeted clinical genetic testing is inadequate for risk stratification. Here, we develop and validate a clinical transcriptome-based assay for stratification of acute myeloid leukemia (AML). Comparison of ribonucleic acid sequencing (RNA-Seq) to whole genome and exome sequencing reveals that a standalone RNA-Seq assay offers the greatest diagnostic return, enabling identification of expressed gene fusions, single nucleotide and short insertion/deletion variants, and whole-transcriptome expression information. Expression data from 154 AML patients are used to develop a novel AML prognostic score, which is strongly associated with patient outcomes across 620 patients from three independent cohorts, and 42 patients from a prospective cohort. When combined with molecular risk guidelines, the risk score allows for the re-stratification of 22.1 to 25.3% of AML patients from three independent cohorts into correct risk groups. Within the adverse-risk subgroup, we identify a subset of patients characterized by dysregulated integrin signaling and RUNX1 or TP53 mutation. We show that these patients may benefit from therapy with inhibitors of focal adhesion kinase, encoded by PTK2, demonstrating additional utility of transcriptome-based testing for therapy selection in myeloid malignancy.

Fig. 1. Experimental overview and short nucleotide variant/indel analysis.

A Overview of datasets used at each stage of the analysis, colored by sequencing project, with sample size for each data set indicated in brackets. For the AML PMP exploratory batch, MDS samples were used for profiling the SNV and gene fusion pipelines, but not for expression-based analyses, and validation RNA-Seq libraries were prepared in triplicate. B–C Matched variant allele frequency (VAF) for variants between whole-genome sequencing (WGS) and RNA-Seq (B) and whole exome sequencing (WES) and RNA-Seq (C), by concordance status, variant type, and coverage status (sites with ≤10x coverage are indicated as ‘Low WGS/WES Depth’). Selected variants discussed in the text are labeled with Human Genome Variation Society (HGVS) nomenclature. D Variant observations in CEBPA. Potentially disruptive mutations for each patient (y-axis) are indicated by their chromosomal coordinate (x-axis), with predicted impact and variant concordance. E Sensitivity and positive predictive value (PPV) for each variant caller in the validation cohort (with 95% confidence intervals).

Fig. 2. Structural variation analysis.

A Number of structural variation events detected by cytogenetics analysis vs. RNA-Seq fusion detection. B $z$-scaled gene expression of GATA2 and MECOM in the AML PMP exploratory patient cohort. MECOM- (red) and KMT2A-related (blue) structural variants are colored. The dashed horizontal line indicates the threshold of $z$ ≥ 3 for MECOM outlier expression.

Fig. 3. Training and validation of the AML prognostic score (APS) gene expression signature.

A APS model coefficients. The y-axis indicates the genes making up the APS set, with the x-axis indicating the model coefficients. B–F Survival plots for the AML PMP (B, n = 154), TCGA LAML (C, n = 173), Beat AML (D, n = 293), AML PMP Prospective (E, n = 42), and TARGET pediatric AML (F, n = 156) cohorts, for above-median and below-median values of the APS value within each cohort. Dashed lines indicate time to median survival. Log-rank p values are indicated for each cohort.

Fig. 4. Univariate survival analysis.

A–C Univariate survival analysis for selected clinical, gene fusion, mutation, stratification, and expression-based predictors of survival for the AML PMP (A), TCGA LAML (B), and Beat AML (C) cohorts. In each panel, the x-axis indicates the hazard ratio, and the y-axis indicates the FDR-adjusted p-value for each tested variable for univariate Cox proportional hazards models. Point sizes are scaled to reflect the proportion of patients affected by the relevant marker, and colored as clinical variables—red, gene fusions—green, SNV/indels—purple, expression—blue. D Distribution of APS values for AML PMP patients, colored by ELN-RNA risk status (favorable —blue, intermediate—green, adverse—red). Horizontal dashed lines indicate the bottom, intermediate, and top terciles of the APS value distribution. Comparison bars show two-sided Wilcoxon test p-values for each pairwise comparison.

Fig. 5. Patient stratification.

A Comparison of patient stratification by ELN-RNA and ELN-RNA-APS models for the AML PMP cohort. The width of each link is scaled to the number of patients shared between each set of risk categories. B–C Survival curves for the ELN-RNA and ELN-RNA-APS models for the AML PMP cohort. D–I Stratification and survival curves for the TCGA LAML and Beat AML patient cohorts. Patient groups are colored as in Fig. 4D (favorable—blue, intermediate—green, adverse—red), with log-rank p values indicated.

Fig. 6. Pathway enrichment analysis.

A Pathway enrichment analysis of IPA canonical pathways overrepresented in differentially expressed genes across all three cohorts. Pathways are ranked by mean activation $z$-score, with point size scaled by Fisher’s exact test p values. B Geneset enrichment analysis (GSEA) against the Reactome pathway database. Pathways are ranked by the number of cohorts the pathway was enriched in, then by the mean FDR-adjusted p value across cohorts. GSEA p values are derived from permutation testing, and corrected for multiple testing using the FDR method. In both panels, gene sets containing PTK2 and/or ITGB3 are highlighted in bold and italic.

Fig. 7. Differential gene expression for patients with first-tercile vs. third-tercile APS values.

A–C Volcano plots of differentially expressed genes for the AML PMP (A), TCGA LAML (B), and Beat AML (C) cohorts. Genes with absolute log₂fold-change of 1 and FDR-adjusted p ≤ 0.1 (based on Wald tests as implemented in DESeq2) are highlighted in red (over-expressed in third-tercile APS) and blue (over-expressed in first-tercile APS). D Comparison of differentially expressed genes between cohorts. All genes which were differentially expressed in any single cohort are displayed with their log2FoldChange values across all cohorts. The top 40 over-expressed and bottom ten under-expressed genes (by mean fold-change across cohorts) are labeled.

Fig. 8. Correlation of PTK2 expression with specific mutations.

A Samples from each cohort ranked by $z$-scaled expression of PTK2. The upper annotation plot for each cohort indicates the presence of selected mutations. Indicated p values are for Kruskal-Wallis tests comparing continuous PTK2 expression against presence or absence of specific mutations. B Top highly correlated proteins with FAK. $r^2$ values indicate Pearson correlation coefficients. C Relative FAK protein expression in patients with and without mutant RUNX1, TP53, or FLT3 in the AML Proteome Atlas. Two-sided $t$-test p-values are indicated for each comparison. D Western blot assessment of FAK and RUNX1 expression in KG1a and THP-1 cell line derivatives. Cell lines were treated with either (C) control shRNA, (1) shRUNX1-1, or (2) shRUNX1-2, with molecular weights quantified in kilodaltons. Quantifications scaled to controls are shown as bar plots. Each blot represents a single experiment, with the exception of RUNX1 in THP-1, which was performed twice with similar results. E, F Colony-forming cell assay for AML cell lines with short-hairpin RNAs against RUNX1, and treated with FAK inhibitors VS-4718 (0.5 μm for KG1a, 1.5 μm for THP-1) or Defactinib (1 μm for both cell lines). The indicated p values correspond to two-sided $t$-tests for each comparison.