Splicing regulatory dynamics for precision analysis and treatment of heterogeneous leukemias

Abstract

The role of splicing dysregulation in cancer is underscored by splicing factor mutations; however, its impact in the absence of such rare mutations remains poorly understood. Prompted by the finding that splicing uniquely resolved genetic subtypes of cancer, we developed an unsupervised computational workflow called OncoSplice to comprehensively define tumor molecular landscapes. In adult and pediatric acute myeloid leukemia (AML), OncoSplice identified the spectrum of driver genetics from splicing profiles alone, defined more than a dozen previously unreported molecular subtypes recurrent across AML cohorts, and discovered a dominant splicing subtype that partially phenocopies U2AF1-mutant splicing. Although pediatric leukemias lack splicing factor mutations, this U2AF1-like subtype similarly spanned pediatric and adult AML genetics and consistently predicted poor prognosis. Using long-read single-cell RNA sequencing, we confirmed that discovered U2AF1-like splicing was shared across cell states, co-opted a healthy circadian gene program, was stable through relapse, and induced a leukemic stem cell program. Pharmacological inhibition of an implicated U2AF1-like splicing regulator, PRMT5, rescued leukemia missplicing and inhibited leukemic cell growth. Finally, genetic deletion of IRAK4, a common target of U2AF1-like and PRMT5 treatment, blocked leukemia development in xenograft models and induced differentiation. This work suggests that broad splicing dysregulation, in the absence of select mutations, is a therapeutic target in heterogeneous leukemias.

Figure 1.. Mutation-defined splicing is largely obscured in leukemia.

(A,B) Heatmap of top marker splicing events (A) and differentially expressed genes (B) in AML (Leucegene RNA-Seq) for a subset of patients and common splicing factor mutations/fusions (n=142, training). Columns=patients, Rows=events/genes. ΔPSI=relative difference in Percent Spliced In (PSI) values (median normalized) . C) Percentage of patients with AML that were accurately classified according to the indicated genetics based on splicing versus gene expression (n=200), based on 3-fold cross-validation (SVM, one versus rest). D) Heatmap of alternative splicing-patterns identified in Leucegene RNA-Seq, identified using a single-cell analysis clustering algorithm (ICGS). ΔPSI values are median normalized per event. E) Cartoon of the OncoSplice computational workflow to define new splicing subtypes and mechanisms of gene regulation from RNA-Seq. These steps consist of: splicing quantification, unsupervised subtype discovery, supervised subtype identification (genetics, multi-factor splicing event correlation) and RNA-regulatory splicing-subtype prediction based on RNA binding protein (RBP) expression, binding sites and CLIP-Seq data.

Figure 2.. OncoSplice uncovers genetically heterogenous subtypes of AML.

A) OncoSplice-defined AML subtypes with coincident cancer genomic variants in 367 adult AML samples (yellow=subtype assigned patient) (table S1). Previously defined AML subtypes (top panel) and new OncoSplice-defined subtypes (bottom panel), annotated for RNA-Seq-detected genomic variants, oncofusions, deletions or structural rearrangements (bold=splice-ICGS reported). For the top panel, final subtypes were revised according to known patient genetics. B) Heatmap of concordant splicing events between OncoSplice-defined subtypes. Hierarchical clustering of the percentage of overlapping splicing events between all pairs of AML subtypes (regulated in the same direction) are shown (black=high percentage, white=low). Clusters of samples dominated by U2AF1 (U2AF1-like or U2AF1 mutation), SRSF2 (SRSF2-like or SRSF2 mutation) or NPM1/FLT3-ITD are denoted by a box (right). C) For major OncoSplice-defined subtypes the number of differentially expressed genes (DEGs) and unique alternative-splicing events (AS) are shown. Subtypes are grouped into those principally defined by AS (left), AS and DEGs (middle) or DEGs (right). D) Splicing example: SashimiPlot of CASP9 splicing in a U2AF1-like and an SRSF2-like patient sample (top). SashimiPlot lines between exons indicate junctions and numbers indicate junction-read counts. The alternative splice event results in predicted CASP9 protein isoforms (bottom) including the pro-apoptotic long CASP9a isoform and the short CASP9b isoform, which lacks the peptidase domain (ExonPlot view AltAnalyze). E) Annotation of the frequency of MultiPath-PSI-defined splice-event types (defined below) associated with each AML subtype (denoted to the left). F) Annotation of the AltAnalyze-predicted impact of splice events on protein domain and protein length in each AML subtype (denoted to the left).

Figure 3.. U2AF1-like splicing partially phenocopies mutation engendered splicing dysfunction.

A) splice-ICGS (U2AF1-like signature) reveals broadly deregulated splicing in most patients with AML. The white boxes indicate 1) RNA-Seq samples with U2AF1-S34 mutations and U2AF1-like splicing and 2) samples with SRSF2-P95 mutations and SRSF2-like splicing. B). Heatmap showing significant splicing events (P-value <0.05, FDR adjusted and ΔPSI =0.1) in adult AML with splicing factor mutations (U2AF1-S34, SRSF2-P95, SF3B1, U2AF1-Q157). The coincidence of U2AF1-S34 and SRSF2-P95 splicing events with U2AF1-like and SRSF2-like, respectively, is illustrated (white boxes). C) Venn diagram displaying AML-subtype-associated splicing events in U2AF1-S34 and U2AF1-like or SRSF2-P95 and SRSF2-like. D) Weblogo analysis of U2AF1 binding-site preferences at the e-3 splice-site position for selective cassette-exon splicing events from U2AF1-S34 mutation or U2AF1-like in panel C. E) The number of cassette exon events included and excluded for selective U2AF1-S34 and all U2AF1-like events are shown for each binding site preference. F) Kaplan-Meier curves for overall survival in patients from TCGA AML (top) and TARGET AML (bottom) with associated coxph p-values (left: all splice-ICGS stringently classified U2AF1-like versus all other considered AMLs). Analysis of TCGA was restricted to cytogenetically normal AMLs with no RNA binding protein (RBP) mutations and under 60 years of age. G) Distribution of AML cell-line aggregate CRISPR-screen scores (CSS) of 287 genes corresponding to U2AF1-like splicing events (common pediatric and adult) that are associated with poor overall survival compared to all CSS genes. CSS is displayed on the y-axis and the density of CSS scores U2AF1-like splicing event genes versus all genes is displayed on the x-axis. A Wilcoxon rank sum p-value (two-sided) was computed for the comparison of CSS between these two groups. Exemplar splicing event genes which are required for leukemic growth in >6 cell lines are indicated. H,I) Select poor-survival associated in U2AF1-like and mutation-associated splicing events. Example sample SashimiPlots (H) and violin plots of PSI values for all patients (I). J) Violin plot displaying the PSI distribution for previously identified U2AF1-S34 or SRSF2-P95 splicing events in U2AF1- and SRSF2-like samples. in=inclusion exon, ex=exclusion exon.

Figure 4.. Single-cell long-read sequencing defines an intrinsic U2AF1-like splicing program in AML.

A) U2AF1-like splicing status in a prior relapse cohort (n=19 patients) of matched AML samples at diagnosis (blue dot) and relapse (pink triangle) (57). The U2AF1-like-score is calculated from the aggregate of 364 poor survival-associated U2AF1-like inclusion versus exclusion ΔPSI splicing event values for each patient. Classification of AML samples as U2AF1-like, SRSF2-like, or Other are based on the range of scores produced for patients from Leucegene assigned by splice-ICGS. Prior annotated epigenetic subtypes (eloci) are indicated (left). B) U2AF1-like splicing status in n=42 patients with AML from multi-timepoint sampling in the BEAT-AML RNA-Seq cohort. C) Experimental design of single-cell(sc) long-read (LR) RNA-Seq of AML samples with defined xenograft LSC activities. D) UMAP of scRNA-Seq profiles from primary AML samples projected into a comprehensive cell atlas of healthy bone marrow cell states using the software cellHarmony. E) Cell frequency of AML bone marrow from the LR-scRNA-Seq. F) Assignment of U2AF1-like splicing scores to each of the most frequently detected bone marrow cell-population clusters for each of the 3 patients with AML. G) Interquartile range (IQR) of protein length distributions comparing all monocytic progenitors from a U2AF1-like (AML-7) and SRSF2-like (AML-12) for unique differential mRNA isoforms protein products (ΔPSI>0.4). H) IQR of pairwise inferred protein isoforms length differences in U2AF1-like versus SRSF2-like monocytic progenitors for each differential splicing event (ΔPSI>0.4). IQR statistical differences are shown as Mann-WhitneyU test p-values.

Figure 5.. U2AF1-like splicing is mediated by PRMT5 and WDR77 expression.

A) Top-associated U2AF1-like splicing factors, comparing gene expression and splicing (Pearson correlation). Correlation values are shown in the upper right quadrants. Scatter plots (lower left quadrants) illustrate the pairwise expression value of the indicated factors. B) AML cell-line CRISPR-screen scores (CSS) for differentially expressed and knockdown (KD) of inferred regulatory RBPs in U2AF1-like or SRSF2-like samples compared to all RBPs. C) The extent of splicing concordance (similarity index) between Leucegene AML U2AF1-like splice events to RBP knockdown (KD) or over-expression (OE) (n=77) in the indicated cell lines. Reference differential splicing events were obtained from hundreds of collected RBP perturbation RNA-Seq datasets with matching controls. D,E) Heatmap of all AML U2AF1-like significant splicing events overlapping between (D) shRNA KD of U2AF1 in K562 cells (ENCODE) or (E) shRNA of WDR77 in MDA-MB-231. F) Concordance between Leucegene U2AF1-like splice events with cell-type and hematological malignancy specific splicing programs. G) Heatmap of Leucegene U2AF1-like significant splicing events overlapping with PRMT5 inhibitor treated MDS-L cell RNA-Seq exonic splicing events. H) SashimiPlot of the IRAK4 gene locus in PRMT5 inhibitor treated and control MDS-L cell RNA-Seq (prior annotated Exon 4 (denoted E6 in AltAnalyze)).

Figure 6:. Treatment with PRMT5 inhibitor PRT543 by targeting of IRAK4 leads to increased myeloid differentiation and impaired MDS/AML progenitor cell function.

A) Western blot of IRAK4-L (long isoform) in MDSL cells following 72 hours of treatment with the PRMT5 inhibitor (PRT543) at 30nM, 100nM and 300nM. B) NF-kB reporter activity in THP1-Blue NF-kB cells with increasing doses of the PRMT5 inhibitor (PRT543) for 24 hours. *** P – value = 5.8×10⁻⁵ (300nM), P – value = 1.0×10⁻⁵ (500nM), P – value = 5.3×10⁻⁶ (1000nM). C) Cell viability in MDSL cells with increasing concentrations of PRT543 50nM, 100nM and 300nM as compared to control. P – value = 4.1×10⁻⁵ (50nM D7), P – value = 4.1×10⁻⁷ (100nM D7), P – value = 3.8×10⁻⁸ (300nM D7), P – value = 4.1×10⁻³ (50nM D21), P – value = 3.2×10⁻⁹ (100nM D21), P – value = 1.5×10⁻¹⁰ (300nM D21). D) MDSL were treated with PRT543 300nM every 72 hours as compared to control and myeloid differentiation assessed on days 21. P – value = 7.2×10⁻³. E) Representative images of Giemsa stained MDSL cells alone and those treated with PRMT5 inhibitor (PRT543) at 300nM treatment for 72 hours. F) MDS patient samples (MDS1, MDS4 and MDS5) treated with PRMT5i and control for 14 days in clonogenic assays and then assessed for myeloid and erythroid differentiation by FACS. G) Colony formation in isogenic WT and IRAK4^KO THP1 and MDSL cell lines (two independent experiments). P – value=0.037 (left panel), P – value=0.00088 (right panel). H) Kaplan Meier survival analysis of NSGS (NOD/SCID mice with IL2RG knockout and transgenic cytokine expression) (n = 5 mice/group) engrafted with WT or IRAK4^KO THP1 or MDSL cells (Data represent one of two independent experiments with similar trends).