Integrative multi-omics and drug response profiling of childhood acute lymphoblastic leukemia cell lines

Abstract

Acute lymphoblastic leukemia (ALL) is the most common childhood cancer. Although standard-of-care chemotherapeutics are sufficient for most ALL cases, there are subsets of patients with poor response who relapse in disease. The biology underlying differences between subtypes and their response to therapy has only partially been explained by genetic and transcriptomic profiling. Here, we perform comprehensive multi-omic analyses of 49 readily available childhood ALL cell lines, using proteomics, transcriptomics, and pharmacoproteomic characterization. We connect the molecular phenotypes with drug responses to 528 oncology drugs, identifying drug correlations as well as lineage-dependent correlations. We also identify the diacylglycerol-analog bryostatin-1 as a therapeutic candidate in the MEF2D-HNRNPUL1 fusion high-risk subtype, for which this drug activates pro-apoptotic ERK signaling associated with molecular mediators of pre-B cell negative selection. Our data is the foundation for the interactive online Functional Omics Resource of ALL (FORALL) with navigable proteomics, transcriptomics, and drug sensitivity profiles at https://proteomics.se/forall .

Fig. 1. Multi-omics profiling of childhood ALL cell lines.

a Schematic workflow of childhood ALL cell lines data generation and analyses conducted in this study. Cells were profiled for protein level using mass spectrometry and mRNA expression using RNA-seq. The majority of the cell lines, excluding derivative or sister cell lines (n = 43) were scanned for their sensitivity to 528 investigational and oncology drugs. b An overview of the most prominent molecular features of the 51 cell lines in this panel. BM bone marrow, PB peripheral blood, PE peripheral effusion, DSRT drug sensitivity and resistance testing. c Principal component analysis (PCA) of the highly variable proteins in the 51 cell lines showed distinct clusters corresponding to BCP-ALL (green), immature B-ALL (blue), T-ALL (pink), and EBV-transformed B cells (gray). d Cell lines grouped by the most prominent features of the described cytogenetic subtypes using proteomics data. BR biological replicate. Source data are provided as a Source Data File 1.

Fig. 2. Synergies and discordances of mRNA and protein levels in childhood ALL cell lines.

a Density plot of Spearman rank order of gene symbol-wise correlation of the mRNA-protein levels (n = 10,738 RNA-protein pairs), with RNA or protein markers stratified by gene symbol. The mRNA-protein levels were positively correlated for 97% of the pairs across the entire matching panel (n = 64 samples), with 78% of the pairs showing significant positive correlation (FDR ≤0.01) assessed by t-distribution adjusted P value. b GSEA of the mRNA-protein level correlations for KEGG pathways (n = 162 ranked pathways) showing that certain specialized signal transduction pathways were enriched among the highly correlating pairs, while pathways associated with the spliceosome, proteasome, and ribosomes belonged to the lowest mRNA-protein correlating pairs. Selected pathways are annotated with blue circles. c Comparison of random pairwise correlations (blue) of quantitative proteomics data and mRNA levels across the 49 ALL cell lines and their replicates to known (yellow) interaction pairs (n = 31,842) from the CORUM database demonstrating higher correlations at the protein level. The statistical difference was assessed using an unpaired two-sided t-test, with P values (p) indicated. The center line in the violin plots represents the median, the vertical bars represent the first and third quartiles and the curve represents the density curve. d Relative log2 ratio protein levels versus log2 TPM (Transcripts Per Million) mRNA levels of Go-Ichi-Ni-San (GINS) complex members across 49 cell lines and their biological replicates. GINS mutations (amino acid substitution) for GINS-complex members are shown for annotated cell lines MOLT-13, KARPAS-45, and REH. Fitted lines were plotted using a linear fit. Dotted lines represent the difference between the actual protein levels and the mRNA-predicted levels. Correlation P values are derived from two-sided t-distribution. e Ranked mRNA–protein correlations (n = 10,738 ranked mRNA-protein pairs) with selected highlighted markers (blue circles) that are frequently involved in or associated with leukemogenesis. Source data are provided as a Source Data File 2.

Fig. 3. MS-based quantification and clustering of the proteome of childhood ALL cell lines.

a Heatmap showing the hierarchical clustering of the total overlap of 9100 identified and quantified proteins across 51 cell lines. b GSEA of the differentially expressed genes in CLC3 versus the remaining T-ALL cell lines using the proteomics and transcriptomic data showing significant enrichment of the spliceosome in the proteomics data (NES = 2.1, adjusted P value q = 5.5e-4) but not in the transcriptomic data (NES = 1.06, adjusted P value q = 0.73). Only one biological replicate from the proteomics data was used to have matching numbers and an unbiased comparison to the transcriptomics data. The significance of NES was assessed by Kolmogorov–Smirnov statistics. c Dendrogram of the consensus clustering of the B-lineage and T-lineage separately showing five CLC for B-lineage (B)-CLC) and seven CLC for T-lineage (T-CLC) cell lines. d Sankey diagram, showing connection nodes between proteomics CLCs and mRNA-based CLCs, by cytogenetic subtype. The approximately unbiased bootstrapping probability (pr) are indicated for each cluster. e Volcano plot of differentially abundant proteins in the KMT2A-AFF1 cell lines compared to the KMT2A-MLLT1 cell lines (top panel) using DEqMS and the TP53 levels in KMT2A-AFF1 cell lines compared to the KMT2A-MLLT1 cell lines (bottom panel) using two-sided, unpaired t-test, P value (p) is shown in the plot. Source data are provided as a Source Data File 3.

Fig. 4. Drug sensitivity and drug target correlations of childhood ALL cell lines.

a Heatmap depicting the sDSS of the 528 drugs from the DSRT across the 43 tested ALL cell lines. The x-axis is the cell lines ordered by the rank from hierarchical clustering (Pearson ward.D2) and the colors are the individual sDSS for each drug. The legends indicate the CLC and cytogenetic subtype of the cell lines. Higher sDSS indicates a more potent and effective drug relative to normal bone marrow cells. b Scatter plot of the median versus variance of the sDSS for each drug across the tested cell lines. Selected drugs are highlighted and colored according to their drug class. c The correlation between the sDSS for the glucocorticoids in the DSRT and relative protein level of the glucocorticoid receptor, NR3C1 for the 43 tested ALL cell lines. The x-axis shows the relative log2 protein level of NR3C1 and the y-axis shows the sDSS for each respective cell line (n = 43). The Pearson correlation coefficient (R) and two-sided t-distribution P values (P) for the comparison are shown in each plot. The linear regression trendline (black) and its 95% confidence interval (shaded gray area) are shown in the graph. The colors indicate the cytogenetic subtype of the cell lines, as annotated in 4e. d Scatter plot of the Spearman correlation versus rank for the sDSS of each drug and the relative log2 protein level of its putative protein target. The drugs are colored by drug classes and selected drug–drug target pairs are highlighted. e Ranked Pearson correlations of tacrolimus sDSS and highlighted proteins of the FKBP family (purple circles) and a scatter plot of tacrolimus sDSS and FKBP10 levels for each of the 29 tested ALL cell lines where FKBP10 was quantified. The x-axis shows the relative log2 protein level of the protein and the y-axis shows the sDSS for each respective cell line (n = 29). The Pearson correlation coefficient (R) and two-sided t-distribution P values (P) for the comparison are shown in each plot. The linear regression trendline (black) and its 95% confidence interval (shaded gray area) are shown in the graph. The colors indicate the cytogenetic subtype of the cell lines. Source data are provided as a Source Data File 4.

Fig. 5. Pharmacoproteomics identifies drug mechanisms of action in childhood ALL cell lines.

a UMAP of the distribution of class-annotated drugs as a two-dimensional reduction. Drugs were correlated to gene symbol-based protein quantification data based on sDSS, and Pearson correlation scores from this analysis were used to generate the UMAP. Each drug is annotated by its class. b UMAP dimensional reduction of each drug, calculated using Pearson correlations between sDSS and protein abundance. Each drug is annotated by a specific target family. c sDSS of the tacedinaline (top panel) and panobinostat (bottom panel) across the tested cell lines. The tacedinaline resistant TCF3-PBX1 fusion cell lines are highlighted with red circles. The red dashed line indicates the selected threshold of sDSS = 8. d CETSA derived T_agg (aggregation temperature) curves for HDAC1 in KOPN-8 and RCH-ACV cells in the presence of DMSO and 20 μM tacedinaline after 3 h incubation. All data were normalized to the response observed for each treatment condition at the lowest test temperature (top panel). Dose-response CETSA of tacedinaline incubated for 3 h at seven different concentrations ranging from 100 μM to 24 nM at 53 °C based on raw chemiluminescence data from the western blot analysis for HDAC1 in KOPN-8 and RCH-ACV cell lines (bottom panel). Data were provided as two individual data points from two independent experiments (n = 2) for each cell line. e Volcano plots of differential drug sensitivity between B-lineage and T-lineage ALL cell lines. The cut-off was set at adjusted t-test P value q ≤ 0.02 and the delta mean sDSS between the B-lineage and T-lineage cell lines set at ≥5. Selected drugs are annotated as blue circles. f The correlation between the sDSS for venetoclax in the DSRT and relative protein level of BCL2 (top panel) and sDSS for A-1155463 and relative protein level of BCL2L1 (bottom panel) for the 43 tested ALL cell lines. The x-axis shows the relative log2 protein level of the protein and the y-axis shows the sDSS for each respective cell line (n = 43). The Pearson correlation coefficient (R) and two-sided t-distribution P values (P) for the comparison are shown in each plot. The linear regression trendline (black) and its 95% confidence interval (shaded gray area) are shown in the graph. The colors indicate the cytogenetic subtype of the cell lines. g Ranked OSU-03012 sDSS and protein Pearson correlations, highlighting members of the TCP1 containing chaperonin complex, HSP90AB1, and PDPK1 interacting proteins from the STRING network database. h Scatter plot of OSU-03012 sDSS and relative log2 protein level of HSP90AB1 and TCP1 for each of the BCP-ALL cell lines (n = 25). The x-axis shows the relative log2 protein level of the protein and the y-axis shows the sDSS for each respective cell line (n = 25). The Pearson correlation coefficient (R) and two-sided t-distribution P values (P) for the comparison are shown in each plot. The linear regression trendline (black) and its 95% confidence interval (shaded gray area) are shown in the graph. Points are labeled by cytogenetic subtype of each represented BCP-ALL cell line. Source data are provided as a Source Data File 5.

Fig. 6. Selective sensitivity of MEF2D-HNRNPUL1 cell lines to bryostatin-1.

a Fold changes of differentially expressed mRNA using edgeR in 20 clinical patient samples with different MEF2D-rearrangements (MEF2D-BCL9 n = 13, MEF2D-HNRNPUL1 n = 4, MEF2D-CSF1R n = 1, MEF2D-SS18 n = 1, and MEF2D-DAZAP1 n = 1) plotted against differentially abundant proteins in the MEF2D-HNRNPUL1 cell lines (n = 3) using DEqMS. Only differentially expressed mRNA that passed a P value ≤0.01 (two-sided t-test) are shown in the scatter plot with selected highlighted genes. N indicates the number of genes passing this criterion. R Pearson correlation coefficient. P two-sided t-distribution P value. b Volcano plot of the differential protein levels in the MEF2D-HNRNPUL1 cell lines and their replicates (n = 6) versus the rest of the BCP-ALL cell lines and replicates (n = 37) using DEqMS (spectra counts adjusted t-test). The cut-off was set at P ≤ 0.05 and the fold change was set to log2(1.5). Selected proteins are highlighted with blue circles. c The sDSS of bryostatin-1 in the MEF2D-HNRNPUL1 cell lines (n = 4) compared to the sDSS for the remaining BCP-ALL cell lines (n = 24). The two-sided t-test P value is indicated as (P). d Ranked bryostatin-1 sDSS and protein level Pearson correlations with selected highlighted proteins in purple. e Western blot of ERK1/2 and pERK1/2 (Thr202 and Tyr204) of MEF2D-HNRNPUL1 fusion cell lines LC4-1, P30-OHKUBO (P30), KASUMI-7 (Kas-7), and KASUMI-9 (Kas-9) after 2 h treatment with 100 nM bryostatin-1. f Viable cell quantification normalized to corresponding mean DMSO viable cell count of the four MEF2D-HNRNPUL1 fusion cell lines, treated with 100 nM of bryostatin-1 or 25 nM PMA. An equal volume of DMSO was used as a control (n = 3). BCP-ALL cell lines ALL-PO, REH, RCH-ACV, NALM-6, SUP-B15, COG-LL-355h, COG-LL-394h, and MHH-CALL-2, lacking the MEF2D-HNRNPUL1 fusion, were subjected to the same treatments. Viable cells were quantified by flow cytometry, excluding zombie aqua dyed non-viable cells. Results are merged from five independent experiments. P values were obtained from unpaired two-sided t-tests. g Viable cell quantification normalized to corresponding mean DMSO viable cell count of the four MEF2D-HNRNPUL1 fusion cell lines, treated with 100 nM of bryostatin-1 alone or in combination with 1 uM MEK inhibitors UO126, trametinib, or selumetinib. Alternatively, to block ERK directly, 1 uM ERK inhibitor ERK 11e was used, and equal volume of DMSO was used as a control (n = 3). Viable cells were quantified by flow cytometry, excluding zombie aqua dyed non-viable cells. Results are merged from three independent experiments. P values were obtained from unpaired two-sided t-tests. Source data are provided as a Source Data File 6.