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.

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 shuffles all cells in the input data. Cells lose all identity information (sample, batch, and cell type).

Example: a cell from mouse 3 batch 1 may be reassigned to any sample in the data (mouse 4 batch 1 or mouse 6 batch 2).

This negative-control method randomly shuffles cells within each batch. Cells lose their sample and cell type identity but remain in their original batch.

Example: A cell from mouse 3 batch 1 may be reassigned to mouse 4 batch 1, but not to another batch (e.g., mouse 6 batch 2).

This negative-control method randomly shuffles cells within each cell type. Cells retain their cell type identity but lose their sample and batch identity.

Example: A B cell remains a B cell but if may be reassigned from mouse 3 batch 1 to mouse 4 batch 1 or mouse 6 batch 2.

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.

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 and cells from split 2 are mapped onto it.
2. A FlowSOM tree is created using cells in the sample from split 2 and cells from split 1 are mapped onto it.
3. For each direction, a similarity measure is computed by comparing cell type proportions of the two splits 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 donors and both mapping directions is computed and reported as the final metric value. Unlabelled cells are excluded from the analysis.

Additionally, for debugging and inspection, the component outputs the per-donor cluster-by-cell-type absolute difference matrices in the output AnnData uns.

Measures how much batch integration shifts the effect size of biological group differences in functional marker expression.

Cohen's d = (mean_A − mean_B) / pooled SD is computed per (functional marker, cell type) pair using per-sample mean expressions, restricted to the pairs in sig_unintegrated (i.e. those found significant in both unintegrated batches by the Wilcoxon metric above). This covers every pair with a genuine baseline group difference, including ones where that difference no longer reached significance after integration — that loss of significance is itself a form of effect size drift this metric is meant to capture.

Group A and group B are determined once from the unintegrated data's obs["group"] values, sorted alphabetically, and reused for every batch and split. Since Cohen's d is signed, this is what keeps a positive value meaning "group A higher than group B" consistently across batches and splits, rather than the sign flipping depending on which group happened to be labelled A in a given batch or split.

Ground truth (unintegrated):

• Cohen's d is computed per batch, then averaged across the two batches.

Integrated:

• Cohen's d is computed separately for split 1 and split 2 (batch-agnostic).
• The absolute deviation from the unintegrated ground truth is computed per split, then averaged: mean(|d_s1 − d_unintegrated|, |d_s2 − d_unintegrated|).
• Taking abs() before averaging prevents opposite-direction errors in the two splits from cancelling each other out.

Score = mean of the above across all (marker, cell type) pairs. A lower score indicates that integration better preserved the original effect sizes.

Assesses whether batch integration preserves biologically meaningful differences in functional marker expression between biological groups (e.g. WT vs KO).

Per-sample mean expression is computed for each (functional marker, cell type) pair. A two-sided Wilcoxon rank-sum test compares the two biological groups. Raw p-values are used without FDR correction, as the minimum achievable p-value with n=3 samples per group is 0.10, making FDR adjustment counterproductive.

Group A and group B are determined once from the unintegrated data's obs["group"] values, sorted alphabetically, and reused for every batch and split. This keeps "group A" and "group B" referring to the same biological group throughout, rather than being re-derived independently per batch/split.

Baseline (unintegrated data):

• Tests are run separately within each batch.
• A pair enters the baseline only if significant (p ≤ 0.1) in BOTH batches.

Preservation (integrated data):

• The same test is run per split, ignoring batch labels.
• A pair is counted as preserved if significant in BOTH split 1 and split 2.

Score = number of preserved pairs / number of baseline pairs. A higher score indicates more biologically relevant group differences were maintained.

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 (cell type) marker expression values of the two splits are first standardised together, using the pooled mean and standard deviation of the two splits, so that the peak calling thresholds mean the same thing for every marker. They are then smoothed using kernel density estimation (KDE) (scipy.stats.gaussian_kde, evaluated on a grid of 100 points), and peaks are identified using the scipy.signal.find_peaks function. For peak calling, the height parameter is set to 0.1 and the prominence parameter is set to 0.01.

Case Definitions:

• Case 1: Consistent peaks in both unintegrated AND integrated data (ideal outcome). Example: A marker shows 2 peaks for both splits in unintegrated AND integrated data.
• Case 2: Consistent peaks in unintegrated data but inconsistent between integrated data and unintegrated data in each split (problematic - batch integration introduced inconsistency). Batch correction broke the consistency. Example 1: A marker shows 2 peaks in both unintegrated splits, but after integration shows 2 peaks in split 1 and 1 peak in split 2. Example 2: A marker shows 2 peaks in both unintegrated splits, but after integration shows 1 peak in both splits
• Case NGT: Inconsistent peaks in unintegrated data, consistent in integrated data (excluded from metric) OR inconsistent peaks in both unintegrated AND integrated data (excluded from metric). These cases are excluded from calculation because there are differences in unintegrated that cannot be accounted for. Example: A marker is inconsistent in unintegrated data (2 and 3 peaks in split 1 and split 2) and is either inconsistent after integration (1 and 2 peaks in split 1 and split 2) or consistent after integration (2 peaks in both split 1 and split 2).

Ratio of consistent peaks is defined as the number of Case 1 occurrences divided by the total of Case 1 and Case 2 occurrences. Cases NGT (where unintegrated data has inconsistent peaks) are excluded from the denominator, but is reported in the final anndata output for diagnostic. A higher score indicates better performance, meaning there are fewer cases with inconsistent peaks after batch integration.

For edge cases where methods return only zero values for a given marker/donor/cell type, but is not the case in the unintegrated data, it is automatically assigned as Case 2.

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.

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).