To create a functional and accessible workflow for cytometry datasets, we were interested in building a pipeline optimized to meet criteria that address challenges in analyses of this nature ( Figure 1A ). This included 1) a framework optimized for interoperability, versatility, and ease of use, 2) designed to receive standard inputs from pre-processing steps including batch correction and 3) provide a set of outputs that can easily serve as inputs to downstream analysis and visualization tools. The pipeline required both the selection of a scalable and accurate clustering algorithm robust to downsampling and the ability to tune and evaluate clustering resolution to best meet the specifics of the biological system and scientific inquiry. Due to the wealth of existing single-cell technologies and the need for multiple measurement types to fully understand complex biological systems, we were additionally interested in developing a pipeline that would accommodate a range of single-cell cytometric data modalities including imaging.

A critical component of the pipeline was to select a clustering algorithm that could 1) scale to a large number of cells and parameters while maintaining reasonable runtimes and memory usage and 2) recapitulate populations identified via expert manual gating. We therefore evaluated a subset of available clustering algorithms, focusing on popular algorithms shown to perform well in the literature (23, 24). These included PhenoGraph (11), a graph-based community detection algorithm; CLARA (25), an extension of the partitioning around medoids algorithm; FlowSOM (9), which uses self-organizing map clustering; and PARC (10), a recent combinatorial graph-based algorithm optimized for scalability. We applied these algorithms to a CyTOF dataset measuring the expression of 42 proteins at the single-cell level on peripheral blood mononuclear cells (PBMCs) from 17 individuals (Methods, Supplemental Table 1 ). We also subset the dataset to various sizes ranging from 1,000 to 10,000 to 50,000 cells per individual to evaluate both speed and memory usage on a range of dataset sizes ( Figure 1B ). While memory usage was similar across algorithms (except for PhenoGraph performance on the full dataset, which did not conclude in a reasonable amount of time), FlowSOM and CLARA utilized the least runtime across dataset sizes. We therefore proceeded to evaluate those two clustering algorithms on the full dataset. Notably, increasing the requested cluster number identified by the algorithm across the variously sized datasets strongly increased runtime for CLARA as compared to FlowSOM ( Figure 1C ), which appeared more time efficient even with larger cluster counts regardless of cell numbers. This indicated that CLARA takes more time than FlowSOM to perform a high-resolution clustering to identify cell subsets.

While our final pipeline does focus on R-based implementation for clustering (FlowSOM and CLARA) and permits maximum interoperability with other R-based normalization, debarcoding, and batch correction packages, early investigations explored annotations from both R and python packages; thus, we benchmarked FlowSOM and PARC cluster annotations against manual gating ( Supplemental Figure 1 ) to evaluate their performance in identifying cell populations. To compare the accuracy of annotating clusters based on CyTOF panel expression to gating ground truth, we developed a custom cell barcode scheme (see Methods—FCS modifications) to identify cells in the FCS and facilitate comparison of a cell’s ground truth gating-based assignment versus how the cell was annotated by two expert immunologists. To select a resolution/grid for annotation and evaluation, we calculated the Davies–Bouldin index (DBI) (19), a within-dataset similarity metric used to evaluate cluster resolution, across a variety of cluster grid sizes ( Figure 1D , Supplemental Figure 2A ). Thirty-six clusters (6 × 6 grid, FlowSOM) or 37 clusters (Resolution 1.3, PARC) were selected for independent annotation by two expert immunologists to compare to “ground truth” manual gates from a third expert immunologist. Both ground truth gating and cluster annotations were performed for major immune cell populations (“coarse-level”, e.g., all CD4⁺ T-cell subsets including memory or naïve are annotated as CD4⁺ T) ( Figures 1E, F ) and for more fine-grained sub-populations (“fine-level”, Supplemental Figures 2B–E ); these annotations were then visualized in UMAP space. Both FlowSOM and PARC coarse annotations performed well in the recapitulation of manual ground truth gating as captured by four evaluation metrics, including accuracy, adjusted Rand index (26), Fowlkes–Mallows index (27), and mutual information (28) ( Figures 1G, H ). As expected, while performance was not as strong on fine annotations, metrics showed reasonably accurate clustering for both clustering algorithms ( Supplemental Figure 2C ). Some sources of error included the fact that global clustering did not isolate small subsets including cDC1s or antibody-secreting cells (ASCs) as their own cluster at the selected resolution due to its low abundance. In addition, the accuracy of some subsets of CD4⁺ and CD8⁺ T cells varied due to their definition based on markers that have a continuum of expression rather than clear positive and negative expression (e.g., CD45RA). However, the unsupervised nature of the clustering facilitated the identification of cell subsets not included in our manual gating, such as CD4⁻CD8⁻ T cells and intermediate monocytes ( Figure 1G ), and B-cell subsets based on the expression of CXCR5 expression or CD38 ( Supplemental 2B ). Notably, while the majority of clusters were easily recognizable and annotated, a small number of clusters were observed to be dispersed across the UMAP visualization ( Supplemental Figure 2F ), rather than having a more homogenous phenotype ( Supplemental Figure 2G ), which illuminated the utility of visualizing each cluster’s dispersion as an output for a developed pipeline. Taken together, these data demonstrate that while PARC and FlowSOM perform similarly compared to expert gating and have similar performance on low numbers of cells, FlowSOM outperforms when considering scalability to higher numbers of cells.

Based on our observations regarding runtime, memory usage, cluster evaluation methods, and performance on manual gating recapitulation, we developed the Cyclone pipeline in R ( Figure 2 ). We designed Cyclone to produce a series of outputs that users can access to select cluster resolution, evaluate cluster quality, annotate clusters, and utilize resulting cluster statistics (phenotypes and abundance) for downstream analyses. To start Cyclone, the user provides FCS or matrix data files as well as files to specify markers and file metadata. Cyclone expects any normalization and batch correction steps to be performed prior to the use of the pipeline. For CyTOF data, this can be performed by established R packages, including premessa (github.com/ParkerICI/premessa) for bead-based normalization and CytoNorm (21) or cyCombine (22) for batch correction. After reading in and arcsinh transforming the data, Cyclone calculates UMAP dimensionality reduction. Cyclone works with either FlowSOM or CLARA for clustering; selecting FlowSOM enables DBI-based cluster resolution optimization prior to user grid selection, while with CLARA, the user selects a single resolution for clustering. If FlowSOM is selected, Cyclone then performs iterative optimization of clustering across a variety of cluster grid sizes, which can then be compared using DBI. After the user selection of desired grid, Cyclone performs clustering, generates summary matrices, proceeds through an optional SCAFFoLD step, and then generates output files and visualizations. While not providing batch correction, Cyclone does provide a means for assessing batch or any other input file metadata via UMAP (split by batch) visualization ( Figure 3A ) as well as clustered heatmaps of cluster frequency with batch or other metadata information overlaid as rugs ( Figure 3B ). In addition to a DBI plot for cluster resolution selection, Cyclone outputs UMAPs colored by cluster ( Figure 3C ) as well as heatmaps of marker expression per cluster for ease of cluster annotation ( Figure 3D ). Additional UMAPs and histogram plots are provided per cluster, showing the cluster’s distribution across the UMAP for an evaluation of cluster dispersion as an indication of cluster quality ( Figure 3E ). Feature UMAP plots of marker expression in all cells are also exported by default ( Figure 3F ). Taken together, the Cyclone pipeline provides a means of clustering, as well as evaluating and annotating these clusters, to then be readily used in downstream analyses.

We developed Cyclone to be easily leveraged by researchers with various backgrounds, especially those with limited coding knowledge and minimal computational resources. We also prioritized interoperability with both upstream and downstream tools for cytometry data processing and analysis ( Figure 4A ). In addition to significant documentation (including vignettes to get users started with the pipeline), we evaluated Cyclone on downsampled datasets to determine whether downsampling to fewer cells could provide an alternative for users to run Cyclone locally rather than needing additional compute resources for analyses. Thus, we downsampled our evaluation dataset to 50,000 cells per sample and coarsely annotated the resulting optimized clusters ( Figure 4B , Supplemental Figures 3A, B ). We then compared the accuracy of these coarse downsampled annotations to the coarse annotations on the full dataset and encouragingly found strong concordance between the full and downsampled datasets ( Figure 4C ). When doing fine annotations ( Supplemental Figure 3C ), metrics again showed reasonably accurate annotations when compared to the full dataset fine annotations ( Supplemental Figure 3D ). Thus, data downsampling presents a practical option for dataset clustering with Cyclone should computing resources be a challenge for some users.

To optimize Cyclone interoperability with upstream and downstream processing, we considered batch correction as a primary upstream target and SCAFFoLD map analysis and accessible visualization as primary downstream targets. For batch correction, we ensured a pathway exists to enter the Cyclone workflow with cytometry data in either a matrix format [as output by cyCombine (22)] or adjusted-FCS files [as output by CytoNorm (21)] to accommodate interoperability with outputs from both commonly used batch correction methods. CytoNorm’s adjusted-FCS files can enter into the pipeline as normal, but Cyclone’s prepare_checkpoint1() provides an entry point for cyCombine’s matrix format. The function takes the same primary inputs as a run without a batch correction step, except for accepting a matrix (with cells in rows and markers in columns) instead of the FCS directory. The function then performs all the same optional steps (arcsinh transformation, control sample removal, and subsampling) before outputting a Checkpoint1.Rdata file. Afterward, users can run Cyclone as normal to continue from the second step of UMAP calculation.

To aid in cluster annotation, SCAFFoLD connects clusters to unique landmark populations based on the cosine similarity in the feature expression space (15). For powering SCAFFoLD map downstream analysis, we added an optional step directly into the Cyclone pipeline; if landmark population FCS files are provided, an output *.scaffold file can be used to generate a SCAFFoLD map (15) via the “scaffold::scaffold.run()” command ( Figure 4D ). For powering accessible visualization and other downstream follow-ups, we ensured that the colors of Cyclone plot outputs use color blindness-accessible color palettes and that data outputs can be used with dittoSeq–a color blindness-friendly visualization tool (16). Although dittoSeq was designed for single-cell RNAseq data, it proves generalizable to other data modalities including high-dimensional cytometry. To directly power dittoSeq integration and downstream statistical analyses, we provide an additional vignette showing how to 1) transform Cyclone data objects into a SingleCellExperiment [(29), Chapter four] object compatible with dittoSeq, 2) generate useful visualizations, 3) add cell-type annotations for each cluster as well as further annotations such as spatial information in the case of imaging data, and 4) run statistics on differences in cluster or cell-type frequencies between samples. Shown as examples are boxplots of how cluster frequencies compare between male and female subjects ( Figure 4E ) and the percent composition of each cluster in terms of subject sex ( Figure 4F ); these functions make probing metadata categories of interest easy to code and assist in producing publication-ready visualization. These accessible (in both the color palate and easily leveraged graphing functions) visualizations are readily created with dittoSeq from Cyclone outputs.

While optimization and development of Cyclone were initially for CyTOF datasets, we realized the need for such analysis techniques on a broader set of similar data modalities, including other cytometry platforms such as spectral flow cytometry. To evaluate Cyclone’s utility for this analysis, we applied the Cyclone pipeline to a spectral flow cytometry dataset of mouse liver samples from a transgenic mouse model for hepatitis B virus (HBV) response (30). Liver leukocytes collected from eight mice on day 8, the peak of the immune response in this model, were analyzed using a myeloid-focused spectral flow cytometry panel ( Supplemental Table 2 ). In a repeat experiment, day 8 liver leukocytes from nine mice were analyzed using a lymphoid-focused panel ( Supplemental Table 3 ). The Cyclone pipeline could be used as built, but we used config files to specify the co-factor for the arcsinh transformation to a value more typical of flow cytometry datasets (see Spectral flow data generation and analysis—Cyclone analysis in Methods). For both experiments, CD45⁺ live cells were provided to the Cyclone pipeline. After selecting the local minimum DBI, we identified and annotated 21 clusters in the T Cell panel ( Figure 5A ) and 17 clusters from the Myeloid panel ( Figure 5B ). Both panels were able to distinguish several T-cell and NK-cell subsets ( Supplemental Figure 4A ) as well as liver resident macrophages (Kupffer cells) and monocyte-derived macrophages ( Supplemental Figure 4B ). As was previously observed for CyTOF data, unsupervised clustering using Cyclone largely recapitulated the populations identified by manual expert gating at a coarse-level ( Figures 5C, D ). Cyclone further enabled the identification of cell subsets across immune cell lineages ( Figures 5C, D ), which could have been missed by manual gating, such as NK-cell subsets with various expression of KLRG1 or CD62L ( Figure 5A , Supplemental Figure 4A ). To evaluate the biological information contained within the unsupervised clustering, we took advantage of the presence of tetramer staining to identify the properties of the HBV-specific T cells. While the tetramer staining was not used as a clustering parameter, we found that the vast majority of CD8⁺ ( Figure 5E ) and CD4⁺ ( Figure 5F ) antigen-specific T cells were enriched in one or two effector T-cell clusters with high expression of activation markers identified by Cyclone for each T-cell subset ( Figure 5G ). Taken together, we observed that Cyclone performed well on spectral flow data and enabled the unsupervised identification of cell phenotypes that are associated with distinct biological features.

Understanding how the phenotype of individual cells relates to the function of multicellular compartments within tissues requires the ability to identify cellular phenotypes with multiple proteins while simultaneously quantifying the spatial distribution and interactions of these cells across large regions of tissue. Along with its applications demonstrated already, Cyclone provides a unique opportunity to analyze the spatial distribution of individual cell phenotypes as well as cellular neighborhoods in an unsupervised manner. To that end, we tested how the Cyclone pipeline compared to a prototypical imaging data analysis pipeline containing image visualization. First, we created a 7-plex immunofluorescent staining panel to be used on a colorectal tumor tissue (CRC1) consisting of a tumor marker (EPCAM), T-cell markers (CD3, CD4, and CD8), and myeloid markers (CD163, HLA-DR, and XCR1) ( Figure 6A ). Next, we utilized DeepCell (31) segmentation software to demarcate individual cells by inputting nuclei (DAPI⁺) and membrane markers (CD3⁺, CD4⁺, CD8⁺, and CD163⁺) and labeled each cell-type annotation based on rational gating parameters (e.g., CD4⁺ T cell = CD3⁺CD4⁺CD8⁻XCR1⁻) ( Supplemental Figure 5A , table). We found that determining the manual thresholding on certain markers such as CD163⁺ expression was not visually clear in designating a suitable cutoff as compared to CD3⁺ expression and that this could ultimately lead to variability in the frequency of cell types annotated ( Supplemental Figure 5A , histogram plots). For example, the threshold cutoff at 11 observed CD163⁺ expression on low-background, segmented cells, while the threshold cutoff at 13 missed CD163⁺ cells. Thus, we opted for a middle ground by choosing threshold 12 for the downstream comparison with our Cyclone pipeline. To run the Cyclone pipeline on this multiplexed immunofluorescence data, we followed the following procedure. Once raw expression values of each marker for every cell were curated and assigned to each cell identified by DeepCell, we obtained a cell per protein expression matrix used to run the Cyclone pipeline (see Methods). After DBI evaluation, we chose a grid size of 2 × 4, which had the lowest DBI value ( Figure 6B ). This resolution generated eight unique clusters spanning immune and non-immune cells ( Figure 6C ). These clusters comprised tumor cells expressing different levels of HLA-DR (Clusters 2 and 4), CD4⁺ and CD8⁺ T cells (Clusters 3 and 8), mononuclear phagocytes (MNPs) (Clusters 5, 6, and 7), and one cluster with low expression for all markers (Cluster 1) ( Figure 6D ). Interestingly, among the different MNPs, Cluster 7 is characterized by the high expression of XCR1⁺ and CD163⁺ and low HLA-DR expression. Since the combination of the Cyclone and DeepCell segmentation outputs provides cluster identities (Cyclone) and x,y coordinates (DeepCell) for each cell within the same data frame, we leveraged this to evaluate where these XCR1⁺CD163⁺ cells were located within the tissue. When overlaying Cluster 7 onto the image of CRC1, cells depicted as Cluster 7 (green arrows), we were able to confirm that they had both XCR1⁺ and CD163⁺ expression ( Figure 6E ). Additionally, these cells had an enriched spatial distribution within the stromal compartment as opposed to cDC1s and other clusters enriched in the tumor compartment ( Figures 6F, G , Supplemental Figure 5B ). We next validated that XCR1⁺CD163⁺ cell types could be found in other colorectal (CRC2) and kidney (KID1) tumor samples subjected to Cyclone ( Supplemental Figure 5C ). However, these samples observed different stromal and tumor enrichments for CD4⁺ T cells, cDC1s, and XCR1⁺CD163⁺ cells compared to CRC1.

Notably, while conventional manual gating for XCR1⁺CD163⁺ cells is possible, the thresholding strategy of comparing independent markers with strict cutoffs likely misses out on what Cyclone identifies when leveraging the entire set of markers to identify cell types ( Supplemental Figure 5D ). Taken together, we highlight the successful application of the Cyclone pipeline on multiplexed imaging data in tumor tissue to discriminate two different XCR1⁺ subsets and their spatial features within the tumor microenvironment.

To scope clustering algorithms for our pipeline, we tested four Python- or R-based widely used clustering tools for cytometry data: CLARA, FlowSOM, PhenoGraph, and PARC. CyTOF analysis allows users to capture hundreds of thousands of cells, and clustering such large datasets requires runtime- and memory-efficient tools that do not compromise clustering performance. One aim of our work was to design a scalable pipeline for large datasets (containing many samples, each with hundreds of thousands of cells). Therefore, we decided to compare the runtime and memory usage of the selected clustering tools. Different parameters affect the runtime and memory usage in different tools. We observed the parameters affecting the number of clusters were the ones that controlled algorithm runtime. To perform a fair comparison between the tools, we identified these parameters affecting cluster counts and used values that produced a similar number of clusters across all four tools. We used k = 24 in CLARA, xdim = 6 and ydim = 6 in FlowSOM, k = 25 in PhenoGraph, and resolution = 1.3 in PARC. These parameters resulted in 30–37 clusters called by each tool. We measured the time for clustering and the memory usage on CentOS nodes of a high-performance computer cluster.

After clustering the data using FlowSOM or PARC, we performed manual annotation of the resulting clusters from both tools based on the median expression of markers. Since each individual may annotate clusters differently, we attempted to account for human-to-human variations in the manual annotation of clusters by having two immunology experts independently annotate the clusters from FlowSOM and PARC. We established “coarse” and “fine”-levels of annotations. Coarse annotations describe cells of different compartment groups (e.g., CD4⁺ T cells, CD8⁺ T cells, and B cells). Fine annotations further parse cell types into subtypes according to their phenotype (e.g., CD8⁺ T cells are split further into naïve, central memory, effector memory, and effector memory re-expressing CD45RA). We then calculated similarities between cluster annotations of single cells and the single-cell annotation based on a third immunology expert gating of the CyTOF data (ground truth annotation). We performed this comparison for subsets of data multiple times and calculated averages. Specifically, we subsampled cluster annotation and manual gating annotation of randomly selected 10,000 cells, calculated accuracy, and adjusted the Rand index, Fowlkes–Mallows index, and mutual information. We repeated this process 10 times with different random seeds and calculated the mean across the iterations for each similarity metric. We performed this analysis for FlowSOM ( Figure 1H , Supplemental Figure 2C ) and PARC annotations (data not shown) from both immunology experts.

Cyclone outputs (checkpoint1.Rdata and checkpoint8.Rdata) were used to create a SingleCellExperiment object (29) in R containing the arcsinh transformed expression matrix, UMAP embeddings, clustering, and cell and sample-level metadata. dittoSeq (16) functions dittoFreqPlot and dittoBarPlot were then used to create boxplots of cluster frequencies per sample and stacked bar plots of batch composition per cluster, respectively.