An R-based reproducible and user-friendly preprocessing pipeline for CyTOF data

Data description

The data used in this pipeline were generated in the context of the Tumor Profiler project, a multi-center observational study investigating the relevance of different innovative technologies, including CyTOF, imaging mass cytometry, single-cell DNA and RNA sequencing, as well as ex vivo drug testing to improve the diagnostic of advanced cancer patients¹⁵.

The samples of interest included blood samples and tumor biopsies collected at the University Hospital Zurich in spring 2020. These samples were assessed by mass cytometry in the context of a set of references including commercially available cell lines, PBMCs from healthy donors and PHA activated PBMCs. PBMCs from patients and healthy donors were collected based on a ficoll gradient¹⁶, and tumor samples were dissociated as previously described¹⁷. Once in single-cell suspension, all samples were stained for 5 min on ice with Cell-ID TM Cisplatin-194Pt (#201194, Fluidigm) to identify dead cells and subsequently fixed with PFA 1.6% (#15710, Electron Microscopy Sciences). Samples were stored as dry pellet at −80° C until CyTOF measurement.

The dataset used in this study was obtained from a single CyTOF run containing nine references, three blood samples and three tumor samples barcoded with a 20-well barcoding plate¹⁷. Reference samples were selected to contain positive and negative populations for each marker included in the study’s antibody panel. This design was chosen to enable for quality control and batch correction across independent experiments based on quantile scaling as described in Schuyler et al.¹¹. After barcoding, pooled cells were stained at 4 degree for 1h with a 40-Ab panel designed to perform an in-depth characterization of the samples’ immune compartment. DNA intercalation was performed with a 1h room temperature incubation in Cell-ID Intercalator-Ir (#201192B, Fluidigm). Finally, the cell suspension was diluted 1:10 in Maxpar® Cell Acquisition Solution (#201240, Fluidigm) and 10% of EQ Four Element Calibration Beads (#201078, Fluidigm), and acquired on a Helios^TM upgraded CyTOF 2 system at a flow rate of 150 events per second.

Throughout this workflow, we will make use of a set of metadata for standard preprocessing steps (normalization, debarcoding and compensation), as well as various quality controls. An overview of the metadata used is given in Table 1. All data used are provided as underlying data¹⁸.

Data organization

Most data used and returned throughout this workflow are kept in an object of Bioconductor’s SingleCellExperiment (SCE) class from the SingleCellExperiment package¹⁹. This data structure can store all single-cell related data (measurement data and transformations thereof; cell, feature and experiment-wide metadata; dimensionality reductions), allowing for synchronized and thus less error-prone data manipulation.

The key component of SCEs are matrix-like assays, where rows are features (targets) and columns are observations (cells), that store the measurement data and any data derived thereof. Metadata associated with cells are stored under colData, feature metadata under rowData, and any experiment-wide metadata may be stored in the metadata slot. Lastly, the SCE can store an arbitrary number of dimensionality reductions under reducedDims. For a more detailed description of usage and structure of SCEs, we refer to the SingleCellExperiment package’s documentation.

Constructing a SingleCellExperiment

By default, flowCore’s read.FCS() function, which underlies read.flowSet() for reading in a set of FCS files, transforms channel intensities and removes events with extreme values. To omit this behavior, we recommend reading in files with arguments transformation = FALSE and truncate_max_range = FALSE; by default, files will be read in by CATALYST’s prepData() function with these settings.

As described above, the SCE class allows the keeping of multiple data transformations in a single object. Thus, when applying a transformation to arrive at expression-like data, we can store the transformed data in a separate assay without overwriting the raw ion count data. In this way, any data generated and used throughout preprocessing (e.g., normalized, compensated or batch-corrected counts and their arcsinh-transformed counterparts) can be in principal retained, and written to intermediate FCS files for backup or quality control outside of R. However, it is worth noting that this procedure could lead to a shortage of memory for large datasets, in which case overwriting the data at each step is advisable; if not specified otherwise, CATALYST overwrites by default.

A SCE can be constructed using CATALYST’s prepData() function, which accepts a path to a directory with one or many FCS files, a character vector of FCS filenames, a single or list of flowFrame(s), or a https://bioconductor.org/packages/3.11/bioc/html/flowCore.html (flowCore package²⁴). By default (transform = TRUE), an arcsinh-transformation with cofactor = 5 is applied to the input (count) data, and the resulting expression matrix is stored in the exprs assay slot of the output SCE:

# construct ’SingleCellExperiment’
fcs <- list.files("data", "acquisition", full.names = TRUE)
(sce <- prepData(fcs, transform = TRUE, cofactor = 5))

## class: SingleCellExperiment
## dim: 63 368152
## metadata(1): experiment_info
## assays(2): counts exprs
## rownames(63): 75As CD15 ... 208Pb CD45
## rowData names(4): channel_name marker_name marker_class use_channel
## colnames: NULL
## colData names(1): sample_id
## reducedDimNames(0):
## altExpNames(0):

Initially, our SCE has two assays containing dual ion counts (assay counts) and cofactor-5 arcsinh-transformed counts (assay exprs). The cofactor used for transformation is stored inside the object’s internal metadata (int_metadata(sce)$cofactor), and the FCS file of origin for each cell under cell metadata column sample_id (accessible via colData(sce)$sample_id or, equivalently, sce$sample_id). In our dataset, FCS files correspond to acquisitions rather than biological samples. Thus, we rename the cell metadata variable sample_id to file_id to avoid ambiguity:

i <- match("sample_id", names(colData(sce)))
names(colData(sce))[i] <- "file_id"

The total number of cells across all acquisitions corresponds to the number of columns in the SCE (ncol(sce)). We can summarize the number of cells in each file by tabulating the file_ids:

data.frame(
     file_id = levels(sce$file_id),
     n_cells = tabulate(sce$file_id))

##                   file_id n_cells
## 1 CyTOF_acquisition_1.FCS   48675
## 2 CyTOF_acquisition_2.FCS  125607
## 3 CyTOF_acquisition_3.FCS  193870

In both mass and flow cytometry, each feature has both a channel and target associated with it. As can be seen from printing the sce variable above, prepData() defaults to using targets as rownames (when available). We can retrieve each feature’s measurement channel using the channels() accessor, and use channel metals and masses to extract the indices of features that are relevant to different preprocessing steps. Namely, we assign channels measuring DNA to the variable dna, and channels for live gating (one DNA channel and cisplatin) to live:

# store character vector or channels names
chs <- channels(sce)
# store DNA & live channel indices
dna <- grep("^Ir", chs)
live <- grep("191|194", chs)

Filtering for stable signal

High quality data generation requires a stable signal throughout the acquisition. Various issues can lead to signal change over time, including unstable flow rate, introduction of air or introduction of metal contamination in the system. These changes in signal intensity can vary in terms of duration and intensity, and can affect all or only a subsets of channels simultaneously. In order to detect regions of the acquisition affected by signal instability, we display the signal for selected channels as a function of time in a scatter plot (Figure 1).

# plot channels of interest vs. time
coi <- chs[c(dna[1], which(rowData(sce)$use_channel))]
plotScatter(sce, chs = c("Time", coi), label = "both") +
     labs(y = "expression") +
     scale_x_continuous(
          expression("Time ("*10^6~"ms)"),
          labels = function(u) u/1e6) +
     theme_bw(base_size = 8) + theme(
          aspect.ratio = 2/3,
          panel.grid = element_blank(),
          axis.text = element_text(color = "black"),
          strip.background = element_rect(fill = NA))

Figure 1. Scatter plots of DNA channel Ir191 and the 41 channels measuring antigens against time.

Bins are colored by cell density; y-axis corresponds to cofactor-5 arcsinh-transformed dual counts.

In this particular experiment, we do not observe time-related signal instability. In case part of the acquisition should be excluded, this could be done by manually gating on the region with stable signal, and subsequent subsetting to only retain cells that fall within the gate’s boundaries (argument pop = "+"). Vice versa, it is possible to select a region with unstable signal, and remove it from the SCE object (pop = "-"). For the sake of completeness, we include how a region of unstable signal could be excluded via manual gating:

# construct ’GatingSet’
ff <- sce2fcs(sce[dna, ], assay = "exprs")
gs <- GatingSet(flowSet(ff))

# apply rectangular gate to exclude unstable signal
min_t <- ...
max_t <- ...
gs_add_gating_method(
    gs, alias = "stable",
     pop = "-", parent = "root",
     dims = paste0("Time,", chs[dna[1]]),
     gating_method = "boundary",
     gating_args = sprintf("min=c(%s,0),max=c(%s,10)", min_t, max_t))

# plot scatter of DNA vs. Time
ggcyto(gs,
     aes_string("Time", chs[dna[1]])) +
     geom_hex(bins = 128) +
     geom_gate("stable") +
     facet_null() + theme_bw() +
     ggtitle(NULL) + theme(
         legend.position = "none",
         panel.grid = element_blank())

# subset selected events
sce <- sce[, gh_pop_get_indices(gs, "stable")]

Normalization

In the case of mass cytometry, signal drift during acquisition due to a progressive loss of sensitivity must be accounted and normalized for. A widely established strategy is to mix samples with polystyrene beads embedded with metal lanthanides, allowing monitoring of instrument performance throughout data acquisition⁸. These beads are in turn used to estimate and correct for the signal’s time drift. When independent experiments have to be analyzed in the same context, variation due to changes in instrument performance over time combined with intervals between scheduled maintenance have to be taken into account as well. In this case, the bead signal should be normalized to a set of reference beads from an earlier acquisition. This ensures that different experiments are normalized to the same level, independent of the CyTOF’s sensitivity.

A MATLAB tool to perform normalization outside of R was available until recently at nolanlab/beadnormalization; current R implementations are available through CATALYST and premessa. CATALYST provides an extension of bead-based normalization as described by Finck et al.⁸, with automated identification of bead singlets (used for normalization), as well as of bead-bead and cell-bead doublets (to be removed), thus eliminating the need for manual gating. This is implemented as follows:

The above procedure is carried out by a single function, normCytof(), which takes as input a SCE and a set of arguments that control the normalization parameters and output format. Here, we specify beads = "dvs", corresponding to DVS Science beads (lanthanides Ce140, Eu151, Eu153, Ho165, Lu175). Secondly, we provide the path to a set of reference beads that are used to compute baseline intensities for normalization (argument norm_to). Lastly, we set overwrite = FALSE to retain both raw and normalized data, and remove_beads = TRUE to exclude bead and doublet events:

# specify path to reference beads
ref_beads <- file.path("data", "normalization_beads.fcs")
# apply bead-based normalization
system.time(res <- normCytof(sce,
     beads = "dvs", norm_to = ref_beads,
     remove_beads = TRUE, overwrite = FALSE))

##    user  system elapsed
##  18.567   1.155  19.930

When remove_beads = TRUE (the default), normCytof() will return a list of three SCEs containing filtered, bead and remove events, respectively, as well as two ggplot objects:

names(res)

## [1] "data"    "beads"   "removed" "scatter" "lines"

The first SCE (res$data) contains the filtered data with additional assay slot "normed" housing normalized expressions. The remaining two SCEs are data subsets that contain any events identified as beads (slot beads) and all removed events (including beads, bead-bead and bead-cell doublets; slot removed), respectively; thus, the beads themselves are a subset of the removed events. Here, we compare the number and percentage of cells contained in each subset:

# view no. of remaining, bead & removed events
ns <- sapply(res[1:3], ncol)
ps <- sprintf("%1.2f", ns/ncol(sce)*100)
data.frame(t(cbind("# events" = ns, "% of total" = ps)))

##              data beads removed
## # events   337525 27303   30627
## % of total  91.68  7.42    8.32

As a first quality control plot, res$scatter (Figure 2) renders scatter plots of bead channels (x-axis) versus DNA (y-axis), where events identified as beads as well as their expression range are highlighted in color; bead events should have low DNA intensity (since they are not cells) and high intensities across all bead channels.

Figure 2. Scatter plots for bead channels vs. DNA.

Events identified as beads are colored in blue; for each bead channel, expression ranges across all bead events are indicated as rectangular gates. Events were downsampled to at most 10,000 for visualization.

Secondly, res$lines (Figure 3) displays smoothed median bead intensities before and after normalization; these typically decrease with time prior to normalization, and should be approximately constant and centered around the baseline after normalization. In our dataset, normalization is performed based on previously acquired reference beads. Thus, baseline values correspond to the reference bead’s mean bead channel intensities. As shown in Figure 3, the bead channel levels are considerably lower after normalization, indicating higher sensitivity in the current run. Importantly, the slight decrease in signal over time is no longer present after normalization.

Figure 3. Running-median smoothed bead intensities vs. time before and after normalization; colored by bead channel.

In order to assess the sensitivity of the CyTOF during acquisition and identify potential issues that would have remained undetected during the tuning of the instrument, we compute the mean bead channel counts across events identified as beads (res$beads subset). A logical vector of which channels correspond to beads is stored under rowData column bead_ch, which we can use to subset the counts assay to include bead channels only.

# compute mean bead channel counts for current run
is_bead <- rowData(res$beads)$bead_ch   # get bead channels
bead_cs <- counts(res$beads)[is_bead, ] # subset counts
rownames(bead_cs) <- chs[is_bead]       # use channels as names
(bead_ms <- rowMeans(bead_cs))           # compute means

##  Ce140Di  Eu151Di  Eu153Di  Ho165Di  Lu175Di
## 2842.343 2111.072 2660.292 2537.831 2323.313

To assess the measurement sensitivity during the current run, we compare the mean bead channel counts computed above to those obtained from 12 previously acquired runs available in metadata table ref_bead_counts.csv. The resulting boxplot (Figure 4) shows that the current run’s sensitivity is relatively high, but well in the range of previous runs.

# read in reference mean bead channel counts
ref <- read.csv(file.path("data", "ref_bead_counts.csv"))

# join into single tidy data.frame
df <- bind_rows(ref, bead_ms, .id = "group")
df <- melt(df, id.var = "group")

# boxplot of reference vs. current run’s mean bead channel counts
ggplot(df, aes(variable, value)) +
    geom_boxplot(data = df[df$group == 1, ]) +
    geom_point(data = df[df$group == 2, ],
        col = "red", pch = 4, stroke = 1) +
    labs(x = "bead channel", y = "mean count") +
    qc_theme + ggtitle(
        "QC on bead channel counts",
        "[-] = reference | x = current run")

Figure 4. Bead channel count quality control.

Boxplot comparing the mean dual ion counts (y-axis) across bead channels (x-axis) obtained for the current run (red crosses) to those from 12 previsouly acquired reference runs (boxes).

After normalization, we overwrite the input dataset with the filtered subset that no longer includes bead events, or bead-bead and bead-cell doublets:

sce <- res$data

Debarcoding

In mass cytometry, samples are often labeled with unique sample-specific barcodes and pooled together for processing and measurement, an approach termed multiplexing²⁶. The most widely used barcoding scheme is based on Zunder et al.¹⁰, and relies on binary palladium-based mass-tag cell barcoding. Here, each sample i = 1, ..., n is either positive or negative for each of m palladium isotopes, resulting in an m-choose-k barcoding scheme, where k is the number of positive barcodes. For example, labeling of three out of six palladium isotopes will result in n-choose-k = 6-choose-3 = 20 unique barcodes. In order to recover the individual samples for further analysis, the pooled dataset is debarcoded (or deconvoluted) computationally.

The single cell debarcoding (SCD) algorithm first sorts each cell’s barcode intensities to assign preliminary barcode IDs such that a cell is assigned to the barcode population for which its barcode intensities are highest. Next, intensities within each barcode population are scaled using the 95th expression quantiles, and thereby brought to a comparable scale. Finally, events whose separation between highest negative and lowest positive barcode intensity is below a threshold value (separation cutoff) are left unassigned.

In the initial SCD algorithm, sample yields are determined by a single global cutoff on the separation between positive and negative barcode populations. Naturally, this procedure is suboptimal when yields as a function of the applied cutoff do not decline simultaneously. To optimize cell yields in such cases, CATALYST provides an option to automatically estimate or specify sample-specific separation cutoffs.

The SCD algorithm is implemented in CATALYST as a three-step procedure: i) preliminary barcode assignment (assignPrelim()); ii) automated estimation of sample-specific separation cutoffs (estCutoffs()); and, iii) application of cutoffs to arrive at final barcode assignments (applyCutoffs()).

Preliminary barcode assigment

For our dataset, a 6-choose-3 = 20 barcoding scheme was used (Figure 5). Five barcodes were unused (empty_1-5), resulting in 15 samples (9 references, 6 samples of interest). We first read the corresponding debarcoding_scheme.csv into R:

Figure 5. 6-choose-3 palladium isotope debarcoding scheme.

Rows correspond to palladium isotopes (barcode channels), columns to barcode identifiers (samples). Each sample is negative (white) or positive (grey) for 3 out of 6 barcode channels, resulting in 20 unique barcode combinations.

# read in debarcoding scheme
fn <- file.path("data", "debarcoding_scheme.csv")
bc_key <- read.csv(fn, row.names = 1, check.names = FALSE)

# all barcodes are positive for exactly 3 barcoding channels
all(rowSums(bc_key) == 3)

## [1] TRUE

During this first debarcoding step, each event is preliminarily assigned to a barcode according to its top-k expressed barcode channels. Here, events whose expression is highest for a combination of barcode channels that does not appear in the debarcoding scheme (bc_key) will be given barcode ID 0 (for “unassigned”). Thus, we can remove empty barcodes from the bc_key variable such that events assigned to these barcodes are left unassigned from the start. Alternatively, one could perform debarcoding using the non-subsetted key, and filter out empty barcodes downstream.

# remove empty barcodes from debarcoding scheme
is_empty <- grepl("empty", rownames(bc_key))
bc_key <- bc_key[!is_empty, ]
bc_ids <- rownames(bc_key)

For preliminary barcode assignment, we use CATLAYST’s assignPrelim() function, providing the input data (sce) and debarcoding scheme (bc_key). If not specified otherwise, assignPrelim() will default to using the exprs assay slot (argument assay). Because we ran normCytof() with overwrite = FALSE, this assay contains arcsinh-transformed raw counts; we set assay = "normexprs" in order to use the normalized values instead:

# do preliminary barcode assignments
system.time(sce <- assignPrelim(sce, bc_key, assay = "normexprs"))

##    user  system elapsed
##  15.150   0.394  15.685

In the returned SCE, feature metadata (rowData) column is_bc indicates whether or not a channel corresponds to a barcode channel:

# view barcode channels
channels(sce)[rowData(sce)$is_bc]

##      MCB1      MCB2      MCB3      MCB4      MCB5      MCB6
## "Pd102Di" "Pd104Di" "Pd105Di" "Pd106Di" "Pd108Di" "Pd110Di"

For each event, barcode identifiers are stored in colData column bc_id. After this preliminary round of assignment, 57980/337525 events (17.18%) have been left unassigned:

# tabulate number of (un)assigned events
table(sce$bc_id == 0)

##
##  FALSE   TRUE
## 279545  57980

Furthermore, for each cell, the barcode channel expressions are scaled relative to the 95th expression percentiles of its respective barcode population. The scaled data is stored in assay slot scaled. Based on these scaled barcode channel intensities, a separation value is computed as the distance between highest negative and lowest positive barcode channel; separations are stored in colData column delta.

Estimation of separation cutoffs

To decide on separation cutoffs, we consider yields upon debarcoding as a function of the applied cutoff (Figure 6). Commonly, this function will be characterized by an initial weak decline, where doublets are excluded, and subsequent rapid decline in yields to zero. In-between, low numbers of counts with intermediate barcode separation give rise to a plateau. Ideally, the applied separation cutoffs should provide a balance between high cell yield and low assignment uncertainty, marking the approximate midpoint of the yield function’s plateau region.

plotYields(sce, which = "0")

Figure 6. Yields plot for a 6-choose-3 debarcoding scheme.

Shown is the distribution of barcode population separations (histogram) and cell yields by sample (lines) as a function of the applied separation cutoff. Left axis corresponds to cell yield in percent; right axis to the total number of cells.

Instead of a single global cutoff, we estimate a sample-specific cutoff to account for barcode population yields that decline in an asynchronous fashion. To this end, we fit both a linear and a three-parameter log-logistic model to each yield function. For the linear fit, we estimate the cutoff as the value for which yields have declined to 50%. For the log-logistic fit, we compute the cutoff as the value for which there is minimal yield decline by minimizing each yield function’s 1st derivative. For each barcode, the final cutoff estimate is computed as the mean of both estimates, weighted with the goodness (residual sum of squares) of each fit (see Methods for details). Thus, the choice of thresholds for the distance between negative and positive barcode populations is: i) automated and ii) independent for each barcode. Nevertheless, reviewing barcode-specific yield plots and, in rare cases, refining the estimated separation cutoffs is advisable (see Figure 7).

Cutoff estimation is performed by CATALYST’s estCutoffs() function, which takes as input a SCE as returned by assignPrelim(); that is, preliminary barcode assignments are required to estimate separation cutoffs. estCutoffs() will store sample-specific cutoff estimates under metadata slot sep_cutoffs, but will leave barcode assignments unchanged.

sce <- estCutoffs(sce)
metadata(sce)$sep_cutoffs

## CellLine_R1 CellLine_R2 CellLine_R3 CellLine_R4     PBMC_R1     PBMC_R2
##  0.13829757  0.13688667  0.09161242  0.12436944  0.13041729  0.18048328
##     PBMC_R3    Tumor_R1    Tumor_R2     PBMC_S1     PBMC_S2     PBMC_S3
##  0.26517571  0.21012223  0.20543857  0.10440146  0.12902187  0.24858758
##    Tumor_S1    Tumor_S2    Tumor_S3
##  0.18444497  0.14691641  0.20818002

We can visually inspect the estimated cutoffs using plotYields() with argument which specifying the barcode ID of interest (Figure 7). In our example, the cutoff estimate nicely marks the midpoint of the yield function’s plateau or, equivalently, the valley between peaks of cell yields. To decrease the stringency of the applied cutoff, and thus increase the resulting cell yield, we could set the sample’s cutoff to e.g. 0.1. Vice versa, a more stringent cutoff of e.g. 0.2 would decrease the cell yield but yield a purer population.

As an alternative to inspecting the cutoff estimate for each sample in R, we could specify which = bc_ids to obtain a list of yield plots for all barcodes; the generated set of plots may be written to a single PDF file via providing plotYields() with an out_path to allow for easy reviewing of the separating cutoffs currently stored within the object.

plotYields(sce, which = "PBMC_R1")

Figure 7. Yields plot for an exemplary sample, including the estimated separation cutoff.

Shown is the distribution of barcode population separations (histogram) and cell yields (line) for the sample as a function of the applied separation cutoff. Left axis corresponds to cell yield in percent; right axis to the total number of cells.

Besides a cutoff on the separation between positive and negative barcode populations, to trim outliers, the SCD algorithms applies an additional cutoff on the Mahalanobis distance (argument mhl_cutoff), a metric that quantifies the distance of a given event to the expression distribution of the barcode population it has been assigned to.

In Figure 6, we can observe that population yields decline synchronously with increasing separation cutoffs, and that we might consider applying a global separation cutoff, e.g., at ∼ 0.15. For this data, yields are in fact similar, independent of whether we apply sample-specific cutoffs or a single global one. Nevertheless, applying sample-specific cutoffs is recommended in order to maximize cell yields while minimizing uncertainty in barcode assignments.

# store preliminary barcode IDs
bc_ids0 <- sce$bc_id

# apply global & sample-specific separation cutoff(s)
sce_glob <- applyCutoffs(sce, sep_cutoffs = 0.15, mhl_cutoff = 30)
sce_spec <- applyCutoffs(sce, mhl_cutoff = 30)

# compare cell yields for both cutoff strategies
c(global = mean(sce_glob$bc_id == 0),
  specific = mean(sce_spec$bc_id == 0))

##    global  specific
## 0.3573839 0.3584949

After debarcoding, we compare the number of events assigned to each barcode population before and after applying separation cutoffs, and filter out events that have been left unassigned (barcode ID 0). As shown in Figure 8, after applying the separation cutoffs, the number of unassigned cells (0) increases, while the number of cells assigned to each barcoding well decreases. We also observe a higher decrease in assigned cells for tumor samples, which underwent a dissociation process and contain more debris. Conversely, highly viable cell lines and PBMCs have a higher recovery yield.

# proceed with sample-specific filtering
sce <- sce_spec

# compute number of events per population
# before vs. after applying separation cutoffs
barplot(rbind(table(bc_ids0), table(sce$bc_id)),
     beside = TRUE, ylab = "cell count",
     las = 2, cex.axis = 0.5, cex.names = 0.5)
legend("topright", fill = c("black", "grey"),
     legend = c("before filtering", "after filtering"))

# remove unassigned events
sce <- sce[, sce$bc_id != 0]

Figure 8. Barplot of cell counts before (black) and after (grey) apply separation cutoffs.

Compensation

Mass cytometry utilizes heavy metals (usually from the lanthanide series) as reporters to label antibodies. As a result, channel crosstalk originating from spectral overlap and autofluorescence is significantly less pronounced in mass cytometry compared to flow cytometry. Yet, spillover due to abundance sensitivity, isotopic impurities, and oxide formation still exists, giving rise to artefactual signal that can impede data interpretability.

A combined experimental-computational pipeline to correct for spillover in mass cytometry data has been proposed by Chevrier et al.⁹ and is implemented in the CATALYST package. In brief, compensation is achieved via the following three-step approach outlined here (see for details).

We will apply a pre-acquired spillover matrix (metadata file spillover_matrix.csv). Thus, we enter at step 3, which involves only compensating the input dataset using CATALYST’s compCytof() function. By default, compCytof() will reuse the cofactor stored in int_metadata(sce)$cofactor for computing arcsinh-transformed data from the compensated counts, thus applying the same transformation as during data preparation and normalization:

# read in pre-computed spillover matrix
sm <- file.path("data", "spillover_matrix.csv")
sm <- read.csv(sm, row.names = 1)
# apply NNLS compensation
system.time(
     sce <- compCytof(sce, sm, method = "nnls",
          assay = "normcounts", overwrite = FALSE))

##    user  system elapsed
##  67.525   7.318  75.337

To visually inspect how compensation affects signal intensities, we can generate scatter plots pre- and post-compensation; an exemplary pair of channels is shown in Figure 9. In such a plot, we can observe a slight positive association between the signals of spill-affected channels, which should be removed upon compensation.

i <- grep("173|174", chs, value = TRUE)
p1 <- plotScatter(sce,
  chs = i,
  label = "channel",
  assay = "normexprs") +
  ggtitle("Uncompensated")
p2 <- plotScatter(sce,
  chs = i,
  label = "channel",
  assay = "compexprs") +
  ggtitle("Compensated") +
  ylab(NULL)
wrap_plots(p1, p2)

Figure 9. Scatter plots for two exemplary channels before (left) and after correction for spillover (right).

Gating

Many events acquired in mass cytometry may in fact be debris, doublets or dead cells, and should be filtered out through a gating step. Here, we suggest a strategy that first applies an elliptical gate on cell events, defined as double positive for the DNA channels Ir191/Ir193. This allows the exclusion of debris and doublets. As a second step, we discard cells that are positive for the dead cell marker Pt194.

These two steps are performed using the openCyto R package²⁰, and the resulting gates are visualized on scatter plots of the channels subjected to gating using ggcyto²³. For consistent visualization, we again define a common plotting theme for scatter plots of channels chs that include the gating boundaries for the specified gate_id:

.scatter <- function(gs, chs, gate_id = NULL,
    subset = ifelse(is.null(gate_id), "root", "_parent_")) {
    p <- ggcyto(gs, max_nrow_to_plot = 1e5,
        aes_string(chs[1], chs[2]), subset) +
        geom_hex(bins = 100) + facet_wrap(~ name, ncol = 5) +
        (if (is.null(gate_id)) list() else geom_gate(gate_id)) +
        ggtitle(NULL) + theme_bw(base_size = 8) + theme(
             aspect.ratio = 1,
             legend.position = "none",
             panel.grid.minor = element_blank(),
             strip.background = element_rect(fill = NA),
             axis.text = element_text(color = "black"),
             axis.text.x = element_text(angle = 45, hjust = 1))
    suppressMessages(p + coord_equal(expand = FALSE,
        xlim = c(-1, 11), ylim = c(-1, 11)))
}

Gating on cells

In order to apply sample-specific gates, we first convert the SCE into a flowSet with a separate frame for each sample (argument split_by = "bc_id"). As gating should be performed on expression-like data (not ion counts), we further specify assay = "exprs" to retain the arcsinh-transformed assay slot. Thirdly, since conversion from SCE to flowCore data structures requires matrix transposition (rows correspond to targets in the SCE, but to events in flowFrame/Sets), we retain only those channels that are relevant when gating of (live) cells: DNA and dead channels, whose indices are stored in variables dna and live.

# subset DNA & live channels
sub <- sce[union(dna, live), ]

# add metadata variable ’i’ to track cell indices
colData(sub) <- DataFrame(
     bc_id = sub$bc_id,
     i = seq_len(ncol(sce)))

# split SCE by sample
fs <- sce2fcs(sub,
     assay = "compexprs",
     split_by = "bc_id",
     keep_cd = TRUE)

# construct ’GatingSet’
gs <- GatingSet(fs)

We apply an elliptical gate (gating_method = "flowClust.2d") to exclude the two lowest density percentiles (quantile = 0.98). Because the input gating set contains a separate frame for each barcode, the gate will be computed separately for each sample.

# apply elliptical gate on DNA channels
gs_add_gating_method(gs,
     alias = "cells",
     pop = "+", parent = "root",
     dims = paste(chs[dna], collapse = ","),
     gating_method = "flowClust.2d",
     gating_args = "K=1,quantile=0.98,target=c(5,5)")

We use ggcyto to produce scatter plots of the DNA channels, with geom_gate("cells") to visualize the gates computed above (Figure 10):

# plot scatter of DNA channels split by sample
.scatter(gs, chs[dna], "cells")

Figure 10. Scatter plots of DNA channels, split by sample and including elliptical cell gates.

Gating on live cells

The wrapper function .live_gate() defines a polygonal gate comprised of a line and a bivariate standard normal density Z, such that cells pass gating when i) their expression is within the qth quantile of Z; and, ii) their expression falls below a line parameterized by intercept i and slope s. In this way, the gate is centered around the expression peak, while excluding cells whose dead channel intensities increases with DNA content.

# define live cell gate plug-in
# x = expression matrix, q = quantile, i = intercept, s = slope
.live_gate <- function(x, q = 0.99, i = 1, s = 0.5) {
    # specifying gating function
    .gating_fun <- function(fr, pp_res, channels = NA, id = "", ...) {
         # subset channels of interest
         x <- exprs(fr[, channels])
         # scale data for comparison w/ ’qnorm()’
         x0 <- scale(x)
         # set boundary level as q-th quantile of standard normal
         z <- qnorm(q)
         # find p(x) for that level
         pd <- dmvnorm(c(z, z))[1]
         px <- dmvnorm(x0)
         # find points above boundary level
         keep1 <- px > pd
         # find points below line y = a + b * x
         keep2 <- (i + s * x0[, 1]) > x0[, 2]
         # intersection of points below line & above threshold level
         pts <- x[keep1 & keep2, ]
         # get boundary points (convex hull)
         pts <- pts[chull(pts), ]
         # return gate
         polygonGate(.gate = pts, filterId = id)
    }
    # register gate
    suppressMessages(
      foo <- register_plugins(
        fun = .gating_fun,
        methodName = "liveGate",
        dep = "mvtnorm",
        "gating"))
}

In contrast to the cell gates above, we apply live gates with sample-specific gating parameters. To this end, we specify a list l containing quantiles q, intercepts i and slopes s for each sample. These parameters are updated iteratively to remove dead cells while retaining cell yields as high as possible. After manual adjustments, we arrive at the following sample-specific gating parameters (Figure 11).

# set default parameters for all samples
l <- lapply(c(q = 0.99, i = 0.9, s = 0.4), function(u)
    setNames(rep(u, length(gs)), sampleNames(gs)))

# adjust parameters for specific samples
l$i[["PBMC_R2"]] <- 1.2
l$i[["PBMC_R3"]] <- 1.2
l$i[["PBMC_S1"]] <- 1.2
l$s[["PBMC_S2"]] <- 0.2
l$i[["PBMC_S2"]] <- 0.6
l$i[["PBMC_S3"]] <- 1.8
l$s[["Tumor_S2"]] <- 0.3
l$i[["Tumor_S2"]] <- 0.6
l$s[["Tumor_S3"]] <- 0.3
l$i[["Tumor_S3"]] <- 0.4

for (i in sampleNames(gs)) {
    # register & apply live gate with sample-specific parameters
    .live_gate(x, q = l$q[i], i = l$i[i], s = l$s[i])
    gs_add_gating_method(gs[i],
         alias = "live",
         pop = "+",
         parent = "cells",
         dims = paste(chs[live], collapse = ","),
         gating_method = "liveGate")
}
.scatter(gs, chs[live], "live")

Figure 11. Scatter plots of DNA and dead cell channels, split by sample and including the dead cell polygon gates.

We display the yield of "cell" and "live" gates on each samples to quickly assess the cell losses occurring at the two gating steps (Figure 12). As expected the "cell" gate leads to a systematic loss of around 1% of cells across all the samples. The "live" gate leads to a stronger reduction of cell yield in the tumor samples, consistent with the fact that those samples, which underwent enzymatic dissociation, contain more dead cells.

# extract gating frequencies
df <- gs_pop_get_stats(gs,
    type = "percent",
    nodes = c("cells", "live"))
df <- rename(df, gate_id = "pop")

# barplot of cell yields after cell/live gating
ggplot(df, aes(sample, percent, fill = gate_id)) +
     geom_bar(width = 2/3, stat = "identity", position = "dodge") +
     scale_x_discrete(limits = bc_ids, expand = c(0, 2/3)) +
     scale_y_continuous(labels = seq(0, 100, 25),
          limits = c(0, 1), expand = c(0, 0)) +
     labs(y = "cell yield (%)") + qc_theme

Figure 12. Barplot of cell and live gating yields.

For each barcode ID (x-axis), frequencies are relative to the total number of cells in the population before gating; bars are colored by gate ID.

We extract a logical vector indicating whether a given event is included in or excluded by the "live" gate applied above by applying gh_pop_get_indices to each sample in gs. Secondly, we extract the cell indices from gs and subset the SCE to keep only cells that passed the "live" gate.

fs <- gs_pop_get_data(gs, "live") # get data from ’GatingSet’
es <- lapply(fs, exprs)            # get expression matrices
es <- do.call("rbind", es)         # join into single data.frame
sce <- sce[, es[, "i"]]            # subset retained cells

Finally, we can again visualize scatter plots of dead channels against DNA as a quality control for the retained subset of cells (Figure 13).

Figure 13. Scatter plots of dead cell channels against DNA, including the subset of cells remaining after live cell gating.

Quality control

Having completed the standard preprocessing steps, we proceed to investigate how the current run compares to prior acquisitions in terms of the number of cells in each reference and sample, and the expression levels of each target. Large parts of the metadata generated by now may no longer be needed, and unnecessarily increases output file sizes for large-scale datasets. Therefore, we will retain only two key cell metadata variables: file_ids containing the FCS filename each cell originates from, and bc_ids containing the barcode population assignments. We secondly rename these variable to make the following quality control steps more intuitive.

# drop all cell metadata except file of origin & barcode IDs
colData(sce) <- colData(sce)[c("file_id", "bc_id")]

# rename cell metadata variable
i <- match("bc_id", names(colData(sce)))
names(colData(sce))[i] <- "sample"

In the debarcoding scheme used for deconvolution of the multiplexed samples (Section Debarcoding), barcode identifiers were chosen to contain all information relevant for each sample. This setup allows us to extract sample metadata directly from the bc_ids. Alternatively, and especially for more complex experimental designs, this information could be stored in a separate metadata table. Such a table could then be used to match the bc_ids with the listed samples, and add arbitrary metadata information (e.g., batch, patient ID, treatment).

In our example, barcode identifiers include each sample’s type (CellLine, PBMC or Tumor), group (R for reference or S for sample of interest), and replicate number; and follow a consistent naming scheme: . We can easily extract these components and store them in the SCE’s cell metadata (colData):

sce$type <- gsub("_.*", "", sce$sample)
sce$group <- gsub("[^R|S]", "", sce$sample)
head(colData(sce))

## DataFrame with 6 rows and 4 columns
##                   file_id      sample        type       group
##                  <factor> <character> <character> <character>
## 1 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R
## 2 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R
## 3 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R
## 4 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R
## 5 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R
## 6 CyTOF_acquisition_1.FCS CellLine_R1    CellLine           R

Quality control (QC) on reference cell counts

As a first quality control, we compare the cell counts of each reference sample (R) to those obtained from 7 previous runs (Figure 14). Since the references are obtained from pre-barcoded aliquots of cells, the number of reference cells acquired gives direct information regarding the cell yield throughout the whole experiment: From cell barcoding to acquisition on the CyTOF. As shown in Figure 10, the current run has a lower yield compared to average run. This observation is consistent across all references measured.

# boxplot of current vs. reference cell counts
ref <- read.csv(file.path("data", "ref_cell_counts.csv"))
run <- table(sce$sample[sce$group == "R"])

# join into single tidy data.frame
df <- bind_rows(ref, run, .id = "group")
df <- melt(df, id.var = "group")

ggplot(df, aes(variable, value)) +
     geom_boxplot(data = df[df$group == 1, ]) +
     geom_point(data = df[df$group == 2, ],
         col = "red", pch = 4, stroke = 1) +
     labs(x = "sample", y = "cell count") +
     qc_theme + ggtitle(
         "QC on reference cell counts",
         "[-] = reference | x = current run")

Figure 14. Reference cell count quality control.

Boxplot comparing the reference cell counts obtained for the current run (red crosses) to those from 7 previsouly acquired runs.

QC on sample cell counts

Secondly, we compare the cell counts for the 4 samples of interest (2 PBMC, 2 tumor samples) to the number of cells recorded for 14 tumor and PBMC samples each (28 reference samples in total) acquired in previous runs (Figure 15). This step provides a first quality assessment of the samples of interest. Here, samples with too few cells will be less reliable, and potentially less representative of the original tissue, making conclusions from downstream analyses more difficult to draw.

ref <- read.csv(file.path("data", "sample_cell_counts.csv"))
run <- table(sce$sample[sce$group == "S"], dnn = "sample")
run <- as.data.frame(run, responseName = "count")
run$type <- sce$type[match(run$sample, sce$sample)]
df <- bind_rows(ref, run, .id = "group")

ggplot(df, aes(type, count)) +
     geom_boxplot(data = df[df$group == 1, ]) +
     geom_point(data = df[df$group == 2, ],
         col = "red", pch = 4, stroke = 1) +
     labs(x = "type", y = "cell count") +
     qc_theme + ggtitle(
         "QC on sample cell counts",
        "[-] = reference | x = current run")

Figure 15. Sample cell count quality control.

Boxplot comparing the sample cell counts obtained for the current run (red crosses) to those from 7 previsouly acquired runs.

QC on mean marker intensities

As the third and final quality control, we compare the 98th expression quantiles across all targets of interest over the pooled references to those obtained from 8 previously acquired runs available in metadata table ref_marker_levels.csv (Figure 15). We chose the 98th percentile to account for the fact that some populations are rare, and we are particularly interested in assessing signal stability for positive cells rather than the median of the population. Since the pooled references are identical from one run to another, this gives a direct indication of the current run’s staining efficacy and enables early identification of antibody degradation over time.

# read in reference data
ref <- file.path("data", "ref_marker_levels.csv")
ref <- read.csv(ref, check.names = FALSE)

# compute 98th expression quantiles
# for reference samples in current run
es <- assay(sce, "compexprs")
es <- es[names(ref), sce$group == "R"]
run <- rowQuantiles(es, probs = 0.98)

# join into single tidy data.frame
df <- bind_rows(ref, run, .id = "group")
df <- melt(df, id.var = "group")

ggplot(df, aes(variable, value)) +
     geom_boxplot(data = df[df$group == 1, ]) +
     geom_point(data = df[df$group == 2, ],
         col = "red", stroke = 0.5) +
     labs(x = "target", y = "98th expression quantile") +
     qc_theme + ggtitle(
         "QC on marker levels",
         "[-] = referece | o = current run")

Figure 16. Mean marker expression quality control.

Boxplot comparing the mean marker expression obtained for the current run (shown in red) to those from 8 previsouly acquired runs.

Batch alignment

Each CyTOF run contains the same set of references. Similar to the approach used by Schuyler et al.¹¹, we use these references as anchors to calculate a channel-specific correction factor by dividing the 98th percentile measured in the current run by the average 98th percentile obtained across the first seven runs of the project. The signal observed in each channel for the samples of interest is then divided by these correction factors derived from the reference samples (Figure 16).

# compute 98th count quantiles via back-transformation
# (using same cofactor as always) & average across replicates
cf <- int_metadata(sce)$cofactor
ref <- data.frame(
  target = colnames(ref),
  mean_count = colMeans(sinh(ref)*cf))

# initialize correction factor of 1 for all channels
cfs <- setNames(rep(1, nrow(sce)), rownames(sce))

# compute batch correction factors for relevant channels
cs <- assay(sce, "compcounts")
csR <- cs[ref$target, sce$group == "R"]
run <- rowQuantiles(csR, probs = 0.98)
cfs[ref$target] <- run / ref$mean_count

# apply marker-specific batch correction (bc)
cs <- sweep(cs, 1, cfs, "/")
assay(sce, "bccounts") <- cs

# apply arcsinh-transformation
assay(sce, "bcexprs") <- asinh(cs/cf)

Normalization

Identification of bead events. Commonly, bead events are identified by manual gating on scatter plots of DNA vs. bead channels where DNA show be low, and expression should be high across all bead channels. Instead, we propose a programmatic way to identify beads that includes detection of bead-bead and cell-bead doublets.

Our normalization strategy leverages the already established SCD algorithm for preliminary tagging of events as beads. In this context, the debarcoding scheme is a 2×(2+m) matrix (Figure 17). Here, columns correspond to the two DNA channels and m barcode channels; rows correspond to barcodes 0 (no bead) and 1 (is bead), where non-bead events are positive for DNA channels only (barcode 11000. . . ), while bead events are negative for DNA and positive for all bead channels (barcode 00111. . . ):

Figure 17. Schematic of the debarcoding scheme used by ‘CATALYST‘’s ‘normCytof()‘ function to identify bead events.

Rows correspond to barcodes, columns to DNA and bead channels. Each barcode is either positive (grey) or negative (blank) for a given channel; cells (barcode 11000...) are positive for DNA and negative for bead channels, bead events (barcode 00111...) are negative for DNA and positive for bead channels.

Upon initially assignment of bead events, we apply a median ± n median absolute deviation (MAD) rule to remove low- and high-signal events from the bead population used for estimating normalization factors. As n decreases, bead populations become more narrow and bead-bead doublets are excluded. The extent to which bead populations are trimmed can be adjusted via argument trim (default 5).

Notably, slight over-trimming does not affect normalization. It is therefore recommended to choose a trim value that is small enough to assure removal of doublets at the cost of reduced bead population sizes.

Correcting for signal-decrease over time. To correct for the effect of event acquisition time on signal intensity, we follow the method proposed by Finck et al.⁸. In essence, bead intensities are smoothed using a median sliding-window with width k (default 500 bead events). At each timepoint, the slope of a line with intercept zero is computed my minimizing the squared error between smoothed bead and mean bead intensities (= baseline). Alternatively, a reference set of beads from which to compute the baseline can be provided. Slopes for non-bead timepoints are obtained via constant interpolation of these slopes. Here, large slopes correspond to significant deviation from the baseline, and small slopes indicate that the signal is already similar to the baseline. Thus, raw bead counts are normalized by multiplication with the fitted slopes at each timepoint.

Debarcoding

Preliminary barcode assignment. The debarcoding process commences by assigning each event a preliminary barcode ID. This requires specification of a binary barcoding scheme (or debarcoding key)

where i = 1, ..., n is the barcode index, j = 1, ..., m a barcode channel, and n, m denote the number of unique barcodes and barcoding channels, respectively. Further, let k_i denote the number of positive barcoding channels for barcode i: $k_i={\displaystyle {\sum }_{j=1}^m{b}_{ij}.}$

If k_i = k ∀ i = 1, ..., n (i.e., every barcode has the same number of positive barcoding channels), the k channels with the highest signal in a given event are considered to be positive, the remaining m − k to be negative. The separation δ of positive and negative events is then defined as the difference between the kth highest and (m − k)th lowest scaled intensity for that event.

Seperation cutoff estimation. When the separation between positive and negative barcoding channels is low, we cannot be confident in the barcode assignment.

For the estimation of cutoff parameters, we consider yields upon debarcoding as a function of the applied cutoffs. Commonly, this function will be characterized by an initial weak decline, where doublets are excluded, and subsequent rapid decline in yields to zero. In between, low numbers of counts with intermediate barcode separation give rise to a plateau. To facilitate robust estimation, we fit a linear and a three-parameter log-logistic function³² to the yields function with drc’s LL.R function³³ (Figure 18). As an adequate cutoff estimate, we target a point that marks the end of the plateau regime and on-set of yield decline to appropriately balance confidence in barcode assignment and cell yield.

We define the linear model cutoff estimate c_LM as the value for which the cell yield Y has declined to half of the initial Yield β₀:

where β₀, β₁ are the intercept and slope of the linear model fit, respectively.

We define the log-logistic model cutoff estimate c_LLM as the value for which the log-logistic function’s decline is minimized relative to its value:

The final cutoff estimate c is defined as the weighted mean between these estimates:

where w is the goodness of the linear fit relative to the log-logistic fit:

Figure 18. Schematic description of automated separation cutoff estimation.

Bar graphs represent the distribution of cells relative to the barcode distance, dotted line scorresponds to yield upon debarcoding as a function of the applied separation cutoff. The yield curve is fitted with a linear regression (blue) and a three parameter log-logistic function (red). The final cutoff estimate (black dashed line) is defined as the mean of estimates derived from both fits, weighted with the goodness of the respective fit.

Compensation

Retrieval of real signal. As in conventional flow cytometry, we model spillover linearly, with the channel stained for as predictor, and spill-effected channels as response. Thus, the intensity observed in a given channel j are a linear combination of its real signal and contributions of other channels that spill into it. Let I denote the (unknown) real and J the observed signal. Further, let s_ij be the proportion of channel j signal that is due to channel i, and w_j the set of channels that spill into channel j. Then

In matrix notation, measurement intensities may be viewed as the convolution of real intensities and a spillover matrix $SM=(s_{ij})\in {ℝ}_{+}^{n\times p},$ where n denotes the number of samples (cells) and p the number of features (channels): J = I · SM. The real signal I can then be retrieved via:

where SM⁻¹ is termed compensation matrix (CM).

While mathematically exact, the solution to this equation will yield negative values, and does not account for the fact that ion counts are strictly non-negative. A computationally efficient way to adress this is to instead use non-negative linear least squares (NNLS), which optimizes the least squares criterion under the constraint of non-negativity:

To solve for I, we apply the Lawson-Hanson algorithm^34,35 for NNLS implemented in the nnls package.

Spillover estimation. Because any signal not in a single stain experiment’s primary channel j results from channel crosstalk, each spill entry s_ij can be approximated by the slope of a linear regression with channel j signal as the response, and channel i signals as the predictors, where i ∈ w_j. computeSpillmat() offers two alternative ways for spillover estimation (Figure 19).

The default method approximates this slope with the following single-cell derived estimate: Let i⁺ denote the set of cells that are possitive in channel i, and $s_{ij}^c$ be the channel i to j spill computed for a cell c that has been assigned to this population. We approximate $s_{ij}^c$ as the ratio between the signal in unstained spillover receiving and stained spillover emitting channel, I_j and I_i, respectively. The expected background in these channels, $m_j^{-}$ and $m_i^{-}$, is computed as the median signal of events that are i) negative in the channels for which spill is estimated (i and j); ii) not assigned to potentionally interacting channels; and, iii) not unassigned, and subtracted from all measurements:

Each entry s_ij in SM is then computed as the median spillover across all cells c ∈ i⁺:

In a population-based fashion, as done in conventional flow cytometry, s_ij is computed as the slope of a line through the medians (or trimmed means) of stained and unstained populations, $m_j^{+}$ and $m_i^{+}$, respectively. Background signal is computed as above and subtracted, according to:

Figure 19. Population versus single-cell based spillover estimation.

In a population-based setting (left), spillover is estimated as the slope of a line through the centers of positive (red) and negative (blue) populations. In a single-cell based setting (right), slopes are computed independently for each cell in the positive population, and spillover is estimated as their median.

On the basis of their additive nature, spill values are estimated independently for every pair of interacting channels. Hereby, we take into account only interactions that are sensible from a chemical and physical point of view: M ± 1 channels (abundance sensitivity), M + 16 channels (oxide formation), and channels measuring isotopes (impurities; Figure 20).

Figure 20. Heatmap of channel interactions expected to exhibit spillover.

Included are only interactions that are sensible from a chemical and physical point of view: adjacent mass channels (abundance sensitivity), +16 mass channels (oxidation), and channels measuring isotopes (impurities).

Alternatively, interactions = "all" will compute a spill estimate for all n · (n − 1) possible interactions, where n denotes the number of measurement parameters. Estimates falling below the threshold specified by th will be set to zero. Lastly, note that diagonal entries s_ii = 1 for all i ∈ 1, ..., n, so that spill is relative to the total signal measured in a given channel.

Underlying data

The CyTOF data as well as all metadata required to run the full pipeline presented herein are available from Figshare as well as the Tumor Profiler website at https://tu-pro.ch/download/catalyst.

Figshare: An R-based reproducible and user-friendly preprocessing pipeline for CyTOF data. https://doi.org/10.6084/m9.figshare.c.5063984.v1¹⁸

This project contains the following underlying data:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).