1 Introduction

Advances in flow and mass cytometry now enable the simultaneous measurement of more than 50 parameters per cell in millions of cells per sample. This high-dimensional data poses new demands on existing analysis techniques. Conventional gating, until now the gold standard for cytometry data analysis, faces challenges in scaling to high-dimensional cytometry data. Issues with conventional gating strategies include inefficiency, subjectivity, and a significant risk of missing unknown or minor cell populations. Hence, computational methods for the unbiased analysis of high-dimensional cytometry data are required.
Several computational tools for cytometry data analysis are available, primarily designed for the R programming language. Examples include CATALYST / CyTOF workflow, FlowSOM, flowCore, flowAI, and diffcyt. These workflows address specific aspects of high-dimensional cytometry data analysis, including handling flow cytometry standard (FCS) files, preprocessing, dimensionality reduction, clustering based on self-organizing maps, and differential abundance analysis. However, interoperability between these packages is often limited and inefficient, and methods for advanced data analysis such as data integration and trajectory analysis are scarce.
Whereas high-dimensional cytometry data analysis methods are less established, algorithms for analyzing single-cell RNA sequencing (scRNAseq) data are highly developed and offer a host of analysis possibilities. The most commonly used framework for scRNAseq data analysis is Seurat, which covers every aspect of data analysis from preprocessing to a range of downstream analyses, including dimensionality reduction and clustering. A notable advantage of scRNAseq packages is their interoperability, since most third-party tools that offer additional analysis options, such as batch correction or trajectory analysis, integrate seamlessly into the Seurat workflow. However, these scRNAseq analysis algorithms are not readily accessible for the analysis of cytometry data.
With these challenges and possibilities in mind, we developed the R package Seumetry, a flow and mass cytometry analysis framework that offers state-of-the-art analysis options across the whole data life cycle. Seumetry includes a broad range of cytometry-specific algorithms for data handling, while it seamlessly integrates with Seurat, providing access to the latest scRNAseq analysis options.

2 Description of data

To demonstrate the functionality and workflow of Seumetry, we generated a high-dimensional spectral flow cytometry dataset comprising 39 parameters of human intestinal immune cells. The test dataset consists of immune cells isolated from the intestinal mucosa of 7 adult human donors, including two anatomical layers: epithelium and lamina propria. The data was manually pre-gated on single live cells. For spectral flow cytometry, the raw recordings are spectrally unmixed in the acquisition software before export. Seumetry works on these unmixed intensities and does not perform the unmixing itself.

3 Setup

To run the vignette, these tools need to be installed additionally.

devtools::install_github("saeyslab/CytoNorm")
devtools::install_github("biosurf/cyCombine")
BiocManager::install("EnhancedVolcano")
install.packages("pheatmap")

Load libraries

library(Seumetry)
library(pheatmap)
library(EnhancedVolcano)
library(ggplot2)
library(dplyr)

Download test dataset

# zenodo url
zenodo <- "https://zenodo.org/records/11935872/"
# create data directory
dir.create("data/fcs", recursive = TRUE, showWarnings = FALSE)
# download panel and metadata
download.file(paste0(zenodo, "files/metadata.csv?download=1"), "data/metadata.csv",
    quiet = TRUE)
download.file(paste0(zenodo, "files/panel.csv?download=1"), "data/panel.csv", quiet = TRUE)
# download fcs files
for (file in read.csv("data/metadata.csv")$file_name) download.file(paste0(zenodo,
    "files/", file, "?download=1"), paste0("data/fcs/", file), quiet = TRUE)

Set global ggplot2 options and define theme to make beautiful plots.

# adjust fill and color globally using options
discrete_colors <- c("#5988B2", "#C9635E", "#67976B", "#C88F67", "#A583B0", "#CDB86A",
    "#ABCCE2", "#7F007F", "#C8A5A3", "#B9E1B8", "#C5807B", "#E0CDB1", "#917C6F",
    "#5988B2", "#D5E8AE", "#CC7E78", "#B4C9E3", "#AB8BAF", "#97C496", "#5F6E9D",
    "#B3CC96", "#AF8B99", "#C88F67", "#C16D6B", "#8E8E8E", "#5E93C2", "#FF8FBF",
    "#AFD192", "#8E8E8E", "#8BCF92", "#5F6E9D", "#8E6A8E", "#A88169", "#9D7BAE",
    "#6D8D7E", "#6E6E6E", "#9D4949", "#F2E89C", "#A9FF93", "#93A9C5", "#B76CA8",
    "#D4857B", "#4C0000", "#BDBA88", "#6E6E6E", "#FFDF7F", "#B2B2B2", "#366E5E",
    "#C6C6C6", "#DDBE8F", "#C6C6C6", "#89603A", "#9C7A8D", "#7CCCBE", "#AFFF9F")
options(ggplot2.discrete.colour = discrete_colors)
options(ggplot2.discrete.fill = discrete_colors)
# define new gradient colors
gradient_colors <- colorRampPalette(c("steelblue", "slategray2", "white", "tan1",
    "firebrick3"))(100)
# define a theme to match plots from other packages with Seumetry plots
set_theme <- list(theme_linedraw(), theme(aspect.ratio = 1, panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()))

4 Metadata and panel

Metadata and panel are data.frames required by Seumetry for loading FCS files, creating a Seurat object, and preprocessing the data.

metadata <- read.csv("data/metadata.csv")
panel <- read.csv("data/panel.csv")

Metadata: a data.frame with at least 2 columns: “file_name” and “sample_id”. All additional columns are used as metadata columns and are added to the Seurat object.

head(metadata)
##                 file_name   sample_id patient_id tissue layer
## 1 1_Ileum_IEL_CD45pos.fcs 1_Ileum_IEL          1  ileum   IEL
## 2 1_Ileum_LPL_CD45pos.fcs 1_Ileum_LPL          1  ileum   LPL
## 3 2_Ileum_IEL_CD45pos.fcs 2_Ileum_IEL          2  ileum   IEL
## 4 2_Ileum_LPL_CD45pos.fcs 2_Ileum_LPL          2  ileum   LPL
## 5 3_Ileum_IEL_CD45pos.fcs 3_Ileum_IEL          3  ileum   IEL
## 6 3_Ileum_LPL_CD45pos.fcs 3_Ileum_LPL          3  ileum   LPL

Panel: a data.frame with at least 2 columns: “fcs_colname” and “antigen”. The fcs_colname are the names of the channels in the FCS file. Antigen is the desired name for downstream analysis. Additionally, columns for transformation are required. See below (transformation) for details.

head(panel)
##         fcs_colname antigen arcsinh_cofactor biexp_neg biexp_pos biexp_width
## 1             APC.A  IL15RA             6000         0       4.0       -1000
## 2         APC.Cy7.A    CD27             6000         0       4.5       -1000
## 3    APC.Fire.810.A    CCR7             6000         0       4.0       -1000
## 4        APC.R700.A   CD127             6000         0       4.5       -1000
## 5 Alexa.Fluor.647.A    CD1C             6000         1       4.5       -1000
## 6           BB515.A   CD141             6000         1       6.0       -1000

5 Loading FCS files

FCS files are loaded based on *.fcs files in the folder and on “file_name” column in metadata data.frame.
Each cell will receive a unique cell ID based on sample_id provided in metadata. Under the hood, create_flowset() wraps flowCore::read.flowSet().

fcs_fs <- create_flowset("data/fcs", metadata)
fcs_fs
## A flowSet with 14 experiments.
## 
## column names(50): FSC.A FSC.H ... eFluor.506.A Time

It is highly recommended to save the “fcs_fs” object.

6 Create Seurat object

All FCS files will be merged and metadata will be added based on sample_id and metadata data.frame. Raw data is saved in Seurat assay “fcs”. Only channels present in the panel data.frame are kept, and these channels are renamed from “fcs_colname” to “antigen”. Channels that were not indicated in the panel data.frame are stored in Seurat assay “unused”. All panel data are stored in slot for easy access.

seu <- create_seurat(fcs_fs, panel, metadata, derandomize = FALSE)
seu
## An object of class Seurat 
## 50 features across 1151826 samples within 2 assays 
## Active assay: fcs (39 features, 0 variable features)
##  1 layer present: counts
##  1 other assay present: unused

Derandomization: For mass cytometry data, it might be desired to derandomize data, which will round intensity values up to the nearest whole number (see https://biosurf.org/cytof_data_scientist.html#34_Data_transformations). This is implemented in the create_seurat function by setting derandomize to TRUE.

# plot cellnumbers per sample
plot_cellnumber(seu)

7 Preprocessing

7.1 Bead normalization

For mass cytometry data, bead normalization may be required. This vignette uses a flow cytometry test dataset, so bead normalization is not applied. Please refer to section “Other features” for a description of bead normalization implementation.

7.2 Compensation

Compensation is based on spillover matrices that can be provided within the FCS file or as a data.frame. For spectral flow cytometry, the data is already unmixed by the acquisition software. An additional compensation can still be applied in Seumetry to correct residual spillover between detectors, which is what we do below. This is not always necessary, so the step can be skipped when the unmixed data already looks clean.

There are different options how to supply the matrix. For details see compensate_data function documentation. Here, we use option 3.
- Option 1) Use external spillover matrix directly: same compensation for all files.
- Option 2) Use spillover matrix saved in FCS files: same compensation for all files.
- Option 3) Use spillover matrix saved in FCS files: compensation matrix used from individual FCS files, thus can be different for each file.

The compensate_data function will use flowCore::compensate() to compensate the raw values and write a new assay into the Seurat object called “comp”. The compensated data are saved in “counts” layer.
Warning: if multiple spillover matrices are present in FCS files, use the correct one!

Here, we use option 3 to compensate the data.

# Check different matrices using:
names(flowCore::spillover(fcs_fs[[1]]))
## [1] "SPILL"      "spillover"  "$SPILLOVER"
# In this case, spillover matrix in column 3 is the correct one.
seu <- compensate_data(fcs_fs, seu, fcs_matrix = 3)
# Check that compensation worked.
plot1 <- plot_cyto(seu, x = "IgD", y = "CD3", assay = "fcs", layer = "counts", scale = "log",
    rasterize = TRUE) + ggtitle("Uncompensated")
plot2 <- plot_cyto(seu, x = "IgD", y = "CD3", assay = "comp", layer = "counts", scale = "log",
    rasterize = TRUE) + ggtitle("Compensated")
plot1 + plot2

7.3 Transformation

The following transformations are possible: arcsinh and biexp. Both stabilize the variance and compress high-intensity events so that populations across the full dynamic range become well separated. The arcsinh transformation is computed directly as asinh(x / cofactor), while the biexponential transformation wraps flowWorkspace::flowjo_biexp().

Transformation is performed on “counts” data of the DefaultAssay of the Seurat object. The DefaultAssay is either “fcs”, “beadnorm”, or “comp” depending on whether bead normalization or compensation was performed or not.

The transform_data function will return a Seurat object with transformed data written into the “data” layer.

Both arcsinh and biexp transformation can and should be used with custom parameters. The arcsinh cofactor (and the biexponential parameters) have a strong effect on how well populations separate and on the downstream clustering, so they should be set per panel rather than left at the default (arcsinh cofactor = 5). Per-marker values are supplied as columns in the panel data.frame (arcsinh_cofactor, or biexp_pos, biexp_neg and biexp_width), which is stored in when the Seurat object is created. Markers without an entry use the default. It is good practice to check the separation visually on biaxial plots after transforming (see below). If you prefer to set the cofactors automatically, you can estimate them per marker by variance stabilization, for example with the immunoClust package, and enter the resulting values in the panel data.frame.

seu <- transform_data(seu, "arcsinh")
plot1 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", layer = "counts", scale = "log",
    rasterize = TRUE) + ggtitle("Untransformed")
plot2 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", layer = "data", rasterize = TRUE) +
    ggtitle("Arcsinh transformation")
plot1 + plot2

7.4 Downsampling

Flow- and Mass cytometry data can contain millions of cells. If the number of cells differs largely between samples, it is advisable to downsample so overall differences are not driven by individual samples with high number of cells. Furthermore, downsampling can decrease computational time, if quick data exploration is desired.
It is highly recommended to save a non-downsampled Seurat object.

# set seed for downsampling for reproducibility
set.seed(42)
# downsample using Seurats subset function
Idents(seu) <- seu$sample_id
seu <- subset(seu, downsample = 20000)
seu
## An object of class Seurat 
## 89 features across 225882 samples within 3 assays 
## Active assay: comp (39 features, 0 variable features)
##  2 layers present: counts, data
##  2 other assays present: fcs, unused
# plot cellnumbers per sample
plot_cellnumber(seu)

8 Quality control

Cytometry is an inherently noisy technology. There will be events close to the axis or artefacts due to antibody aggregation. Seumetry provides multiple quality control tools to achieve clean data for downstream analysis.

8.1 Removal of outliers

It can occur with cytometry, that some events have extremely negative or positive values. Here, we can remove these outlier events 1) by setting a manual threshold for each channel or 2) by using an automatic removal algorithm based on isolation forest.

8.1.1 Manual removal

Use a named vector to indicate threshold for each channel. The function can take positive or negative thresholds.

# example: check a channel pair to determine if it contains outlier events
plot_cyto(seu, "CD4", "CD3", style = "2d_density", assay = "comp", layer = "data")

# manually set thresholds based on cyto plots
thresholds_neg <- c(CD3 = -4, CD4 = -2)
thresholds_pos <- c(CD3 = 4, CD4 = 4.5)
# remove values above or below thresholds
seu_manual <- remove_outliers_manual(seu, thresholds_neg)
seu_manual <- remove_outliers_manual(seu_manual, thresholds_pos, negative = FALSE)
plot_cyto(seu_manual, "CD4", "CD3", style = "2d_density", assay = "comp", layer = "data")

8.1.2 Automatic removal

Outliers are detected using an isolation forest. In ungated flow cytometry data, this algorithm mainly removes axis-near events, doublets, and debris. Even when the data has been pre-gated on single live cells, as in this dataset, some dying or membrane-compromised cells and other artefacts can slip through that gating, and these are picked up here. The score above which an event is regarded as an outlier can be supplied; it refers to the approximate depth it takes to isolate an observation. The threshold is usually between 0.6 and 0.8, with a default threshold of 0.7. The lower the threshold, the more cells are labeled as outliers.

The function will save two features for each cell into meta.data of the Seurat Object: 1) the isolation forest score (outlier_score) and 2) whether an event passed the threshold or not (outlier).

# identify outliers (default threshold: 0.7)
seu <- detect_outliers(seu, score_threshold = 0.7)
plot1 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", layer = "data") + ggtitle("Before removal")
plot2 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", layer = "data", style = "point",
    color = "outlier_score") + ggtitle("Outlier score")
plot3 <- plot_cyto(seu, x = "CD4", y = "CD3", assay = "comp", layer = "data", style = "point",
    color = "outlier") + ggtitle("Outlier")
plot4 <- plot_cyto(subset(seu, subset = outlier == FALSE), x = "CD4", y = "CD3",
    assay = "comp", layer = "data") + ggtitle("After removal")
plot1 + plot2 + plot3 + plot4

Remove outliers from Seurat object.

# remove outliers from Seurat object
seu <- subset(seu, subset = outlier == FALSE)

8.2 Removal of aggregates

With high-dimensional flow cytometry, antibody aggregates can occur. These are usually characterized by highly co-linear events that form diagonal structures in XY plots. This algorithm can identify potentially problematic channel combinations and identify aggregates in these channels.

Here is an example of such aggregates:

plot1 <- plot_cyto(seu, x = "CD28", y = "CXCR5")
plot2 <- plot_cyto(seu, x = "CXCR3", y = "CD1C")
plot1 + plot2

The first step is to identify channel combinations that potentially contain aggregates. This is done by selecting only double positive events (default: events > 1) and running a pearson correlation. The default threshold of labeling a channel potentially containing aggregates is pearson R > 0.7. The lower the threshold, the more channels are labeled as potentially containing aggregates.

The Pearson step is intentionally permissive. Its only job is to pick candidate channel pairs to inspect in the next (RANSAC) step, and it removes no cells, so the exact threshold is not critical. The relevant parameters of detect_aggregate_channels() are: - threshold (Pearson R, default 0.7): lower it if a channel pair you suspect contains aggregates is not being flagged; raise it to inspect fewer pairs. - pos_threshold (default 1): the cutoff separating positive from negative events, since aggregates are only detectable among double-positive events. Raise it for channels with high background. - percent_cells (default 0.5%): the minimum fraction of double-positive cells required before a pair is flagged, which guards against a high Pearson R driven by only a handful of cells.

The function “detect_aggregate_channels” returns a list of 2 matrices and 1 data.frame.
- Matrix 1: pearson correlation matrix.
- Matrix 2: binary correlation matrix based on threshold of pearson R (default threshold = 0.7).
- Data.frame: contains the channel combinations that passed the pearson R threshold.

problem_channels <- detect_aggregate_channels(seu, threshold = 0.7)

These correlation matrices can be plotted, for example, using a heatmap.

pheatmap::pheatmap(problem_channels[[1]], main = "Pearson R", border_color = NA)

pheatmap::pheatmap(problem_channels[[2]], main = "Pearson R>0.7", border_color = NA)

Next, we can plot the channel combinations that potentially contain aggregates and manually assess if these channels are really problematic.

plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
    "Channel_1"], y = problem_channels[[3]][i, "Channel_2"])
do.call(gridExtra::grid.arrange, c(plots, ncol = 4))

rm(plots)

If a channel combination does not look like it contains aggregates, remove that row from the data.frame (problem_channels[[3]]).

Next, we can detect aggregates in the problematic channels using a modified RANSAC algorithm. Each cell will get an aggregate_score, which is the number of channel combinations in which it was labeled as an aggregate. By default, potential aggregates are only labeled real aggregates if they occur in 2 or more channel combinations. The aggregate_score and whether a cell has passed the aggregate_score threshold (default >= 2) is stored in meta.data of the Seurat object.

This RANSAC step is where cells are actually flagged. A pair picked up by the permissive Pearson screen is not removed at that stage, it is only fitted here. Aggregates show up as events lying on a near x = y line through the origin, with the slope constrained to roughly 0.75 to 1.25, while real cells form the surrounding cloud. This is also why genuine biological co-expression, for example in NK cells or activated T cells, does not get systematically removed: co-expression gives an elevated but diffuse correlation, whereas aggregates fall on a tight, almost proportional line, and only events on that line are flagged. The key parameters of detect_aggregates() are: - score_threshold (default >= 2): how many channel pairs must flag a cell before it is called an aggregate, and the main precision/recall dial. Raise it for fewer false positives (more conservative); lower it to remove aggregates more aggressively. - fit_threshold: the residual cutoff for counting an event as lying on the aggregate axis; tighten it (smaller value) to require stricter collinearity. - min_to_fit / max_iteration (passed to the RANSAC fit): if many channel pairs are skipped because no model could be fitted, the function suggests reducing min_to_fit or increasing max_iteration.

seu <- detect_aggregates(seu, problem_channels[[3]])
# Check the aggregate removal performance
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
    "Channel_1"], y = problem_channels[[3]][i, "Channel_2"], style = "point", color = "aggregate_score")
do.call(gridExtra::grid.arrange, c(plots, ncol = 4))

rm(plots)
plots <- list()
for (i in 1:nrow(problem_channels[[3]])) plots[[i]] <- plot_cyto(seu, x = problem_channels[[3]][i,
    "Channel_1"], y = problem_channels[[3]][i, "Channel_2"], style = "point", color = "aggregate")
do.call(gridExtra::grid.arrange, c(plots, ncol = 4))

rm(plots)

If the result is satisfactory, the aggregates can be removed.

seu <- subset(seu, subset = aggregate == FALSE)

It is possible that some aggregates still remain. In that case it is advisable to fine tune the parameters described above. The most common adjustments are lowering the Pearson threshold in detect_aggregate_channels() to screen more channel pairs, or lowering score_threshold in detect_aggregates() to label aggregates more aggressively. Alternatively, remaining artefacts may be removed by detecting outliers again using the isolation forest.

9 Dimensionality reduction & clustering

Dimensionality reduction, clustering & visualizations are accessible via the Seurat workflow. We also integrated a more classic clustering approach for cytometry data: FlowSOM.

9.1 Dimensionality reduction

Multiple reductions are possible via the Seurat workflow: PCA, UMAP or tSNE. If UMAPs and clustering do not separate celltypes well, it can help to only select features for PCA that distinguish expected cell subsets, e.g. CD3, CD4, and CD8 if working with T cells.

Furthermore, input data and data processing can be adjusted:
- No centering (data resembles raw intensities more closely)
- No scaling (preserves relationship of absolute marker intensities)
- Use (selected) features directly as input

Here, all features are used for scaling, centering and PCA. For UMAP, 10 PCs are used.

# scale data
seu <- Seurat::ScaleData(seu, features = row.names(seu), do.scale = TRUE, do.center = TRUE)
# run PCA for all features
seu <- Seurat::RunPCA(seu, features = row.names(seu), approx = FALSE)
# run UMAP using 10 PCs
seu <- Seurat::RunUMAP(seu, dims = 1:10)

Integration of multiple datasets or batch correction methods based on PCA such as CCA integration or harmony can also be used (not done in this tutorial).

All Seurat visualizations are available. For more detail, see https://satijalab.org/seurat/

plot1 <- Seurat::DimPlot(seu, group.by = "sample_id") + set_theme
plot2 <- Seurat::DimPlot(seu, group.by = "layer") + set_theme
plot1 + plot2

Seurat::VlnPlot(seu, features = c("CD8", "CD4"), pt.size = 0) & set_theme & RotatedAxis()

It is recommended to use a custom colorpalette when using UMAP to plot expression of markers.

Seurat::FeaturePlot(seu, features = c("CD8", "CD4")) & scale_color_gradientn(colors = gradient_colors) &
    set_theme

9.2 Clustering

For clustering, multiple methods such as Louvain algorithm are accessible via Seurat:

# run clustering using 10 PCs
seu <- Seurat::FindNeighbors(seu, dims = 1:10)
seu <- Seurat::FindClusters(seu, resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 222369
## Number of edges: 5162098
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8732
## Number of communities: 15
## Elapsed time: 100 seconds
Seurat::DimPlot(seu, group.by = "seurat_clusters") + set_theme

The resolution argument of Seurat::FindClusters controls how fine-grained the clustering is, with higher values giving more clusters. There is no single correct value. We recommend clustering at a somewhat higher resolution than you think you need, so that small or rare populations are resolved instead of being absorbed into larger ones. Clusters that split a single biological population can then be merged afterwards based on their marker expression (see the cluster-marker step below), for example by relabeling them into one cell type. This overclustering-then-merging strategy is also recommended by FlowSOM metaclustering and the CATALYST / CyTOF workflow (cytofWorkflow, via CATALYST::mergeClusters).

We also implemented FlowSOM and metaclustering: run_FlowSOM() wraps FlowSOM::BuildSOM() and FlowSOM::metaClustering_consensus().

# run flowSOM and metaclustering
seu <- run_FlowSOM(seu, metaclusters = 10, xdim = 10, ydim = 10)
Seurat::DimPlot(seu, group.by = "SOM_cl") + set_theme

Seurat::DimPlot(seu, group.by = "SOM_metacl") + set_theme

9.3 Cluster markers

Not only visualizations, but also other Seurat features are fully functional. For example, FindAllMarkers can be used to help annotate clusters.

markers <- Seurat::FindAllMarkers(seu, only.pos = TRUE)
top5 <- markers %>%
    group_by(cluster) %>%
    slice_max(avg_log2FC, n = 5)
Seurat::DotPlot(seu, unique(top5$gene), assay = "comp") + scale_color_gradientn(colors = gradient_colors) +
    theme_linedraw() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
    RotatedAxis()

10 Differential abundance

Differential abundance (DA) analysis is provided by the Seumetry function differential_abundance(), which wraps the edgeR package. DA is done using a generalized linear model (GLM, edgeR::glmFit()) to model number of cells per group given attribute/feature x condition (e.g. seurat_clusters x treatment) and a likelihood ratio test (LRT, edgeR::glmLRT()) to make a contrast between all conditions.

Using the function with “check_coeff = TRUE” will return all coefficients to design the proper contrast for the DA analysis.

# first check the contrasts
differential_abundance(seu, attribute = "seurat_clusters", group_by = "sample_id",
    formula = as.formula("~0+layer"), check_coeff = TRUE)
## [1] "layerIEL" "layerLPL"

Here, we compare IEL vs LPL using a contrast of c(1, -1).

# calculate differential abundance
da_res <- differential_abundance(seu, attribute = "seurat_clusters", group_by = "sample_id",
    formula = as.formula("~0+layer"), contrast = c(1, -1))
head(da_res)
##        logFC    logCPM        LR       PValue         FDR
## 12 -5.428742 12.773777 11.797022 0.0005932555 0.008898833
## 13  3.623494 11.592334  8.342532 0.0038727591 0.029045693
## 14  3.026072  8.603269  4.191187 0.0406346314 0.203173157
## 9  -1.502019 14.242988  2.688074 0.1011019542 0.274754636
## 4   1.233421 16.568219  2.583875 0.1079570620 0.274754636
## 7  -1.091062 14.614803  2.498226 0.1139746511 0.274754636

The results can be exported as a table and/or plotted, for example, using EnhancedVolcano.

EnhancedVolcano::EnhancedVolcano(da_res, lab = rownames(da_res), drawConnectors = TRUE,
    x = "logFC", y = "FDR", title = "IEL vs LPL", pCutoff = 0.1, FCcutoff = 0.1,
    ylim = c(0, 2.5), xlim = c(-6, 6)) + set_theme

11 Differential expression

In addition to Seurat “FindMarkers” function, Seumetry also includes a feature to find differential expression (DE) of markers based on a pseudobulk analysis. To this end, the median fluorescent intensity is calculated using the Seumetry function “median_expression”. The DE analysis is provided by the Seumetry function DE_pseudobulk(), which wraps the limma package: a linear model is fit (limma::lmFit()) followed by empirical Bayes statistics for differential expression (limma::eBayes()).

de_res <- DE_pseudobulk(seu, fixed_vars = "layer", contrast = "layerIEL-layerLPL")
head(de_res)
##              logFC     AveExpr         t    P.Value adj.P.Val         B
## HLA-DR -0.20004267 0.198944790 -2.768178 0.01548676 0.4149983 -3.002072
## CCR6    0.15600921 0.002797601  2.603752 0.02128196 0.4149983 -3.204036
## PD1    -0.07599857 0.165656719 -2.211019 0.04479843 0.5621678 -3.678805
## TCRgd  -0.09593926 0.267605577 -2.074175 0.05765824 0.5621678 -3.839342
## CD28   -0.10582602 0.241387962 -1.607892 0.13092954 0.8506270 -4.352114
## CXCR3   0.24238317 0.642482201  1.505640 0.15513876 0.8506270 -4.455014

The results can be exported as a table and/or plotted, for example, using EnhancedVolcano.

EnhancedVolcano::EnhancedVolcano(de_res, lab = rownames(de_res), drawConnectors = TRUE,
    x = "logFC", y = "P.Value", title = "IEL vs LPL", pCutoff = 0.1, FCcutoff = 0.1,
    xlim = c(-0.5, 0.5), ylim = c(0, 2)) + set_theme

12 Additional visualizations

12.1 Cytometry-like plots

Sometimes it is useful to create cytometry-like plots such as density plots. Seumetry includes various plotting functionalities using “plot_cyto”.

plot_cyto(seu, x = "CD3", y = "CD4")

plot_cyto(seu, x = "CD3", y = "CD4", style = "point", color = "sample_id")

plot_cyto(seu, x = "CD3", style = "density", color = "sample_id")

12.2 Sample-level PCA

A principal component analysis (PCA) based on median marker expression per sample can give a good impression of differences between conditions or replicates. Seumetry includes the “plot_pca” function to compute a sample-level PCA. The PCA can also be grouped by other factors than sample (based on Seurat Object meta.data).

plot_pca(seu)

plot_pca(seu, group_by = "sample_id", color = "layer")

12.3 Frequency plots

Seumetry includes a function to visualize the frequency or proportion of cells based on metadata columns. For example, it can be useful to assess whether the frequency of clusters is stable across different samples or differs between conditions.

plot_frequency(seu, "seurat_clusters", "sample_id")

plot_frequency(seu, "seurat_clusters", "layer")

12.4 Cell numbers

plot_cellnumber(seu, "sample_id")

13 Other features

13.1 Bead normalization

The test dataset used throughout this vignette is a flow cytometry dataset, and bead normalization is only available for mass cytometry data. Here, we implemented bead normalization in the Seumetry function bead_norm(), which wraps CATALYST::normCytof() from the CATALYST package. To showcase the bead normalization implementation in Seumetry, we use the raw_data provided by CATALYST.

library(CATALYST)
data("raw_data")

Now we need to prepare a panel and add cell IDs to load data into Seumetry.

# prepare panel dataframe
panel_cat <- data.frame(fcs_colname = unname(flowCore::parameters(raw_data[[1]])@data[["name"]]),
    antigen = unname(flowCore::parameters(raw_data[[1]])@data[["desc"]]))
panel_cat <- panel_cat[!is.na(panel_cat$antigen), ]
# give cells an ID (required for Seurat Object creation)
row.names(raw_data[[1]]@exprs) <- paste0("raw_data_1_", 1:nrow(raw_data[[1]]@exprs))
row.names(raw_data[[2]]@exprs) <- paste0("raw_data_2_", 1:nrow(raw_data[[2]]@exprs))
# create seurat + transform (cofactor 5 default)
seu_cat <- create_seurat(raw_data, panel_cat)
# to identify beads, normCytof() requires transformed data
seu_cat <- transform_data(seu_cat, "arcsinh")

Finally we can perform bead normalization. It will use the DefaultAssay of the Seurat Object and save normalized raw intensities in “beadnorm” assay. Transformation should be done afterwards. Based on CATALYST documentation, parameter k refers to “median window used for bead smoothing” and affects visualizations only.

seu_cat <- bead_norm(seu_cat, beads = "dvs", plot_res = TRUE, k = 50)

Furthermore, the returned Seurat Object contains two new metadata columns:
- beads: events that likely represent beads.
- bead_doublet: events that likely represent beads OR bead-cell doublets.
To remove beads and bead doublets from the object, you can use subset.

seu_cat <- subset(seu_cat, subset = bead_doublet == FALSE)
seu_cat
## An object of class Seurat 
## 123 features across 4851 samples within 3 assays 
## Active assay: beadnorm (56 features, 0 variable features)
##  2 layers present: counts, data
##  2 other assays present: fcs, unused

13.2 Integration of third-party tools

13.2.1 Conversion of Seurat Object

To allow integration with a wide variety of single-cell and cytometry libraries, we implemented a function (convert_seurat) to convert Seurat Objects to SingleCellExperiment, flowFrame, and flowSet objects. These can be used, for example, with batch correction methods (CytoNorm, cyCombine, gaussNorm) or other third-party libraries (diffcyt, CATALYST).

fs <- convert_seurat(seu, to = "FS", split_by = "sample_id")
fs
## A flowSet with 14 experiments.
## 
## column names(39): IL15RA CD27 ... CD11C IgD
ff <- convert_seurat(seu, to = "FF")
ff
## flowFrame object 'FlowFrame from Seurat'
## with 222369 cells and 39 observables:
##        name               desc     range  minRange  maxRange
## $P1  IL15RA              APC.A         6  -3.19209         5
## $P2    CD27          APC.Cy7.A         6  -2.07227         5
## $P3    CCR7     APC.Fire.810.A         5  -2.35942         4
## $P4   CD127         APC.R700.A         6  -1.83344         5
## $P5    CD1C  Alexa.Fluor.647.A         7  -2.81388         6
## ...     ...                ...       ...       ...       ...
## $P35   CD19    Spark.NIR.685.A         5  -3.30207         4
## $P36  CD123 Super.Bright.436.A         7  -3.30280         6
## $P37    CD4     cFluor.YG584.A         5  -3.11489         4
## $P38  CD11C       eFluor.450.A         6  -2.91153         5
## $P39    IgD       eFluor.506.A         6  -3.23538         5
## 290 keywords are stored in the 'description' slot
sce <- convert_seurat(seu, to = "SCE")
sce
## class: SingleCellExperiment 
## dim: 39 222369 
## metadata(0):
## assays(1): counts
## rownames(39): IL15RA CD27 ... CD11C IgD
## rowData names(6): fcs_colname antigen ... biexp_pos biexp_width
## colnames(222369): 1_Ileum_IEL_0 1_Ileum_IEL_6 ... 7_Ileum_LPL_8968
##   7_Ileum_LPL_8969
## colData names(20): orig.ident nCount_fcs ... SOM_cl SOM_metacl
## reducedDimNames(2): pca umap
## mainExpName: NULL
## altExpNames(0):

13.2.2 Example integration: CytoNorm

Here, we show an example of integrating a third-party tool. We use CytoNorm as an example, as it is widely used as a cytometry batch correction method.

CytoNorm is used for batch correction, but just to showcase the integration with Seumetry, we will use it as a donor- or sample normalization method. To this end, we will pool a subset of cells from each sample and use it as a reference sample. In a second step, we will normalize all samples based on the model trained with the reference sample.

First we create a flowSet of compensated and normalized intensities from a subset of cells from each sample. We split the samples by mucosal layer. This will be used as a reference for CytoNorm. We also create a flowSet of all samples.

# subset Seurat object
seu_sub <- subset(seu, downsample = 1000)
# create reference FlowSet
fs_ref <- convert_seurat(seu_sub, to = "FS", split_by = "layer", assay = "comp",
    layer = "data")
# create FlowSet of all samples
fs <- convert_seurat(seu, to = "FS", split_by = "sample_id", assay = "comp", layer = "data")

Next, we train the CytoNorm model.

# train model
model <- CytoNorm::CytoNorm.train(files = fs_ref, labels = c("IEL", "LPL"), channels = row.names(seu),
    transformList = NULL, plot = TRUE, seed = 42)

Finally, we can normalize all samples based on this model. CytoNorm will write FCS files for each normalized sample.

# match each sample's layer in the order of the flowSet frames
labels <- metadata$layer[match(flowCore::sampleNames(fs), metadata$sample_id)]
# run normalization
fs_norm <- CytoNorm::CytoNorm.normalize(model = model, files = fs, labels = labels,
    transformList = NULL, prefix = "", transformList.reverse = NULL, outputDir = "data/cytonorm")

The resulting normalized FCS files can be loaded, and the intensities can be added to a new assay of the Seurat object and visualized.

# name of FCS channels was changed during conversion; adjust panel accordingly
panel_new <- panel
panel_new$fcs_colname <- panel_new$antigen
# create Seurat object from flowSet of normalized data
seu_norm <- create_seurat(fs_norm, panel_new, metadata, derandomize = FALSE)
# add normalized intensities to other Seurat object
seu[["cytonorm"]] <- CreateAssayObject(data = GetAssayData(seu_norm))
plot1 <- plot_cyto(seu, x = "CD4", assay = "comp", layer = "data", style = "density",
    color = "sample_id") + ggtitle("Unnormalized")
plot2 <- plot_cyto(seu, x = "CD4", assay = "cytonorm", layer = "data", style = "density",
    color = "sample_id") + ggtitle("Normalized")
plot1 + plot2

13.2.3 Example integration: cyCombine

Another widely used cytometry batch correction method is cyCombine, which groups cells with a self-organizing map and corrects each group with ComBat. cyCombine’s batch_correct() function can operate directly on a Seurat object, so a Seumetry object can be passed without any conversion. It only requires a “batch” column in the meta.data and corrects the transformed intensities in the “data” layer of the active assay. We pass it the compensated assay (“comp”) on its own (via DietSeurat), so the correction is applied to those markers.

As with the CytoNorm example above, our test dataset contains no technical batches, so to showcase the integration we treat each donor (patient_id) as a “batch” while preserving the biological mucosal layer (IEL vs LPL) as a covariate.

# cyCombine requires a 'batch' column in the meta.data (factor for plotting)
seu$batch <- factor(seu$patient_id)
# hand cyCombine a clean object with only the compensated, arcsinh-transformed
# assay ('comp')
seu_input <- DietSeurat(seu, assays = "comp")
seu_cc <- cyCombine::batch_correct(seu_input, covar = "layer", markers = row.names(seu_input),
    norm_method = "scale", seed = 42)
# add corrected intensities to a new assay of the original object
seu[["cycombine"]] <- seu_cc[["cyCombine"]]

The corrected intensities can be visualized like any other assay. Here, the donor-to-donor spread in CD4 is reduced after correction.

plot1 <- plot_cyto(seu, x = "CD4", assay = "comp", layer = "data", style = "density",
    color = "batch") + ggtitle("Uncorrected")
plot2 <- plot_cyto(seu, x = "CD4", assay = "cycombine", layer = "data", style = "density",
    color = "batch") + ggtitle("Batch corrected")
plot1 + plot2

13.3 Median fluorescent intensity

The Seumetry function “median_expression” can be used similar to Seurat “AverageExpression” function, but uses the median instead of the average. This function can be useful, for example, to plot median expression per cluster or condition.

sample_mfi <- median_expression(seu, group_by = "sample_id")
head(sample_mfi)
##        1_Ileum_IEL 1_Ileum_LPL 2_Ileum_IEL 2_Ileum_LPL  3_Ileum_IEL 3_Ileum_LPL
## IL15RA -0.17184436 -0.19478872  0.08925699  0.11194346 -0.067248831 -0.09096448
## CD27   -0.02169249 -0.05417867  0.01792439 -0.03766002  0.005206601 -0.01126040
## CCR7   -0.06621804 -0.12664701 -0.03386669 -0.04040400 -0.041614615 -0.03376832
## CD127   0.88408459  1.30156459  0.61199466  0.87384912  0.627081433  0.73614117
## CD1C    0.09165838  0.10414324 -0.18360783 -0.18490663 -0.046887197 -0.03263464
## CD141  -0.27540470 -0.66448338 -0.05679416 -0.13388351  0.117899077  0.15485369
##        4_Ileum_IEL  4_Ileum_LPL 5_Ileum_IEL 5_Ileum_LPL  6_Ileum_IEL
## IL15RA -0.10841902 -0.171803202 -0.10551499 -0.15142469 -0.127147375
## CD27    0.02086014  0.002260578  0.01017423  0.09883527 -0.008690137
## CCR7   -0.04067647 -0.049189390 -0.08166193 -0.02394237 -0.041863951
## CD127   0.63469098  0.674176088  0.45213541  0.16694058  0.080616511
## CD1C    0.01590664  0.095685465  0.02145164  0.27190240  0.315171889
## CD141  -0.23024512 -0.247446724 -0.64072933 -0.51548617 -0.375830419
##        6_Ileum_LPL 7_Ileum_IEL 7_Ileum_LPL
## IL15RA -0.13363610 -0.12398204 -0.21331916
## CD27   -0.01000960 -0.01981501 -0.02549146
## CCR7   -0.07634711 -0.06969665 -0.12100583
## CD127   0.90268127  0.28923809  0.17169467
## CD1C    0.27682553  0.20881542  0.32878685
## CD141  -0.42677419 -0.38622042 -0.48161871

13.4 Export of FCS files

Furthermore, sometimes it can be useful to use external software that relies on FCS files. For this purpose, Seumetry includes the export_fcs() function, which wraps flowCore::write.FCS(). This can be done, for example, to visualize cells of a specific cluster with external cytometry software.

# subset Seurat object to specific cluster
seu_sub <- subset(seu, subset = seurat_clusters == 1)
# export FCS file containing cells from this cluster
export_fcs(seu_sub, filename = "cluster_1.fcs", assay = "fcs", layer = "counts")