Cyto Batch Integration

CytofBatchadjust corrects batch effects across cytometry data by aligning signal intensity peaks for each channel across batches. The algorithm uses technical replicates included in each barcode set as reference points to anchor each batch. Adjustments factors are calibrated using anchor samples representing each barcode set, then applied to all samples composing a batch.

This implementation uses samples from each group as control samples.

CytofBatchadjust corrects batch effects across cytometry data by aligning signal intensity peaks for each channel across batches. The algorithm uses technical replicates included in each barcode set as reference points to anchor each batch. Adjustments factors are calibrated using anchor samples representing each barcode set, then applied to all samples composing a batch.

This implementation uses samples only from one group as control samples.

Corrects for batch effects by fitting linear models, gains statistical power via an EB framework where information is borrowed across genes. This uses the implementation combat.py

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine using all control samples (replicates, in cyCombine terms), with a square SOM grid and correct the batches to batch 1. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine using all control samples (replicates, in cyCombine terms), with a square SOM grid and correct the batches to a midpoint derived from all batches. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine without any control samples (replicates, in cyCombine terms), with a square SOM grid and correct the batches to batch 1. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine without any control samples (replicates, in cyCombine terms), with a square SOM grid and correct the batches to a midpoint derived from all batches. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine using control samples only from one biological group (replicates, in cyCombine terms), with a square SOM grid and correct the batches to batch 1. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

cyCombine perform batch integration using self-organizing maps and ComBat. It first uses self-organizing maps (SOM) to group similar cells into clusters, then applies a ComBat-based method to correct batch effects within each cluster. Here, we run cyCombine using control samples only from one biological group (replicates, in cyCombine terms), with a square SOM grid and correct the batches to a midpoint derived from all batches. The size of the SOM grid is varied linearly between value of 6-16, with default set to 8 following the default value in create_som function in cyCombine. Rlen, which is the number of times data is presented to the SOM network and can impact clustering quality. It is varied linearly between value of 10-20 with default set to 10.

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

Here, we run CytoNorm using all control samples available, aligning the batches to batch 1.

The parameter nQ, which specifies the number of quantiles used when computing the splines is varied linearly between value of 80-120, with default set to 99 following the default value provided by CytoNorm. Clustering was performed by FlowSOM. The number of cells clustered by FlowSOM is set to be number of cells in the smallest control samples or 1,000,000, whichever is the smaller, multiplied by how many control samples there are in the data. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm.

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

Here, we run CytoNorm using all control samples available, aligning to a midpoint derived from all batches.

The parameter nQ, which specifies the number of quantiles used when computing the splines is varied linearly between value of 80-120, with default set to 99 following the default value provided by CytoNorm. Clustering was performed by FlowSOM. The number of cells clustered by FlowSOM is set to be number of cells in the smallest control samples or 1,000,000, whichever is the smaller, multiplied by how many control samples there are in the data. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm.

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

In this CytoNorm version, an aggregate of each batch is created and subsequently used as a proxy for the control samples, aligning the batches to batch 1.

The number of cells used to create an aggregate is set as the number of cells in the smallest sample or 1,000,000, whichever is the smaller, multiplied by how many samples there are in the batch. Clustering was performed by FlowSOM. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm. The number of cells clustered by FlowSOM is set as the number of cells in the smallest aggregate or 1,000,000, whichever is the smaller, multiplied by how many batches there are in the data (because 1 aggregate = 1 batch).

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

In this CytoNorm version, an aggregate of each batch is created and subsequently used as a proxy for the control samples, aligning to a midpoint derived from all batches.

The number of cells used to create an aggregate is set as the number of cells in the smallest sample or 1,000,000, whichever is the smaller, multiplied by how many samples there are in the batch. Clustering was performed by FlowSOM. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm. The number of cells clustered by FlowSOM is set as the number of cells in the smallest aggregate or 1,000,000, whichever is the smaller, multiplied by how many batches there are in the data (because 1 aggregate = 1 batch).

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

Here, we run CytoNorm using only control samples from one donor, aligning the batches to batch 1.

The parameter nQ, which specifies the number of quantiles used when computing the splines is varied linearly between value of 80-120, with default set to 99 following the default value provided by CytoNorm. Clustering was performed by FlowSOM. The number of cells clustered by FlowSOM is set to be number of cells in the smallest control samples or 1,000,000, whichever is the smaller, multiplied by how many control samples there are in the data. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm.

CytoNorm corrects batch effects by using reference control samples (aliquots of one sample, technical replicates) included with each batch. It clusters cells using FlowSOM, then trains a model on the control samples to learn how marker expression distributions differ across batches for each population. It then uses splines to align these distributions to a common reference (either a midpoint derived from all batches or to a single batch).

Here, we run CytoNorm using only control samples from one donor, aligning to a midpoint derived from all batches.

The parameter nQ, which specifies the number of quantiles used when computing the splines is varied linearly between value of 80-120, with default set to 99 following the default value provided by CytoNorm. Clustering was performed by FlowSOM. The number of cells clustered by FlowSOM is set to be number of cells in the smallest control samples or 1,000,000, whichever is the smaller, multiplied by how many control samples there are in the data. The size of the SOM grid is varied linearly between value of 6-16, with default set to 15 following the default value provided by CytoNorm. The number of metaclusters is varied linearly between value of 8-20, with default set to 10 following the default value provided by CytoNorm.

CytoVI is a deep generative model that utilizes antibody-based single-cell profiles to learn a biologically meaningful latent representation of each cell. It is part of the scvi-tools framework and is built upon the variational autoencoder (VAE) architecture.

This method batch-normalizes a set of cytometry data samples by identifying and aligning the high density regions (landmarks or peaks) for each channel. The data of each channel is shifted in such a way that the identified high density regions are moved to fixed locations called base landmarks. Normalization is achieved in three phases:

1. identifying high-density regions (landmarks) for each flowFrame in the flowSet for a single channel
2. computing the best matching between the landmarks and a set of fixed reference landmarks for each channel called base landmarks
3. manipulating the data of each channel in such a way that each landmark is moved to its matching base landmark. Please note that this normalization is on a channel-by-channel basis

NOTE: The default implementation uses max.lms=2, although for some channels it is not possible to compute 2 landmarks, resulting in an error. In order to fully automate the batch normalization process, this implementation checks whether it is possible to compute 2 landmarks, and if not, it sets max.lms=1 for that channel.

Harmony is a general-purpose R package with an efficient algorithm for integrating multiple data sets. It is especially useful for large single-cell datasets such as single-cell RNA-seq.

Limma removeBatchEffect is a method that uses a linear model and matrix decomposition to remove batch effects from a dataset. It first fits a linear model to the data, then decomposes the model matrix into a set of orthogonal components. The batch effect is then removed by subtracting the component corresponding to the batch effect from the data.

Seurat RPCA performs batch integration by projecting each query dataset into the PCA space of a goal batch, and identifying anchors using reciprocal PCA (RPCA). RPCA identifies mutual nearest neighbors between the goal batch and the remaining batches in their shared low-dimensional space (PCA), which are used to align and integrate the batches into the goal batch.

We ran Seurat RPCA implemented in Seurat v4.4.0, available on their github. This is because subsequent version (>= v5) does not support getting corrected count matrix, only corrected PC space. See: https://github.com/satijalab/seurat/issues/8551.

We varied the number of PCs and nearest neighbours considered when running RPCA.

Seurat RPCA performs batch integration by identifying mutual nearest neighbors (anchors) between all batches using reciprocal PCA (RPCA). Instead of merging batches into a single PCA space, each batch is projected into the PCA space of the others, and anchors are found where cells are mutual nearest neighbors across these projections. These anchors are then used to compute batch correction vectors, aligning the batches into a shared corrected space.

We ran Seurat RPCA implemented in Seurat v4.4.0, available on their github. This is because subsequent version (>= v5) does not support getting corrected count matrix, only corrected PC space. See: https://github.com/satijalab/seurat/issues/8551

We varied the number of PCs and nearest neighbours considered when running RPCA.

The component works by reading and writing back the 'unintegrated' data without performing any operation.

The method returns only samples from one batch (batch 1), for both the splits. This act as a positive control method for the metrics that compare technical replicates between runs ('horizontal' metrics), as the difference between a sample compared to itself should be zero. It also acts as a positive control for 'vertical' metrics, as the samples being returned are from the same batch.

This negative control randomly permutes cell-to-sample (hence batch) assignments while keeping each cell's measured markers unchanged. This destroys any biological and batch specific structure but preserves marker expression.

Purpose:

• Provide a baseline to verify that integration methods outperform random assignment of cells to batches.

Example:

• A cell from a KO sample in batch 1 may be reassigned to any sample in the whole data. It may be reassigned to another KO sample in batch 1 or 2, or to a WT sample in batch 1 or 2.

This negative-control method randomly permutes cell-to-cell type assignments. Cells remain assigned to their original batch (batch effects preserved). Within each batch, cells are reassigned to random samples, destroying biological/sample-specific structure (e.g., KO vs WT differences).

Purpose:

• Evaluate whether an integration method preserves differences between samples and biological groups while removing batch effects.

Example:

• A cell from a KO sample in batch 1 may be reassigned to any sample in batch 1 (KO or WT), but it will never be moved to batch 2.

This negative-control method randomly permutes cell-to-cell type assignments. Cells will be assigned to any cell types, regardless of their original cell type or sample of origin or batch of origin.

Purpose:

• Evaluate whether an integration method preserves differences between cell types while removing batch effects.

Example:

• A Neutrophil from a KO sample in batch 1 may be reassigned to any cell type (B cell, T cell, Monocyte, etc.) from any sample in any batch.

First, a simple linear model is fitted for each paired sample and marker to determine the fraction of variance ($R^{2}$) explained by the batch covariate $B$. The average batch R-squared is then computed as the average of the $R^{2}$ values across all paired samples, markers, cell types. As a result, $\overline{R^2_B}_{cell\ type}$ quantifies how much of the total variability in the data is driven by batch effects. Consequently, lower values are desirable.

$\overline{R^2_B}{cell\ type} = \frac{1}{NCM}\sum{\substack{(x_{\mathrm{split1}},,x_{\mathrm{split2}})\ \textit{paired samples}}}^{N} \sum_{j=1}^{C} \sum_{i=1}^{M},R^2!\bigl(\mathrm{marker}_i \mid B\bigr)$

Where:

• $N$ is the number of paired samples, where $x_{\mathrm{split1}}$ and $x_{\mathrm{split2}}$ are the two technical replicates that have been batch-corrected. Technical replicates belong to different batches.
• $M$ is the number of markers
• $C$ is the number of cell types
• $B$ is the batch covariate

A high value of $\overline{R^2_B}{global}$ or $\overline{R^2_B}{cell\ type}$ indicates that the batch variable explains a large portion of the variance in the data, which indicates a higher level of batch effects. A good performance on $\overline{R^2_B}{global}$ but not on $\overline{R^2_B}{cell\ type}$ might indicate that the batch effect correction is not addressing cell type specific batch effects.

First, a simple linear model is fitted for each paired sample and marker to determine the fraction of variance ($R^{2}$) explained by the batch covariate $B$. The average batch R-squared is then computed as the average of the $R^{2}$ values across all paired samples, markers. As a result, $\overline{R^2_B}_{global}$ quantifies how much of the total variability in the data is driven by batch effects. Consequently, lower values are desirable.

$\overline{R^2_B}{global} = \frac{1}{N*M}\sum{\substack{(x_{\mathrm{split1}},,x_{\mathrm{split2}})\ \textit{paired samples}}}^{N} \sum_{i=1}^{M} ,R^2!\bigl(\mathrm{marker}_i \mid B\bigr)$

Where:

• $N$ is the number of paired samples, where $x_{\mathrm{split1}}$ and $x_{\mathrm{split2}}$ are the two technical replicates that have been batch-corrected. Technical replicates belong to different batches.
• $M$ is the number of markers
• $B$ is the batch covariate

A higher value of $\overline{R^2_B}_{global}$ indicates that the batch variable explains more of the variance in the data, which indicates a higher level of batch effects.

Compute the cell type local inverse simpson index (cLISI) for each cell, then return the median as the final score.

Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference between two probability distributions.

Here, EMD is used to compare marker expression distributions between paired samples from the same donor quantified across two different batches. For each paired sample, cell type, and marker, the marker expression values are first converted into probability distributions. This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1. The wasserstein_distance function from SciPy is then used to calculate the EMD between the two probability distributions belonging to the same cell type, marker, and a given paired samples. This is then repeated for every cell type, marker, and paired sample. Finally, the average of all these EMD values is computed and reported as the metric score.

A high score indicates large overall differences in the distributions of marker expressions between the paired samples, suggesting poor batch integration. A low score means the small differences in marker expression distributions between batches, indicating good batch integration.

Earth Mover Distance (EMD), also known as the Wasserstein metric, measures the difference between two probability distributions.

Here, EMD is used to compare marker expression distributions between all integrated samples from the same group. For each pair of samples, cell type, and marker, the marker expression values are first converted into probability distributions. This is done by binning the expression values into a range from -100 to 100 with a bin width of 0.1. The wasserstein_distance function from SciPy is then used to calculate the EMD between the two probability distributions belonging to the same cell type, marker, and a given paired samples. This is then repeated for every cell type, marker, and paired sample. Finally, the average of all these EMD values is computed and reported as the metric score.

A high score indicates overall, there is a large difference in distribution of marker expression after batch integration. A low score means that overall, the samples are well integrated.

The metric is based on the FlowSOM algorithm, a popular method which uses self-organizing maps for the viasualization/interpretation/clustering of cytometry data. The FlowSOM algorithm creates a tree structure that represents the relationships between different cell populations in the data.

For each paired sample of technical replicates (where 'split1' are integrated technical replicates from split 1 and 'split2' are integrated technical replicates from split 2)

1. A FlowSOM tree is created using cells in the sample from split 1.
2. Cells from the sample in split 2 are mapped onto the FlowSOM tree created in step 1.
3. A similarity measure is computed by comparing cell type proportions of 'split2' and 'split1' in each flowsom cluster.

Ideally, the proportions of cell types in the flowsom clusters of the paired samples should be similar, as they are technical replicates.

The FlowSOM mapping similarity measure can be expressed as follows: $\text{FlowSOM mapping similarity} = 100 - \text{FlowSOM mapping dissimilarity}$

The $\text{FlowSOM mapping dissimilarity}$ is:

$\text{FlowSOM mapping dissimilarity} = \sum_{k=1}^{K}w_{k}\sum_{c=1}^{C}\abs{P^{split1}{k,c} - P^{split2}{k,c}}$

Where:

• $K$ is the number of flowsom clusters
• $C$ is the number of cell types
• $w_{k}$ is the weight of flowsom cluster $k$. It refers to the number of cells in flowsom cluster $k$, divided by the total number of cells in the flowsom tree. Note: the flowsom tree contains all the cells of a sample pair (both integrated and validation).
• $P^{split1}_{k,c}$ is the percentage of cell type $c$ in flowsom cluster $k$ of the split 1 sample
• $P^{split2}_{k,c}$ is the percentage of cell type $c$ in flowsom cluster $k$ of the split 2 sample

The average FlowSOM mapping similarity among all paired samples is computed and reported as the final metric value. Unlabelled cells are excluded from the analysis.

Compute the integration local inverse simpson index (iLISI) for each cell, then return the median as the final score.

The metric compares the number of cell type specific marker expression peaks between unintegrated and batch integrated data. The number of peaks is calculated using the scipy.signal.find_peaks function. The metric is calculated as the absolute difference between the number of peaks in the unintegrated and batch-normalized data. The (cell type) marker expression profiles are first smoothed using kernel density estimation (KDE) (scipy.stats.gaussian_kde), and then peaks are then identified using the scipy.signal.find_peaks function. For peak calling, the prominence parameter is set to 0.1 and the height parameter is set to 0.05*max_density. Ratio of inconsistent peaks is defined as number of cases where the number of peaks differ between the two splits in the batch normalized data divided by the total number of cases. Cases where there are different number of peaks between the two splits in the unintegrated data are ignored from the denominator. A lower score indicates better performance, means there are less cases with inconsistent peaks after batch integration.

CyTOF data (whole blood) from 2 human donors, each with unstimulated and stimulated samples. The stimulated samples were treated with Interferonα (IFNα) and Lipopolysaccharide (LPS). The mass cytometry cytometry panel includes 21 surface markers for cell phenotyping and 14 markers for functional analysis of the signaling responses. The data has been arcsinh transformed with cofactor 5 and pregated to only include events which fall either in the CD66-CD45+ or CD66+ gates.

Mass cytometry data of PBMC from CLL patients and healthy donors. Data was originally used in the CytofRUV paper. Data has been preprocessed (bead normalised, arcsinh with cofactor 5, pregated on live single cells).

Flow cytometry data of spleens from 4 WT (IKK2 fl/fl CD11c-cre +/+) and 4 KO (IKK2 fl/fl CD11c-cre Tg/+) B6 mice, measured with a 22-color panel and 2 different instrument settings. Data has been preprocessed (compensated with a batch-specific compensation matrix, logicle transformed, cleaned with PeacoQC and pregated on live single CD45+ cells).