Proteogenetic drug response profiling elucidates targetable vulnerabilities of myelofibrosis

Abstract

Myelofibrosis is a hematopoietic stem cell disorder belonging to the myeloproliferative neoplasms. Myelofibrosis patients frequently carry driver mutations in either JAK2 or Calreticulin (CALR) and have limited therapeutic options. Here, we integrate ex vivo drug response and proteotype analyses across myelofibrosis patient cohorts to discover targetable vulnerabilities and associated therapeutic strategies. Drug sensitivities of mutated and progenitor cells were measured in patient blood using high-content imaging and single-cell deep learning-based analyses. Integration with matched molecular profiling revealed three targetable vulnerabilities. First, CALR mutations drive BET and HDAC inhibitor sensitivity, particularly in the absence of high Ras pathway protein levels. Second, an MCM complex-high proliferative signature corresponds to advanced disease and sensitivity to drugs targeting pro-survival signaling and DNA replication. Third, homozygous CALR mutations result in high endoplasmic reticulum (ER) stress, responding to ER stressors and unfolded protein response inhibition. Overall, our integrated analyses provide a molecularly motivated roadmap for individualized myelofibrosis patient treatment.

Fig. 1. Single-cell analysis of mutant CALR protein expression.

A Staining of a cell line panel with the CALR mutant-specific antibody CALRm. Representative images of the included cell lines are shown (left panels; blue: DAPI; red: CALRm). Boxplots (right) represent single-cell CALRm intensities for the panel as quantified by intracellular flow cytometry analysis (right panel). p-Values indicate Student’s t-test significance: **** = p < 0.0001. Exact p-values are reported in Source Data. Box plots indicate the median (horizontal line) and 25% and 75% ranges (box), and whiskers indicate the 1.5× interquartile range above or below the box. B Single-cell immunofluorescence imaging of MF PBMCs stained with CALRm. Representative images of different MF PBMC samples are shown (left panels; blue: DAPI; red: CALRm). Single-cell mean CALRm intensities of four replicate wells with 25 images per well are shown as a violin plot with percentages of thresholded CALRm-positive cells annotated as text (right panel). Representative results of two independent repeats are shown. Violin colors indicate the CALR VAF of the corresponding MF patient. See also Supplementary Fig. 1.

Fig. 2. Clinical and molecular determinants of cellular heterogeneity in MF.

A Outline of the integrative pharmacoscopy and clinical proteotyping workflow of MF patients as performed in this study. B Circos plot of clinical annotations of the MF PBMC cohort. Legend indicates color codes for the included variables. For discrete variables, patient numbers are included, and for continuous variables, the range and median are represented. Proteotyping outliers are annotated by H and T for HSPC and T-cell proteotypes, respectively. See also Supplementary Data 1. C Cell population sizes as defined by the CNN-based cell classifier across the DMSO-treated conditions of the MF PBMC cohort. Annotations indicate the mutation status and PB blast counts of the respective patients. D t-SNE embedding of the CNN class probabilities for DMSO-treated cells across the cohort of 43 MF patients (n = 41,286 patient cells; up to 250 cells randomly selected from HSPC, T-cell, monocyte, and other cells for each MF patient). Top left panel: CNN-based cell types. Top right panels: single-cell CALRm (upper) and pSTAT5 (lower) intensities, scaled by z-score normalization. Bottom panels: single-cell CD3/CD14 and CD34 intensities and nuclear area. E CALRm/pSTAT5 intensities per cell type and driver mutation. Intensities are scaled by z-score normalization per patient. Boxplots indicate the range of scaled intensities across MF CALR and JAK2 patients, and dots represent individual patient values (MF CALR, n = 14; MF JAK2, n = 26). Fill colors indicate cell types as in Fig. 1C. F CALRm intensities per cell type and homozygous (VAF > = 75%) or heterozygous MF CALR (VAF < 75%). Intensities are scaled by z-score normalization across all MF CALR cells (MF CALR heterozygous, n = 11; MF CALR homozygous, n = 3). Boxplots indicate the range of scaled intensities, and dots represent individual patient values. Fill colors indicate cell types as in Fig. 1C. Asterisks indicate non-adjusted two-sided Student’s t-test significance: **** = p < 0.0001, *** = p < 0.001, ** = p < 0.01, * = p < 0.05; exact p-values are reported in Source Data. Boxplots as in Fig. 1A. See also Supplementary Fig. 2.

Fig. 3. Pharmacoscopy elucidates drugs specifically depleting the malignant MF clone.

A Summary of pharmacoscopy-based drug responses, split by MF CALR (left; n = 14 patients) and JAK2 (right; n = 26 patients). The x-axis shows the signed significance of HSPC drug responses per patient, averaged across the genetically stratified subcohorts. The y-axis shows the signed significance of oncogenic CALRm and pSTAT5 drug responses per patient for MF CALR and JAK2, respectively, averaged across the genetically stratified subcohorts. p-values indicate Student’s t-tests of summarized drug responses across replicates and concentrations compared to the respective control condition (DMSO, PBS, or isotype control). Values indicate −log10 p-values, signed positively or negatively for on- or off-target response, respectively, averaged across the CALR MF and JAK2 MF cohorts. The dot color indicates drug class; the dot size is the fraction of patients that have a significant on-target effect averaged across the two drug scores. B Association of MF patient drug responses with clinical factors. For every drug and for both drug response readouts, an ANOVA was performed for selected factors. Significant drug response-factor associations (ANOVA p < 0.05) were counted (x-axis) per factor (y-axis) and are displayed as either associating with sensitivity (blue) or resistance (red). C Protein pathway-level associations to MF patient drug responses. For every drug, the HSPC drug responses were correlated (Spearman) to HSPC protein levels across patients. GSEA was performed on the ranked correlations, calculating the enrichment of positive or negative protein–drug response correlations. Significant GSEA enrichment (p < 0.05) was counted (x-axis) per KEGG pathway term (y-axis), and the top 15 positives (blue; left panel) and negative (red; right panel) enriched terms are displayed. See also Supplementary Figs. 3 and 4.

Fig. 4. CALR driver mutations sensitize MF HSPCs to BET inhibition.

A BET and HDAC inhibitor HSPC responses across the cohort split by driver mutation (MF CALR, n = 14 patients; MF JAK2, n = 26 patients). Responder = significant (p < 0.05) on-target response, non-responder = non-significant on-target, resistant = off-target response. B Signed ANOVA significance of BET and HDAC inhibitor HSPC drug responses to clinical parameters. C BET inhibitor drug responses of a co-cultured CALR CRISPRed cell line panel analyzed by pharmacoscopy. A representative experiment of two is shown, with two technical replicates across two conditions with 25 images per well. The y-axis depicts the change in relative fraction per cell line in the drug-treated condition compared to DMSO. p-Values indicate significance as determined by the Student’s t-test of CALR MUT depletion across the different replicates and concentrations compared to DMSO-treated cells. D Analysis of protein processes significantly up- or downregulated in BET (left) and HDAC inhibitor (right) responding patients. t-Tests were performed for each protein, comparing protein levels in HSPCs of responding MF patients compared to the non-responding and resistant patients. GSEA was performed on proteins ranked according to signed t-test significance. GSEA normalized enrichment scores (NES; x-axis) and enrichment significance (y-axis) are shown. E HSPC protein levels of selected Ras signaling pathway members. Patients are grouped according to HSPC response to BET and HDAC inhibitors (Responder, n = 28; resistant/non-responder, n = 22). Protein expression (y-axis) of proteins shown is significantly different for responders for at least 3 BET/HDAC inhibitors. p-Values indicate the Student’s t-test significance. Boxplots as in Fig. 1a. F Analysis of protein process–drug response associations of Bcl-2 inhibitors as in (D). Asterisks indicate non-adjusted two-sided significance: **** = p < 0.0001, *** = p < 0.001, ** = p < 0.01, * = p < 0.05; exact p-values are reported in Source Data. See also Supplementary Figs. 3 and 4.

Fig. 5. Large-scale clinical proteotyping elucidates the MF-specificity of targetable protein alterations.

A Workflow for the proteotyping of myeloproliferative neoplasm (MPN) samples. Proteotypes of granulocytes isolated from healthy donors (HD; n = 31) and a clinical cohort of essential thrombocythemia (ET; n = 42) and MF patients (n = 40) from two different medical institutions were analyzed using DIA MS/MS. HD are age- and gender-matched to the MPN cohort. B Circos plot of clinical annotations for the HD/MPN granulocyte cohort. Legend indicates color codes for the included variables. For discrete variables, patient numbers are included, and for continuous variables, the range and median are represented. Extended clinical annotations can be found in Supplementary Data 1. C t-SNE embedding of granulocyte proteotypes. Dimensionality reduction is based on all quantified proteins (n = 4056). MPN patients and HD are labeled by disease and mutation status. D Proteotype similarities between HD and MPN subcohorts. Boxplots show distributions of Euclidean distances between proteotypes of selected subcohorts (self-comparisons were excluded). Asterisks indicate non-adjusted two-sided Student’s t-test significance: **** = p < 0.0001; exact p-values are reported in Source Data. Boxplots as in Fig. 1A. E Workflow overview of machine learning-based approach to identify the minimal protein signature that classifies HD and MPN patients by disease and mutation status. The dataset is split into disjoint test, train, and validation cohorts, after which a recursive feature elimination (RFE) model is applied to prioritize the proteins most discriminatory for disease and mutation status. Based on the selected proteins, a multi-layered perceptron classifier is trained, of which the scores and selected proteins are reported. The final protein signature is derived from the most frequent scoring proteins. F Protein selection frequencies across 100 RFE runs. The main panel shows protein selection frequencies (y-axis) are shown for all proteins selected at least once (x-axis). Insert shows the same data for the most frequently selected proteins. The gray box indicates the 15-protein signature of proteins selected in at least 85 runs. G t-SNE embedding of the 15-protein signature for MPN patients and HD. MPN patients and HD are labeled by disease and mutation status. See also Supplementary Figs. 4 and 5.

Fig. 6. Protein and pathway-level alterations linked to disease and driver mutation status.

A Volcano plots of comparisons underlying the classification of the RFE-selected protein signature (n = 15). Volcanoes indicate log2FC (x-axis) and significance (y-axis) of Student’s t-test results comparing granulocyte protein levels for the indicated groups. Members of the 15-protein signature are highlighted and colored by reaching significant differences only in the HD vs. MPN comparison (MPN general alterations), in the MF vs. ET comparison (Disease-specific alterations), or also in the MF CALR vs. MF JAK2 comparison (Mutation-specific alterations). B Association network of protein and pathway-level alterations linked to disease and driver mutation status. For every protein quantified in the granulocyte proteotypes, expression levels were associated with general disease (MPN vs. HD), diagnosis (MF vs. ET), and mutation status (CALR vs. JAK2) by ANOVA analysis. For all clinical factors, the top 20 significantly associated proteins with an adjusted ANOVA p-value < 0.001 are shown. The edges connecting these proteins to factors represent the corresponding signed ANOVA significance. For a pathway-level view, all signed ANOVA p-values were ranked per factor, after which GSEA enrichment analyses were performed. Significantly enriched KEGG pathway terms are linked to the respective factors with edges representing signed GSEA p-values (adjusted p < 0.05). For discrete clinical factors, red edges indicate association with the first term and blue with the second term. For VAF, red edges indicate positive and blue negative association. Orange proteins indicate those proteins within the network that are also members of the 15-protein signature. Green edges represent STRING-defined physical protein-protein interactions (physical interaction STRING score > 0.7). Brown and gray backgrounds indicate subnetworks of interest that are experimentally followed up on. See also Supplementary Fig. 5.

Fig. 7. Proliferative MCM-high MF responds to the targeting of survival signaling and DNA replication.

A Correlation of MCM4 and MCM7 protein levels across the granulocyte cohort. Levels normalized by subtracting the median HD expression. B ANOVA significance of averaged MCM4 and MCM7 protein level associations with clinical parameters across the MF cohort. The red color indicates a positive association with continuous parameters (PB blasts and age) or to the first term of discrete factors (mutation, sex, and treatment). C Z-score normalized MCM4 and MCM7 protein levels across MF patients stratified by PB blast counts (<1%, n = 64; 1–5%, n = 39; >5%, n = 6). p-Values indicate the Student’s t-test significance. D Protein–protein interaction network of the top-50 proteins correlating to MCM4 and MCM7 protein levels across T-cell, HSPC, and granulocyte proteotypes. Edges: STRING score > 0.7. Proteins colored by average Pearson correlation. E Pathway enrichment analysis results of the 50-protein signature. Significant positive KEGG term enrichments (p < 0.05) are shown. F Ranking of averaged MCM4 and MCM7 protein levels per cell type and MF patient for the 20 MF patients present in both the granulocyte and MF PBMC cohorts. G Quantification of the percentage of MCM7, Ki-67, and yH2Ax positive PBMCs. Dots summarize the staining quantification of 2–4 patients (brackets) per MCM group averaged across three replicate wells and 25 images per well. Bars indicate the average per group. H Representative images of in vitro HSPC DNA replication fork fibers (left panel) and corresponding quantification (right panel). Dots summarize average fiber lengths of HSPCs isolated from 1 to 6 patient samples (brackets) per MCM group with 93–370 fibers per sample. Bars indicate the average per group. I On-target HSPC drug responses correlated to the average expression of the 50-protein replicative signature. x-axis, Pearson correlations; y-axis, corresponding significance. J Scatterplots of the three highest correlating drug responses of (I). Linear fit (blue line) and 95% confidence intervals (gray area) are shown. Pearson correlation and significance are indicated. K Vosaroxin responses of FACS-isolated CD34 + MF patient HSPCs from three vosaroxin-sensitive patients. Asterisks indicate non-adjusted two-sided significance: **** = p < 0.0001, *** = p < 0.001, ** = p < 0.01, * = p < 0.05; exact p-values are reported in Source Data. Boxplots as in Fig. 1A. See also Supplementary Fig. 6.

Fig. 8. Homozygous CALR mutations lead to ER stress-associated vulnerabilities.

A Gene set enrichment analysis (GSEA) of ranked correlations between granulocyte protein levels and CALR/JAK2 granulocyte VAF. KEGG terms with an adjusted enrichment p < 0.0001 are shown. B Protein–protein interaction network of proteins significantly correlating to CALR VAF. Proteins with significant Pearson correlations are shown (p < 0.01). Edges: physical interaction STRING score > 0.7. Node color: Pearson correlations; Outer rings: selected protein pathway memberships. Blue shading indicates the interaction cluster of CALR and associated ER proteins. C Heatmap of ER stress cluster protein expression levels across MF patient granulocytes. Patient VAF and protein class are indicated. D Volcano plot of MF CALR homozygous versus all other MF CALR proteotypes. Proteins are colored green and annotated if significantly different and previously described to be upregulated in ER-stressed cell lines. E Heatmap of ER stress cluster protein expression levels across cell types and cell lines. Protein levels normalized as log2FC as indicated in the legend. F Single-cell expression levels of CALRm correlated with HSPA5 across PBMCs of two MF CALR patients. Representative images (left panel) and quantification (right panel) of four replicate wells with 25 images per well are shown. Representative results of two independent repeats are shown. Pearson correlation and significance are reported. G Significance of depletion of CALR mutant expressing cells by ER stress inducers and UPR inhibitors across two MF CALR patient samples. The dashed line represents the diagonal. H Comparison of oncogenic drug responses between MF CALR homozygous and either MF CALR heterozygous (x-axis) or MF JAK2 (y-axis). Axes indicate an averaged difference in signed −log10 significances. Selected drugs are labeled and colored by drug class. I Oncogenic drug responses correspond to selected MF CALR homozygous-specific drug responses shown in (H). Boxplots indicate the range of drug responses across the MF cohort, dots represent individual patients. J Carfilzomib responses of FACS-isolated CD34+ MF patient HSPCs. Boxplots represent log2FC differences between DMSO- and drug-treated conditions. Asterisks indicate non-adjusted two-sided Student’s t-test significance: **** = p < 0.0001, ** = p < 0.01, * = p < 0.05; exact p-values are reported in Source Data. Boxplots as in Fig. 1A. See also Supplementary Fig. 7.