Large collections of genome-wide data can facilitate the characterization of disease states and subtypes, permitting pan-cancer analysis of molecular phenotypes and evaluation of disease context for new therapeutic approaches. We analyzed 9,544 transcriptomes from more than 30 hematologic malignancies, normal blood cell types, and cell lines, and showed that disease types could be stratified in a data-driven manner. We then identified cluster-specific pathway activity, new biomarkers, and in silico drug target prioritization through interrogation of drug target databases. Using known vulnerabilities and available drug screens, we highlighted the importance of integrating molecular phenotype with drug target expression for in silico prediction of drug responsiveness. Our analysis implicated BCL2 expression level as an important indicator of venetoclax responsiveness and provided a rationale for its targeting in specific leukemia subtypes and multiple myeloma, linked several polycomb group proteins that could be targeted by small molecules (SFMBT1, CBX7, and EZH1) with chronic lymphocytic leukemia, and supported CDK6 as a disease-specific target in acute myeloid leukemia. Through integration with proteomics data, we characterized target protein expression for pre-B leukemia immunotherapy candidates, including DPEP1. These molecular data can be explored using our publicly available interactive resource, Hemap, for expediting therapeutic innovations in hematologic malignancies.
This study describes a data resource for researching derailed cellular pathways and candidate drug targets across hematologic malignancies.
Gene expression profiles facilitate genome-wide analyses that can stratify patient subtypes and identify the activity patterns of various cellular pathways under different biological conditions (1, 2). Even though a large number of studies have characterized hematologic malignancies and normal blood cell types at genome-wide level since the introduction of microarray technology, most include only tens to hundreds of samples and focus on one disease. Thus, understanding the complete heterogeneity and similarity of diseases states and their subtypes remains an open challenge. Moreover, many hematologic malignancies are rare on the population level, necessitating collecting data across study cohorts.
Hematologic malignancies include acute and chronic leukemias of myeloid and lymphoid lineage, B-, T-, and NK-cell lymphomas, and multiple myeloma, and a number of premalignant conditions such as myelodysplastic syndrome (MDS), and myeloproliferative neoplasms (MPN). These diseases have highly variable genetic features, unique clinical courses, and varying therapeutic approaches. There is also a marked difference in prevalence, genetic background, and prognosis between adult and pediatric blood cancers. In children, acute lymphoblastic leukemia (ALL) is the most common hematologic malignancy, while in adults, non-Hodgkin lymphomas (NHL), followed by multiple myeloma, chronic lymphocytic leukemia (CLL), and acute myeloid leukemia (AML) are the most common. Treatment is moving toward increased utilization of targeted therapies in combination with traditional chemotherapies. Targeted therapies include tyrosine kinase inhibitors such as those developed against BCR-ABL fusion found in CML and some ALL cases, or antibody therapies including CD38 targeting in multiple myeloma, and engineered CAR-T cells recognizing cell surface CD19 or CD22 antigens in relapsed ALL and NHL (3–5). Yet, current therapies to treat hematologic malignancies rely heavily on drugs that target DNA metabolism in actively proliferating cells or intracellular signaling events that are involved in proliferation (6). Although these drugs have markedly improved progression-free survival, redundancy in signaling and the failure to eradicate quiescent cells (7) can facilitate the rapid development of therapy resistance. Testing a wider portfolio of new drug targets or repurposing drugs with established clinical indications represents promising strategies (6, 7). Molecular profile–guided approaches hold promise to improve the efficiency of this process (8).
We present here a resource that organizes samples from patients with cancer, healthy donors, and those at premalignant stages for comparative analysis based on both curated annotations and data-driven clustering of molecular phenotypes. This hematologic pan-cancer analysis permits the identification of clinically relevant molecular features and the exploration of new drug-targeting approaches across the disease hierarchy. The data and analysis tools are made available as an interactive online resource Hemap (http://hemap.uta.fi/) that synthesizes the curated genome-wide data across different disease subtypes.
Materials and Methods
Dataset retrieval and extraction of sample annotation data
Transcriptome datasets for Hemap were retrieved from the NCBI Gene Expression Omnibus (GEO) database (9) and represent samples hybridized to hgu133Plus2 genome-wide microarrays. The metadata were retrieved on the basis of matching disease ontology terms for hematologic malignancies against sample annotations (R/Bioconductor GEOmetadb package, “gsm” and “gse” tables), followed by manual curation, resulting in 10,470 samples. Refer to Supplementary Methods for details. Eight leukemia types, 8 B-cell lymphoma types, 7 T/NK lymphomas, multiple myeloma, and 4 proliferative disorders are represented by primary patient samples, while in total 36 disease types are included considering also their subcategories and cell line data (Supplementary Tables S1 and S2).
Data preprocessing and quality control
Samples with a median of raw probe intensity distribution deviating more than 1.5 in log2 scale from the median of all medians were deemed outliers and filtered out as well as those with an interquartile range deviating more than 0.75 from the median of interquartile ranges. Finally, duplicate samples, as well as all disease types with less than three samples (and samples assigned to these), were removed, resulting in 9,544 samples that were processed using the RMA probe summarization algorithm (10) with probe mapping to Entrez Gene IDs (from BrainArray version 18.0.0, ENTREZG). Finally, we employed the bias-correction method (11) to correct for any remaining technical differences (Supplementary Fig. S1). BeatAML (12) clinical and mutation data were obtained from source data file 41586_2018_623_MOESM3_ESM.xlsx. RNA sequencing (RNA-seq) counts were obtained from the authors. Genes with over 1 cpm expression in over 1% of the samples were kept. Data were normalized using limma voom and quantile normalization.
Dimensionality reduction methods are unsupervised methods that use measures of (dis)similarity and an optimization strategy to return as output sample coordinates in a lower dimension. Metrics of continuity, trustworthiness, and k-NN error were used to assess how well the visualization in 2D preserved their relative placement in the original coordinate space. We tested Gaussian Process Latent Variable Model (GPLVM; ref. 13), Locally Linear Embedding (LLE; ref. 14), Principal Component Analysis (PCA; ref. 15), Probabilistic Principal Component Analysis (PPCA; ref. 16), Sammon Mapping (SM; ref. 17), and t-Distributed Stochastic Neighbor Embedding (t-SNE; ref. 18; see Supplementary Methods for parameters). Comparison of the different methods encouraged the selection of t-SNE maps, specifically the Barnes–Hut approximated version of t-SNE implementation (BH-SNE; ref. 19). In final analysis, 15% genes with highest variance were used in construction of t-SNE maps (see ref. 20 for justification).
Assignment of cluster centers on t-SNE maps
Kernel density–based clustering algorithm (mean-shift clustering with bandwidth parameter set to 2.5, LPCM-package in R), was used to cluster the data following the dimensionality reduction. This method allows the discovery of sample sets that share similar features without having to prespecify the number of clusters. The term “cluster” is used in the text to refer to this computational clustering result, and the term “group” is used in context of visual examination. Pairwise statistical analysis between different sample features and clusters was performed as in ref. 21, based on Spearman correlation and the Bonferroni method for multiple hypothesis testing correction (see Supplementary Methods for details).
Discretizing gene expression with mixture models
Microarray hybridization generates background signal, which we would like to distinguish from real expression signal. The large sample size of Hemap was leveraged for fitting gene-specific models to cluster the gene expression in two components (expressed and not detected, denoted by 1 and −1, respectively). Gaussian finite mixture models were fitted by expectation-maximization algorithm (R package mclust version 4.3). If the uncertainty value from the model was more than 0.1, we assigned a value of 0 to denote low level. In addition, each log2 expression value lower than 4 was assigned a value −1 and values higher than 10 a value of 1. This was done to assure minimal amount of misclassifications of data samples to wrong components. The model was chosen by fitting both equal and variable variance models and ultimately choosing the model that achieved a higher Bayesian information criterion (BIC) to avoid overfitting. For drug target analysis, we utilized an adjustment for genes where background distribution was not found (gene is always expressed), or if over 90% of the samples had uncertain expression based on the model classification. Expressed state was assigned when >60% of the uncertain samples had expression above 6. Not detected status was reevaluated similarly (60% at level below 6).
Gene set analysis
The pathway and gene set enrichment analysis available in the Hemap resource was generated on the basis of gene sets retrieved from MsigDB v5.0 (molecular signatures; ref. 22), Wikipathways (06.2015; ref. 23), Recon 1 (metabolic pathways; ref. 24), Pathway Commons 7 (25), and DSigDB v1.0 (drug targets; ref. 26). Gene sets were limited to contain between 5 and 500 expressed genes (as defined above) per gene set, resulting in 19,680 gene sets that were evaluated across the dataset. The gene set variation analysis (GSVA; ref. 27), GSVA package 1.13.0 in R, was used to calculate gene set enrichment scores (positive for increased and negative for decreased expression) for each sample (parameters mx.diff = F, tau = 0.25, rnaseq = F). Significance was evaluated on the basis of empirical P values calculated using 1,000 random permutations of genes within the gene set, separately for gene set sizes 5–20, 25, 30, 40, 50, 75, 100, 200, 300, 400, and 500 to correct for differences in gene set sizes. Hypergeometric test was used to calculate enrichment of significant scores in a specific cluster.
Data sources used for evaluating drug targeting approaches
Drugs in clinical trials for leukemias, lymphomas, or multiple myeloma were obtained from ClinicalTrials.gov (accessed March 7th, 2018) maintained by the NIH (Bethesda, MD), including ongoing and terminated trials. Leukemia clinical trials were further sorted to those with clinical indication associated with AML, pre-B-ALL, CML, CLL, or multiple leukemia types. Drugs in use based on approved status in Finland were provided by the Finnish Pharmaceutical Information Centre Ltd and drugs approved by the FDA for leukemia, lymphoma, and myeloma were queried from FDA website (fda.gov–Drugs–Information on Drugs; Supplementary Table S3). A table of gene-level details for each drug was obtained from DSigDB (DSigDBv1.0 Detailed.txt; ref. 26) and integrated to Hemap in silico drug screening analysis. The list of drugs targeting epigenetic modifiers is based on the gene list with 124 genes available from ChEMBL_20 Target Classification Hierarchy (Supplementary Table S4; ref. 28). Analysis using TTD (Therapeutics Targets Database; ref. 29), DGIdb3.0 (30) for FDA-approved drugs across a wider disease context (31) as a source database was based on a total of 11,373 unique drugs and 1,270 unique genes. Drugs in use and in clinical trial included high confidence targets that were reported in several databases or had an associated PubMed identifier. A surface marker gene list with total of 996 genes was obtained from Cell Surface Protein Atlas (32) to evaluate putative novel immunotherapy targets.
Drug target in silico analysis in hierarchical manner
A disease hierarchy: (i) All disease samples; (ii) disease combinations; (iii) leukemia, lymphoma, myeloma; (iv) AML, pre-B-ALL, T-ALL, CLL; (v) disease clusters; was used to evaluate disease or subtype-specific drug target expression. Statistical significance of binary feature enrichment (e.g., high expression state) in a particular sample group was first evaluated using the hypergeometric enrichment test, followed by Bonferroni adjustment of P values. If >90% of the samples had high expression for a gene in the disease context, Inf score was assigned instead of −log10 P value (hypergeometric test would not be meaningful if the sample size was close to the whole population). Each significant gene was uniquely assigned to the disease group with the lowest P value. In the case of equal P values, a broader disease category was prioritized using the disease hierarchy. As a second filtering level, the Wilcoxon test was used to test whether the drug target gene is expressed at higher level in cancer compared with normal erythroid, myeloid, B-lymphoid, or T-lymphoid samples. One normal sample group comparison was accepted for downstream analysis (with the respective comparison annotated as failed). In silico drug analysis was benchmarked using two case studies: drugs from ref. 33 and secondly known vulnerabilities (in clinical use/trial). Success rate was reported for drug target gene expression in disease, specificity for disease/subtypes, and higher expression relative to normal cells.
BeatAML drug analysis
Spearman correlation was computed for each drug AUC values and cancer map clusters, drug target genes, or target gene mutations. Furthermore, mutations with at least 5 observations and significant correlation Padj < 0.05 to drug AUC values or significant Fisher exact test Padj < 0.05 in cancer map clusters were added as features that could improve drug sensitivity prediction.
From total of 122 drugs, 47 were excluded on the basis of three criteria. First, 25 drugs with IC50 lower quartile below 10 nmol/L were excluded as these drugs have limited efficacy. Second, 9 drugs with less than 80 samples with measured drug responses were excluded. Third, only drugs with drug target information were kept, resulting in total of 75 drugs. The elastic net implemented in glmnet (34) was trained using 10-fold cross-validation using caret (35) trainControl and repeatedcv method. Caret function train and its functionality tuneGrid were used to optimize alpha parameter denoting the L1 and L2 regularization term proportions for elastic net. Each drug had three categories of features to fit the model: clusters, drug target gene expression, or mutations. To test the importance of each category in model fitting, sample order was randomly shuffled for one category while the original order was preserved for the other categories. Therefore, if the shuffled category features were important for the model fit, model overall fit should decrease as the other features are unchanged. This process was repeated 100 times and median of RMSE and R2 values were computed. Only drugs with good fit when using all the features were kept, having R2 over 0.25 and RMSE less than 0.9.
Drug sensitivity testing using patient and healthy donor samples
Bone marrow aspirates or peripheral blood samples were obtained from patients with AML (n = 52) and healthy donors (n = 15) after informed written consent using protocols approved by a local Institutional Review Board and in accordance with the Declaration of Helsinki. Mononuclear cells (MNC) were isolated by density gradient separation (Ficoll-Paque PREMIUM; GE Healthcare) and immediately used for drug testing. Cells were maintained in Mononuclear Cell Medium (MCM; Promocell) or in a 25% HS-5 conditioned medium plus 75% RPMI1640 medium mix (CM). Palbociclib and idarubicin (from Selleck) were solvated in dimethyl sulfoxide and plated in 5 different concentrations in 10-fold dilutions on 384-well plates using an Echo acoustic dispenser (Labcyte), 1–10,000 nmol/L for palbociclib; 0.1–1,000 nmol/L for idarubicin. A total of 1 × 104 cells were added per well and incubated with the drugs for 3 days at 37°C, 5% CO2. Viability was measured using the CellTiter-Glo reagent (Promega) according to the manufacturer's instructions and using the PHERAstar (BMG LABTECH) or SpectraMax Paradigm (Molecular Devices) plate readers. Sensitivity to the drugs was quantified using a drug sensitivity score (DSS), which is a modified AUC-based metric described previously (36). A selective DSS value was calculated by subtracting the mean DSS of the healthy bone marrow controls from the DSS of individual AML samples.
Anti-DPEP1 antibody (Atlas antibodies, rabbit polyclonal IgG against human renal dipeptidase 1, product number: HPA009426, lot number: A57960) was used with the dilution 1:2,500 to stain formalin-fixed and paraffin-embedded bone marrow trephine biopsy samples from pediatric patients with pre-B-ALL from the Pirkanmaa ERVA area between the years 2000 and 2017. Diagnostic samples (n = 126; including also one Burkitt lymphoma and a T-lymphoblastic leukemia/lymphoma case) were stained with a Ventana Benchmark GX using UltraView Universal DAB Detection Kit. Cytoplasmic and membranous staining was graded negative if less than 20% of the leukemic blasts were stained, positive if 20% or less than 50% of the blasts were positively stained and strong positive if 50% or a greater proportion of the blasts were positive. The analysis was performed by two pathologists without the knowledge of the patient data or the interpretation of the other analyst. The samples and clinical data were studied with the approval of the Tampere University Hospital Ethical committee (#R16054 and #R13109) and in accordance with the Declaration of Helsinki.
Interactive web resource for data analysis
The interactive online resource and the accompanying user guide for the Hemap resource are described in more detail in the Supplementary Methods and available at http://hemap.uta.fi/.
Integrating transcriptomes to characterize molecular states across hematologic malignancies
For the comparative analysis of hematologic malignancies on molecular level, we assembled gene expression profiles from the NCBI GEO database (9), comprising patient samples representing different cancers and proliferative disorders, and including cell lines and normal blood cell types as controls. Sample annotations were curated, and each sample was assigned a disease category. After data quality control and bias correction (see Materials and Methods; Supplementary Fig. S1), 9,544 samples comprise the final dataset (denoted “Hemap” samples) for subsequent analysis, including 7,279 patient samples (mainly diagnostic) from hematologic malignancies (Fig. 1A; Supplementary Tables S1 and S2).
To enable discovery and statistical comparison of previously known and novel molecular phenotypes alongside the annotated disease classes, we utilized a data-driven approach that allows discovery of sample groups and visualizes these for interpretation. First, we compared dimensionality reduction methods that allow visualization of complex data in 2D space. The data representation accuracy were quantitatively assessed using the metrics of continuity, trustworthiness, and k-nearest neighbor (k-NN) classifier error (see Materials and Methods; Supplementary Fig. S2). As a result, t-SNE (18) and its approximation, Barnes–Hut–SNE (BH-SNE; ref. 19), was selected, as it performed robustly (continuity and trustworthiness, 0.9860 and 0.9943, respectively) in two dimensions and still preserved the neighborhood structure (k-NN error 0.0668; Supplementary Fig. S2). The t-SNE map was then utilized for density-based clustering to assign each sample to a cluster (Fig. 1B, see t-SNE Materials and Methods for details) and the results were compared with annotated disease classes (Fig. 1C). We conclude that both quantitative and biological assessments confirm that our approach faithfully organizes the samples in an unsupervised manner based on their molecular phenotype and disease type. We denote the resulting data-driven sample grouping as the Hemap “cancer map” in the following text.
Comparative analysis associates clinical annotations and pathway activity to the molecular disease stratification
The 2D cancer map revealed a clinically relevant substructure (Supplementary Table S2), as exemplified by the different B-cell lymphomas and pre-B-ALL cytogenetic subtypes (colored in Fig. 2A and B, respectively), providing biological validation for separation of relevant phenotypes on the cancer map. A detailed comparison to annotations is presented in Supplementary Table S2. Next, statistically significant associations of clusters with gene expression levels, clinical annotations and pathway enrichment scores across different databases were calculated (see Materials and Methods). These results can be interactively tabulated and visualized using the online Hemap resource. We selected five most significant pathways at disease cluster level, or those matching pre-B-ALL subtype clusters (Fig. 2C) for visualization in a heatmap (see also Supplementary Table S5). In AML, the pathways for hematopoietic stem cell differentiation, pentose phosphate pathway, renin–angiotensin system, IL8/CXCR1-mediated signaling events and C-MYB transcription factor networks were most significantly enriched. These reflect well the known progenitor-like phenotype of AML cells. Pentose phosphate pathway, on the other hand, represents a recently uncovered vulnerability (37, 38) that is important for AML growth. Similar disease-relevant pathways were also uncovered from T-ALL (TCR pathway), CLL (BCR signaling pathway), lymphomas (cell adhesion molecules; ref. 39), and multiple myeloma (N-glycan biosynthesis; refs. 40, 41). In pre-B-ALL clusters, processes related to transcriptional regulation were highly significant (including histone modification, CTCF pathway, and RNA processing). WNT signaling (42, 43) was found as a cluster 29-specific (t1;19) enriched pathway, which matches its known relevance in these TCF3-PBX1 fusion carrying cases. Samples expressing a gene or pathway of interest can be visualized as shown in Fig. 2D, distinguishing the progenitor-like MLL subtype of pre-B-ALL based on the lack of expression of the differentiation marker MME (also known as CD10) that is used in clinical diagnostics (Fig. 2E). Similarly, most significant associations between disease clusters and drug signatures can be examined by e-staining their significance (in red), as illustrated by association of PI3K inhibitor BEZ235 gene set signature from DsigDB to pre-B-ALL (Fig. 2F), which validates a known association between a drug and a disease subtype. Further analysis on the BEZ235 gene set and several case studies on how to generate novel hypothesis are presented in the accompanying User guide to demonstrate different analysis (refer to “Explore” and “e-staining” examples).
Pan-cancer analysis to recognize vulnerabilities across disease contexts
Parallel to molecular stratification, the diversity of patient profiles in Hemap has the potential to support the development of new therapeutic strategies by leveraging the information about the expression profiles across hematologic malignancies. We analyzed the specificity of drug target expression states across patient groups in a hierarchical manner (Materials and Methods), as illustrated in Fig. 3A (see also Supplementary Table S3 for a list of drugs and their targets and Supplementary Table S4 for significant associations listed by disease hierarchy). The corresponding significance ranking for targets of approved drugs is shown as heatmaps in Fig. 3A and B, where the columns represent different disease contexts and gene targets (in rows) are sorted according to their most significant association. The clinical indication for the drug(s) that could be used to target each gene is indicated in the panel on the right, while e-staining results for example drug targets are shown in Fig. 3C (see also Supplementary Fig. S3). Proteasome-targeting drugs bortezomib and carfilzomib are in use for lymphomas and multiple myeloma. Accordingly, 10 of 20 genes encoding the proteasome subunits are associated to this disease hierarchy level, or to the pan-cancer category, with highest significance (Fig. 3A). In comparison, for precision drugs such as the antibody drugs elotuzumab (SLAMF7, P < 1e-315) or daratumumab (CD38, P = 1e-196) approved for multiple myeloma, or rituximab (MS4A1, P < 1e-315 in LY+CLL) used in lymphomas and CLL (Fig. 3A), the specific gene targets can be examined. Among all known vulnerabilities (drugs in clinical use/trial) a gene-level analysis detected 84% of targets expressed and 69% were associated with highest specificity score (−log10) to the respective disease context (see Supplementary Table S3). This is exemplified by the comparison of genes with significant association to lymphoid leukemias (Fig. 3B). BCL2 targeted by venetoclax is shown as an example of an approved target in CLL that our analysis associates with this disease context and with potential for targeting in multiple myeloma. The genes marked with asterisk, including IL2RA, indicate targets of drugs approved for other hematologic malignancies. Our analysis associated these with repurposing potential in CLL and/or ALL. FLT3 is a recently approved target with disease cluster-specific expression in B-lymphoid and myeloid leukemias.
Utility of molecular disease stratification for evaluating drug screen results
Next, we examined leukemias at disease subtype level from two ex vivo drug screening datasets (12, 33). Venetoclax had lower efficacy in T-ALL versus B-ALL and lowest efficacy was in t1;19 samples in the ALL drug screen (33), which agrees with venetoclax target BCL2 gene expression in Hemap (Fig. 3C). Topotecan and dasatinib had the opposite profile, also in agreement with subtype-specific expression of their targets TOP1MT and LCK (Supplementary Fig. S3). Taken together, out of 15 drugs from this ALL screen tested with our hierarchical analysis, 14 (93%) had a candidate target expressed and 12 (80%) received highest target indication in ALL (Supplementary Table S4). Using the larger beatAML dataset (12), we set out to examine in an unbiased manner what matters more in predicting drug responsiveness: target expression, genetic lesions traditionally used to stratify patients, or the molecular phenotype as defined by clustering of transcriptome states. We implemented models using elastic nets, where a model for each drug (75 in total) was fit using these three categories of features. To test their importance for model fitting, sample order was randomly shuffled for one category while the original order was preserved for the other categories. The results for 11 drugs that achieved the best overall model fit (R2 > 0.25) are shown in Fig. 4A, including venetoclax, panobinostat (HDAC inhibitor), palbociclib (CDK4/6 inhibitor), 7 kinase inhibitors (many targeting FLT3), and an ALK inhibitor. The average R2 value from 100 tests is colored in the heatmap and summarized as a boxplot next to it. If the shuffled feature was important for the model fit, a decrease in R2 is expected (shift from darker red to dimmer or blue colors) as the other features are unchanged. For venetoclax, this analysis implicated target gene expression as the main predictor (Supplementary Fig. S4). For FLT3-targeting compounds, FLT3 mutation status was implicated as the top predictive feature (Supplementary Fig. S4). However, overall, the lack of cluster features in the model resulted in lowest predictive power. The disease clusters were the best predictors for palbociclib and panobinostat, whereas mutation status had no effect on their model fit. Panobinostat and palbociclib showed opposite drug responses in clusters 13, 2, 6 compared with cluster 1 (Fig. 4B). Hemap clusters 17, 5, and 6 corresponded to these clusters (Supplementary Fig. S4) and were similarly enriched for NPM1 and FLT3 mutations or PML-RARA fusion in both datasets. Comparison of clinical phenotypes revealed that blast morphology was different between the clusters, linking maturation level to the differential drug response (Fig. 4C; Supplementary Fig. S4).
Classical targets involved in DNA metabolism (TOP2A and B) and clinically interesting targets, including CDK6, BCL2, MDM2 and VEGFR2 from clinical trials, ranked highly in our disease hierarchy analysis, as shown in Fig. 3 and Supplementary Fig. S3. However, when compared with normal cell types, only 7% of the targets had higher expression in disease than in normal cells (Supplementary Table S3). Palbociclib target CDK6 is highly expressed in all acute leukemias compared with normal blood cell types, while TOP2A has high mRNA levels also in normal blood cells (Fig. 4D). To evaluate drug sensitivity that is specific to cancer cells, an experimental ex vivo screening approach is exemplified in Fig. 4D by comparing in AML patient cell responses to the CDK4/6 inhibitor palbociclib and idarubicin targeting TOP2A (see Materials and Methods). Drug sensitivity and selective drug sensitivity scores (DSS and sDSS, respectively; see Materials and Methods; ref. 36) are compared in box plots (Fig. 4E). Overall, the AML patient bone marrow ex vivo cultures were more responsive to idarubicin (refer to Supplementary Fig. S4 for AML cell line data). However, a negative score indicating higher effect on normal bone marrow cell viability was observed for idarubicin in a larger fraction of AML cases compared with palbociclib. This observation of nonspecific response, implied by negative sDSS score, is consistent with our predictions from Hemap data. Therefore, the normal samples included in Hemap could provide valuable additional information for drug target selection. Comparison of BCL2 and BCL2L1 (also known as BCL-XL) levels are presented as another example in Supplementary Fig. S4, relevant to venetoclax versus navitoclax toxicity in targeting apoptosis. The Advanced Use Case in the Hemap User guide extends this analysis using pathway activities and drug chemical screen data.
Evaluating new therapeutic strategies in a pan-hematologic cancer context
Epigenetic regulation has emerged as an important mechanism that can corrupt the gene regulatory network (44), motivating novel therapeutic approaches. Utilizing the disease spectrum in Hemap, we performed a pan-hematologic cancer analysis of epigenetic modifiers (Supplementary Table S5), focusing on genes encoding proteins that are validated targets of small-molecule drugs (available from ChEMBL; ref. 28). We found elevated expression of this set of genes significantly enriched in CLL, T-ALL, and clusters 28 (pre-B-ALL) and 32 [AML; Fig. 5A, hypergeometric test Padj values 0.0003, 0.0074, 0.0127, 0.0174, respectively; see also Supplementary Table S5 for additional mutation frequency information (45) for the genes shown]. The expression state for six most significant genes from CLL are shown on the Hemap cancer map (Fig. 5A) and from independent validation RNA-seq data (Fig. 5B; ref. 46).
A second promising new strategy, immunotherapy, can kill cancer cells by targeting surface proteins with antibodies (47) or chimeric antigen receptors (48). However, side effects due to targeting normal blood cells along with development of resistance occur (49). To provide a rational basis for extending the target repertoire, we used disease hierarchy analysis to rank 996 candidates available in the Cell Surface Protein Atlas (Supplementary Table S6; ref. 32), resulting in broad, disease, and subtype-specific candidates. The top ranked candidate genes in our analysis correspond to those that are uniformly high expressed within the specified disease context. The stem cell antigen CD33, actively pursued for treatment of AML (50), is among highly ranked surface targets in clinical trials shown in Supplementary Fig. S3. Next, we obtained proteomics profiles from 19 patients with B-ALL (51) to compare our ranked list for pre-B-ALL (refer to Supplementary Table S6) to protein-level expression. The trend between in silico drug target rank and protein detection rate is plotted in Fig. 5C. Validation rate for top candidates was above 75%. The highly ranked surface targets CLEC14A, DPEP1, CELSR2, MME, SDK2, INSR, GPM6B, ELFN2, FLT3, SLC22A16, FLT4, and APCDD1 correspond to those with higher expression in patients with pre-B-ALL compared with normal blood cells (see also Supplementary Fig. S5). The high gene expression state of DPEP1 (Fig. 5D) in pre-B-ALL was further validated at protein level based on IHC of diagnostic bone marrow biopsies. The grading from 117 ALL bone marrows and 16 samples representing other lymphoid malignancies or normal lymphoid tissues is presented in Table 1 and illustrated in Fig. 5E.
|.||DPEP1 IHC .||.|
|.||Negative .||Positive .||Strong positive .||Total .|
|.||DPEP1 IHC .||.|
|.||Negative .||Positive .||Strong positive .||Total .|
Abbreviations: BL, Burkitt lymphoma; CHL, classical Hodgkin lymphoma; MCL, mantle cell lymphoma; NSCHL, nodular sclerosis classical Hodgkin lymphoma; PTCL, peripheral T-cell lymphoma.
To further facilitate the utilization of the data, precalculated results are accessible via our interactive web resource (http://hemap.uta.fi/) including the expression state for 4,277 drug target gene sets and 1,094 drug response signatures, which can be further investigated in the context of the 12,433 pathways and molecular signatures (see Materials and Methods and User guide examples). Disease hierarchy analysis for the curated list of drug to target gene associations (11,373 drugs; 1,182 genes) from the Therapeutic Target Database (TTD; ref. 29), DGIdb (30), and targets of FDA-approved drugs across disease (31) is available in Supplementary Table S4. In this manner, in silico drug target selection based on Hemap can leverage gene and pathway expression, as evaluated across cancer types and normal blood cell types.
The integration of available genome-wide data from patients allows uncovering shared disease mechanisms and new therapeutic options. Recent work has highlighted that molecular and genetic data that helps stratify patients can dramatically increase the likelihood of success during clinical development (8, 52). However, in several cancer types, including those of hematopoietic and lymphoid tissues, the majority of data have been collected by separate studies concentrating on certain cancer types, which hinders the identification of cancer type–specific features. We present an interactive online resource, Hemap (http://hemap.uta.fi/), for analysis across multicenter gene expression datasets to investigate disease subgroups and compare molecular phenotypes across 9,544 samples from hematologic malignancies.
In practice, the samples included to Hemap are inaccessible to most clinical researchers. The Hemap resource serves to repurpose data from public repositories for clinical interpretation in an intuitive manner that does not require data analysis expertise. In future versions of Hemap, we plan to include also RNA-seq studies. Currently, the resource already contains the The Cancer Genome Atlas AML dataset and the User Guide includes examples using this data. Alongside curated disease assignment, we present a data-driven approach that organizes and integrates heterogeneous sample collections in an unbiased manner. To facilitate this, we demonstrated how unsupervised clustering and dimensionality reduction methods, here by the t-SNE method, can be used for organizing the molecular profiles for further downstream analysis. The high level of performance of t-SNE has been shown in context of various data types (18, 53, 54). In this manner, genes characterizing the patient clusters can be identified for further delineation of their functional role. In CLL, our analysis implicates high expression of several polycomb group proteins (SFMBT1, CBX7, and EZH1) in CLL that could be targeted by small molecules, in line with chromatin state data (46), and their mutation (45) frequencies, highlighting the importance to consider the spectrum of genetic and epigenetic changes in these malignancies. Earlier studies have implicated epigenetic plasticity as a key driver of CLL evolution during treatment (55). Specifically, CLL cases had little to no genetic subclonal evolution, while significant recurrent DNA methylation changes were enriched for regions near Polycomb targets (55). To further elucidate the mechanisms, inclusion of posttreatment data and integrating methylation, chromatin marker, and mutation profiles represent important future directions in developing the Hemap resource.
From a therapeutic perspective, approaches for the development of treatment strategies with a broad disease focus and molecular subtype resolution are urgently needed. We used Hemap to provide a roadmap for candidate drug therapies that allows prioritizing new candidates based on disease specificity. Our analysis recapitulated known vulnerabilities, providing additional confirmation for targets in current clinical trials: several compounds targeting Bcl2 have been developed and have shown promise in treating both CLL and non-Hodgkin's lymphoma (56, 57). However, navitoclax that also targets BCL2L1 (also known as BCL-XL) displays platelet toxicity. This potential for off-target effects was visible as high gene expression level in the erythroid lineage, supporting the choice of the more selective venetoclax. The prevalent high expression also in multiple myeloma and pre-B-ALL found in our study provides a rationale for the expansion of the testing of these compounds in lymphoid malignancies. This suggestion is additionally supported by a recent study showing that these compounds have promise in MLL-rearranged ALL (58), a pre-B-ALL subtype corresponding to cluster 29 in our dataset. However, Hemap analysis predicts insensitivity in T-ALL and t1;19 subtype, matching recent ALL drug screen data (33). Similarly, the elevated expression of the p53-regulating MDM2 in pre-B-ALL fits with recent data on successful application of antagonists in clinical trials (59), and mechanisms for its high expression ETV6-RUNX1–positive leukemias (60).
Currently, no drug screens have been carried out in primary patient cells across the spectrum of hematologic malignancies in Hemap. The utility of Hemap for drug repurposing was demonstrated in our recent study that identified dasatinib as a targeted therapy for a subgroup of patients with T-ALL (61). Here, we examined drug screen datasets to examine how differential drug responsiveness could be linked to disease subclusters and drug targets identified from the cancer maps. Using the beatAML dataset, we systematically compared the importance of mutations, clusters, and drug target gene expression in predicting drug responses. For drugs with a good overall model fit, clusters were the best predictors. However, the importance of each predictor was largely influenced by drug type. Best predictor for venetoclax response was BCL2 expression level. FLT3 mutation status and other mutations were the best predictors for kinase inhibitors. In contrast, disease clusters were the best predictors for palbociclib and panobinostat responses to which mutation status had no effect on model fit. Comparison of clusters in which panobinostat and palbociclib showed opposite drug responses revealed that blast morphology was different, linking maturation level to differential drug response. Furthermore, their drug targets were differentially expressed in these clusters, pointing out the importance of integrating context and drug target expression for in silico drug screening. Surprisingly, the HDAC expression pattern revealed cytosolic members (HDAC6 and HDAC10; ref. 62) as resistance markers, while nuclear HDAC4 and HDAC9 correlated with sensitivity. Our analysis also supported CDK4/6 as disease-specific targets that are known to act as critical activators of normal and leukemic hematopoietic stem cells (63). Here, palbociclib compared favorably to idarubicin regarding patient blast sensitivity against normal bone marrow cells, reflected in mRNA data from Hemap. The selectivity over normal cells may improve further using combination therapy (63) that allows decreasing the dose. However, additional parameters such as drug target protein level, drug metabolism, and cell proliferation rate further contribute to sensitivity and therefore not all patients matching a molecular subtype or expressing the target mRNA can be expected to respond favorably.
Cancer cells display remarkable plasticity: resistance to recently approved CD19-targeting CAR-T therapy has been shown to occur via mutations or splicing defects at the CD19 locus or lineage switching (49). To combat the diversity of resistance mechanisms, there is a demand to diversify the target repertoire. In pre-B-ALL, we identified promising surface protein candidates, prioritizing targets with consistently high levels within the given disease context, and low levels in normal blood cell types. Over 75% of the highly ranked candidates were confirmed using proteomics (51), and additional literature confirmation was found for five candidates. Moreover, we validated DPEP1 as a potential surface target in pre-B-ALL by IHC staining of diagnostic bone marrow biopsies. Positive staining was found in each subtype for majority of cases, except in MLL where both the Hemap gene expression data and protein staining indicated low or undetectable levels. The validation cohort consisted of pediatric cases, while Hemap analysis included also adult samples. DPEP1 is a zinc-dependent metalloproteinase that is expressed aberrantly in several cancers, and has been implicated as a potential therapeutic target in colon and pancreatic cancers (64, 65). In future, increased availability of protein-level data from different hematologic malignancies will allow confirming additional targets.
In conclusion, the interactive Hemap resource facilitates comparative analyses across multiple hematologic malignancies. We envision that the mechanistic insight gained by concomitant identification of molecular subtypes, genetic aberrations, and derailed cellular pathways will expedite therapeutic innovations and clinical utility.
Disclosure of Potential Conflicts of Interest
No potential conflicts of interest were disclosed.
Conception and design: P. Pölönen, K. Porkka, M. Nykter, M. Heinäniemi
Development of methodology: P. Pölönen, J. Mehtonen, J. Lin, S. Häyrynen, M. Nykter
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): A. Mäkinen, D. Malani, C.A. Heckman, O. Lohi
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P. Pölönen, J. Mehtonen, J. Lin, T. Liuksiala, S. Häyrynen, A. Mäkinen, A. Kumar, V. Pohjolainen, K. Porkka, K.J. Granberg, O. Lohi, M. Nykter, M. Heinäniemi
Writing, review, and/or revision of the manuscript: P. Pölönen, J. Mehtonen, J. Lin, S. Häyrynen, S. Teppo, A. Mäkinen, K. Porkka, C.A. Heckman, P. May, K.J. Granberg, O. Lohi, M. Nykter, M. Heinäniemi
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Teppo, A. Mäkinen, P. May, K.J. Granberg
Study supervision: V. Hautamäki, K.J. Granberg, O. Lohi, M. Nykter, M. Heinäniemi
We would like to thank the Finnish Pharmaceutical Information Centre Ltd for information provided regarding approved drugs in hematologic malignancies, Dr. Aik Choon Tan for the advice with drug signature gene sets, staff of the FIMM Technology Center High-Throughput Biomedicine Unit and Sequencing Lab, members of the T. Aittokallio, C.A. Heckman, O. Kallioniemi, S. Mustjoki, K. Porkka, and K. Wennerberg groups for their help to analyze drug response profiles, and Prof. I. Shmulevich, Dr. Sheila Reynolds, and Roger Kramer for assistance with feature matrix analysis. The work was supported by grants from the Academy of Finland (project no. 312043 to M. Nykter; 310829 to M. Nykter; 259038 to K.J. Granberg; 276634 to M. Heinäniemi 277816 to O. Lohi), The Finnish Cultural Foundation (Interdisciplinary Science Workshops, to M. Heinäniemi), Sigrid Juselius Foundation (to M. Nykter), and Cancer Society of Finland (to M. Nykter, M. Heinäniemi, O. Lohi, K.J. Granberg), Paulo Foundation (to O. Lohi), Foundation for Pediatric Research (to O. Lohi), Jane and Aatos Erkko Foundation (to O. Lohi), Nokia Foundation (to V. Hautamäki), and University of Eastern Finland (to M. Heinäniemi). We thank CSC–IT Center for Science and UEF bioinformatics center for providing computational resources.
The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.