Background fluorescence and spreading error are major contributors of variability in high‐dimensional flow cytometry data visualization by t‐distributed stochastic neighboring embedding

Abstract Multidimensional single‐cell analysis requires approaches to visualize complex data in intuitive 2D graphs. In this regard, t‐distributed stochastic neighboring embedding (tSNE) is the most popular algorithm for single‐cell RNA sequencing and cytometry by time‐of‐flight (CyTOF), but its application to polychromatic flow cytometry, including the recently developed 30‐parameter platform, is still under investigation. We identified differential distribution of background values between samples, generated by either background calculation or spreading error (SE), as a major source of variability in polychromatic flow cytometry data representation by tSNE, ultimately resulting in the identification of erroneous heterogeneity among cell populations. Biexponential transformation of raw data and limiting SE during panel development dramatically improved data visualization. These aspects must be taken into consideration when using computational approaches as discovery tools in large sets of samples from independent experiments or immunomonitoring in clinical trials.

In the last years, much effort has been dedicated to harmonize flow cytometry assays for basic science and clinical trials, especially in regard to machine standardization (3)(4)(5)(6), reagents (7), protocols (8), and data analysis (9,10). The technology is extremely robust and reproducible in identifying rare cell populations in independent experiments performed by different operators (11), rather much variability is introduced when independent users analyze the same data due to subjective gating (9). As a consequence, unsupervised algorithms for the identification of cell populations in highdimensional data sets would provide an objective approach for data analysis (9). As these approaches make the use of every single cell in the input, technical variability between independent experiments or measurement errors may easily generate artefacts and lead to the misidentification of cell populations. Here, we tested the applicability of tSNE to 27-parameter single-cell flow cytometry data.

Sample Collection
All experiments using human buffy coats were approved by the Humanitas Research Hospital IRB. Peripheral blood mononuclear cells (PBMCs) were isolated from buffy coats by Ficoll separation and frozen in liquid nitrogen according to standard procedures, as described in Ref. (11).

Flow Cytometer Quality Control (QC)
Machine QC was performed according to the protocol developed by Perfetto et al. (3). Briefly, optimal voltage settings of PMTs were determined in two steps: first, by defining the separation of dimly stained Cytocal beads (Thermo Fisher Scientific, Waltham, MA, USA) over unstained compensation beads (Compbeads; BD Biosciences, San Jose, CA, USA) at 50 V intervals (range 350-800 V), and second, by defining the separation of quantum simply cellular beads (QSCB; previously stained with fluorescently conjugated monoclonal antibodies) over unstained Compbeads at AE25 V of the target voltage defined in the first step. Inter-experiment QC was performed by running single-peak Rainbow beads (Spherotec, Lake Forest, IL, USA) and unstained Compbeads and by adjusting PMT voltages to reach target values. Laser delays were adjusted manually.

Polychromatic Flow Cytometry for the Detection of Surface and Intracellular Antigens
All data were acquired between January and December 2017 on a FACS Symphony A5 flow cytometer (BD Biosciences) equipped with five lasers (UV, 350 nm; violet, 405 nm; blue, 488; yellow/green, 561 nm; red, 640 nm; all tuned at 100 mW, except UV tuned at 60 mW) or with a FACSAria cell sorter equipped with four lasers (violet, 405 nm; blue, 488; yellow/green, 561 nm; red, 640 nm; all tuned at 50 mW). Flow cytometry data were compensated in FlowJo by using single-stained controls (Compbeads incubated with fluorescently conjugated antibodies), as described in Ref. 12. Cells previously frozen in liquid nitrogen were used in all experiments. Vials of PBMCs from three different donors were thawed in five independent, replicate experiments and stained with a 25-color flow cytometry panel of antibodies as described below. Cells were added to RPMI1640 supplemented with 10% fetal bovine serum and 1% penicillin-streptomycin, 1% L-glutamine, and 20 μg/ml deoxyribonuclease I from bovine pancreas (Sigma-Aldrich, St. Louis, MO, USA). Subsequently, cells were extensively washed and stained immediately with the Zombie Aqua Fixable Viability Kit (BioLegend, San Diego, CA, USA) in phosphate-buffered saline for 15 min at room temperature (RT). Cells were then washed and stained with the combination of antibodies, purchased from BD Biosciences, BioLegend, or eBioscience, as indicated in Supporting Information Table S1. Antibodies were previously titrated to define optimal concentration, as described in Refs. 12,13. Staining of CXCR5, CCR7, and CXCR3 was performed for 20 min at 37 C, while all other surface markers (except CD3) were stained for 20 min at RT. Ki-67, CD3, and transcription factors were detected intracellularly following the fixation of cells with the FoxP3/transcription factor staining buffer set (eBioscience, Waltham, MA, USA) according to manufacturer's instructions and by incubating with specific antibodies for 30 min at 4 C.

High-Dimensional Single-Cell Data Preprocessing and Analysis by tSNE
27-Parameter Flow Cytometry Standard (FCS) 3.0 files were imported into FlowJo software version 9 or 10 and left untreated or biexponentially transformed (the same transformation for all files, performed in version 10) prior to tSNE analysis. Alternatively, background fluorescence events below the arbitrary threshold of 100, the maximum common value for all fluorescence parameters that is below the threshold of positive antigen expression in this dataset, were transformed by a process referred to as "background randomization." Briefly, all fluorescence events comprised in the interval between the lowest value of the log distribution (as obtained from the machine) and 100 were randomly distributed evenly across all channels. Replicate samples are expected to contain a similar number of events in this interval, although with different distributions across fluorescence channels, and for this reason are forced to have the same distribution by background randomization.
Following transformation, 27-parameter samples were assigned with a computational barcode for their unique identification, concatenated (3,000 CD8 + T cells per sample for a total of 45,000 CD8 + T cells) and visualized with tSNE (Barnes-Hut implementation) in FlowJo. The following parameters were tuned in preliminary experiments (14), then used at default values (unless otherwise specified) as we observed no substantial changes in the final distribution: iterations, 1000; perplexity, 40; initialization, deterministic; theta, 0.5; eta, 200. In some experiments, perplexity was changed to determine variability in tSNE cluster distribution and data display/representation ( Figure 3a). All parameters except for lineage markers (CD3, CD4, and CD8), Aqua dead cell marker and Eomes transcription factor, were included in the analysis.
To test the influence of background fluorescence on the final tSNE representation independently of FlowJo software, PBMC samples from two different healthy donors were stained with a nine-color panel (Supporting Information  Table S1) within the same experiment and acquired on a FACSAria cell sorter. Markers were selected in a way that CD8 + T N cells are defined only by negative expression, that is, CD56 − , HLA-DR − , CD4 − , and CD45RO − . Cells were further gated as CCR7 + , although CCR7 was not included in the analysis. In this way, only background, and not positive fluorescence (i.e., protein abundance on the cell surface) defines CD8 + T N cells in the tSNE display. Compensated FCS 3.0 files (containing 20,000 cells/sample) were exported as .txt, concatenated in R analyzed using the Rtsne function (from the Rtsne package, using default options with the exception of PCA parameter that was set equal to "FALSE") to calculate the coordinates for the 2D tSNE plot. Resulting data were further converted in FCS 3.0 format and analyzed with FlowJo 10.

Influence of Spreading Error (SE) on tSNE Map Calculation
To test the effect of SE (15) on the tSNE map calculation, we designed a polychromatic panel where SE is contributed only by CD57 BV605 spreading into the YG610 channel (in this experiment detecting anti-CCR7 PE-CF594; Figure 2a). A PBMC sample was stained as above with a backbone panel including Zombie Aqua and antibodies directed to recognizing the following antigens: CD3, CD4, CD8, CD27, CD45RA, CD45RO, CCR7, and CD56. After wash, cells were split into four replicate tubes, three of which were subsequently stained with increasing amounts of anti-CD57 conjugated to BV605. The fourth tube was left unstained as a control. Then, 20,000 cells from each of the four samples (80,000 in total) were concatenated and visualized with tSNE. The following markers were given as input: CD27, CD45RA, CD45RO, CCR7, and CD56. To make sure that the final result was dependent only on SE and not by the abundance of CD57 expression, the CD57 parameter was excluded from tSNE calculation.

RESULTS
We evaluated the usefulness of tSNE to determine the heterogeneity of multidimensional single-cell data generated by 27-parameter polychromatic flow cytometry. For this purpose, PBMC samples from three healthy donors were thawed in five independent, replicate experiments and stained with 24 different fluorescently conjugated monoclonal antibodies directed to recognize antigens expressed by CD8 + T cells. A viability dye was used to exclude dead cells. The full gating strategy used to identify CD8 + T cells is depicted in Supporting Information Figure S1a. A single concatenated file containing 45,000 events from the different experiments was subjected to tSNE analysis (see "Methods" section for details on the procedure). Following debarcoding of replicate samples, we noticed that tSNE maps were substantially different in independent experiments, specifically runs 1-2, and runs 3-4 were mostly similar to each other (on the basis of the abundance of cells identified in the arbitrary gate) and that those from run 5 were the most different compared to runs 1-4 ( Figure 1a). We thus aimed to investigate the reason at the basis of such variability, which characterized both tSNE1 and tSNE2 dimensions (Figure 1b). First of all, we tested whether samples from different runs harbored differential abundance of immune populations, for example due to the loss of certain subsets as a consequence of the freeze/thaw procedure. Therefore, we quantified the expression of single markers, one at a time, by standard gating and identified similar percentages of antigen expression (Supporting Information Figure S1b). These data also confirm the reproducibility of the flow cytometry technology. Overall, more variability between measurements could be observed for those markers not having a clear separation over negative cells, in our case T-bet and Eomes. In particular, Eomes displayed the highest variability among replicate stainings and was thus excluded from further tSNE analysis.
Differential distribution of events in the tSNE map among independent, replicate samples could be due to real biological heterogeneity (e.g., subsets composition) or to technical artifacts in the single-cell data. To test this, we reasoned that a phenotypically homogeneous T-cell population would be identified by a single island in the map. Therefore, we selected highly purified CD8 + naive T (T N ) cells on the basis of a combination of six markers that are known to be associated with naivety (i.e., CD45RA + CCR7 + CD27 + CD95 − Tbet − CD73 + ; Figure 1c) (16,17) and plotted them on top of the tSNE map. Surprisingly, tSNE classified T N cells into multiple putative subpopulations (Figure 1d), despite the relative homogeneity of T N cell immunophenotype (either positive antigen expression or negative antigen expression; Figure 1e). Since data points (cells) are distributed in a stochastic way in the final tSNE map, we tested whether the final representation is dependent on the initial seed. To this end, we ran the algorithm several times with random seeds, and obtained very similar results (data not shown). We then thought of additional factors other than true biological heterogeneity that could drive the "clustering" of highdimensional flow cytometry events in the tSNE space, including differences in (1) compensation, (2) daily instrument performance, (3) positive signals, as a consequence of point 2 but also possibly caused by some sort of batch effect, as well as differences in background values (i.e., the biological negative fraction not containing fluorochrome-conjugated antibodies), and (4) SE generated by errors in the measurement of photons at red/far red wavelengths (15).
The impact of fluorescence compensation on the final tSNE representation is intuitive, as incorrect compensation may result in newly generated populations of cells in the multidimensional space. Despite subtle differences in the final compensated data may occur between samples, correct compensation was generated by running experiment-specific, single-stained controls and by carefully examining single pairwise combinations (12)   background fluorescence values in different experiments (Supporting Information Figure S3). We observed slight changes in background levels only in run 3, which can hardly explain differential tSNE distributions across the five runs. Positive signals from reference beads were highly comparable. As for beads, we observed reproducible measurements of antigens that are expressed on CD8 + T N cells, that is, CXCR3, CD73, CD45RA, CD27, CCR7, and CD98 (Figure 1e), indicating that a possible batch effect due to sample preparation is negligible. Instead, we noticed substantial variability in the distribution of the background for some parameters such as CD25, CD57, CD69, and Ki67 (Figure 1e), which ultimately affected the final tSNE output by erroneously identifying virtual subpopulations of cells (Figure 1f, left). Such putative heterogeneity identified by tSNE among CD8 + T N cells has no biological meaning (indeed CD8 + T N cells express little, if any CD69, CD25, or CD57; Figure 1e) and seems to be driven, at least in part, by cells with very low fluorescence values that are piled on the low end of the log scale, an effect generally referred to as the "log artifact" (18). By plotting the differential fluorescence intensity of background events on top of the tSNE map depicting T N cells, as identified from Figure 1d, it was possible to appreciate that T N "clusters" indeed harbor different background values of CD25, CD57, and CD69 (Figure 1f, right). Randomization of background, that is, the redistribution of negative events evenly across all channels below an arbitrary threshold (set at 100 for all parameters), improved the reproducibility of the tSNE plots (Supporting Information Figures S4a and S4b), thereby demonstrating that background values are indeed a source of variability.
Despite the reproducibility of the flow cytometry assay, a residual batch effect may still be present due to small variability in positive fluorescence across independent experiments. To overcome this, we analyzed two samples that were stained within the same experiment with the same analytical workflow. We also validated our results independently of the flow cytometer and the analytical platform. Also in this case, tSNE displays of CD8 + T N cells from the two samples were distinct (Supporting Information Figure S5a) and correlated with different levels of CD45RO, but not of CD56, HLA-DR, and CD4 background levels (Supporting Information Figures S5b  and S5c). Specifically, low background levels of CD45RO were associated with high tSNE1 axis and low tSNE2 axis values, and vice versa. These data confirm our previous observation that the influence of background fluorescence on tSNE display is a generalized phenomenon.
We next tested the contribution of SE to variability in tSNE representation. SE is directly proportional to the square root of the intensity of the staining, may vary between samples as a result of differential antigen abundance, and can potentially generate new cell clusters as a consequence of fluorescence spreading in the "negative" and "positive" regions of secondary dimensions (Figure 2a). Despite having exactly the same immunophenotypes and virtually the same biological functions, one could hypothesize that cells in the spread interval could be inaccurately mapped to different tSNE regions because defined by very different values in the multidimensional space. To test the individual contribution of SE, we mapped "negative" and "positive" SE on the top of a control sample that did not contain CD57 BV605 (Figure 2b). In this way, we noticed that while cells with positive SE values overlap to tSNE islands already occupied by other cells (Figure 2c; indeed, cells double positive for CD57 and CCR7 may exist in the sample), cells with negative SE values fill empty regions of the tSNE plot (Figure 2d), thus further contributing to generate new islands with irrelevant biological meaning. In conclusion, different background measurements and SE are major sources of variability in multidimensional single-cell polychromatic flow cytomertry data representation by tSNE.
SE is an intrinsic physical characteristic of the dyes and cannot be corrected following data acquisition (4). However, SE can be minimized by mounting more powerful lasers so to increase photon emission, and by optimizing filter/mirror combinations or during panel development, if required in specific flow cytometry panels so to increase detection of dimly expressed antigens (12). In contrast, background fluorescence representation can be modified by using computational approaches, the most common of which is biexponential transformation (18). While positive data remain untouched, this approach circumvents the log artifact by normalizing data around zero, thus decreasing background variability between samples (Supporting Information Figure S6). Figure 3a shows that biexponential transformation dramatically improved tSNE representation of five replicate 27-parameter flow cytometry analyses. Improvement could be also observed by correlating tSNE axes values from two experiments with the highest variability (i.e., run 1 and run 5; Figure 3b). Finally, CD8 + T N cells (identified as in Figure 1c) appeared relatively homogeneous following biexponential transformation when compared with nontransformed data (Figures 3c). Similar results were obtained with a different panel (Supporting Information Figure S7).

DISCUSSION
Identification of trends in multidimensional single-cell data that are associated with a certain pathology or immune status has been pursued since the development of flow cytometers capable of measuring more than 7-8 colors (19). The advent of new reagents and technologies such as CyTOF and scRNAseq extended this capability to dozens and thousands of parameters, respectively, thus making the development of new computational approaches for the visualization of such data a fundamental need (20). tSNE is currently the most popular approach in this regard as it displays single cells on an intuitive 2D graph and where relative antigen or gene expression can be visualized by the color tone (2). Despite limited in its analysis capability compared to the abovementioned technologies, polychromatic flow cytometry, now capable to measure up to 30 parameters simultaneously, is still the most popular and versatile single-cell technology that is available in thousands of laboratories around the world. It is therefore anticipated that approaches such as tSNE will be largely utilized to display such data. While flow cytometry and CyTOF share multiple aspects in terms of protocols, reagents, and concept, the nature of their single-cell data is different. Since flow cytometry relies on fluorescence measurements, scientists have to inevitably deal with background fluorescence, compensation, and SE that otherwise are poorly present in CyTOF.
Here, we have identified background fluorescence and SE to impact substantially the visualization of multidimensional single-cell data by tSNE, specifically resulting in the generation of virtual "islands" that do not reflect real heterogeneity at the phenotypic or biological level. Differences in background, but not positive fluorescence, could be observed among replicate samples run in independent experiments. This was not due to machine performance, which confirmed to be very stable over time, rather variability in background seems to be sample dependent (Supporting Information Figure S4). These low (negative) fluorescence signals are irrelevant when analyzing data by standard gating, as no biological information beyond the number of cells contained in the negative fraction should be considered. Nevertheless, fluorescence "heterogeneity" in the negative interval may lead to substantial confusion when displaying multidimensional datasets with approaches such as tSNE, where every single fluorescence value, being it positive or negative, is considered.
Randomizing negative events below an arbitrary threshold evenly across all channels improved visualization but did not solve the issue in full, probably because of the lack of an algorithm capable to introduce channel-specific thresholds (in fact, threshold of positivity is different in every single channel; Supporting Information Figure S4c), rather biexponential transformation, which normalizes negative data around zero without impacting the positive ones, substantially reduces the variability among independent experiments, thus enabling to obtain comparable tSNE plots of replicate, independent measurements. However, it was not possible to obtain complete overlap of tSNE plots between the biologically identical samples (Figure 3a). We attribute this to small variations in background or positive fluorescence, compensation or SE that ultimately affect multiple combinations in the multidimensional space and that are captured by the tSNE algorithm. Indeed, we further identified SE as a major contributor of variation in multidimensional flow data representation. This study describes a situation where a single pairwise fluorochrome combination that is dramatically affected by SE is sufficient to generate cell clusters with irrelevant biological meaning. As SE can affect multiple combinations in the 30-parameter space (4) (Lugli and Roederer, unpublished observation), it is anticipated that data representation by tSNE would be progressively worsened when increasing the number of parameters.
A wrong assumption is that cells occupying different areas of the tSNE display harbor intrinsic biological heterogeneity; therefore, it is not rare to see the use of tSNE maps to identify putative subsets in a given population or set of samples. On the basis of our results, gating on the tSNE map followed by phenotypic analysis would lead to overfragmentation of immunophenotypes that complicates, rather than simplifies, multidimensional data analysis. Instead, clustering algorithms should be run in parallel to identify relevant cell subpopulations, as it is now done for scRNAseq. Along with this, it is of foremost importance to validate computational results with raw data, i.e., once the clusters are identified, to evaluate the relative abundance of antigen expression in that specific subpopulation either by percentage or MFI, if not both, followed by critical evaluation of immunophenotypes.
Moreover, beyond machine, reagent, and computational standardization, immunophenotyping of replicate samples must be included in independent experiments, so to estimate the overall variability of the assay. Integration of classical gating and novel computational approaches will be pivotal to achieve the discovery of novel biological insights in high-dimensional polychromatic flow cytometry data.