Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Methods. Author manuscript; available in PMC 2016 July 1.
Published in final edited form as:
PMCID: PMC4468028

Methods for discovery and characterization of cell subsets in high dimensional mass cytometry data


The flood of high-dimensional data resulting from mass cytometry experiments that measure more than 40 features of individual cells has stimulated creation of new single cell computational biology tools. These tools draw on advances in the field of machine learning to capture multi-parametric relationships and reveal cells that are easily overlooked in traditional analysis. Here, we introduce a workflow for high dimensional mass cytometry data that emphasizes unsupervised approaches and visualizes data in both single cell and population level views. This workflow includes three central components that are common across mass cytometry analysis approaches: 1) distinguishing initial populations, 2) revealing cell subsets, and 3) characterizing subset features. In the implementation described here, viSNE, SPADE, and heatmaps were used sequentially to comprehensively characterize and compare healthy and malignant human tissue samples. The use of multiple methods helps provide a comprehensive view of results, and the largely unsupervised workflow facilitates automation and helps researchers avoid missing cell populations with unusual or unexpected phenotypes. Together, these methods develop a framework for future machine learning of cell identity.

1. Introduction

1.1. High dimensional single cell biology

Single cell biology is transforming our understanding of the biological mechanisms driving human diseases and healthy tissue development [1]. Mass cytometry is a recently developed technology that enables simultaneous detection of more than 40 features on individual cells [2, 3]. High dimensional mass cytometry measurements are single cell, quantitative, and well-suited to unsupervised computational analysis. New analysis tools have been created to take advantage of the massive amounts of data that result from high content single cell techniques like mass cytometry. Variations of many of these tools have been developed and applied for gene expression analysis, a field facing similar problems with data dimensionality. These tools draw on advances in machine learning and statistics that are not yet widely applied in biological studies. Many of these tools are complementary and address different aspects of data analysis, and it can be challenging for biologists to know when and how to use these tools to get the most out of their data. Advances have also been made in automating and standardizing the flow cytometry data analysis workflow [46]. Here, we present a modular workflow focused on high dimensional single cell analysis that combines multiple tools to provide a comprehensive view of both cells and populations. Rather than making the workflow fully automated, the goal here was to combine the complementary benefits of expert analysis and machine learning. This approach maintains single cell views, provides automatic population assignment for each cell, and facilitates statistical comparison of the key cellular features that characterized each population. This semi-supervised workflow facilitates comparison of populations discovered by different computational approaches, in different clinical samples, or using different biological features (e.g. RNA expression, cell surface protein expression, and cell signaling).

An advantage of traditional analysis in flow cytometry is the reliance on identification of known, prominent populations with strong supporting biology in the literature. Given the typical panel size for fluorescent experiments, this type of supervised analysis is fast and usually adequate. Unfortunately, expert manual gating has been shown to be particularly prone to inter-operator variability [7] and a tendency to overlook cell populations [810]. Recent efforts have developed new tools for high dimensional cytometry data that bring in elements of machine learning and statistical analysis, including clustering [1114], dimensionality reduction [8], variance maximization [15], mixture modeling [6, 1618], spectral clustering [19], neural networks [20], and density-based automated gating [21]. Here, we highlight use of these tools in a sequential single cell bioinformatics workflow (Table 1). In particular, different tools address aspects of data visualization, dimensionality reduction, population discovery, and feature comparison. It can be valuable to apply multiple tools in order to view data in different ways and fully extract biological meaning at the single cell level (Figure 1) and the population level (Figure 2 and Figure 3). After identifying cell subsets with the aid of computational tools, measured features, such as protein expression in the examples here, can be compared between and within the subsets. Traditional statistics used include medians, variance, and fold changes. Other statistical methods such as histogram statistics and probability binning have also been used to compare distributions in flow cytometry data [2224].

Fig. 1
Distinguishing initial populations with viSNE analysis of per-cell protein expression and expert gating
Figure 2
Revealing cell subsets with SPADE analysis of population hierarchy, cell abundance, and median protein expression
Figure 3
Characterizing cell subsets with a heatmap analysis of median protein expression and hierarchical clustering of proteins and populations
Table 1
A modular machine learning workflow for semi-supervised high-dimensional single cell data analysis

1.2. Overview of the analysis workflow

The workflow presented here was applied to a CyTOF dataset from the analysis of healthy human bone marrow and a diagnostic sample of blood from a patient with acute myeloid leukemia (AML). The annotated FCS files and a step-by-step guide are available online from Cytobank ( [25] and FlowRepository ( [26]. This workflow was developed for use with high-dimensional mass cytometry data. However, it can also be applied to fluorescent flow cytometry data. The main steps presented consist of event restriction, population discovery, and population characterization. Each of these aspects of data analysis can be achieved with a variety of techniques (Table 1), and some tools address multiple steps. By sequentially combining three different techniques, this workflow draws on the strengths of specific tools, keeps biologists in touch with single cell views, and enables analysis of data from different studies and single cell platforms.

In the case of the example dataset here, the overall biological goal was to identify and compare three populations of cells: leukemia cells (AML blasts) and non-malignant cells (non-blasts) in the blood of a leukemia patient, and bone marrow cells from a healthy donor. In the analysis workflow, cell events were first manually gated based on event length and DNA content to include intact, single cells (Figure 1) [11]. Next, visualization of stochastic neighbor embedding (viSNE) was used to identify and gate major subsets (Figure 1). Gated cells from healthy bone marrow and AML were then analyzed by spanning-tree progression analysis of density-normalized events (SPADE) to discover and compare cell subsets (Figure 2). Finally, the cell subsets identified by SPADE were further characterized using complete linkage hierarchical clustering and a heatmap in R (Figure 3). The details of mass cytometry data collection and processing prior to initial cell selection (gating) are not covered in detail here. These early steps include experiment design, collection of data at the instrument (and instrument setup), any normalization, and transformation of the data to an appropriate scale (Table 1).

The initial event restriction step that begins the workflow focuses the analysis on populations of cells. The goal at this step is to remove events that do not contribute useful information while making minimal changes to the data and not over-focusing. Event restriction is traditionally performed using biaxial gating (Table 1), but given the high dimensionality of mass cytometry data, use of viSNE (Figure 1) can simplify the process of distinguishing initial populations and avoid overlooking cells with unusual or unanticipated phenotypes. The second step, cell subset identification, is also traditionally performed by expert gating (Table 1). However, clustering tools such as SPADE [12] (Figure 2), Misty Mountain [13], and Citrus [14], among others, can be used to automatically assign cells to groups or clusters in high dimensional data. In the workflow here, the goal is to find all the phenotypic clusters of cells in healthy bone marrow, AML blasts, and non-blast cells from AML blood (Figure 2). As the final step, characterization of discovered cell subsets takes place downstream of manual gating or automated discovery tool implementation, and generally consists of feature expression comparison with heatmaps, violin plots, and histogram overlays for visualization, as well as data modeling and other statistical analysis. This workflow emphasizes integration of automated, unsupervised approaches with minimal human gating and processing. This type of semi-supervised cell population discovery and characterization can decrease human bias and variability and identify phenotypically unusual or rare cell subpopulations.

1.3. Advantages of machine learning tools: dimensionality reduction, clustering, and modeling

Not all tools perform the same analysis functions. Three functions that are useful for high-content single cell analysis include dimensionality reduction, clustering of cells into populations, and modeling. SPADE and viSNE both include dimensionality reduction steps that project multi-dimensional data into a lower dimensional space for visualization and further interpretation. These algorithms aim to preserve key high-dimensional phenotypic relationships between cells when visualizing and comparing them in 2D space. Depending on the structure of the data, other dimensionality reduction tools might be used (Table 1). Locally linear embedding (LLE) and isometric mapping (ISOMAP) are designed for the types of continuous phenotypic distributions seen in developmental progressions. ISOMAP accounts for geodesic distance in addition to local linear distances between high dimensional data points in order to reduce the dimensions of continuous and non-linear data [27, 28]. A similar principle is applied with LLE, where locally linear embedding of similar data points in high dimensional space is preserved while allowing for a non-linear global embedding of the data during projection into low dimensional space [29]. In contrast, multidimensional scaling (MDS) and principal component analysis (PCA) preserve linear, multi-dimensional variance. One of the advantages of PCA and other techniques, such as joint clustering and modeling [30], is the creation of a model that can be applied to newly analyzed samples. In addition to the unsupervised tools discussed here, population analysis techniques that include some supervision can be particularly useful for mapping features across known developmental progressions [31, 32].

Notably, dimensionality reduction alone does not assign cells to groups. Here, dimensionality reduction with viSNE is used to aid expert interpretation of cluster identity. In this example, cells are projected onto a biaxial plot space by viSNE and then gated. Thus, viSNE is being used to see the phenotypic relationships of the cells according to all 27 protein features. This can help researchers visualize high dimensional data without losing rare populations that are best observed in single cell views. Following t-SNE or viSNE analysis, a human expert can look for cell clusters or major populations, as is the case here (Figure 1), or a computational tool can identify cell clusters (Table 1), as with t-SNE + DensVM analysis [33]. As the workflow becomes increasingly unsupervised, it is especially important to include a single cell view early in the analysis so that expert can perform quality checks and get a sense of the overall biological results.

2. Data collection, processing, and initial population identification

2.1. Data collection

In mass cytometry, as with fluorescent flow cytometry, single cell suspensions are stained with metal-conjugated antibodies specific to molecules of interest. At the mass cytometer, cells are aerosolized and streamed single-file into argon plasma where they are atomized and ionized. The resulting ion cloud passes through a quadrupole to exclude low mass ions and enrich for reporter ions whose abundance is proportional to cellular features. These reporter ions are quantified by time of flight mass spectrometry [34] and recorded in an IMD format file. These data are typically parsed into single cell events and converted to a flow cytometry standard (FCS) file for analysis [35]. Many software programs can handle FCS files, including Cytobank (, FlowJo (, R/Bioconductor, MATLAB, Cytoscape, and GenePattern ( [36]. Text files containing the expression matrix (where rows are cells and features are columns, and there is a median intensity value for each cell) can also be extracted directly from the IMD file from the cytometer or from the FCS file for manual analysis outside of flow cytometry analysis software. In Cytobank, export of text files with intensity values is available from the FCS file details page. An expression matrix can also be extracted from the FCS file in R and MATLAB using FCS file parsing functions. In R, the package “flowCore” can be used to extract the intensity values from the FCS file using the exprs() function [37]. In MATLAB, the tool “FCS data reader” includes the function fca_readfcs() to extract the intensity values of FCS files [38].

Here, the healthy human bone marrow sample analyzed was obtained as a de-identified sample left over from diagnostic analysis of non-cancerous tissue in the Vanderbilt Immunopathology core. Acute myeloid leukemia peripheral blood samples were collected from consented patients. In all cases, samples are collected in accordance with the Declaration of Helsinki following protocols approved by Vanderbilt University Medical Center (VUMC) Institutional Review Board. The patient blood sample evaluated here was collected at the time of diagnosis following initial evaluation and prior to any treatment.

2.2. Data processing and scale transformation

In order to prepare data for dimensionality reduction and analysis, initial processing steps aim to ensure the quality of cell events and perform appropriate scale transformations. Quality control varies by user and is especially important when conducting studies across time or using data from different instruments. Data normalization using internal bead controls can be applied as part of this data processing [39]. In this case, the two samples were collected sequentially on the same instrument and no signal normalization was required. Efforts are underway to facilitate comparison of data among groups and centers and to report elements of panel design, instrument settings, data processing, and normalization. MIFlowCyt is a data standard set by International Society for Advancement of Cytometry (ISAC) that specifies the minimum amount of information that must be included in an FCS file to ensure reproducibility and transparency [40]. ISAC has also established a file format for classification results from flow cytometry data (CLR) [41] that handles cell classification from manual or automated identification and compliments the Gating-ML file format that was developed for sharing biaxial gate classifications [42]. Additionally, there have been efforts to standardize and compare computational flow cytometry data analysis tools. The FlowCAP project compares automated tools for cytometry data analysis using standardized datasets [5]. EuroFlow is a consortium of research groups that optimize flow cytometry protocols and analysis methods and set standards for the field of immunology and hematological studies [43]. Reporting of optimized antibody panels has also been standardized in the form of Optimized Multicolor Immunofluorescence Panels (OMIPs) [44]. Cytobank ( and FlowRepository ( provide online access to annotated cytometry data files, including mass cytometry datasets [25, 26].

Because cytometry data are log-normal, a log or log-like scale is typically used to visualize and interpret the data. Commonly used scales include inverse hyperbolic sine (arcsinh), logarithmic, and logicle (also referred to as “bi-exponential”) scales [45]. Logicle or log-like scales more accurately represent the spread of data around 0 than logarithmic scales, given that modern cytometers can produce negative and zero values that cannot be transformed using logarithmic scales. The implementation of the arcsinh scale here was first used for fluorescent flow cytometry [46] and is now standard for mass cytometry. Typically, a cofactor is included as part of the arcsinh scale transformation as a way of setting a channel specific minimum significance threshold. The cofactor is set to an intensity value below which differences are not significant. For mass cytometry, cofactors typically range from 3 to 15 and depend on background and signal to noise with the detection channel and antibody-metal conjugate. In fluorescent flow cytometry, cofactors generally range from 25 to 2000 and are especially useful in correcting for channel specific differences in spreading of negative events that depend on fluorophore selection, compensation, and instrument setup. For fluorescent flow cytometry data, appropriate compensation must also be applied prior to analysis in order to correct for any spillover between channels. Algorithms have been developed for fluorescent cytometry to automatically determine scale transformations [47, 48]. Applying an appropriate scale transformation prior to computational analysis is critical because it impacts quantification of distance between cells in the same way that it affects visualization of distance in biaxial plots.

2.3. Initial population identification and quality assessment

Beginning data analysis with a single cell view reveals the quality of the data and allows experts to spot rare cell subsets or artifacts that can be obscured in aggregate analysis. It is valuable to review the single cell data to verify computational analysis results, and it is vital in publications to provide representative single cell views of findings. Here, intact single cells were gated by human analysis of event length and iridium intercalator uptake (Figure 1). This initial gating might be accomplished various ways, such as use of cisplatin exclusion to identify live cells [49]. Event length is generally higher for the mass cytometry equivalent of ‘cell doublets’ that can occur when the signal from two cells is not well separated in time. Intercalator uptake helps mark all cells and is proportional to nucleic acid content [11, 34]. A biaxial view of each channel was then used to evaluate data quality prior to computational analysis. If no intercalator positive events are seen in this view, it suggests that there were no cells in the sample or there was an error in DNA intercalator staining. Once intact, single cells have been identified (Figure 1A), a quick check using traditional biaxial plots or histograms can be used to ensure there is no clear overstaining. Severe overstaining results in errors while collecting data on the cytometer because event length is too great and individual cell events cannot be distinguished. Additionally, checks could be included at this step for contaminant signals. Atomic mass contaminants, such as barium and lead, can be found in water, buffers or glassware. Collecting data for the corresponding channels (137 and 138 for Ba, 208 and 209 for Pb) can be used to track these contaminants. In summary, intact single cells are first gated by a human expert. This step may be automated, but it represents an opportunity for quality assessment and initial familiarization with the data prior to computational analysis.

3. Unsupervised machine learning tools

3.1. viSNE

viSNE is a cytometry analysis tool that employs t-stochastic neighbor embedding (t-SNE) in mapping individual cells in a two or three-dimensional map that is based on their high dimensional relationships [8, 50]. viSNE can be used to provide a human readable two-dimensional (2D) view of cells that are arranged in a way that approximates high-dimensional phenotypic similarity. viSNE is implemented in MATLAB and Cytobank [25], and the Cytobank implementation of viSNE is shown here (Figure 1). viSNE can be run using a single population of cells or multiple populations drawn from one or more files. However, cell features selected for analysis must have been measured on all cell populations in a comparable way and features must be measured on comparable scales. It is sometimes helpful to subsample cell events from populations to speed the analysis or test robustness. Sampling can be ‘equal’ with respect to the starting populations, in order to ensure that each cell population is represented on the viSNE map by the same number of cells, or ‘proportional’, so that each population is represented by a number of cells proportional to its abundance. When data are thought to contain rare cell subsets, subsampling should be avoided to preserve rare cells. Initial gating can be used to focus the analysis on a population of interest and increase its relative abundance. Here, equal numbers of cells were selected from the AML PBMC and healthy marrow files for the viSNE analysis.

The cell features selected for viSNE mapping affect the structure of the viSNE map. Markers that vary highly between cell subsets will polarize subsets, placing them farther apart in tightly grouped islands. Markers with low variance on subsets will cause those cells to be placed closer together on the map. Thus, including markers that are not expressed on any cells can result in compression of islands on the map and loss of subset polarization. Features that might contribute to clustering can be selected in an unsupervised manner based on variance. For example, features that vary more in disease than in healthy controls might be particularly useful in stratifying cells associated with distinct molecular subgroups [51]. Here, all 27 markers in the panel were included in viSNE mapping because all were expressed and variable on the cells in the samples. The displayed viSNE map shows cells from the AML patient file only (Figure 1). The resulting viSNE map showed a broad distribution of heterogeneous CD45lo AML cells and several distinct islands of non-blast cells (Figure 1B). Relative protein expression as heat intensity can be viewed for each marker in the panel and are shown here for the 27 markers on the panel (Figure 1C). The two main populations of AML blast and non-blast cells were then manually gated from the viSNE map and exported as separate FCS files for further comparison to healthy bone marrow cells using SPADE and heatmaps. All healthy marrow cells were exported from the viSNE analysis as no additional gating was required to identify major populations. Depending on the sample and biological question, it may be useful to gate initial major populations from several or all files in this step of the analysis.

The MATLAB implementation of viSNE can be accessed through the freely downloadable cyt tool ( Cyt employs a user interface that allows for selection of features for mapping, selection of files or gates to be mapped, an interface for visualizing parameter intensity on a heat scale, and a tool within the interface for manually gating populations resulting viSNE map.


SPADE is an algorithm that includes dimensionality reduction, clustering of cells into populations (also referred to as ‘nodes’), and visualization using a 2D minimum spanning tree. Data must be appropriately scaled and intact cell events gated prior to SPADE analysis as described above. Here, this is done prior to viSNE gating. SPADE has been implemented in Cytobank, R, Cytoscape (, and MATLAB. In R, the package “spade” includes functions to implement individual steps of SPADE and to execute a comprehensive SPADE analysis [52]. CytoSPADE is a plugin available for use in Cytoscape that provides a user interface with the R implementation ( The MATLAB implementation of SPADE requires the SPADE V2.0 MATLAB tool that is freely downloadable ( Here, the Cytobank implementation of SPADE was used to compare populations identified in viSNE guided gating. User-defined parameters for SPADE analysis include downsampling, feature selection, and a target number of nodes. Target downsampling, which can be indicated as either a percentage of cells or an absolute number, specifies how much weight to give clusters of varying density. A lower downsampling percentage increases the likelihood that sparse regions of density will be given their own clusters rather than being grouped into clusters with regions of higher density. When a sample is thought to contain rare subsets of cells, entering a lower downsampling value can help distinguish these cells as a separate population [11, 12]. Feature selection in SPADE can also be based on selecting highly variable or biologically relevant markers, as described above for viSNE. The number of nodes indicates the target number of clusters (i.e. cell subsets) that the algorithm should produce, and 200 nodes is a good default for standard mass cytometry datasets containing ~105 to 107 total intact single cells. Including more nodes in the analysis helps to assign rare subsets to their own clusters. These clusters can be easily combined in a process called “bubbling”, in which a human expert manually refines the cluster identity. A table of basic statistics, such as median intensity of each feature, is generated for each population of cells identified by SPADE and can be downloaded as a text file. Cell subsets identified by SPADE can additionally be exported as individual FCS files for further analysis, as in the heatmap analysis shown here (Figure 3).

In the example here, three populations were analyzed. The two populations of AML blast and non-blast cells identified by viSNE (Figure 1) were compared with the population of healthy bone marrow cells stained with the same mass cytometry antibody panel. Here, a concatenated file containing all three populations was also included to allow visualization of all cells simultaneously on one tree (Figure 2C). SPADE can initialize with a fixed or random seed and is random in the Cytobank implementation. The same random seed can be set from run to run in the MATLAB implementation. However, when new files are added to the analysis, a different tree can still stem from the same seed, which necessitates re-running the analysis to include any additional files. For this analysis, the downsampling percentage was set to 1%, the target number of nodes was 100, and the features used for clustering were all 27 measured markers in the panel. The resulting SPADE trees are shown in Figure 2.

Including SPADE in this analysis workflow has several advantages. First, SPADE produces a visualization of population abundances by altering the sizes of each node depending on how many cells it encompasses. For example, it can be seen in the SPADE tree that the non-blast AML cells fall almost exclusively into one node, reflecting their relative homogeneity and the lack of normal immune cell populations in the AML patient’s blood (Figure 2). Clustering with SPADE also assigns each cell to a discrete group, which minimizes analysis variability and prevents loss of cells that are outside of gated regions in manual biaxial gating. In a standard SPADE analysis, the algorithm is asked to “over cluster”, producing hundreds of relatively small clusters rather than grouping cells into fewer, larger groups. This over clustering gives high resolution to improve rare subset identification and allows for a thorough annotation and characterization of all potentially discrete biological populations in the heatmap analysis.

4. Characterizing and visualizing populations

4.1. Population heatmaps

With some algorithms it is not straightforward to compare the results of an analysis of one set of samples with the results from another set of samples. For example, with SPADE it is not straightforward to map a new sample onto an existing minimum spanning tree defined using different samples. Instead, a new SPADE analysis is generally run that includes both the new and old samples. In contrast, a heatmap can be used to compare populations identified in different analysis runs of SPADE or populations identified by different clustering techniques. Heatmaps also provide a compact view that facilitates comparing many populations according to a large variety of measured features. In heatmaps, different types of biological and clinical information can also be used to group populations or assessed for association with resulting groups [46, 53]. While population heatmaps provide an intuitive, high-level view of the results, they can obscure variation within subsets [25]. To address this, statistics other than median expression can be shown in the heatmap, such as variance or the 95th percentile of expression [1, 54].

For the last step in the workflow here, tables of statistics for the hundreds of cell subsets identified in the three starting populations (Figure 2) were exported from SPADE as text files listing median expression of each feature for each cell subset. Cell subsets were excluded from further analysis if they contained less than 1% of the cells in the starting population. This arbitrary threshold was set in order to exclude sparse clusters where low cell number could potentially increase the error of reported medians. Here, the 1% threshold resulted in exclusion of approximately 5% of the total cell events from heatmap characterization. The table of statistics was then imported into R using the “read.table” function from the R.utils package [55] and visualized as a hierarchically clustered heatmap using the “heatmap.2” function in the gplots package (Figure 3A) [56]. Output of a hierarchical dendrogram as part of the heatmap can be specified as one of the input parameters of the heatmap.2 function. The R package “stats” also offers a function called “heatmap” that performs the same function as heatmap.2 with slight differences in visualization options. After the clustered heatmap was generated, expert analysis was used to assign biological classifications to each group of populations in the hierarchical clustering, and included the same populations seen in viSNE (Figure 1B) and SPADE (Figure 2B): dendritic cells (DCs), monocytes, natural killer cells (NKs), CD8+ T cells, CD4+ T cells, B cells, immature myeloid cells (Imm. myel.), four subsets of AML blast cells (AML1 through AML4), and erythroid blast cells (Ery. bl.) (Figure 3B).

Use of a clustered heatmap in the workflow allows for simultaneous visualization of several markers for the same clusters (population of cells) from multiple files. Furthermore, nodes are hierarchically clustered, and this clustering can be pruned at various levels by the user to further group the nodes into biological populations. It is also important to note that the distance between nodes has quantitative meaning in the clustered heatmap dendrogram, as opposed to the distances on the SPADE tree that are for visualization purposes and not quantitative. Heatmap analysis therefore compliments the SPADE visualization by facilitating simultaneous visualization of nodes from multiple files and by quantifying phenotypic distances between the nodes.

4.2. Other packages and flowCore

There are many R packages designed for statistical and visual analysis of flow cytometry data, including flowCore [37], flowViz [57], flowStats [58], and flowClust [58], among others. These packages include functions for producing heat maps, histograms, bar plots, biaxial density plots, and are part of efforts to automate and standardize computational analysis of cytometry data [5, 6]. Apart from the R packages designed for flow cytometry data analysis, other analysis and visualization packages can be applied to single cell data. For example, box and whisker plots or violin plots can be produced to show median, range, and the distribution of the feature in each subset.

5.0. Other Considerations for automated flow cytometry data analysis

5.1. Algorithm Selection

Three major considerations when choosing tools or algorithms for flow cytometry data analysis include 1) linear vs. non-linear measurement, 2) supervised or unsupervised approaches, and 3) need for modeling. The first consideration is whether a linear or non-linear method of dimensionality reduction is best for the data. Phenotypic relationships between cells may follow a ‘creode’, or necessary path, that is non-linear with respect to protein expression (i.e. co-expression or co-variance of molecules is not linearly correlated with important progressions in cellular identity or trajectories in data space) [10]. In this case, nonlinear dimensionality reduction tools may better preserve the high dimensional phenotypic relationships between cells compared to tools that assume a linear relationship between variables. The second consideration is whether an unsupervised or supervised method is needed. In an exploratory analysis where novel populations are anticipated, unsupervised approaches will minimizes the risk of overlooking the populations. Lastly, a consideration is whether or not the goal of analysis is to build a model. Mixture modeling tools can be implemented for analysis of flow cytometry data that will produce a model as output for downstream analysis. Additional issues to consider include 1) selection of features, which is generally initiated by hypotheses and pragmatic concerns and then narrowed to include those features with biologically meaningful variation [51], and 2) aspects of statistical power, including sample size, cluster density, and false discovery rate (FDR). It is vital to calculate FDR or a related statistic, such as the f-measure, in cases where a truth is known [5].

5.2. Scalability of Workflow

Biomedical studies that employ flow and mass cytometry often accrue large numbers of samples over long periods of time. This and similar workflows can be adapted to accommodate data from these large studies. In order to account for experimental or instrument variability, normalization is necessary in these cases in order for compare samples run at different times or from different instruments. Bead normalization has been optimized for use with mass cytometry to control for machine variability between runs [39, 59, 60]. Polystyrene beads embedded with heavy metal isotopes are run with every sample as a standard that can be used to correct MI values for each event based on technical variability. When samples accrue over a long period of time, a key consideration is that new results may not be easily mapped back to the original viSNE map or SPADE tree without re-analysis. This is one advantage of heatmaps, which compare samples according to a simple ‘model’ of the data, such as median expression of selected features.

This workflow as presented includes manual intervention that could be prohibitive when analyzing many data files simultaneously. While all steps of this analysis could generally be batched and automated, human review of single cell data is advantageous at workflow breakpoints to verify computational results and spot artifacts. Cytobank and other flow cytometry data analysis software allow for rapid, simultaneous viewing and pre-processing of multiple files, including scale transformation and gating. viSNE analysis can currently be run on up to 800,000 cells in Cytobank, and this limit is pragmatic, not theoretical. Many files can be run simultaneously by subsampling cells equally or proportionally from each file prior to the viSNE run. SPADE can also be run on many files simultaneously, and data files with cluster information can be quickly downloaded in a compressed folder.

Import of text files into R and selection of nodes based on the number of cells they contain can be automated and batched for highly scalable and rapid heatmap analysis. However, a potential limitation of large-scale analyses is the visualization of all nodes simultaneously on the heatmap. It may be useful in these cases to segment the SPADE tree into major populations by “bubbling” and then building separate heatmaps from each bubble rather than for the whole tree. Depending on the expected prevalence of rare cells in the dataset, the user can request fewer nodes in the SPADE run in order to decrease the final number of clusters to be analyzed and visualized on the heatmap.

5. Conclusions

Data analysis in cytometry remains largely manual, supervised, and focused on large changes in magnitude of expression. As new tools are developed to assist in gating, reduce dimensionality, and automate analysis, it is important to show biologists the value of these tools and to integrate them into workflows that can become routine. The workflow presented here blends supervised and unsupervised analysis tools so that biologists can visualize results at the single cell level while still getting an accurate view of the big picture. Combining tools also allows the analyst to visualize data in multiple ways, which can be useful to extract the most meaning from a data set. Existing tools allow for identification of populations based on single cell expression profiles and characterization of these subsets using standard statistics, including expression magnitude, marker variance, and subset abundance. Going forward, tools that quantify cellular heterogeneity, identify critical population features, and assign biological identity to machine-identified subsets will be particularly useful in filling out the toolkit.


  • Bioinformatics strategy for analysis of high-dimensional single cell data
  • Unsupervised approaches reveal and characterize cells with unexpected phenotypes
  • Modular approach facilitates development and testing of new tools
  • Sequential use of machine learning tools combines complementary strengths of each

Supplementary Material


This study was supported by R25 CA136440-04 (K.E.D.), NIH/NCI K12 CA090625 (P.B.F), R00 CA143231-03 (J.M.I.), the Vanderbilt-Ingram Cancer Center (VICC, P30 CA68485), and VICC Young Ambassadors and VICC Hematology Helping Hands awards. Thanks to Mikael Roussel for helpful discussions of myeloid cell identity markers.


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Conflicts of Interest: JMI has a financial interest as co-founder and board member in Cytobank Inc., a software company for single cell data analysis. No other conflicts.


1. Irish JM, Doxie DB. Current topics in microbiology and immunology. 2014;377:1–21. [PMC free article] [PubMed]
2. Bandura DR, Baranov VI, Ornatsky OI, Antonov A, Kinach R, Lou X, Pavlov S, Vorobiev S, Dick JE, Tanner SD. Analytical chemistry. 2009;81:6813–6822. [PubMed]
3. Ornatsky O, Bandura D, Baranov V, Nitz M, Winnik MA, Tanner S. Journal of immunological methods. 2010;361:1–20. [PubMed]
4. Finak G, Frelinger J, Jiang W, Newell EW, Ramey J, Davis MM, Kalams SA, De Rosa SC, Gottardo R. PLoS computational biology. 2014;10:e1003806. [PMC free article] [PubMed]
5. Aghaeepour N, Finak G, Flow CAPC, Consortium D, Hoos H, Mosmann TR, Brinkman R, Gottardo R, Scheuermann RH. Nature methods. 2013;10:228–238. [PMC free article] [PubMed]
6. Pyne S, Hu X, Wang K, Rossin E, Lin TI, Maier LM, Baecher-Allan C, McLachlan GJ, Tamayo P, Hafler DA, De Jager PL, Mesirov JP. Proceedings of the National Academy of Sciences of the United States of America. 2009;106:8519–8524. [PubMed]
7. Maecker HT, Rinfret A, D'Souza P, Darden J, Roig E, Landry C, Hayes P, Birungi J, Anzala O, Garcia M, Harari A, Frank I, Baydo R, Baker M, Holbrook J, Ottinger J, Lamoreaux L, Epling CL, Sinclair E, Suni MA, Punt K, Calarota S, El-Bahi S, Alter G, Maila H, Kuta E, Cox J, Gray C, Altfeld M, Nougarede N, Boyer J, Tussey L, Tobery T, Bredt B, Roederer M, Koup R, Maino VC, Weinhold K, Pantaleo G, Gilmour J, Horton H, Sekaly RP. BMC immunology. 2005;6:13. [PMC free article] [PubMed]
8. Amir el AD, Davis KL, Tadmor MD, Simonds EF, Levine JH, Bendall SC, Shenfeld DK, Krishnaswamy S, Nolan GP, Pe'er D. Nature biotechnology. 2013;31:545–552. [PMC free article] [PubMed]
9. Krutzik PO, Clutter MR, Nolan GP. J Immunol. 2005;175:2357–2365. [PubMed]
10. Irish JM. Nature immunology. 2014;15:1095–1097. [PubMed]
11. Bendall SC, Simonds EF, Qiu P, Amir el AD, Krutzik PO, Finck R, Bruggner RV, Melamed R, Trejo A, Ornatsky OI, Balderas RS, Plevritis SK, Sachs K, Pe'er D, Tanner SD, Nolan GP. Science. 2011;332:687–696. [PMC free article] [PubMed]
12. Qiu P, Simonds EF, Bendall SC, Gibbs KD, Jr., Bruggner RV, Linderman MD, Sachs K, Nolan GP, Plevritis SK. Nature biotechnology. 2011;29:886–891. [PMC free article] [PubMed]
13. Sugar IP, Sealfon SC. BMC bioinformatics. 2010;11:502. [PMC free article] [PubMed]
14. Bruggner RV, Bodenmiller B, Dill DL, Tibshirani RJ, Nolan GP. Proc Natl Acad Sci U S A. 2014;111:E2770–E2777. [PubMed]
15. Newell EW, Sigal N, Bendall SC, Nolan GP, Davis MM. Immunity. 2012;36:142–152. [PMC free article] [PubMed]
16. Mosmann TR, Naim I, Rebhahn J, Datta S, Cavenaugh JS, Weaver JM, Sharma G. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2014 [PMC free article] [PubMed]
17. Naim I, Datta S, Rebhahn J, Cavenaugh JS, Mosmann TR, Sharma G. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2014 [PMC free article] [PubMed]
18. Chen X, Hasan M, Libri V, Urrutia A, Beitz B, Rouilly V, Duffy D, Patin E, Chalmond B, Rogge L, Quintana-Murci L, Albert ML, Schwikowski B, Milieu Interieur C. Clinical immunology. 2015 [PubMed]
19. Zare H, Shooshtari P, Gupta A, Brinkman RR. BMC bioinformatics. 2010;11:403. [PMC free article] [PubMed]
20. Tong DL, Ball GR, Pockley AG. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2015 [PubMed]
21. Qian Y, Wei C, Eun-Hyung Lee F, Campbell J, Halliley J, Lee JA, Cai J, Kong YM, Sadat E, Thomson E, Dunn P, Seegmiller AC, Karandikar NJ, Tipton CM, Mosmann T, Sanz I, Scheuermann RH. Cytometry. Part B, Clinical cytometry. 2010;78(Suppl 1):S69–S82. [PMC free article] [PubMed]
22. Roederer M, Moore W, Treister A, Hardy RR, Herzenberg LA. Cytometry. 2001;45:47–55. [PubMed]
23. Bagwell CB, Hudson JL, Irvin GL., 3rd The journal of histochemistry and cytochemistry : official journal of the Histochemistry Society. 1979;27:293–296. [PubMed]
24. Overton WR. Cytometry. 1988;9:619–626. [PubMed]
25. Kotecha N, Krutzik PO, Irish JM. Current protocols in cytometry / editorial board J. Paul Robinson, managing editor … [et al.] 2010:17. Chapter 10 Unit 10.
26. Spidlen J, Breuer K, Brinkman R. Current protocols in cytometry / editorial board, J. Paul Robinson, managing editor … [et al.] 2012:18. Chapter 10 Unit 10.
27. Tenenbaum JB, de Silva V, Langford JC. Science. 2000;290:2319–2323. [PubMed]
28. Becher B, Schlitzer A, Chen J, Mair F, Sumatoh HR, Teng KWW, Low D, Ruedl C, Riccardi-Castagnoli P, Poidinger M, Greter M, Ginhoux F, Newell EW. Nature immunology. 2014:1181–1189. [PubMed]
29. Roweis ST, Saul LK. Science. 2000;290:2323–2326. [PubMed]
30. Pyne S, Lee SX, Wang K, Irish J, Tamayo P, Nazaire MD, Duong T, Ng SK, Hafler D, Levy R, Nolan GP, Mesirov J, McLachlan GJ. PloS one. 2014;9:e100334. [PMC free article] [PubMed]
31. Bendall SC, Davis KL, Amir el AD, Tadmor MD, Simonds EF, Chen TJ, Shenfeld DK, Nolan GP, Pe'er D. Cell. 2014;157:714–725. [PMC free article] [PubMed]
32. Inokuma MS, Maino VC, Bagwell CB. Journal of immunological methods. 2013;397:8–17. [PubMed]
33. Becher B, Schlitzer A, Chen J, Mair F, Sumatoh HR, Teng KW, Low D, Ruedl C, Riccardi-Castagnoli P, Poidinger M, Greter M, Ginhoux F, Newell EW. Nature immunology. 2014;15:1181–1189. [PubMed]
34. Ornatsky O, Baranov VI, Bandura DR, Tanner SD, Dick J. Journal of immunological methods. 2006;308:68–76. [PubMed]
35. Spidlen J, Moore W, Parks D, Goldberg M, Bray C, Bierre P, Gorombey P, Hyun B, Hubbard M, Lange S, Lefebvre R, Leif R, Novo D, Ostruszka L, Treister A, Wood J, Murphy RF, Roederer M, Sudar D, Zigon R, Brinkman RR. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2010;77:97–100. [PMC free article] [PubMed]
36. Spidlen J, Barsky A, Breuer K, Carr P, Nazaire MD, Hill BA, Qian Y, Liefeld T, Reich M, Mesirov JP, Wilkinson P, Scheuermann RH, Sekaly RP, Brinkman RR. Source Code Biol Med. 2013;8:14. [PMC free article] [PubMed]
37. Ellis B HP, Hahne F, Meur NL, Gopalakrishnan N, Spidlen J. R package version 1.32.2.
38. Balkay L. MATLAB Central File Exchange. 2014 Jun; Retrieved.
39. Finck R, Simonds EF, Jager A, Krishnaswamy S, Sachs K, Fantl W, Pe'er D, Nolan GP, Bendall SC. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2013;83:483–494. [PMC free article] [PubMed]
40. Lee JA, Spidlen J, Boyce K, Cai J, Crosbie N, Dalphin M, Furlong J, Gasparetto M, Goldberg M, Goralczyk EM, Hyun B, Jansen K, Kollmann T, Kong M, Leif R, McWeeney S, Moloshok TD, Moore W, Nolan G, Nolan J, Nikolich-Zugich J, Parrish D, Purcell B, Qian Y, Selvaraj B, Smith C, Tchuvatkina O, Wertheimer A, Wilkinson P, Wilson C, Wood J, Zigon R, F. International Society for Advancement of Cytometry Data Standards Task. Scheuermann RH, Brinkman RR. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2008;73:926–930. [PMC free article] [PubMed]
41. Spidlen J, Bray C, I.D.S.T. Force. Brinkman RR. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2015;87:86–88. [PMC free article] [PubMed]
42. Spidlen J, Leif RC, Moore W, Roederer M, F. International Society for the Advancement of Cytometry Data Standards Task. Brinkman RR. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2008;73A:1151–1157. [PMC free article] [PubMed]
43. Kalina T, Flores-Montero J, van der Velden VH, Martin-Ayuso M, Bottcher S, Ritgen M, Almeida J, Lhermitte L, Asnafi V, Mendonca A, de Tute R, Cullen M, Sedek L, Vidriales MB, Perez JJ, te Marvelde JG, Mejstrikova E, Hrusak O, Szczepanski T, van Dongen JJ, Orfao A, EuroFlow C. Leukemia. 2012;26:1986–2010. [PMC free article] [PubMed]
44. Mahnke Y, Chattopadhyay P, Roederer M. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2010;77:814–818. [PubMed]
45. Herzenberg LA, Tung J, Moore WA, Herzenberg LA, Parks DR. Nature immunology. 2006;7:681–685. [PubMed]
46. Irish JM, Myklebust JH, Alizadeh AA, Houot R, Sharman JP, Czerwinski DK, Nolan GP, Levy R. Proc Natl Acad Sci U S A. 2010;107:12747–12754. [PubMed]
47. Moore WA, Parks DR. Cytom Part A. 2012;81A:273–277. [PMC free article] [PubMed]
48. Parks D, Roederer M, Moore WA. Cytom Part A. 2004;59A:87–87.
49. Fienberg HG, Simonds EF, Fantl WJ, Nolan GP, Bodenmiller B. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2012;81:467–475. [PMC free article] [PubMed]
50. van der Maaten L, Hinton G. J Mach Learn Res. 2008;9:2579–2605.
51. Irish J, Hovland R, Krutzik P, Perez O, Bruserud O, Gjertsen B, Nolan G. Cell. 2004;118:217–228. [PubMed]
52. Q.P. Linderman M, Simonds E, Bjornson Z. R package version 1.14.0.
53. Irish JM, Hovland R, Krutzik PO, Perez OD, Bruserud O, Gjertsen BT, Nolan GP. Cell. 2004;118:217–228. [PubMed]
54. Kotecha N, Flores NJ, Irish JM, Simonds EF, Sakai DS, Archambeault S, Diaz-Flores E, Coram M, Shannon KM, Nolan GP, Loh ML. Cancer cell. 2008;14:335–343. [PMC free article] [PubMed]
56. Gregory R. Warnes BB, Bonebakker Lodewijk, Gentleman Robert, Liaw Wolfgang Huber Andy, Lumley Thomas, Maechler Martin, Magnusson Arni, Moeller Steffen, Schwartz Marc, Venables Bill. 2015
57. G.R. Ellis B, Hahne F, Meur NL, Sarkar D. R package version 1.30.0.
58. H.F. Lo K, Brinkman R, Gottardo R. BMC bioinformatics. 2009:10. [PMC free article] [PubMed]
59. Abdelrahman AI, Dai S, Thickett SC, Ornatsky O, Bandura D, Baranov V, Winnik MA. Journal of the American Chemical Society. 2009;131:15276–15283. [PMC free article] [PubMed]
60. Abdelrahman AI, Ornatsky O, Bandura D, Baranov V, Kinach R, Dai S, Thickett SC, Tanner S, Winnik MA. J Anal At Spectrom. 2010;25:260–268. [PMC free article] [PubMed]
61. Meehan S, Walther G, Moore W, Orlova D, Meehan C, Parks D, Ghosn E, Philips M, Mitsunaga E, Waters J, Kantor A, Okamura R, Owumi S, Yang Y, Herzenberg LA, Herzenberg LA. Immunologic research. 2014;58:218–223. [PMC free article] [PubMed]
62. Irish JM, Myklebust JH, Alizadeh AA, Houot R, Sharman JP, Czerwinski DK, Nolan GP, Levy R. Proceedings of the National Academy of Sciences of the United States of America. 2010;107:12747–12754. [PubMed]
63. Geoffrey Hinton SR. Advances in neural information processing systems. 2002
64. Hahne F, LeMeur N, Brinkman RR, Ellis B, Haaland P, Sarkar D, Spidlen J, Strain E, Gentleman R. BMC bioinformatics. 2009;10:106. [PMC free article] [PubMed]
65. Van Gassen S, Callebaut B, Van Helden MJ, Lambrecht BN, Demeester P, Dhaene T, Saeys Y. Cytometry. Part A : the journal of the International Society for Analytical Cytology. 2015
66. Shekhar K, Brodin P, Davis MM, Chakraborty AK. Proceedings of the National Academy of Sciences of the United States of America. 2014;111:202–207. [PubMed]