Comprehensive and unbiased multiparameter high-throughput screening by compaRe finds effective and subtle drug responses in AML models

  1. Morteza Chalabi Hajkarim
  2. Ella Karjalainen
  3. Mikhail Osipovitch
  4. Konstantinos Dimopoulos
  5. Sandra L Gordon
  6. Francesca Ambri
  7. Kasper Dindler Rasmussen
  8. Kirsten Grønbæk
  9. Kristian Helin
  10. Krister Wennerberg  Is a corresponding author
  11. Kyoung-Jae Won  Is a corresponding author
  1. Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Denmark
  2. Institute for Molecular Medicine Finland (FIMM), Helsinki Institute of Life Science, University of Helsinki, Finland
  3. Rigshospitalet, Denmark
  4. Centre for Gene Regulation and Expression, School of Life Sciences, University of Dundee, United Kingdom
  5. Cell Biology Program and Center for Epigenetics Research, Memorial Sloan Kettering Cancer Center (MSKCC), United States

Abstract

Large-scale multiparameter screening has become increasingly feasible and straightforward to perform thanks to developments in technologies such as high-content microscopy and high-throughput flow cytometry. The automated toolkits for analyzing similarities and differences between large numbers of tested conditions have not kept pace with these technological developments. Thus, effective analysis of multiparameter screening datasets becomes a bottleneck and a limiting factor in unbiased interpretation of results. Here we introduce compaRe, a toolkit for large-scale multiparameter data analysis, which integrates quality control, data bias correction, and data visualization methods with a mass-aware gridding algorithm-based similarity analysis providing a much faster and more robust analyses than existing methods. Using mass and flow cytometry data from acute myeloid leukemia and myelodysplastic syndrome patients, we show that compaRe can reveal interpatient heterogeneity and recognizable phenotypic profiles. By applying compaRe to high-throughput flow cytometry drug response data in AML models, we robustly identified multiple types of both deep and subtle phenotypic response patterns, highlighting how this analysis could be used for therapeutic discoveries. In conclusion, compaRe is a toolkit that uniquely allows for automated, rapid, and precise comparisons of large-scale multiparameter datasets, including high-throughput screens.

Editor's evaluation

This paper aims to address the current gap in the efficient analysis of large-scale multiparameter flow cytometry and other datasets. The authors offer a software toolkit with an efficient algorithm for comparing numerous samples at once. The study is well presented and is relevant to single cell analysis research.

https://doi.org/10.7554/eLife.73760.sa0

eLife digest

Biology has seen huge advances in technology in recent years. This has led to state-of-the-art techniques which can test hundreds of conditions simultaneously, such as how cancer cells respond to different drugs. In addition to this, each of the tens of thousands of cells studied can be screened for multiple variables, such as certain proteins or genes. This generates massive datasets with large numbers of parameters, which researchers can use to find similarities and differences between the tested conditions.

Analyzing these ‘high-throughput’ experiments, however, is no easy task, as the data is often contaminated with meaningless information, or ‘background noise’, as well as sources of bias, such as non-biological variations between experiments. As a result, most analysis methods can only probe one parameter at a time, or are unautomated and require manual interpretation of the data.

Here, Chalabi Hajkarim et al. have developed a new toolkit that can analyze multiparameter datasets faster and more robustly than current methods. The kit, which was named ‘compaRe’, combines a range of computational tools that automatically ‘clean’ the data of background noise or bias: the different conditions are then compared and any similarities are visually displayed using a graphical interface that is easy to explore.

Chalabi Hajkarim et al. used their new method to study data from patients with acute myeloid leukemia (AML) and myelodysplastic syndrome, two forms of cancer that disrupt the production of functional immune cells. The toolkit was able to identify subtle differences between the patients and categorize them into groups based on the proteins present on immune cells.

Chalabi Hajkarim et al. also applied compaRe to high-throughput data on cells from patients and mouse models with AML that had been treated with large numbers of specific drugs. This revealed that different cell types in the samples responded to the treatments in distinct ways.

These findings suggest that the toolkit created by Chalabi Hajkarim et al. can automatically, rapidly and precisely compare large multiparameter datasets collected using high-throughput screens. In the future, compaRe could be used to identify drugs that illicit a specific response, or to predict how newly developed treatments impact different cell types in the body.

Introduction

Technological developments have accelerated the generation of large-scale multiparameter screening data through methodologies such as high-content microscopy and high-throughput flow cytometry (Boutros et al., 2015; Saeys et al., 2016; Caraus et al., 2015). These technologies can test hundreds of samples (such as drug treatments) each with tens of thousands of events (e.g. cells) labeled for numerous biomarkers (such as cytoplasmic or membrane markers). However, analyzing this massive multiparameter data to provide an overview of similarities and differences between hundreds of samples is still a challenge (Boutros et al., 2015; Saeys et al., 2016; Caraus et al., 2015). This analytical challenge is further complicated by various sources of bias and noise often existing in the data, such as batch effect and signal drift (a gradual shift in the marker intensity across a multi-well plate) (Boutros et al., 2015; Saeys et al., 2016; Caraus et al., 2015).

There have been efforts to cluster samples from large-scale multiparameter (multidimensional) screening data. A simple approach is to use a representative value for each cell marker such as median fluorescence intensity (MFI) for clustering samples (Cossarizza et al., 2019). However, using a single representative value can easily lead to loss of information about biologically relevant variance within and between cell subpopulations. Meta-clustering with single-cell clustering algorithms has been suggested to cluster samples based on the similarity of the centroids of cell subpopulations identified in the individual samples (Qiu et al., 2011; Levine et al., 2015; Van Gassen et al., 2015; Ogishi et al., 2021). While these algorithms are widely used in single-cell data analysis for clustering cells, they are not efficient for clustering of samples. This is because centroid-based analysis can be misleading when subclusters are not sufficiently distinct or the number of sub-clusters varies. Additionally, the heavy computing cost of meta-clustering makes it poorly suited for analyses of large datasets with many samples. Manual gating and machine learning based on prior knowledge have been used to cluster samples (Amir et al., 2019; Bruggner et al., 2014), but using prior knowledge for subpopulation identification can both lead to biased interpretations and failure to make de novo discoveries. Dimension reduction methods (Lvd and Hinton, 2008; Amir et al., 2013; McInnes et al., 2018) coupled with the Jensen-Shannon divergence (JSD) metric have also been used to cluster multidimensional samples (Amir et al., 2013). These algorithms including factor analysis and principal component analysis (PCA) still require excessive computing costs with an inherent information loss. It is also important to note that none of the methodologies developed so far efficiently correct for sources of bias and noise in large-scale multiparameter screening data.

Available computational toolkits (BioScience E, 2020; Potdar et al., 2020; Boutros et al., 2006) mostly allow for single-parameter or unautomated analyses of large-scale screening data using the aforementioned methods. In these toolkits, each well should be first represented by a single parameter such as cell counts or centroids or they require manual intervention. To provide a useful toolkit for precise and effective interpretation of small- to large-scale multiparameter screening data, we developed compaRe. This toolkit has several unique modules for quality control, bias correction, pairwise comparisons, clustering, and data visualization. The quality control and bias correction modules can effectively reveal and remove various sources of bias in the screening data. compaRe clusters samples by measuring the similarity between them using a dynamic mass-aware gridding algorithm. This algorithm increases the robustness of the toolkit to the size of data and signal shift (a technical term referring to batch effect and signal drift), while guaranteeing fast clustering, as it does not bear the computing cost of dimension reduction and subsampling. The toolkit is available both as a command-line version and a graphical user interface (GUI) version that provides various visualizations to help with the interpretation of its readouts.

compaRe performed robustly in the presence of background noise and batch effects even where these input data artifacts could not be corrected. compaRe analyses of multiparameter mass and flow cytometric data from acute myeloid leukemia (AML) and myelodysplastic syndrome (MDS) patient samples revealed interpatient heterogeneity and recognizable phenotypic profiles. When applied to high-throughput flow cytometry of the dose response of AML samples treated with various drugs, compaRe successfully corrected for various sources of bias and clustered the samples based on their response to treatment, allowing for detection of both drastic and subtle phenotypic responses.

Results

compaRe is a comprehensive toolkit for multiparameter screening data

compaRe is designed to analyze the data from small to large-scale multiparameter screening assays such as high-throughput flow cytometry, high-content microscopy, mass cytometry, and standard flow cytometry. The toolkit comprises several modules for quality control, bias correction, clustering, and visualization. Figure 1 shows the modules for a high-throughput flow cytometry of AML samples taken from a mouse model treated with various drugs. During quality control, several sources of bias such as autofluorescence, bioluminescence, carryover effect, edge effect, signal drift, and cell viability drift (drift in the number of live cells across the plate) were identified. The bias correction module could effectively correct for signal and cell viability drifts (two main sources of bias in high-throughput screening with fluorescent markers) using regression analysis (Figure 1, Materials and methods).

compaRe is a comprehensive suite for multiparameter screening data.

High-throughput flow cytometry generates massive multidimensional data from hundreds of samples. compaRe’s quality control (QC) module reveals several sources of bias in the assay such as signal (intensity difference between the top left and bottom right corners) and cell viability drifts. These two are corrected for in the bias correction modules within and between the plates. compaRe performs a pairwise similarity calculation between the samples using dynamic gridding and forming hypercubes (represented by distinct colors). The portions of the data within individual hypercubes are used to calculate similarity. Clustering is performed based on similarity. The graphical user interface (GUI) provides several ways to thoroughly explore and visualize the read-outs.

At the core of the compaRe toolkit is a module for pairwise comparisons of samples. It measures the similarity between two samples using a dynamic mass-aware gridding algorithm (Figure 1, Materials and methods). Given two samples, the algorithm divides the higher dimensional space (formed by, for example, cell surface markers) of the samples individually into several spatial units called hypercubes. The average difference between proportions of data points present in corresponding hypercubes across the samples is used to represent similarity. In this setting, the module becomes robust to signal shift and data size difference between the two samples (Appendix 1). This module generates a similarity (affinity) matrix for the clustering module.

The clustering module uses a graphical algorithm (Figure 1, Materials and methods). Initially, all nodes (samples) are connected forming a complete weighted graph wherein weights represent similarity values. The graph is then pruned to remove potential false positive edges using a threshold inferred from negative controls (untreated samples). After constructing a linked graph, clustering is tantamount to finding maximal cliques (complete subgraphs that cannot be extended), each containing samples with similar responses. compaRe benefits from parallel computing and modular design. Its modular design allows the modules to run independently; thus, the similarity and clustering modules of compaRe can be potentially applied to any problem space.

compaRe is ultra-fast and robust to background noise and batch effect

To evaluate the robustness of compaRe’s comparison module to noise and batch effect, we benchmarked it against JSD with UMAP (for simplicity just JSD) and meta-clustering with PhenoGraph (for simplicity just meta-clustering) (Levine et al., 2015). We analyzed the publicly available mass cytometry data of a total of 21 bone marrow aspirate samples collected from 16 pediatric AML patients and five healthy adult donors labeled for detection of 16 cell surface markers (Levine et al., 2015). We introduced random noise with Gaussian distribution to the 16 parameters of each sample to simulate a batch effect. In this setting, although the added noise undermines similarity, the overall cell population configuration remains intact, and consequently the simulated samples will still have the highest similarity with their original samples.

Even with the added noise, the comparison module correctly identified similar samples (Figure 2a). Conversely, the batch effect seriously compromised the performance of both meta-clustering and JSD, showing several maximum similarities other than the originals (Figure 2b and C). In additional comparison with FlowSOM and SPADE, other commonly used tools for flow cytometry, compaRe’s performance far exceeded their performance (Appendix 1—figure 1). This result demonstrates the advantage of using dynamic gridding for comparison of samples in the presence of noise or batch effect.

compaRe robustly measures the similarity between samples in the presence of batch effect.

Similarity matrix generated by compaRe is shown in (a). Size and color of dots represent the level of similarity. Self-comparisons were removed. Noise was added (marked with *) to the original 21 mass cytometry samples of bone marrow aspirates from 16 pediatric AML patients (S) and five healthy adult donors (H). Similarity matrices using JSD with UMAP and meta-clustering with PhenoGraph are shown in (b) and (c), respectively. The run time of comparing each sample to itself is shown in (d). Samples were sorted based on their size.

Notably, compaRe took only 25 min to analyze the 21 samples (210 pairwise comparisons), without subsampling or dimension reduction. Meanwhile, meta-clustering and JSD took 39 hr and 10 hr respectively. For the feasibility of JSD, we subsampled each sample to 100,000 cells (default value suggested in Amir et al., 2013). When we fixed this limit to 60% of each sample, the computing time of JSD increased to 3 days. To investigate the relation between run time and sample size, we compared each sample to itself and sorted measured times based on sample size (Figure 2d). The run time increased steeply for both meta-clustering and JSD as the sample size increased, while the increase for compaRe was almost unnoticeable.

To further show that compaRe can identify phenotypic changes from a high-dimensional dataset, we used a subset of the data with three healthy and two AML samples stained with 29 (15 membrane and 14 intracellular signaling) markers (Appendix 1—figure 2). Taking H1 as reference, we gradually removed 25%, 50%, 75%, and 100% of cells from a target cluster identified by PhenoGraph. The gradual removal can be regarded as a phenotypic change and the 75% reduction can potentially resemble a rare cell population (a small cluster of cells). As shown in the UMAP projections, the similarity decreased concurrently and more drastically after 100% reduction when phenotypic changes were detected, indicating compaRe is sensitive to phenotypic changes and the existence of rare cell populations.

compaRe reveals interpatient similarity

Non-AML myeloid neoplasias such as MDS can evolve to become AML. Over time, about one-third of all MDS cases develop into AML (DeVita and Lawrence, 2015; Niederhuber et al., 2020). The risk of developing AML largely depends on the MDS subtype at the time of diagnosis, with high-risk MDS developing into AML more often than the lower-risk MDS subtypes (Greenberg et al., 2012). As many immunophenotypic abnormalities are not unique to MDS, several diagnostic flow cytometric antibody panels have been proposed (van Dongen et al., 2012; Alhan et al., 2016). The EuroFlow AML/MDS antibody panel (van Dongen et al., 2012) aims at the parallel identification and categorization of AML and MDS. Both diseases are heterogeneous, affecting multiple cell lineages and multiple maturation stages. Therefore, this panel concerns major myeloid lineages (neutrophilic, monocytic and erythroid) and the detection of abnormal lymphoid maturation profiles in four tubes. The panel uses four backbone markers to identify myeloblasts and an additional set of 15 markers devoted to the characterization of myeloid lineages (Supplementary files 1 and 2).

Unlike the backbone markers, the characterization markers are divided into each tube exclusively. This design was made so that characterization markers from different tubes can be inferred on the same backbone marker subpopulations, but the design makes it impossible to form a multiparameter dataset which is required for clustering methods. However, as compaRe’s comparison module can compare cell population morphologies even in subspaces, we were able to use it to measure similarities between patient samples.

We analyzed 25 bone marrow mononuclear cell samples collected from 16 MDS patients and 9 AML patients (Supplementary file 3). The comparison module provided a detailed overview of similarities of samples. As expected, the AML samples exhibited a great amount of interpatient heterogeneity compared to the MDS samples (Figure 3a and b) with all MDS samples clustered together, and the AML samples spread over three clusters. To verify the performance of the module, we visualized the pairwise comparisons using UMAP projection (Figure 3c and Appendix 1—figures 326). The measured similarities perfectly matched the projections so that from top left to bottom right, as the similarity decreases, the degree of overlap decreases, and the number of exclusive cell populations increases.

compaRe highlights immunophenotypic similarities.

(a) The similarity band plot visualizes the similarity between a sample specified by its row (band) and other samples measured by compaRe (H: higher-risk MDS, L: lower-risk MDS and A: AML). Each band was independently transformed by an exponential function to emphasize the highest and the lowest similarity values. (b) A graphical representation of the similarities. The graph nodes (samples) were clustered by a random walk. (c) The UMAP projection of A1 sample against the other patient samples is provided as an example. The other projections are given in Appendix 1—figures 326. The projections were sorted based on similarity.

We further investigated how different the three groups of the AML samples were (Figure 4 and Appendix 1—figure 27). AML samples 1 and 9 of the blue cluster were confirmed to have a high degree of monocytic differentiation with marked expression of the monocytic maturation markers CD14, CD35, CD64, and CD300e. The AML samples of the green cluster, on the other hand, represented a cluster of poorly differentiated AML cases with low expression of differentiation markers and high expression of the stem cell/progenitor markers CD34 and CD117. Unlike the blue cluster with high monocytic differentiation, and the green cluster with poor monocytic differentiation, the AML samples 2 and 5 of the red cluster included both positive and negative populations of CD11b which is a common granulocytic and monocytic maturity marker, a feature observed in all MDS samples as well (Appendix 1—figure 27).

Immunophenotypic profiles of two different groups of AML patients.

Each row shows the UMAP projection of AML samples 1 and 9 (red and orange) vs AML samples 3, 4, 6-8 (blue) of the green cluster of Figure 3b stained by the markers available in each tube.

In conclusion, compaRe’s comparison module can be used to optimize true cytometric n-dimensional immunophenotypic characterization of patient samples. Interpretation can then be performed in a conventional manner assisted by lower-dimensional projection tools such as PCA and UMAP that promptly provide a phenotypic profile of the patient samples.

Identifying cell-subtype-specific drug responses in mouse AML cells

We applied compaRe to high-throughput flow cytometry data to identify cell subtype-specific responses evoked by antineoplastic agents in leukemic spleen cells from an AML mouse model. Splenic cells were sorted for c-Kit cell surface expression, allowing for the enrichment of stem/progenitor-type leukemic cells. On ex vivo expansion, these cells continuously expand and differentiate in a similar way as in vivo with a clear stem cell/progenitor population and partial differentiation towards CD11b/Gr-1 or CD16/CD32-expressing myeloid cells. After ex vivo expansion, the leukemic cells were plated onto multi-well plates containing a library of 116 antineoplastic agents including surface and nuclear receptor inhibitors and activators, enzyme inhibitors and, cytotoxic chemotherapy in a five-point concentration range, as well as 20 negative control wells (Supplementary file 4). After 72 hr of drug exposure, we stained the cells with fluorescently labeled antibodies against three cell surface markers (CD16/32, Gr-1 and CD11b) and quantified cell surface marker expression using a high-throughput flow cytometer.

compaRe corrected the intraplate signal drift, sources of bias in cell numbers, as well as inter-plate sources of bias (Appendix 1—figure 28). After clustering and clique analysis, we obtained 134 cliques, each sharing similar drug responses (Supplementary file 5).

To get an overview of the assay, we generated a dispersion map of the clusters (Figure 5a and b and Materials and methods). We identified a distinct response group characterized with decreased Gr-1 and concomitant increase of CD16/CD32 as compared to control (Group one in Figure 5a). Most of the cliques included in this response group consisted of drugs in high concentrations with cytotoxic/cytostatic effects. However, some drugs in this group had a milder effect on live cell numbers, and these were enriched for mitogen-activated protein kinase (MAPK) pathway-associated inhibitors (Figure 5c, Supplementary file 6). For instance, trametinib (2.5 nM) in clique 23 (C23) showed a marked decrease of Gr-1 and increase of CD16/CD32, further confirming the results of compaRe (Figure 5d). The MAPK pathway is a regulator of diverse cellular processes such as proliferation, survival, differentiation, and motility (Dhillon et al., 2007). Our findings suggest that MAPK signaling controls the differentiation and/or proliferation towards Gr-1-/CD16+ cells.

compaRe analysis identifies several distinct cell subtype-specific responses in a high-throughput flow cytometry screening of mouse AML cells.

(a) A UMAP plot of cliques identified by compaRe. Cliques are colored by Gr-1 and CD16/CD32 MFIs. Group one is characterized with reduced Gr-1 and increased CD16/CD32 as compared with control. Group two has increased Gr-1 expression compared with control. (b) Heatmap of marker MFIs. Values are normalized between 0 and 1 per marker to make cross-comparisons possible. Cliques containing control, trametinib (2.5 nM) (C23), molibresib and birabresib (C100 and C110), and vincristine (C80) are marked. (c) Waterfall plot of compounds belonging to response group 1, showing live cell count as a percentage of control treatment (DMSO). (d) Density scatter plots for Control (DMSO), C23, C100, and C80.

In high concentration, molibresib and birabresib, inhibitors of BET proteins BRD2, BRD3, and BRD4, caused a reduction in live cell counts but also a reduction of MFI in all the measured markers, which corresponds to the loss of differentiation marker positive cells (Gr-1+, CD11b+, CD16/CD32 high) (Figure 5b: C100, C110, Figure 5d). The BRD2/3/4 proteins regulate transcription via recognition of acetylated lysines on histones and concomitant recruitment of other transcription and chromatin remodeling factors to enhance transcriptional activity (Ferri et al., 2016). The enrichment of undifferentiated cells could therefore be due to an early block in differentiation or that inhibition of BRD2/3/4 has led to a general decrease of cell surface protein transcription.

In this cell model, the leukemic stem-like cells are expected to be present within the differentiation marker negative population. These cells are potential targets for treatments against leukemia. We observed response group 2 (Figure 5a) had a higher MFI in marker Gr-1 as compared to control, the increase was very slight and seemed to be linked to toxic drug concentrations. However, three drugs, vincristine (C80), tazemetostat, and tretinoin clearly reduced the proportion of differentiation marker negative cells (Figure 5d). Interestingly, these three drugs have distinct modes of action: vincristine is a microtubule polymerization inhibitor, tazemetostat inhibits the histone methyltransferase EZH2, and tretinoin is a retinoic acid receptor agonist (Supplementary file 6).

Taken together, compaRe analysis of the high-throughput flow cytometry screening data allowed rapid identification of several distinct phenotypic responses in this mouse AML model, as well as the cellular signals that drive them. Drugs of different mechanism of action can still cluster together if the cellular processes they affect converge in a specific model. Drug response in association with genetic alterations can be one of the applications of compaRe. The genetic alteration could be visualized in the clusters that compaRe identifies.

Identifying highly selective signal transduction inhibitors in human AML cells

We further applied compaRe to the drug screening data from an AML patient sample. Primary AML bone marrow mononuclear cells were dispensed into a 384-multiwell plate containing a library of 40 drugs and drug combinations in seven-point concentration ranges (Supplementary file 7). After 72 hr of drug exposure, the cells were stained with fluorescently labeled antibodies against a panel of AML-related cell surface markers (CD45, CD34, CD38, CD117, HLA-DR, CD45-RA, CD3 and a mix of myeloid differentiation-related markers). A high-throughput flow cytometer was used to quantify cell surface marker expression.

compaRe analysis identified several distinct response groups (Figure 6a, Supplementary file 8). Response group one had notably higher MFIs in the CD34 and CD38 channels compared to controls. Interestingly, the increase in MFIs was due to a drug concentration-dependent appearance of a CD34+/CD38+ cell population that was barely detectable in the DMSO control samples (Figure 6b). The appearance of this CD34+/CD38+ population was also concomitant to a general increase in live cell count (Figure 6c). Altogether, seven different drugs had the same effect (Figure 6d), most of them being selective signal transduction inhibitors such as trametinib (MEK inhibitor), copanlisib (PI3K inhibitor), and PIM447 (PIM kinase inhibitor).

Identification of drugs that induce expansion of CD34+/CD38+ cells in an AML patient sample.

(a) UMAP of cliques identified by compaRe. Cliques are colored by CD34 and CD38 MFIs. Response groups of interest are indicated using a dashed line. (b) Example of response group 1: density scatter plot of markers CD34 and CD38 in different concentrations of PIM kinase inhibitor PIM447. (c) Count of live cells after 72 hr exposure to different concentrations of PIM kinase inhibitor PIM447. (d) Table of drugs that induced expansion of the CD34+/CD38+ cell population.

Response group 2 consisted of two drugs: birabresib and lenalidomide in different concentrations. These induced a decrease in the MFI of CD45-RA and CD45 channels (Appendix 1—figure 29a). In the case of lenalidomide, this response was likely due to cell toxicity and/or growth inhibition (Appendix 1—figure 29b). Interestingly, the birabresib response was very pronounced without the loss of live cell numbers, (Appendix 1—figure 29b) but with a decrease in the MFI in the cell differentiation marker mix channel (Appendix 1—figure 29c).

compaRe also detected response group three as distinct from the controls. This group includes treatment with tretinoin (several concentrations), navitoclax, and mitoxantrone (low dose). Further validation showed the phenotypic response in group three is subtle but with a distinct increase in CD34+ cells (Appendix 1—figure 29d). This result highlights compaRe analysis is sensitive enough to identify these subtle changes.

Discussion

Technological advancements in multiparameter high-throughput screening have enabled testing thousands of biological conditions in a short amount of time. This requires algorithmic development to analyze the large amount of data generated by such technologies. We developed an automated comprehensive toolkit, compaRe, for robust analysis of small- to large-scale multidimensional screening data with several modules for quality control, bias correction, comparison, clustering, and visualization.

The toolkit is unique in many ways. Its quality control and bias correction modules can correct for signal and cell viability drifts in large-scale fluorescence-based screening assays using regression analysis. Its comparison module utilizes a dynamic mass-aware gridding algorithm, which substantially reduces the computing cost and provides robustness to signal shift (batch effect and signal drift). Alternative approaches such as meta-clustering and JSD require both sub-sampling of the data, with the possible loss of valuable subpopulations, and considerably more computing time.

We tested the robustness of the comparison module to batch effect and noise through simulation. The module effectively circumvented the batch effect while JSD and meta-clustering significantly suffered from it. The poor accuracy of meta-clustering demonstrates the drawback of using cluster centroids for similarity comparison across samples while the poor performance of JSD indicates that this approach can work well only in the absence of signal shift. It is of particular note that compaRe does not need subsampling or dimension reduction of the input data.

Multiparameter cytometric analysis of immunophenotypes of AML and MDS patient samples by the comparison module coupled with the EuroFlow AML/MDS antibody panel revealed interpatient heterogeneity and recognizable phenotypic profiles. Even though EuroFlow markers are divided into several discrete tubes, compaRe’s comparison module can compare the cell population distribution to measure similarities between patient samples.

We investigated several types of responses evoked by different doses of antineoplastic agents in two high-throughput flow cytometry screening assays of an AML mouse model and an AML human patient. We could identify subtle but distinct phenotypic drug-induced changes. We also identified drugs with different mechanism of action but similar responses. In general, we showed that drugs will cluster together if the cellular processes they affect converge in a specific model.

The quality control and bias correction modules could successfully correct for signal and cell viability drifts in these studies. In our explored assays, signal drift was obviously associated with the order in which wells were read. It was caused by the time differences in antibody incubation across the plate as the high-throughput flow cytometer requires more than one hour to sample all wells in a 384-well plate. For high-density assay plate formats with large numbers of wells, this can cause gradual incremental influences in intensity and cell viability. Therefore, when aligning wells along the order that the flow cytometer sampled the wells, we found a linear trend in MFIs. We benefited from regression analysis to remove the effect of signal shifts.

During the analyses, the compaRe toolkit made it easy to explore and compare highly complex datasets in a substantially reduced timeline. It is equipped with multithreading and can run through command-line interface on a computer server or GUI on a desktop. The GUI provides the investigator with numerous interactive visualization tools including cell staining, graphical representation, and gating. In sum, it provides a total package for fast, accurate, and readily interpretable multiparameter screening data analysis.

Materials and methods

Mass cytometry of healthy and pediatric AML bone marrow aspirates

Request a detailed protocol

Mass cytometry dataset for 21 samples labeled with 16 surface markers collected from 16 pediatric AML patients obtained at diagnosis and five healthy adult donors (Levine et al., 2015) were downloaded from Cytobank Community with the experiment ID 44185. There are 378 FCS files in this experiment with one FCS file for each of 21 patients for each of 17 conditions (two basal replicates and 16 perturbations). All FCS files from a single patient had been pooled then clustered with the PhenoGraph algorithm. Each file includes a column named PhenoGraph that specifies the PhenoGraph cluster to which each event was assigned as an integer. A value of 0 indicates no cluster was assigned because the cells were identified as outliers during some stage of analysis. Using the PhenoGraph column, we determined centroids of cell clusters, and used PhenoGraph to meta-cluster them as described in Levine et al., 2015 To generate the similarity matrix, we adapted an approach similar to that of compaRe such that each meta-cluster as a spatial unit was treated like a hypercube. We set compaRe’s n to four for this assay (Materials and methods and Appendix 1).

High-throughput flow cytometry of AML mouse model

Request a detailed protocol

AML primary splenic cells from Npm1+/cA (Vassiliou et al., 2011); Flt3+/ITD (Lee et al., 2007); Dnmt3a+/- (Kaneda et al., 2004) Mx1-Cre+ (Kühn et al., 1995) moribund mice were sorted for c-Kit positivity and expanded ex vivo. AML cells were treated with a library of 116 chemotherapy and immunotherapy antineoplastic agents in a five-point concentration range (Supplementary file 4). Treated samples were stained with three informative cell surface antibodies (Supplementary file 9) and fluorescence was detected using a high-throughput flow cytometer iQue Screener Plus (Intellicyt). We set compaRe’s n to five for this assay.

High-throughput flow cytometry of an AML human patient sample

Request a detailed protocol

Mononuclear cells were isolated from a donated human bone marrow aspirate from an AML patient (Danish National Ethical committee/National Videnskabsetisk Komité permit 1705391). The cells were treated with a library of 40 chemotherapy and targeted antineoplastic agents in a seven-point concentration range (Supplementary file 7) for 72 hr. Cells were subsequently incubated with fluorescently labeled antibodies targeting 11 informative cell surface proteins in eight fluorescence channels (Supplementary file 10). Samples were read using a high-throughput flow cytometer (iQue Screener Plus, Intellicyt). We set compaRe’s n to three for this assay.

Flow cytometry of AML and MDS patients

Request a detailed protocol

Clinical flow cytometry data using a slightly modified AML panel as described by the Euroflow Consortium (van Dongen et al., 2012) from 25 bone marrow aspirates from MDS and AML patients from Rigshospitalet (Copenhagen, DK) were used for analysis. Each sample was analyzed using a total of four tubes (Euroflow AML panel tubes 1–4) with eight antibodies in each tube (Supplementary files 1 and 2). Acquisition of data was performed on a FACS Canto (Becton Dickinson Immunocytometry Systems), and data analysis was done in the Infinicyt software (Cytognos, Salamanca, Spain). We set compaRe’s n to five for this assay.

Quality control (QC)

Request a detailed protocol

Multiwell plate heatmaps of medians come in handy in QC to reveal issues such as signal and cell viability drifts occurring during screening. However, as a typical heatmap has an equally spaced color palette, small but significant differences between wells may be obscured and not visible. Therefore, we normalized the color palette by the distribution of the medians. Also, before clustering, we removed outliers in the negative controls that were different from the others in terms of similarity values measured by compaRe.

Correcting signal and cell viability drifts

Request a detailed protocol

Depending on the protocol by which wells are processed, time may become a major concern so that some specific wells may have lower or higher values than expected. To correct for these sources of bias, we employed a two-step correction: intra-plate shift (signal drift) correction and inter-plate shift (batch effect) correction. For a given plate, we first fit a linear regression model and then vertically translate points (well values) with respect to the learned line as it rotates to the slope zero. After correcting for the intra-plate bias, the inter-plate bias is corrected by aligning medians of the plates, that is, translating to a common baseline.

Similarity calculation using dynamic gridding

Request a detailed protocol

To measure the similarity between two datasets, compaRe divides each dimension into n subsets for each dataset individually so that a dataset with d dimensions (markers) will be gridded into at most nd spatial units called hypercubes. compaRe grids only the part of the space encompassing data points, avoiding empty regions. It then measures the proportion of data points for either dataset within each of the corresponding hypercubes. The difference between the two proportions is indicative of the similarity within that relative spatial position represented by each hypercube. The similarity in the exclusive hypercubes is considered 0. We employed local outlier factor (Breunig et al., 2000) for anomaly detection and removing noise cells. Averaging these differences across all the hypercubes indicates the amount of similarity between the two datasets.

compaRe captures the configuration of data enabling it to measure similarity even without correcting for signal drift or batch effect (Appendix 1). This way, two technical replicates analyzed by two different instruments or configurations suffering from signal shift will still have the highest similarity. To generate a similarity matrix of multiple input samples, compaRe runs in parallel. The similarity matrix could then be used for identifying clusters of samples such as drugs with similar dose responses.

Graphical clustering of samples

Request a detailed protocol

To cluster samples, we developed a graphical clustering algorithm in which initially all nodes (samples) are connected forming a weighted complete graph wherein edges represent similarity between nodes. This graph is then pruned to remove potential false positive edges for a given cutoff inferred from negative controls. The optimal cutoff turns out to be the minimum weight in the maximum spanning tree of negative control nodes. After pruning, some samples may end up being connected to the negative controls (biologically inactive agents) and some disconnected (active agents). After constructing this graph, clustering is tantamount to finding maximal cliques among potent agents. In addition to maximal cliques, it also reports communities (a clique is a subset of a community). Communities can be seen as loose clusters. In a community, unlike a clique, similarity is not necessarily transitive meaning that if A is similar to B and B is similar to C, A is not necessarily similar to C. If these were three drugs within a community, concluding they had an equal response was not necessarily right unless they would form a clique.

Dispersion graph and Dispersion map

Request a detailed protocol

compaRe visualizes the similarity of samples in the form of a dispersion graph by constructing their maximum spanning tree (Appendix 1, Appendix 1—figure 30). compaRe also uses UMAP to represent a dispersion map of clusters. The map is constructed using the centroid (median) of each clique. An informative map shows different groups by coloring the centroids according to their value. These groups are mostly the identified communities the cliques come from.

Availability of data

Request a detailed protocol

Mass cytometry datasets were downloaded from Cytobank Community with the experiment ID 44185. AML mouse and human high-throughput flow cytometry data have been deposited in FLOWRepository with the repository IDs FR-FCM-Z357 and FR-FCM-Z3DP respectively. Flow cytometry data of AML and MDS patients have been deposited in FLOWRepository with the repository ID FR-FCM-Z3ET. Acquisition, installation and more technical details are available in compaRe’s online tutorial on (https://github.com/morchalabi/COMPARE-suite, swh:1:rev:df2feaf6aa982e0f6f077eb85f26acce6bb61063, Chalabi, 2022b). Similarity measurement and clustering modules as stand-alone tools have been merged into a separate R package and are available for download at (https://github.com/morchalabi/compaRe, swh:1:rev:594106b1e34c17b405064f1a0f9fb39975a4ec79, Chalabi, 2022a).

Appendix 1

High-throughput flow cytometry of AML mouse model

Leukemic spleen cells were sorted for c-Kit positivity from Npm1+/cA; Flt3+/ITD; Dnmt3a+/-; Mx1-Cre+ moribund mice. Shortly, c-Kit+ splenic cells were expanded for two passages in StemPro-34 SFM media (Gibco) with 100 µM 2-Mercaptoethanol (Gibco), 20 ng/ml murine SCF, 10 ng/ml murine IL-3 and 10 ng/ml IL-6 added (Peprotech), with complete media change every two/three days. Aliquots of one million cells were frozen down in 90% media 10% DMSO. Frozen aliquots were taken up and expanded for one week before drug screening. 5000 cells in 25 µl of media per well was seeded into 384-well plates (Greiner) containing a library of 116 compounds (Supplementary file 4) in a five-point concentration range. After 72 h incubation at 37 °C, 15 µl of medium was aspirated from each well and antibodies (Supplementary file 9) were added to drug plates using acoustic dispensing. Plates were incubated 40 min at RT, covered from light. Next, dead cell dye 7-AAD (BD) was added, and samples were read using a high-throughput flow cytometer iQue Screener Plus (Intellicyt). To remove noise from the data by excluding the most broadly toxic treatments, doublets and dead cells were omitted (Appendix 1—figure 31) and only samples with at least 1,000 live cells were selected for further analyses (selected 465 wells out of 600).

High-throughput flow cytometry of human AML

Donated MNCs from human bone marrow aspirates (Danish National Ethical committee/National Videnskabsetisk Komité permit 1705391) were thawed and allowed to rest overnight in assay media: StemSpan II-SFEM (StemCell), 100 U/ml penicillin/streptomycin (Thermo), including the following human recombinant cytokines from Preprotech (unless otherwise stated), 50 ng/ml Flt3 ligand (StemCell), 10 ng/ml IL3, 10 ng/ml IL-1beta, 20 ng/ml IL6, 20 ng/ml G-CSF, 20 ng/ml GM-CSF, and 10 ng/ml SCF, and the following compounds diluted in DMSO (Merck) 1 µM UM729 (Selleckchem) and 500 nM StemRegenin-1 (MedChemExpress). Before being counted and re-suspended in fresh assay media at a density of 5 × 105 cells/ml. A 20 µl/well was plated in 384-well conical bottom plates (Greiner Bio-One) containing 25 nl of compounds (Supplementary file 7) in DMSO. After 72 hr incubation at 37 °C, 95% RH, 5% CO2 antibodies and viability dye were added to the plates using acoustic dispensing (Echo, Labcyte). Plates were incubated for 1.5 hr covered from light at RT. The samples were then run on an iQue Screener Plus (Intellicyt) high-throughput flow cytometer. The data was gated to remove noise, doublets, and dead cells (Appendix 1—figure 30). The antibodies and stains used are described in Supplementary file 10.

Signal and cell viability drifts correction in compaRe

To correct signal drift, we employed a two-step correction: intra-plate correction and inter-plate correction. For a given plate, we first fit a linear regression model and then vertically translate points (MFIs) with respect to the leaned line as it rotates to slope zero. This is because the relative distance between the points must be retained as much as possible, and no point must be translated to x+y quadrant after correction. To make sure the learned line is not affected by outliers, we first removed them using the interquartile range. In this way, a point at y,x is translated to ybmx+b,x after intra-plate correction. The correction coefficient bmx+b derives from the ratio of y-coordinates of any point on the regression line before and after translation: yy=bmx+b where y is translated y, m is the slope and b is the intercept of the line. This ratio holds true for all other points in the xy-plane.

After correcting for intra-plate signal drift, inter-plate signal drift is corrected by aligning MFI medians of the plates, that is, translating to a common baseline. Let b be the baseline, and b be the median of corrected MFIs in a plate, then the inter-plate correction coefficient is given by bb , and a point at y,x is translated to (ybb,x) . The same approach is employed for correcting cell viability bias (Appendix 1—figure 28).

Similarity measurement in compaRe

compaRe can measure the similarity between two datasets with many variables (dimensions) and observations (data points). compaRe divides each dimension into n subsets so that a dataset with d dimensions will be divided into at most nd spatial units called hypercubes. The hypercubes are formed for either dataset individually. It, then measures the proportion of the observations within each of the corresponding hypercubes. The difference between the two proportions is indicative of the similarity within that relative spatial position represented by that hypercube so that for two similar datasets this difference is near zero in the majority of the hypercubes. Averaging these differences across all the hypercubes indicates the amount of similarity between the two datasets.

It is important to compare two samples across their corresponding hypercubes representing the same relative spatial positions. This means a universal numbering rule is required to ensure having corresponding hypercubes for the two samples in the end. This problem can be modeled as a tree that at each level l (dimension) grows nl new branches (divisions) (Appendix 1—figure 32). However, as the number of branches increases exponentially with l, implementing the tree is infeasible. To overcome this problem, we instead employed a dynamic algorithm in which the hypercube number of each observation is dynamically updated at each iteration. In this approach, the child node number must be found from its parent’s, that is, previous iteration.

Rewriting the branch numbers to include more information reveals that if rl-1=n0++nl-2+fl-1+sl-1nl-2 is the parent node’s number, the child node’s number will be rl=n0++nl-2+nl-1+nfl-1+sl-1+slnl-1 where l is the child’s level, fl=0,,nl-1-1 is the number of families behind, and sl=0,...,l-1 is the number of siblings behind. Therefore, to find child node rl, we first need to calculate fl-1 and sl-1 of its parent as follows:

(1) sl-1=rl-1-n++nl-2nl-2
(2) fl-1=rl-1-n0++nl-2-sl-1nl-2

It can be noticed that rl-1 and sl are always known, fl=nfl-1+sl-1 , and n0+...+nl-1-1 is actually the largest node number at the l th level. Therefore, the problem we need to dynamically solve for each child at each dimension as the tree grows is:

(3) rl=n0++nl-1+fl+slnl-1

Since the similarity metric decreases for each exclusive hypercube, it is important to rid the two samples of outliers lying significantly far from the subpopulations of observations. However, at the same time we need to make sure smaller subpopulations (like rare cell subpopulations) are not mistaken for outliers. We employed local outlier factor which is a powerful tool for anomaly detection. Figure 1 shows an actual AML dataset with three surface markers dissected by compaRe wherein each distinct color corresponds to data points within one abstract hypercube.

compaRe captures the morphology of high dimensional data enabling it to measure similarity even in the presence of moderate signal shift. For example, two technical replicates analyzed by two different instruments or configurations suffering from signal shift will still have the highest similarity by compaRe unless the shift is severe or has modified the morphology of the cell populations which practically does not happen as a result of batch effect or signal drift. This strategy helps compaRe circumvent signal drift or batch effect left uncorrected. Considering that any signal drift correction is essentially an approximate method, this feature is an advantage for compaRe, because together with the correction method they create a synergistic effect.

compaRe is a mass-aware approach meaning it forms hypercubes only around concentrations of data points avoiding areas which are devoid of data points. This substantially speeds up the process by saving a lot of CPU time and memory space making it feasible to compare datasets with numerous variables. As an example, dividing each dimension blindly into just three regions yields more than 1.5 billion regions for consideration for a dataset with as few as 19 surface markers. In practice, however, it turns out many of these regions are empty so using a mass-aware gridding instead of blind gridding improves the comparison complexity from Θ(nd) to Ond (asymptotic notations to represent algorithmic complexity). Even if no region is empty, since compaRe benefits from dynamic programming, it can still finish the process quite fast. Changing n tunes the level of smoothing so that a value between 3–5 works for most assays.

Dynamic programming is key for reducing processing power. In general, the goal is to bin/grid data into relative expression groups (hypercubes). Gridding can be implemented by a simple algorithm dividing each dimension in each iteration. However, as pointed out above, after a couple of rounds, this naïve algorithm turns out to be infeasible. Therefore, one need a more efficient algorithm for implementing gridding. Dynamic programming turned out to be quite effective. What makes dynamic programming very effective is its ability to memorize the values computed in the previous iterations avoiding recomputing potentially expensive algebraic operations (Appendix 1-equation 3).

To generate a similarity matrix of multiple input samples, compaRe runs in parallel for the samples in the upper-triangular submatrix using a multithreading approach. The similarity matrix could then be used for identifying clusters of samples such as drugs with similar dose responses like predicting the mechanism of action of drugs in development.

Appendix 1—figure 1
Performance of meta-clustering with SPADE FlowSOM in the presence of batch effect.

Similarity matrices generated by FlowSOM and SPADE are shown in (a) and (b) respectively. Size and color of dots represent the level of similarity. Self-comparisons were removed. Noise was added (marked with *) to the original 21 mass cytometry samples of bone marrow aspirates from 16 pediatric AML patients (S) and five healthy adult donors (H).

Appendix 1—figure 2
Phenotypic characterization in a high-parameter heterogeneous population of cell types.

Cells from a target cluster (an immunophenotypic cell population) were gradually removed to contort its configuration. We used a dataset of 3 healthy and two pediatric AML bone marrow mononuclear cell samples from the data provided in the 6th reference. Samples were stained with 29 (15 membrane and 14 intracellular signaling) markers. Taking H1 as reference, we gradually removed 25%, 50%, 75% and 100% (phenotypic changes) of cells from the target cluster identified by PhenoGraph. To capture the higher heterogeneity harbored in the AML samples, we set compaRe's n to 4 while we set it to 3 for healthy samples. Each column was scaled individually retaining mutual differences.

Appendix 1—figure 3
UMAP projections of A2 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 4
UMAP projections of A3 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 5
UMAP projections of A4 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 6
UMAP projections of A5 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 7
UMAP projections of A6 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 8
UMAP projections of A7 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 9
UMAP projections of A8 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 10
UMAP projections of A9 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 11
UMAP projections of H1 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 12
UMAP projections of H2 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 13
UMAP projections of H3 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 14
UMAP projections of H4 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 15
UMAP projections of H5 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 16
UMAP projections of L1 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 17
UMAP projections of L2 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 18
UMAP projections of L3 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 19
UMAP projections of L4 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 20
UMAP projections of L5 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 21
UMAP projections of L6 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 22
UMAP projections of L7 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 23
UMAP projections of L8 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 24
UMAP projections of L9 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 25
UMAP projections of L10 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 26
UMAP projections of L11 sample against all other patient samples.

From top left to bottom right, the similarity measured by compaRe decreases as the degree of overlap decreases and the number of exclusive cell populations increases.

Appendix 1—figure 27
Band plots of AML and MDS patient samples.

The immunophenotype of each patient sample is shown in a multiparameter band-dot plot (HrMDS: higher-MDS, LrMDS: lower-MDS). Rectangles gate positive and/or negative populations of monocytic maturation markers as well as the CD11b marker.

Appendix 1—figure 28
Correcting signal and cell viability drift.

(a) Intra- and inter-plate signal drift correction. Accumulation of green-blue tiles in the bottom right corner of the left heatmap shows signal drift in CD11b expression for drugs in plate 1. Sorting median expressions (MFIs) of wells into reading order (column wise left to right) reveals a linear slope. After correction, the slope becomes non-positive (intra-plate correction). Still, there are different baselines between the two plates. Matching median lines of corrected values of all plates correct for this bias (inter-plate correction). (b) Intra- and inter-plate cell viability correction. Accumulation of green tiles in the bottom right corner of the left heatmap shows cell viability drift (7-AAD marker). We follow similar steps with (a) for cell viability correction.

Appendix 1—figure 29
Birabresib treatment leads to loss of CD45 and CD45-RA expression without loss of live cell numbers.

(a) Birabresib response as density scatter plot, CD45 vs CD45-RA. (b) Count of live cells per different concentrations of lenalidomide and birabresib. (c) Heatmap of birabresib response in all marker channels. (d) Example of response group 3: density scatter plots of DMSO-control vs. tretinoin 375 nM in different marker channels.

Appendix 1—figure 30
Dispersion graph.

The (maximum spanning) tree demonstrates the dispersion of tens of potent antineoplastic agents around the control node containing negative controls (DMSO) and impotent agents. The drug library was analyzed by high-throughput flow cytometry coupled with compaRe in an AML human sample. Edge color and label show the amount of similarity between the agents. Impotent drugs are those which were similar enough to negative controls for a cutoff inferred during clustering. As the tree branches and spreads, drugs with stronger potency, usually with higher doses, tend to lie farther from the control node. Using the graph, the investigator can easily pick potent agents such as hits. The graph may also be potentially used to investigate different paths for mechanism of action, leading to different branches.

Appendix 1—figure 31
Removal of noise, dead cells and doublet cells from mouse and human AML sample drug screening data.

(a) AML mouse model drug screening. (b) AML human sample drug screening. Cells were separated from debris using a side scatter height (SSC-H) vs forward scatter height (FSC-H) plot. Singlet cells were determined from FSC-H vs forward scatter area (FSC-A) plot. Live cells were separated from dead cells using a dead-cell-labelling dye, either 7-AAD or DRAQ7.

Appendix 1—figure 32
Demonstration of compaRe algorithm using a 2-dimensional table.

It first forms an abstract square grid (red) encompassing all the data points within the range (1.1, 9.6). At the top level, all the cells (table rows) are in the region number (RN) 0. First iteration divides the first dimension formed by CD1 marker into 3 ( = n) subsets. Assuming a left-first numbering rule, the RN column is dynamically updated (blue column) for each subset using some information such as current RN (grey column), current dimension and possible number of families and siblings behind. For instance, child node 12 has parent node 3, could have two siblings (node 6, node 9) and two families (parent 1, parent 2) behind, although children 11 and 9 were never born as marked with ☓. Final leaves are called hypercubes (HCs). The corresponding grid on the biplot demonstrates that two regions which were devoid of data points have not been assigned any hypercube. For comparing two samples, they are first jointly normalized between a range. The tree graph is just for better visualization and will not be implemented.

Data availability

Mass cytometry datasets were downloaded from Cytobank Community with the experiment ID 44185. AML mouse and human high-throughput flow cytometry data have been deposited in FLOWRepository with the repository IDs FR-FCM-Z357 and FR-FCM-Z3DP respectively. Flow cytometry data of AML and MDS patients have been deposited in FLOWRepository with the repository ID FR-FCM-Z3ET. Acquisition, installation and more technical details are available in compaRe's online tutorial on (https://github.com/morchalabi/COMPARE-suite, (copy archived at swh:1:rev:df2feaf6aa982e0f6f077eb85f26acce6bb61063)). Similarity measurement and clustering modules as stand-alone tools have been merged into a separate R package and are available for download at (https://github.com/morchalabi/compaRe, (copy archived at swh:1:rev:594106b1e34c17b405064f1a0f9fb39975a4ec79)).

The following data sets were generated
    1. Morteza CH
    (2022) GitHub
    ID github.com/morchalabi/compaRe. Comprehensive and unbiased multiparameter high-throughput screening by compaRe finds effective and subtle drug responses in AML models.
The following previously published data sets were used
    1. Levine JH
    2. Simonds EF
    3. Bendall SC
    4. Davis KL
    5. Amir el AD
    6. Tadmor MD
    (2015) Cytobank
    ID 44185. Data-Driven Phenotypic Dissection of AML Reveals Progenitor-like Cells that Correlate with Prognosis.

References

  1. Book
    1. BioScience E
    (2020)
    Eliminating Data Analysis Bottlenecks with IQue Forecyt Software
    Essen BioScience.
    1. Cossarizza A
    2. Chang H-D
    3. Radbruch A
    4. Acs A
    5. Adam D
    6. Adam-Klages S
    7. Agace WW
    8. Aghaeepour N
    9. Akdis M
    10. Allez M
    11. Almeida LN
    12. Alvisi G
    13. Anderson G
    14. Andrä I
    15. Annunziato F
    16. Anselmo A
    17. Bacher P
    18. Baldari CT
    19. Bari S
    20. Barnaba V
    21. Barros-Martins J
    22. Battistini L
    23. Bauer W
    24. Baumgart S
    25. Baumgarth N
    26. Baumjohann D
    27. Baying B
    28. Bebawy M
    29. Becher B
    30. Beisker W
    31. Benes V
    32. Beyaert R
    33. Blanco A
    34. Boardman DA
    35. Bogdan C
    36. Borger JG
    37. Borsellino G
    38. Boulais PE
    39. Bradford JA
    40. Brenner D
    41. Brinkman RR
    42. Brooks AES
    43. Busch DH
    44. Büscher M
    45. Bushnell TP
    46. Calzetti F
    47. Cameron G
    48. Cammarata I
    49. Cao X
    50. Cardell SL
    51. Casola S
    52. Cassatella MA
    53. Cavani A
    54. Celada A
    55. Chatenoud L
    56. Chattopadhyay PK
    57. Chow S
    58. Christakou E
    59. Čičin-Šain L
    60. Clerici M
    61. Colombo FS
    62. Cook L
    63. Cooke A
    64. Cooper AM
    65. Corbett AJ
    66. Cosma A
    67. Cosmi L
    68. Coulie PG
    69. Cumano A
    70. Cvetkovic L
    71. Dang VD
    72. Dang-Heine C
    73. Davey MS
    74. Davies D
    75. De Biasi S
    76. Del Zotto G
    77. Dela Cruz GV
    78. Delacher M
    79. Della Bella S
    80. Dellabona P
    81. Deniz G
    82. Dessing M
    83. Di Santo JP
    84. Diefenbach A
    85. Dieli F
    86. Dolf A
    87. Dörner T
    88. Dress RJ
    89. Dudziak D
    90. Dustin M
    91. Dutertre C-A
    92. Ebner F
    93. Eckle SBG
    94. Edinger M
    95. Eede P
    96. Ehrhardt GRA
    97. Eich M
    98. Engel P
    99. Engelhardt B
    100. Erdei A
    101. Esser C
    102. Everts B
    103. Evrard M
    104. Falk CS
    105. Fehniger TA
    106. Felipo-Benavent M
    107. Ferry H
    108. Feuerer M
    109. Filby A
    110. Filkor K
    111. Fillatreau S
    112. Follo M
    113. Förster I
    114. Foster J
    115. Foulds GA
    116. Frehse B
    117. Frenette PS
    118. Frischbutter S
    119. Fritzsche W
    120. Galbraith DW
    121. Gangaev A
    122. Garbi N
    123. Gaudilliere B
    124. Gazzinelli RT
    125. Geginat J
    126. Gerner W
    127. Gherardin NA
    128. Ghoreschi K
    129. Gibellini L
    130. Ginhoux F
    131. Goda K
    132. Godfrey DI
    133. Goettlinger C
    134. González-Navajas JM
    135. Goodyear CS
    136. Gori A
    137. Grogan JL
    138. Grummitt D
    139. Grützkau A
    140. Haftmann C
    141. Hahn J
    142. Hammad H
    143. Hämmerling G
    144. Hansmann L
    145. Hansson G
    146. Harpur CM
    147. Hartmann S
    148. Hauser A
    149. Hauser AE
    150. Haviland DL
    151. Hedley D
    152. Hernández DC
    153. Herrera G
    154. Herrmann M
    155. Hess C
    156. Höfer T
    157. Hoffmann P
    158. Hogquist K
    159. Holland T
    160. Höllt T
    161. Holmdahl R
    162. Hombrink P
    163. Houston JP
    164. Hoyer BF
    165. Huang B
    166. Huang F-P
    167. Huber JE
    168. Huehn J
    169. Hundemer M
    170. Hunter CA
    171. Hwang WYK
    172. Iannone A
    173. Ingelfinger F
    174. Ivison SM
    175. Jäck H-M
    176. Jani PK
    177. Jávega B
    178. Jonjic S
    179. Kaiser T
    180. Kalina T
    181. Kamradt T
    182. Kaufmann SHE
    183. Keller B
    184. Ketelaars SLC
    185. Khalilnezhad A
    186. Khan S
    187. Kisielow J
    188. Klenerman P
    189. Knopf J
    190. Koay H-F
    191. Kobow K
    192. Kolls JK
    193. Kong WT
    194. Kopf M
    195. Korn T
    196. Kriegsmann K
    197. Kristyanto H
    198. Kroneis T
    199. Krueger A
    200. Kühne J
    201. Kukat C
    202. Kunkel D
    203. Kunze-Schumacher H
    204. Kurosaki T
    205. Kurts C
    206. Kvistborg P
    207. Kwok I
    208. Landry J
    209. Lantz O
    210. Lanuti P
    211. LaRosa F
    212. Lehuen A
    213. LeibundGut-Landmann S
    214. Leipold MD
    215. Leung LYT
    216. Levings MK
    217. Lino AC
    218. Liotta F
    219. Litwin V
    220. Liu Y
    221. Ljunggren H-G
    222. Lohoff M
    223. Lombardi G
    224. Lopez L
    225. López-Botet M
    226. Lovett-Racke AE
    227. Lubberts E
    228. Luche H
    229. Ludewig B
    230. Lugli E
    231. Lunemann S
    232. Maecker HT
    233. Maggi L
    234. Maguire O
    235. Mair F
    236. Mair KH
    237. Mantovani A
    238. Manz RA
    239. Marshall AJ
    240. Martínez-Romero A
    241. Martrus G
    242. Marventano I
    243. Maslinski W
    244. Matarese G
    245. Mattioli AV
    246. Maueröder C
    247. Mazzoni A
    248. McCluskey J
    249. McGrath M
    250. McGuire HM
    251. McInnes IB
    252. Mei HE
    253. Melchers F
    254. Melzer S
    255. Mielenz D
    256. Miller SD
    257. Mills KHG
    258. Minderman H
    259. Mjösberg J
    260. Moore J
    261. Moran B
    262. Moretta L
    263. Mosmann TR
    264. Müller S
    265. Multhoff G
    266. Muñoz LE
    267. Münz C
    268. Nakayama T
    269. Nasi M
    270. Neumann K
    271. Ng LG
    272. Niedobitek A
    273. Nourshargh S
    274. Núñez G
    275. O’Connor J
    276. Ochel A
    277. Oja A
    278. Ordonez D
    279. Orfao A
    280. Orlowski-Oliver E
    281. Ouyang W
    282. Oxenius A
    283. Palankar R
    284. Panse I
    285. Pattanapanyasat K
    286. Paulsen M
    287. Pavlinic D
    288. Penter L
    289. Peterson P
    290. Peth C
    291. Petriz J
    292. Piancone F
    293. Pickl WF
    294. Piconese S
    295. Pinti M
    296. Pockley AG
    297. Podolska MJ
    298. Poon Z
    299. Pracht K
    300. Prinz I
    301. Pucillo CEM
    302. Quataert SA
    303. Quatrini L
    304. Quinn KM
    305. Radbruch H
    306. Radstake TRDJ
    307. Rahmig S
    308. Rahn H-P
    309. Rajwa B
    310. Ravichandran G
    311. Raz Y
    312. Rebhahn JA
    313. Recktenwald D
    314. Reimer D
    315. Reis e Sousa C
    316. Remmerswaal EBM
    317. Richter L
    318. Rico LG
    319. Riddell A
    320. Rieger AM
    321. Robinson JP
    322. Romagnani C
    323. Rubartelli A
    324. Ruland J
    325. Saalmüller A
    326. Saeys Y
    327. Saito T
    328. Sakaguchi S
    329. Sala-de-Oyanguren F
    330. Samstag Y
    331. Sanderson S
    332. Sandrock I
    333. Santoni A
    334. Sanz RB
    335. Saresella M
    336. Sautes-Fridman C
    337. Sawitzki B
    338. Schadt L
    339. Scheffold A
    340. Scherer HU
    341. Schiemann M
    342. Schildberg FA
    343. Schimisky E
    344. Schlitzer A
    345. Schlosser J
    346. Schmid S
    347. Schmitt S
    348. Schober K
    349. Schraivogel D
    350. Schuh W
    351. Schüler T
    352. Schulte R
    353. Schulz AR
    354. Schulz SR
    355. Scottá C
    356. Scott-Algara D
    357. Sester DP
    358. Shankey TV
    359. Silva-Santos B
    360. Simon AK
    361. Sitnik KM
    362. Sozzani S
    363. Speiser DE
    364. Spidlen J
    365. Stahlberg A
    366. Stall AM
    367. Stanley N
    368. Stark R
    369. Stehle C
    370. Steinmetz T
    371. Stockinger H
    372. Takahama Y
    373. Takeda K
    374. Tan L
    375. Tárnok A
    376. Tiegs G
    377. Toldi G
    378. Tornack J
    379. Traggiai E
    380. Trebak M
    381. Tree TIM
    382. Trotter J
    383. Trowsdale J
    384. Tsoumakidou M
    385. Ulrich H
    386. Urbanczyk S
    387. van de Veen W
    388. van den Broek M
    389. van der Pol E
    390. Van Gassen S
    391. Van Isterdael G
    392. van Lier RAW
    393. Veldhoen M
    394. Vento-Asturias S
    395. Vieira P
    396. Voehringer D
    397. Volk H-D
    398. von Borstel A
    399. von Volkmann K
    400. Waisman A
    401. Walker RV
    402. Wallace PK
    403. Wang SA
    404. Wang XM
    405. Ward MD
    406. Ward-Hartstonge KA
    407. Warnatz K
    408. Warnes G
    409. Warth S
    410. Waskow C
    411. Watson JV
    412. Watzl C
    413. Wegener L
    414. Weisenburger T
    415. Wiedemann A
    416. Wienands J
    417. Wilharm A
    418. Wilkinson RJ
    419. Willimsky G
    420. Wing JB
    421. Winkelmann R
    422. Winkler TH
    423. Wirz OF
    424. Wong A
    425. Wurst P
    426. Yang JHM
    427. Yang J
    428. Yazdanbakhsh M
    429. Yu L
    430. Yue A
    431. Zhang H
    432. Zhao Y
    433. Ziegler SM
    434. Zielinski C
    435. Zimmermann J
    436. Zychlinsky A
    (2019) Guidelines for the use of flow cytometry and cell sorting in immunological studies (second edition)
    European Journal of Immunology 49:1457–1973.
    https://doi.org/10.1002/eji.201970107
  2. Book
    1. DeVita VT
    2. Lawrence TS
    (2015)
    Devita, Hellman, and Rosenberg’s Cancer: Principles & Practice of Oncology
    Williams & Wilkins Publishers.
    1. Lvd M
    2. Hinton G
    (2008)
    Visualizing Data using t-SNE
    Journal of Machine Learning Research 9:2579–2605.
  3. Book
    1. Niederhuber JE
    2. Armitage J
    3. Doroshow J
    4. Kastan M
    5. Tepper J
    (2020)
    Abeloff’s Clinical Oncology
    Abeloff Leo.

Decision letter

  1. Aleksandra M Walczak
    Senior and Reviewing Editor; CNRS LPENS, France
  2. Aik Choon Tan
    Reviewer; Moffitt Cancer Center, United States

Our editorial process produces two outputs: i) public reviews designed to be posted alongside the preprint for the benefit of readers; ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "Comprehensive and unbiased multiparameter high-throughput screening by compaRe finds effective and subtle drug responses in AML models" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Aleksandra Walczak as the Senior Editor. The following individual involved in review of your submission has agreed to reveal their identity: Aik Choon Tan (Reviewer #2).

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

Please address all the reviewers comments. The reviewers agree it is a strong paper that is well written but certain improvements can still be made. Especially please give the reader more context in the introduction and compare your method with other existing methods.

Reviewer #1:

Comprehensive and unbiased multiparameter high-throughput screening by compaRe finds effective and subtle drug responses in AML models by Hajkariim et al. introduces a pipeline for pre-processing and analyzing data from multiplex flow cytometry and other technologies. Preprocessing steps include algorithms for correcting common sources of bias in such data. Another key feature is a robust approach to measuring cell similarity across samples. Among the strengths are that the manuscript is well-written, the analysis pipeline is well-motivated, and illustrated with apt examples. The similarity measure is very interesting as well.

There are a few weaknesses as well. It is not completely clear to me how this pipeline agrees and disagrees with common practice in the field. References 1-3, cited to document ongoing analytic challenges, are all at least 5 years old. Comparisons to other approaches, including the use Jensen-Shannon Divergence for similarity, make a convincing case that the proposed method is both effective and computationally efficient, but it is not clear if the comparators represent true standard of practice, or mere straw men. Methodologies are complex and can be difficult to follow, especially the similarity measure.

I would like to see the current state of the field described more clearly in the introduction, as context for the current effort. What makes this unique and important today? This type of clustering is common in single cell sequencing analysis, is it also commonly done in multiplex flow etc, maybe with less computationally demanding tools than JSD? If not, why not?

I think I understand how the similarity measure here works, though its hard to follow the details. Its very interesting, but to be honest, I can't decide if I think its a good solution or needlessly complex. The key, as I understand it, is the binning into "relative" expression groups. This is, after all, how 2-dimensional flow data is commonly interpreted – with plots treated as a visual 2 X 2 table, with row and column boundaries existing largely in the eye of the beholder. The methods need to be clearer throughout, esp. this part. It might help to add another author, a third party who has to understand the method from scratch, and who, having not lived with the details as long, might be expected to provide a more user friendly sense of the overall approach.

Reviewer #2:

In this manuscript, Hajkarim et al. developed compaRe, a user friendly software suite (written in R) for analyzing high-throughput, multi-parameter screening data. There are several modules included in the compaRe toolkit, which can be individually invoked to perform specific tasks, such as quality control, bias correction, pairwise comparisons, clustering and data visualization. All of these modules are available as command-line version and a GUI version for users to use in data analysis, visualization and results interpretation. The authors showed the utility of their toolkit in analyzing multiparameter mass and flow cytometric data from AML and MDS patient samples. Through this analysis using compaRe, the authors showed that they can identify patient heterogeneity and drug response profiles. Overall, this is a well organized and written manuscript describing the development of the new compaRe toolkit. The method is clearly described, and the user manual/tutorial is easy to follow. It seems like compaRe will be a useful toolkit for the research community, which is eager for a one-stop pipeline for analyzing high-throughout multiparameter screening data.

et al.

Strengths:

1. All of these modules are available as command-line version and a GUI version for users to use in data analysis, visualization and results interpretation.

2. The authors showed the utility of their toolkit in analyzing multiparameter mass and flow cytometric data from AML and MDS patient samples. Through this analysis using compaRe, the authors showed that they can identify patient heterogeneity and drug response profiles.

3. Overall, this is a well organized and written manuscript describing the development of the new compaRe toolkit. The method is clearly described, and the user manual/tutorial is easy to follow.

4. It seems like compaRe will be a useful toolkit for the research community, which is eager for a one-stop pipeline for analyzing high-throughout multiparameter screening data.

Weaknesses:

1. The current manuscript didn't compare with some other existing programs/software in analyzing flow and mass cytometry data. It will be important to compare compaRe with existing tools, to show the strengths and weaknesses of compaRe with other tools.

2. The authors could think about adding an additional module to integrate other "omics" data (e.g. such as mutational or gene expression/signatures or pathways), this could be useful for doing the clustering step or to identify patients having the same mutational profiles.

Reviewer #3:

Hajkarim et al. implement an algorithm in their presented toolkit compaRe to compare samples based on the similarities of samples, distinct from the more commonly used meta-clustering approaches, such as PhenoGraph, or dimensional reduction with Jenssen-Shannon Divergence analysis. Similarities among samples are calculated based on the proportions of cells within a sample belonging to an n-dimensional "hypercubes" (or "hypergridding" that is actually mass-aware and not blind) that are stratified by expression levels for n number of markers. The authors demonstrate that this method is much more time-efficient, obviates subsampling, and is robust to batch effects. This method is particularly appropriate for large-scale datasets, facilitating the comparison of numerous samples which would be helpful in screening efforts. The manuscript is written and presented well.

Major strengths:

1. The study demonstrates sufficiently strong support for the toolkit's ability to determine similarity across samples and its computing efficiency with Figure 2, an important advantage of this tool.

2. Compared to other approaches, the method is advantageous for identifying groups of samples that may be similar in a very large-scale dataset. CompaRe does not require (or make use of) manual expert annotation of meta-clusters. The workflow is efficient and unbiased.

Major weakness:

A major weakness of the current presentation of the study is that it has not clearly demonstrated the toolkit's utility in exploring specific phenotypes in-depth within a high-parameter dataset. The following are two examples in which this limitation is relevant, and the authors may address this to strengthen this paper, if in fact detailed phenotyping is considered by the authors as an important feature of the toolkit. If not, the authors should revise the manuscript as such.

First, the authors stated that their approach can be used to "optimize true cytometric n-dimensional immunophenotypic characterization" even in the setting of multi-panel workflow. However, the demonstration was based on samples that seem to have predominant phenotypes that are almost mutually exclusive. It is unlikely that this toolkit would be useful for reliable phenotypic characterization in a largely heterogeneous population of cell types, e.g. even in normal peripheral blood, unless a high number of parameters was concurrently acquired within the dataset. This is an inherent limitation. Nonetheless, a revision to demonstrate how compaRe can evaluate specific clusters of phenotypes with biological significance from a high-parameter dataset (20-30 marker cytometry) would be very helpful.

Second, the authors refer to the method's ability to include all rare cell subsets in the analysis, i.e. the ability to forego any subsampling. The work does not, however, demonstrate clearly how the presence of a rare cell subset in a given sample influences its similarity to other samples. Thus, the toolkit's value of being able to include such a rare subset in the analysis remains unestablished. It would be beneficial to include such a test to see whether the algorithm is sensitive to such changes.

https://doi.org/10.7554/eLife.73760.sa1

Author response

Essential revisions:

Please address all the reviewers’ comments. The reviewers agree it is a strong paper that is well written but certain improvements can still be made. Especially please give the reader more context in the introduction and compare your method with other existing methods.

Reviewer #1:

[…] There are a few weaknesses as well. It is not completely clear to me how this pipeline agrees and disagrees with common practice in the field. References 1-3, cited to document ongoing analytic challenges, are all at least 5 years old. Comparisons to other approaches, including the use Jensen-Shannon Divergence for similarity, make a convincing case that the proposed method is both effective and computationally efficient, but it is not clear if the comparators represent true standard of practice, or mere straw men. Methodologies are complex and can be difficult to follow, especially the similarity measure.

I would like to see the current state of the field described more clearly in the introduction, as context for the current effort. What makes this unique and important today? This type of clustering is common in single cell sequencing analysis, is it also commonly done in multiplex flow etc, maybe with less computationally demanding tools than JSD? If not, why not?

Thank you for your comments. We believe compaRe is a unique tool in that it is designed specifically to cluster samples. Other methods that are commonly used for the analysis of mass/flow cytometry data are primarily designed to cluster cells within samples. The obtained cell clustering information is subsequently used for assessing similarity between samples by meta-clustering. Also, using clustering-free algorithms is often computationally expensive and not as robust as one may wish. For example, JSD was employed to compare samples, but it is not efficient for comparing a large number of samples, and, as we showed in the results, is computationally expensive.

We had previously highlighted these limitations by comparing the outcome of a compaRe analysis to meta-clustering and JSD. We have now provided the comparisons with other well-established and commonly used cell clustering approaches for mass/flow cytometry: FlowSOM and SPADE (main text lines 113-115 and Appendix 1-figure 31). To make these points clearer, we have modified the Introduction. Specifically, we added that:

i. Main text lines 33-35: ‘These technologies can test hundreds of samples (such as drug treatments) each with tens of thousands of events (e.g., cells) labeled for numerous biomarkers (such as cytoplasmic or membrane markers).’

ii. Main text lines 43-48: ‘Meta-clustering with single-cell clustering algorithms has been suggested to cluster samples based on the similarity of the centroids of cell subpopulations identified in the individual samples (5-8). While these algorithms are widely used in single-cell data analysis for clustering cells, they are not efficient for clustering of samples. This is because centroid-based analysis can be misleading when subclusters are not sufficiently distinct or the number of sub-clusters varies.’

iii. Main text 56-58: ‘Available computational toolkits mostly allow for single-parameter or unautomated analyses of large-scale screening data using the aforementioned methods. In these toolkits, each well should be first represented by a single parameter such as cell counts or centroids or they require manual intervention.’

I think I understand how the similarity measure here works, though it’s hard to follow the details. It’s very interesting, but to be honest, I can't decide if I think it’s a good solution or needlessly complex. The key, as I understand it, is the binning into "relative" expression groups. This is, after all, how 2-dimensional flow data is commonly interpreted – with plots treated as a visual 2 X 2 table, with row and column boundaries existing largely in the eye of the beholder. The methods need to be clearer throughout, esp. this part. It might help to add another author, a third party who has to understand the method from scratch, and who, having not lived with the details as long, might be expected to provide a more user friendly sense of the overall approach.

Dynamic programming is key for reducing processing power. We try to explain the concept more clearly as follows (also added to the Appendix 1 lines 114-120). As you pointed out, the goal is to bin/grid data into relative expression groups (hypercubes). Gridding can be implemented by a simple algorithm dividing each dimension in each iteration. However, after a couple of rounds, this naïve algorithm turns out to be infeasible. This can be shown through visualizing the gridding algorithm by a tree graph (Appendix 1-figure 29). The tree clearly shows that the comparison complexity of the gridding algorithm is exponential since the tree grows exponentially. Therefore, we were forced to use a more efficient algorithm for implementing gridding. Dynamic programming (algorithm) turned out to be quite effective. What makes dynamic programming very effective is its ability to memorize the values computed in the previous iterations avoiding recomputing potentially expensive algebraic operations (Appendix 1-equation 3). Moreover, considering that any signal drift correction is essentially an approximate method, the robustness of this algorithm to noise creates a synergistic effect. Therefore, we believe that the current description about similarity measure is required.

Reviewer #2:

[…] 1. The current manuscript didn't compare with some other existing programs/software in analyzing flow and mass cytometry data. It will be important to compare compaRe with existing tools, to show the strengths and weaknesses of compaRe with other tools.

Thank you for your comments. We previously benchmarked compaRe’s similarity module against meta-clustering with several cell clustering methods and selected widely used PhenoGraph which performed better than others (considering several factors). We also used JSD as a clustering-free approach. In response to the request from the reviewer, we have now added the comparisons with SPADE and FlowSOM, two other commonly used algorithms for the analysis of flow and mass cytometry data (main text lines 113-115 and Appendix 1-figure 31). However, please note that available automated and commonly used tools are primarily cell clustering/automatic gating approaches. These include methods such as PhenoGraph, SPADE, and FlowSOM. A key reason why compaRe outperforms these and other methods is that they are designed to cluster cells rather than samples, while compaRe clusters samples.

2. The authors could think about adding an additional module to integrate other "omics" data (e.g. such as mutational or gene expression/signatures or pathways), this could be useful for doing the clustering step or to identify patients having the same mutational profiles.

Drug response in association with genetic alterations is one of the applications of compaRe. The genetic alteration can be visualized in the clusters that compaRe identifies. However, at the moment, we do not have relevant data to apply this type of analysis for visualization. To highlight that the functionality for this type of analyses exists within compRe, we have now explained this in the main text lines 212-213.

Reviewer #3:

[…] A major weakness of the current presentation of the study is that it has not clearly demonstrated the toolkit's utility in exploring specific phenotypes in-depth within a high-parameter dataset. The following are two examples in which this limitation is relevant, and the authors may address this to strengthen this paper, if in fact detailed phenotyping is considered by the authors as an important feature of the toolkit. If not, the authors should revise the manuscript as such.

First, the authors stated that their approach can be used to "optimize true cytometric n-dimensional immunophenotypic characterization" even in the setting of multi-panel worklow. However, the demonstration was based on samples that seem to have predominant phenotypes that are almost mutually exclusive. It is unlikely that this toolkit would be useful for reliable phenotypic characterization in a largely heterogeneous population of cell types, e.g. even in normal peripheral blood, unless a high number of parameters was concurrently acquired within the dataset. This is an inherent limitation. Nonetheless, a revision to demonstrate how compaRe can evaluate specific clusters of phenotypes with biological significance from a high-parameter dataset (20-30 marker cytometry) would be very helpful.

Thank you for your comments. To address this comment, we generated Appendix 1-figure 32 and added the result to the main text lines 124-131. As compaRe measures the similarity between cluster morphologies, we carried out an analysis in which cells from a cluster (an immunophenotypic cell population) were gradually removed to contort its configuration. We used a dataset of 3 healthy and 2 pediatric AML bone marrow mononuclear cell samples from the data provided in the 6th reference. Samples were stained with 29 (15 membrane and 14 intracellular signaling) markers. Taking H1 as reference, we gradually removed 25%, 50%, 75% and 100% (phenotypic changes) of cells from a target cluster identified by PhenoGraph. As the UMAP projections show, the similarity decreased concurrently and more drastically when phenotypic changes were detected.

Second, the authors refer to the method's ability to include all rare cell subsets in the analysis, i.e. the ability to forego any subsampling. The work does not, however, demonstrate clearly how the presence of a rare cell subset in a given sample influences its similarity to other samples. Thus, the toolkit's value of being able to include such a rare subset in the analysis remains unestablished. It would be beneficial to include such a test to see whether the algorithm is sensitive to such changes.

We previously discussed in the manuscript that ‘the poor performance of JSD indicates that this approach can work well only in the absence of signal shift. It is of particular note that compaRe does not need subsampling or dimension reduction of the input data, avoiding the risk of losing rare subpopulations.’ The use of full data rather than subsampled data allows for more robust analyses, in part, simply, owing to using more data and the ability to establish a more refined overall cell cluster morphology of each sample.

Additionally, in the analysis explained above (Appendix 1-figure 32), the 75% sample potentially resembles a rare cell population, a cell cluster with a small number of cells. The analysis demonstrated compaRe was sensitive enough to maintain a certain similarity in the presence of a rare cell population, as the similarity score was reduced when this population was gone in the 100% sample. Nonetheless, since subsampling does not necessarily lose rare subpopulations and to avoid confusion, we removed the phrase ‘avoiding the risk of losing rare subpopulations’ from the discussion.

https://doi.org/10.7554/eLife.73760.sa2

Article and author information

Author details

  1. Morteza Chalabi Hajkarim

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Formal analysis, Methodology, Software, Validation, Visualization, Writing – original draft, Writing – review and editing, conceived, designed, and wrote the study with equal contribution
    Contributed equally with
    Ella Karjalainen
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2039-2676
  2. Ella Karjalainen

    Institute for Molecular Medicine Finland (FIMM), Helsinki Institute of Life Science, University of Helsinki, Helsinki, Finland
    Contribution
    Conceptualization, Formal analysis, Methodology, Validation, Writing – original draft, Writing – review and editing, conceived, designed, and wrote the study with equal contribution
    Contributed equally with
    Morteza Chalabi Hajkarim
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9865-5384
  3. Mikhail Osipovitch

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Software, developed the GUI
    Competing interests
    No competing interests declared
  4. Konstantinos Dimopoulos

    Rigshospitalet, Copenhagen, Denmark
    Contribution
    Data curation, Validation, assembled and annotated clinical cytometry data and assisted in its analysis
    Competing interests
    No competing interests declared
  5. Sandra L Gordon

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Formal analysis, Validation, Writing – original draft, designed and completed the AML patient sample drug screening, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0270-8291
  6. Francesca Ambri

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Data curation, Formal analysis, Validation, designed and completed the AML patient sample drug screening
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1999-9294
  7. Kasper Dindler Rasmussen

    Centre for Gene Regulation and Expression, School of Life Sciences, University of Dundee, Dundee, United Kingdom
    Contribution
    Methodology, generated the mouse cell models, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7344-4177
  8. Kirsten Grønbæk

    1. Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    2. Rigshospitalet, Copenhagen, Denmark
    Contribution
    Funding acquisition, Supervision, assembled and annotated clinical cytometry data and assisted in its analysis, Writing – review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1535-9601
  9. Kristian Helin

    1. Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    2. Cell Biology Program and Center for Epigenetics Research, Memorial Sloan Kettering Cancer Center (MSKCC), New York, United States
    Contribution
    Funding acquisition, Supervision, generated the mouse cell models
    Competing interests
    No competing interests declared
  10. Krister Wennerberg

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing – original draft, conceived, designed, and wrote the study with equal contribution, Writing – review and editing
    For correspondence
    krister.wennerberg@bric.ku.dk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1352-4220
  11. Kyoung-Jae Won

    Biotech Research and Innovation Centre (BRIC) and Novo Nordisk Foundation Center for Stem Cell Biology (DanStem), University of Copenhagen, Copenhagen, Denmark
    Contribution
    Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Supervision, Validation, Writing – original draft, conceived, designed, and wrote the study with equal contribution, Writing – review and editing
    For correspondence
    kyoung.won@bric.ku.dk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2924-9630

Funding

Novo Nordisk Foundation (NNF17CC0027852)

  • Kirsten Grønbæk
  • Kristian Helin
  • Krister Wennerberg
  • Kyoung-Jae Won

Kræftens Bekæmpelse (R223‐A13071)

  • Kirsten Grønbæk
  • Kristian Helin
  • Krister Wennerberg
  • Kyoung-Jae Won

Lundbeck Foundation (R313-2019-421)

  • Kyoung-Jae Won

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported through Novo Nordisk Foundation (Novo Nordisk Foundation Center for Stem Cell Biology, DanStem; Grant Number NNF17CC0027852) and Danish Research Center for Precision Medicine in Blood Cancers funded by the Danish Cancer Society (Grant number R223‐A13071) and Greater Copenhagen Health Science Partners.

Ethics

Human subjects: The informed consent, and consent to publish of patient samples in this study has been approved by the Danish National Science Ethics Committee/National Videnskabsetisk Komite: Målrettet behandling af patienter med blodsygdomme, license no. 1705391.

Senior and Reviewing Editor

  1. Aleksandra M Walczak, CNRS LPENS, France

Reviewer

  1. Aik Choon Tan, Moffitt Cancer Center, United States

Publication history

  1. Preprint posted: January 9, 2021 (view preprint)
  2. Received: September 9, 2021
  3. Accepted: February 14, 2022
  4. Accepted Manuscript published: February 15, 2022 (version 1)
  5. Version of Record published: April 20, 2022 (version 2)
  6. Version of Record updated: May 3, 2022 (version 3)

Copyright

© 2022, Chalabi Hajkarim et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 648
    Page views
  • 115
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Morteza Chalabi Hajkarim
  2. Ella Karjalainen
  3. Mikhail Osipovitch
  4. Konstantinos Dimopoulos
  5. Sandra L Gordon
  6. Francesca Ambri
  7. Kasper Dindler Rasmussen
  8. Kirsten Grønbæk
  9. Kristian Helin
  10. Krister Wennerberg
  11. Kyoung-Jae Won
(2022)
Comprehensive and unbiased multiparameter high-throughput screening by compaRe finds effective and subtle drug responses in AML models
eLife 11:e73760.
https://doi.org/10.7554/eLife.73760

Further reading

    1. Computational and Systems Biology
    2. Medicine
    Riku Klén et al.
    Research Article

    New SARS-CoV-2 variants, breakthrough infections, waning immunity, and sub-optimal vaccination rates account for surges of hospitalizations and deaths. There is an urgent need for clinically valuable and generalizable triage tools assisting the allocation of hospital resources, particularly in resource-limited countries. We developed and validate CODOP, a machine learning-based tool for predicting the clinical outcome of hospitalized COVID-19 patients. CODOP was trained, tested and validated with six cohorts encompassing 29223 COVID-19 patients from more than 150 hospitals in Spain, the USA and Latin America during 2020-22. CODOP uses 12 clinical parameters commonly measured at hospital admission for reaching high discriminative ability up to nine days before clinical resolution (AUROC: 0·90-0·96), it is well calibrated, and it enables an effective dynamic risk stratification during hospitalization. Furthermore, CODOP maintains its predictive ability independently of the virus variant and the vaccination status. To reckon with the fluctuating pressure levels in hospitals during the pandemic, we offer two online CODOP calculators, suited for undertriage or overtriage scenarios, validated with a cohort of patients from 42 hospitals in three Latin American countries (78-100% sensitivity and 89-97% specificity). The performance of CODOP in heterogeneous and geographically disperse patient cohorts and the easiness of use strongly suggest its clinical utility, particularly in resource-limited countries.

    1. Cancer Biology
    2. Computational and Systems Biology
    Gökçe Senger et al.
    Research Article

    Aneuploidy, a state of chromosome imbalance, is a hallmark of human tumors, but its role in cancer still remains to be fully elucidated. To understand the consequences of whole-chromosome-level aneuploidies on the proteome, we integrated aneuploidy, transcriptomic and proteomic data from hundreds of TCGA/CPTAC tumor samples. We found a surprisingly large number of expression changes happened on other, non-aneuploid chromosomes. Moreover, we identified an association between those changes and co-complex members of proteins from aneuploid chromosomes. This co-abundance association is tightly regulated for aggregation-prone aneuploid proteins and those involved in a smaller number of complexes. On the other hand, we observe that complexes of the cellular core machinery are under functional selection to maintain their stoichiometric balance in aneuploid tumors. Ultimately, we provide evidence that those compensatory and functional maintenance mechanisms are established through post-translational control and that the degree of success of a tumor to deal with aneuploidy-induced stoichiometric imbalance impacts the activation of cellular protein degradation programs and patient survival.