Most of the pipeline and visualizations presented herein were adapted from Nowicka et al. (2017)’s “CyTOF workflow: differential discovery in high-throughput high-dimensional cytometry datasets”. For the complete workflow, go here.

# load required packages
library(CATALYST)
library(cowplot)
library(flowCore)
library(diffcyt)
library(scater)
library(SingleCellExperiment)

1 Example data

PBMC_fs:
a flowSet holding PBMCs samples from 4 patients, each containing between 500 and 1000 cells. For each sample, the expression of 10 cell surface and 14 signaling markers was measured before (REF) and upon BCR/FcR-XL stimulation (BCRXL) with B cell receptor/Fc receptor crosslinking for 30’, resulting in a total of 8 samples.
This data set represents a subset of data originating from Bodenmiller et al. (2012) that was also used in the citrus paper (Bruggner et al. 2014).
PBMC_panel:
a data.frame containing each marker’s column name in the FCS file (fcs_colname column), its targeted protein marker (antigen column), and the marker_class (“type” or “state”).
PBMC_md:
a data.frame where rows correspond to samples, and columns specify each sample’s file_name, sample_id, condition, and patient_id.

# load example data
data(PBMC_fs, PBMC_panel, PBMC_md)
PBMC_fs

## A flowSet with 8 experiments.
## 
##   column names:
##   CD3(110:114)Dd CD45(In115)Dd pNFkB(Nd142)Dd pp38(Nd144)Dd CD4(Nd145)Dd CD20(Sm147)Dd CD33(Nd148)Dd pStat5(Nd150)Dd CD123(Eu151)Dd pAkt(Sm152)Dd pStat1(Eu153)Dd pSHP2(Sm154)Dd pZap70(Gd156)Dd pStat3(Gd158)Dd CD14(Gd160)Dd pSlp76(Dy164)Dd pBtk(Er166)Dd pPlcg2(Er167)Dd pErk(Er168)Dd pLat(Er170)Dd IgM(Yb171)Dd pS6(Yb172)Dd HLA-DR(Yb174)Dd CD7(Yb176)Dd

head(PBMC_panel)

##      fcs_colname antigen marker_class
## 1 CD3(110:114)Dd     CD3         type
## 2  CD45(In115)Dd    CD45         type
## 3 pNFkB(Nd142)Dd   pNFkB        state
## 4  pp38(Nd144)Dd    pp38        state
## 5   CD4(Nd145)Dd     CD4         type
## 6  CD20(Sm147)Dd    CD20         type

head(PBMC_md)

##                 file_name sample_id condition patient_id
## 1 PBMC_patient1_BCRXL.fcs    BCRXL1     BCRXL   Patient1
## 2   PBMC_patient1_Ref.fcs      Ref1       Ref   Patient1
## 3 PBMC_patient2_BCRXL.fcs    BCRXL2     BCRXL   Patient2
## 4   PBMC_patient2_Ref.fcs      Ref2       Ref   Patient2
## 5 PBMC_patient3_BCRXL.fcs    BCRXL3     BCRXL   Patient3
## 6   PBMC_patient3_Ref.fcs      Ref3       Ref   Patient3

The code snippet below demonstrates how to construct a flowSet from a set of FCS files. However, we also give the option to directly specify the path to a set of FCS files (see next section).

# download exemplary set of FCS files
url <- "http://imlspenticton.uzh.ch/robinson_lab/cytofWorkflow"
fcs_zip <- "PBMC8_fcs_files.zip"
download.file(paste0(url, "/", fcs_zip), destfile = fcs_zip, mode = "wb")
unzip(fcs_zip)

# read in FCS files as flowSet
fcs_files <- list.files(pattern = ".fcs$")
fs <- read.flowSet(fcs_files, transformation = FALSE, truncate_max_range = FALSE)

2 Data preparation

Data used and returned throughout differential analysis are held in objects of the SingleCellExperiment class. To bring the data into the appropriate format, prepData() requires the following inputs:

x: a flowSet holding the raw measurement data, or a character string that specifies a path to a set of FCS files.
panel: a 2 column data.frame that contains for each marker of interest i) its column name in the raw input data, and ii) its targeted protein marker.
md: a data.frame with columns describing the experimental design.

Optionally, features will specify which columns (channels) to keep from the input data. Here, we keep all measurement parameters (default value features = NULL).

(sce <- prepData(PBMC_fs, PBMC_panel, PBMC_md))

## class: SingleCellExperiment 
## dim: 24 5428 
## metadata(2): experiment_info cofactor
## assays(1): exprs
## rownames(24): CD3 CD45 ... HLA_DR CD7
## rowData names(3): channel_name marker_name marker_class
## colnames: NULL
## colData names(3): sample_id condition patient_id
## reducedDimNames(0):
## spikeNames(0):
## altExpNames(0):

We provide flexibility in the way the panel and metadata table can be set up. Specifically, column names are allowed to differ from the example above, and multiple factors (patient ID, conditions, batch etc.) can be specified. Arguments panel_cols and md_cols should then be used to specify which columns hold the required information. An example is given below:

# alter panel column names
panel2 <- PBMC_panel
colnames(panel2)[1:2] <- c("channel_name", "marker")

# alter metadata column names & add 2nd condition
md2 <- PBMC_md
colnames(md2) <- c("file", "sampleID", "cond1", "patientID")
md2$cond2 <- rep(c("A", "B"), 4)

# construct SCE
prepData(PBMC_fs, panel2, md2, 
    panel_cols = list(channel = "channel_name", antigen = "marker"),
    md_cols = list(file = "file", id = "sampleID", 
        factors = c("cond1", "cond2", "patientID")))

Note that, independent of the input panel and metadata tables, the constructor will fix the names of mandatory slots for latter data accession (sample_id in the rowData, channel_name and marker_name in the colData). The md table will be stored under experiment_info inside the metadata.

3 Diagnostic plots

3.1 `plotCounts`: Number of cells measured per sample

The number of cells measured per sample may be plotted with plotCounts, or directly accessed via n_cells(). This plot should be used as a guide together with other readouts to identify samples where not enough cells were assayed.

n_cells(sce)

## BCRXL1   Ref1 BCRXL2   Ref2 BCRXL3   Ref3 BCRXL4   Ref4 
##    528    881    665    438    563    660    934    759

plotCounts(sce, color_by = "condition")

3.2 `plotMDS`: Multi-dimensional scaling plot

A multi-dimensional scaling (MDS) plot on median expresion values may be rendered with plotMDS. Such a plot will give a sense of similarities between samples in an unsupervised way and of key difference in expression before conducting any formal testing. In our example, we can see a clear separation between reference (REF) and stimulation condition (BCRXL).

plotMDS(sce, color_by = "condition")

3.3 `plotExprHeatmap`: Heatmap of (scaled) median marker expressions

plotExprHeatmap will show a heatmap on median marker intensities with hierarchically clustered rows (samples) and columns (markers). This plot should give an idea of which markers will drive sample clustering, and how similiar samples are in their expression profile. We specify bin_anno = TRUE to display expression values inside each bin, and row_anno = TRUE to include row annotations for each factor in metadata(daf).

plotExprHeatmap(sce, bin_anno = TRUE, row_anno = TRUE)

4 Clustering

4.1 `cluster`: FlowSOM clustering & ConsensusClusterPlus metaclustering

CATALYST provides a simple wrapper to perform high resolution FlowSOM clustering and lower resolution ConsensusClusterPlus metaclustering. By default, the data will be initially clustered into xdim = 10 x ydim = 10 = 100 groups. Secondly, the function will metacluster populations into 2 through maxK (default 20) clusters. To make analyses reproducible, the random seed may be set via seed. By default, if the colData()$marker_class column is specified, the set of markers with marker class “type” will be used for clustering. Alternatively, the markers that should be used for clustering can be specified with argument cols_to_use.

# specify markers to use for clustering
lineage_markers <- c("CD3", "CD45", "CD4", "CD20", 
    "CD33", "CD123", "CD14", "IgM", "HLA_DR", "CD7")
sce <- cluster(sce, features = lineage_markers, 
    xdim = 10, ydim = 10, maxK = 20, verbose = FALSE, seed = 1)

Let K = xdim x ydim be the number of FlowSOM clusters. cluster will add information to the following slots of the input SingleCellExperiment:

rowData:
- cluster_id: cluster ID as inferred by FlowSOM. One of 1, …, K.
colData:
- marker_class: factor "type" or "state". Specifyies whether a marker has been used for clustering or not, respectively.
metadata:
- SOM_codes: a table with dimensions K x (# type markers). Contains the SOM codes.
- cluster_codes: a table with dimensions K x (maxK + 1). Contains the cluster codes for all metaclusterings.
- delta_area: a ggplot object (see below for details).

4.2 Delta area plot

The delta area represents the amount of extra cluster stability gained when clustering into k groups as compared to k-1 groups. It can be expected that high stability of clusters can be reached when clustering into the number of groups that best fits the data. The “natural” number of clusters present in the data should thus corresponds to the value of k where there is no longer a considerable increase in stability (pleateau onset). For more details, the user can refer to the original description of the consensus clustering method (Monti et al. 2003).

# access & render delta area plot
metadata(sce)$delta_area

Differential analysis with CATALYST

2 April 2020

Package

Contents

1 Example data

2 Data preparation

3 Diagnostic plots

3.1 `plotCounts`: Number of cells measured per sample

3.2 `plotMDS`: Multi-dimensional scaling plot

3.3 `plotExprHeatmap`: Heatmap of (scaled) median marker expressions

4 Clustering

4.1 `cluster`: FlowSOM clustering & ConsensusClusterPlus metaclustering

4.2 Delta area plot

Differential analysis with CATALYST

2 April 2020

Package

Contents

1 Example data

2 Data preparation

3 Diagnostic plots

3.1 plotCounts: Number of cells measured per sample

3.2 plotMDS: Multi-dimensional scaling plot

3.3 plotExprHeatmap: Heatmap of (scaled) median marker expressions

4 Clustering

4.1 cluster: FlowSOM clustering & ConsensusClusterPlus metaclustering

4.2 Delta area plot

3.1 `plotCounts`: Number of cells measured per sample

3.2 `plotMDS`: Multi-dimensional scaling plot

3.3 `plotExprHeatmap`: Heatmap of (scaled) median marker expressions

4.1 `cluster`: FlowSOM clustering & ConsensusClusterPlus metaclustering