2.1 Summary

The diffcyt methodology consists of two main components: (i) high-resolution clustering and (ii) empirical Bayes moderated tests adapted from transcriptomics.

We use high-resolution clustering to define a large number of small clusters representing cell populations. By default, we use the FlowSOM clustering algorithm (Van Gassen et al., 2015) to generate the high-resolution clusters, since we previously showed that this clustering algorithm gives excellent clustering performance together with fast runtimes for high-dimensional cytometry data (Weber and Robinson, 2016). However, in principle, other algorithms that can generate high-resolution clusters could also be used.

For the differential analyses, we use methods from the edgeR package (Robinson et al., 2010; McCarthy et al., 2012), limma package (Ritchie et al., 2015), and voom method (Law et al., 2014). These methods are widely used in the transcriptomics field, and have been adapted here for analyzing high-dimensional cytometry data. In addition, we provide alternative methods based on generalized linear mixed models (GLMMs), linear mixed models (LMMs), and linear models (LMs), originally implemented by Nowicka et al. (2017) (available in the CyTOF workflow from Bioconductor).

2.2 Differential abundance (DA) and differential states (DS)

The diffcyt methods can be used to test for differential abundance (DA) of cell populations, and differential states (DS) within cell populations.

To do this, the methodology requires the set of protein markers to be grouped into ‘cell type’ and ‘cell state’ markers. Cell type markers are used to define clusters representing cell populations, which are tested for differential abundance; and median cell state marker signals per cluster are used to test for differential states within populations.

The conceptual split into cell type and cell state markers also facilitates biological interpretability, since it allows the results to be linked back to known cell types or populations of interest.

2.3 Flexible experimental designs and contrasts

The diffcyt model setup enables the user to specify flexible experimental designs, including batch effects, paired designs, and continuous covariates. Linear contrasts are used to specify the comparison of interest.

2.4 More details

A complete description of the statistical methodology and comparisons with existing approaches are provided in our paper introducing the diffcyt framework (Weber et al., 2018).

4.1 Dataset

For the example workflow in this vignette, we use the Bodenmiller_BCR_XL dataset, originally from Bodenmiller et al. (2012).

This is a publicly available mass cytometry (CyTOF) dataset, consisting of paired samples of healthy peripheral blood mononuclear cells (PBMCs), where one sample from each pair was stimulated with B cell receptor / Fc receptor cross-linker (BCR-XL). The dataset contains 16 samples (8 paired samples); a total of 172,791 cells; and a total of 24 protein markers. The markers consist of 10 ‘cell type’ markers (which can be used to define cell populations or clusters), and 14 ‘cell state’ or signaling markers.

This dataset contains known strong differential expression signals for several signaling markers in several cell populations, especially B cells. In particular, the strongest observed differential signal is for the signaling marker phosphorylated S6 (pS6) in B cells (see Nowicka et al., 2017, Figure 29). In this workflow, we will show how to perform differential tests to recover this signal.

4.2 Load data from ‘HDCytoData’ package

The Bodenmiller_BCR_XL dataset can be downloaded and loaded conveniently from the HDCytoData Bioconductor ‘experiment data’ package. It can be loaded in either SummarizedExperiment or flowSet format. Here, we use the flowSet format, which is standard in the flow and mass cytometry community. For some alternative analysis pipelines, the SummarizedExperiment format may be more convenient. For more details, see the help file for this dataset in the HDCytoData package (run library(HDCytoData) and ?Bodenmiller_BCR_XL).

suppressPackageStartupMessages(library(HDCytoData))

## snapshotDate(): 2018-04-27

# Download and load 'Bodenmiller_BCR_XL' dataset in 'flowSet' format
d_flowSet <- Bodenmiller_BCR_XL_flowSet()

## snapshotDate(): 2018-04-27

## see ?HDCytoData and browseVignettes('HDCytoData') for documentation

## downloading 0 resources

## loading from cache 
##     'C:/Users/biocbuild/Documents/AppData/.ExperimentHub/1120'

suppressPackageStartupMessages(library(flowCore))

# check data format
d_flowSet

## A flowSet with 16 experiments.
## 
##   column names:
##   Time Cell_length CD3 CD45 BC1 BC2 pNFkB pp38 CD4 BC3 CD20 CD33 pStat5 CD123 pAkt pStat1 pSHP2 pZap70 pStat3 BC4 CD14 pSlp76 BC5 pBtk pPlcg2 pErk BC6 pLat IgM pS6 HLA-DR BC7 CD7 DNA-1 DNA-2 group_id patient_id sample_id population_id

# sample names
pData(d_flowSet)

##                                                                  name
## PBMC8_30min_patient1_BCR-XL.fcs       PBMC8_30min_patient1_BCR-XL.fcs
## PBMC8_30min_patient1_Reference.fcs PBMC8_30min_patient1_Reference.fcs
## PBMC8_30min_patient2_BCR-XL.fcs       PBMC8_30min_patient2_BCR-XL.fcs
## PBMC8_30min_patient2_Reference.fcs PBMC8_30min_patient2_Reference.fcs
## PBMC8_30min_patient3_BCR-XL.fcs       PBMC8_30min_patient3_BCR-XL.fcs
## PBMC8_30min_patient3_Reference.fcs PBMC8_30min_patient3_Reference.fcs
## PBMC8_30min_patient4_BCR-XL.fcs       PBMC8_30min_patient4_BCR-XL.fcs
## PBMC8_30min_patient4_Reference.fcs PBMC8_30min_patient4_Reference.fcs
## PBMC8_30min_patient5_BCR-XL.fcs       PBMC8_30min_patient5_BCR-XL.fcs
## PBMC8_30min_patient5_Reference.fcs PBMC8_30min_patient5_Reference.fcs
## PBMC8_30min_patient6_BCR-XL.fcs       PBMC8_30min_patient6_BCR-XL.fcs
## PBMC8_30min_patient6_Reference.fcs PBMC8_30min_patient6_Reference.fcs
## PBMC8_30min_patient7_BCR-XL.fcs       PBMC8_30min_patient7_BCR-XL.fcs
## PBMC8_30min_patient7_Reference.fcs PBMC8_30min_patient7_Reference.fcs
## PBMC8_30min_patient8_BCR-XL.fcs       PBMC8_30min_patient8_BCR-XL.fcs
## PBMC8_30min_patient8_Reference.fcs PBMC8_30min_patient8_Reference.fcs

# number of cells
fsApply(d_flowSet, nrow)

##                                     [,1]
## PBMC8_30min_patient1_BCR-XL.fcs     2838
## PBMC8_30min_patient1_Reference.fcs  2739
## PBMC8_30min_patient2_BCR-XL.fcs    16675
## PBMC8_30min_patient2_Reference.fcs 16725
## PBMC8_30min_patient3_BCR-XL.fcs    12252
## PBMC8_30min_patient3_Reference.fcs  9434
## PBMC8_30min_patient4_BCR-XL.fcs     8990
## PBMC8_30min_patient4_Reference.fcs  6906
## PBMC8_30min_patient5_BCR-XL.fcs     8543
## PBMC8_30min_patient5_Reference.fcs 11962
## PBMC8_30min_patient6_BCR-XL.fcs     8622
## PBMC8_30min_patient6_Reference.fcs 11038
## PBMC8_30min_patient7_BCR-XL.fcs    14770
## PBMC8_30min_patient7_Reference.fcs 15974
## PBMC8_30min_patient8_BCR-XL.fcs    11653
## PBMC8_30min_patient8_Reference.fcs 13670

# number of columns
dim(exprs(d_flowSet[[1]]))

## [1] 2838   39

# expression values
exprs(d_flowSet[[1]])[1:6, 1:6]

##       Time Cell_length        CD3     CD45       BC1        BC2
## [1,] 33073          30 120.823265 454.6009  576.8983 10.0057297
## [2,] 36963          35 135.106171 624.6824  564.6299  5.5991135
## [3,] 37892          30  -1.664619 601.0125 3077.2668  1.7105789
## [4,] 41345          58 115.290245 820.7125 6088.5967 22.5641403
## [5,] 42475          35  14.373802 326.6405 4606.6929 -0.6584854
## [6,] 44620          31  37.737877 557.0137 4854.1519 -0.4517288

4.3 Alternatively: load data from ‘.fcs’ files

Alternatively, you can load data directly from a set of .fcs files using the following code. Note that we use the options transformation = FALSE and truncate_max_range = FALSE to disable automatic transformations and data truncation performed by the flowCore package. (The automatic options in the flowCore package are optimized for flow cytometry instead of mass cytometry data, so these options should be disabled for mass cytometry data.)

# Alternatively: load data from '.fcs' files
files <- list.files(
  path = "path/to/files", pattern = "\\.fcs$", full.names = TRUE
)
d_flowSet <- read.flowSet(
  files, transformation = FALSE, truncate_max_range = FALSE
)

4.4 Set up meta-data

Next, we set up the ‘meta-data’ required for the diffcyt pipeline. The meta-data describes the samples and protein markers for this experiment or dataset. The meta-data should be saved in two data frames: experiment_info and marker_info.

The experiment_info data frame contains information about each sample, including sample IDs, group IDs, batch IDs or patient IDs (if relevant), continuous covariates such as age (if relevant), and any other factors or covariates. In many experiments, the main comparison of interest will be between levels of the group IDs factor (which may also be referred to as condition or treatment; e.g. diseased vs. healthy, or treated vs. untreated).

The marker_info data frame contains information about the protein markers, including channel names, marker names, and a vector to identify the class of each marker (cell type or cell state).

Below, we create these data frames manually. Depending on your experiment, it may be more convenient to save the meta-data in spreadsheets in .csv format, which can then be loaded using read.csv.

Extra care should be taken here to ensure that all samples and markers are in the correct order. In the code below, we display the final data frames to check them.

# Meta-data: experiment information

# check sample order
filenames <- as.character(pData(d_flowSet)$name)

# sample information
sample_id <- gsub("^PBMC8_30min_", "", gsub("\\.fcs$", "", filenames))
group_id <- factor(
  gsub("^patient[0-9]+_", "", sample_id), levels = c("Reference", "BCR-XL")
)
patient_id <- factor(gsub("_.*$", "", sample_id))

experiment_info <- data.frame(
  group_id, patient_id, sample_id, stringsAsFactors = FALSE
)
experiment_info

##     group_id patient_id          sample_id
## 1     BCR-XL   patient1    patient1_BCR-XL
## 2  Reference   patient1 patient1_Reference
## 3     BCR-XL   patient2    patient2_BCR-XL
## 4  Reference   patient2 patient2_Reference
## 5     BCR-XL   patient3    patient3_BCR-XL
## 6  Reference   patient3 patient3_Reference
## 7     BCR-XL   patient4    patient4_BCR-XL
## 8  Reference   patient4 patient4_Reference
## 9     BCR-XL   patient5    patient5_BCR-XL
## 10 Reference   patient5 patient5_Reference
## 11    BCR-XL   patient6    patient6_BCR-XL
## 12 Reference   patient6 patient6_Reference
## 13    BCR-XL   patient7    patient7_BCR-XL
## 14 Reference   patient7 patient7_Reference
## 15    BCR-XL   patient8    patient8_BCR-XL
## 16 Reference   patient8 patient8_Reference

# Meta-data: marker information

# source: Bruggner et al. (2014), Table 1

# column indices of all markers, lineage markers, and functional markers
cols_markers <- c(3:4, 7:9, 11:19, 21:22, 24:26, 28:31, 33)
cols_lineage <- c(3:4, 9, 11, 12, 14, 21, 29, 31, 33)
cols_func <- setdiff(cols_markers, cols_lineage)

# channel and marker names
channel_name <- colnames(d_flowSet)
marker_name <- gsub("\\(.*$", "", channel_name)

# marker classes
# note: using lineage markers for 'cell type', and functional markers for 
# 'cell state'
marker_class <- rep("none", ncol(d_flowSet[[1]]))
marker_class[cols_lineage] <- "type"
marker_class[cols_func] <- "state"
marker_class <- factor(marker_class, levels = c("type", "state", "none"))

marker_info <- data.frame(
  channel_name, marker_name, marker_class, stringsAsFactors = FALSE
)
marker_info

##     channel_name   marker_name marker_class
## 1           Time          Time         none
## 2    Cell_length   Cell_length         none
## 3            CD3           CD3         type
## 4           CD45          CD45         type
## 5            BC1           BC1         none
## 6            BC2           BC2         none
## 7          pNFkB         pNFkB        state
## 8           pp38          pp38        state
## 9            CD4           CD4         type
## 10           BC3           BC3         none
## 11          CD20          CD20         type
## 12          CD33          CD33         type
## 13        pStat5        pStat5        state
## 14         CD123         CD123         type
## 15          pAkt          pAkt        state
## 16        pStat1        pStat1        state
## 17         pSHP2         pSHP2        state
## 18        pZap70        pZap70        state
## 19        pStat3        pStat3        state
## 20           BC4           BC4         none
## 21          CD14          CD14         type
## 22        pSlp76        pSlp76        state
## 23           BC5           BC5         none
## 24          pBtk          pBtk        state
## 25        pPlcg2        pPlcg2        state
## 26          pErk          pErk        state
## 27           BC6           BC6         none
## 28          pLat          pLat        state
## 29           IgM           IgM         type
## 30           pS6           pS6        state
## 31        HLA-DR        HLA-DR         type
## 32           BC7           BC7         none
## 33           CD7           CD7         type
## 34         DNA-1         DNA-1         none
## 35         DNA-2         DNA-2         none
## 36      group_id      group_id         none
## 37    patient_id    patient_id         none
## 38     sample_id     sample_id         none
## 39 population_id population_id         none

4.5 Set up design matrix (or model formula)

To calculate differential tests, the diffcyt functions require a design matrix (or model formula) describing the experimental design. (The choice between design matrix and model formula depends on the differential testing method used; see help files for the differential testing methods for details.)

Design matrices can be created in the required format using the function createDesignMatrix(). Design matrices are required for methods diffcyt-DA-edgeR (default method for DA testing), diffcyt-DA-voom, and diffcyt-DS-limma (default method for DS testing).

Similarly, model formulas can be created with the function createFormula(). Model formulas are required for the alternative methods diffcyt-DA-GLMM (DA testing) and diffcyt-DS-LMM (DS testing).

In both cases, flexible experimental designs are possible, including blocking (e.g. batch effects or paired designs) and continuous covariates. See ?createDesignMatrix or ?createFormula for more details and examples.

suppressPackageStartupMessages(library(diffcyt))

# Create design matrix
# note: selecting columns 1 and 2, which contain group IDs and patient IDs
design <- createDesignMatrix(experiment_info, cols_design = 1:2)

4.6 Set up contrast matrix

A contrast matrix is also required in order to calculate differential tests. The contrast matrix specifies the comparison of interest, i.e. the combination of model parameters assumed to equal zero under the null hypothesis.

Contrast matrices can be created in the required format using the function createContrast(). See ?createContrast for more details.

Here, we are interested in comparing condition BCR-XL against Reference, i.e. comparing the BCR-XL level against the Reference level for the group_id factor in the experiment_info data frame. This corresponds to testing whether the coefficient for column group_idBCR-XL in the design matrix design is equal to zero. This contrast can be specified as follows. (Note that there is one value per coefficient, including the intercept term; and rows in the final contrast matrix correspond to columns in the design matrix.)

# Create contrast matrix
contrast <- createContrast(c(0, 1, rep(0, 7)))

# check dimensions
nrow(contrast) == ncol(design)

## [1] TRUE

4.7 Differential testing

The steps above show how to load the data, set up the meta-data, set up the design matrix, and set up the contrast matrix. Now, we can begin calculating differential tests.

Several alternative options are available for running the diffcyt differential testing functions. Which of these is most convenient will depend on the types of analyses or pipeline that you are running. The options are:

• Option 1: Run wrapper function using input data loaded from .fcs files. The input data can be provided as a flowSet, or a list of flowFrames, DataFrames, or data.frames.

• Option 2: Run wrapper function using previously created CATALYST daFrame object.

• Option 3: Run individual functions for the pipeline.

The following sections demonstrate these options using the Bodenmiller_BCR_XL example dataset described above.

# Test for differential abundance (DA) of clusters

# note: using default method 'diffcyt-DA-edgeR' and default parameters
# note: include random seed for reproducible clustering
out_DA <- diffcyt(
  d_input = d_flowSet, 
  experiment_info = experiment_info, 
  marker_info = marker_info, 
  design = design, 
  contrast = contrast, 
  analysis_type = "DA", 
  seed_clustering = 123
)

## preparing data...

## transforming data...

## generating clusters...

## FlowSOM clustering completed in 6.7 seconds

## calculating features...

## calculating DA tests...

# display results for top DA clusters
head(topClusters(out_DA$res))

## DataFrame with 6 rows and 6 columns
##   cluster_id             logFC           logCPM               LR
##     <factor>         <numeric>        <numeric>        <numeric>
## 1         27 -4.65660517834391 12.9427269809787 216.261344152929
## 2         94  -3.6045220849935 13.3478248512472 197.290596122544
## 3          6 -4.96768550524626 13.1754232518767 182.726744338295
## 4         79 -5.14672338106299 12.6959108971985 148.651987281042
## 5         98  -4.5323224133241 12.1767213397505 136.100688454772
## 6         52 -3.05247579631149 13.2659255276742 132.522346367888
##                  p_val                p_adj
##              <numeric>            <numeric>
## 1 5.91440301360392e-49 5.91440301360391e-47
## 2 8.14899746295606e-45 4.07449873147803e-43
## 3 1.23051114578201e-41 4.10170381927338e-40
## 4 3.41673917598719e-34 8.54184793996798e-33
## 5 1.89664707647658e-31 3.79329415295316e-30
## 6 1.15004840686007e-30 1.91674734476678e-29

# calculate number of significant detected DA clusters at 10% false discovery 
# rate (FDR)
threshold <- 0.1
res_DA_all <- topClusters(out_DA$res, all = TRUE)
table(res_DA_all$p_adj <= threshold)

## 
## FALSE  TRUE 
##    23    77

# Test for differential states (DS) within clusters

# note: using default method 'diffcyt-DS-limma' and default parameters
# note: include random seed for reproducible clustering
out_DS <- diffcyt(
  d_input = d_flowSet, 
  experiment_info = experiment_info, 
  marker_info = marker_info, 
  design = design, 
  contrast = contrast, 
  analysis_type = "DS", 
  seed_clustering = 123, 
  plot = FALSE
)

## preparing data...

## transforming data...

## generating clusters...

## FlowSOM clustering completed in 6.3 seconds

## calculating features...

## calculating DS tests...

# display results for top DS cluster-marker combinations
head(topClusters(out_DS$res))

## DataFrame with 6 rows and 9 columns
##   cluster_id   marker          ID            logFC          AveExpr
##     <factor> <factor> <character>        <numeric>        <numeric>
## 1         80      pS6          80 2.90666271134132 2.27979487686086
## 2         99      pS6          99 2.70733271286747 2.35968733272817
## 3         90      pS6          90 2.32371127842426 2.75850533480623
## 4         99   pPlcg2          99 1.82819801669505 1.43253227739984
## 5         89      pS6          89 2.74458436384616 2.38936478394934
## 6         88      pS6          88 3.35978747883548 1.67477703096993
##                  t                p_val                p_adj
##          <numeric>            <numeric>            <numeric>
## 1 28.0982971194917 2.00632614550801e-11 2.80885660371121e-08
## 2  21.896790034549 2.81940312381832e-10 1.97358218667282e-07
## 3 19.5544339191185 9.27623188594239e-10 4.32890821343978e-07
## 4  17.168891891842  3.6178443913468e-09  1.0129964295771e-06
## 5 17.2866326692267 3.36909959252136e-09  1.0129964295771e-06
## 6  14.747677031564 1.74847387209903e-08 4.07977236823107e-06
##                  B
##          <numeric>
## 1 16.5325732055535
## 2 14.1125967440274
## 3 13.0760620282177
## 4 11.7112780918237
## 5  11.730411432573
## 6 10.1205460966827

# calculate number of significant detected DS cluster-marker combinations at 
# 10% false discovery rate (FDR)
threshold <- 0.1
res_DS_all <- topClusters(out_DS$res, all = TRUE)
table(res_DS_all$p_adj <= threshold)

## 
## FALSE  TRUE 
##   581   819

# Prepare data
d_se <- prepareData(d_flowSet, experiment_info, marker_info)

# Transform data
d_se <- transformData(d_se)

# Generate clusters
# note: include random seed for reproducible clustering
d_se <- generateClusters(d_se, seed_clustering = 123)

## FlowSOM clustering completed in 6 seconds

# Calculate cluster cell counts
d_counts <- calcCounts(d_se)

# Calculate cluster medians
d_medians <- calcMedians(d_se)

# Test for differential abundance (DA) of clusters
res_DA <- testDA_edgeR(d_counts, design, contrast)

# display results for top DA clusters
head(topClusters(res_DA))

## DataFrame with 6 rows and 6 columns
##   cluster_id             logFC           logCPM               LR
##     <factor>         <numeric>        <numeric>        <numeric>
## 1         27 -4.65660517834391 12.9427269809787 216.261344152929
## 2         94  -3.6045220849935 13.3478248512472 197.290596122544
## 3          6 -4.96768550524626 13.1754232518767 182.726744338295
## 4         79 -5.14672338106299 12.6959108971985 148.651987281042
## 5         98  -4.5323224133241 12.1767213397505 136.100688454772
## 6         52 -3.05247579631149 13.2659255276742 132.522346367888
##                  p_val                p_adj
##              <numeric>            <numeric>
## 1 5.91440301360392e-49 5.91440301360391e-47
## 2 8.14899746295606e-45 4.07449873147803e-43
## 3 1.23051114578201e-41 4.10170381927338e-40
## 4 3.41673917598719e-34 8.54184793996798e-33
## 5 1.89664707647658e-31 3.79329415295316e-30
## 6 1.15004840686007e-30 1.91674734476678e-29

# calculate number of significant detected DA clusters at 10% false discovery 
# rate (FDR)
threshold <- 0.1
table(topClusters(res_DA, all = TRUE)$p_adj <= threshold)

## 
## FALSE  TRUE 
##    23    77

# Test for differential states (DS) within clusters
res_DS <- testDS_limma(d_counts, d_medians, design, contrast, plot = FALSE)

# display results for top DS cluster-marker combinations
head(topClusters(res_DS))

## DataFrame with 6 rows and 9 columns
##   cluster_id   marker          ID            logFC          AveExpr
##     <factor> <factor> <character>        <numeric>        <numeric>
## 1         80      pS6          80 2.90666271134132 2.27979487686086
## 2         99      pS6          99 2.70733271286747 2.35968733272817
## 3         90      pS6          90 2.32371127842426 2.75850533480623
## 4         99   pPlcg2          99 1.82819801669505 1.43253227739984
## 5         89      pS6          89 2.74458436384616 2.38936478394934
## 6         88      pS6          88 3.35978747883548 1.67477703096993
##                  t                p_val                p_adj
##          <numeric>            <numeric>            <numeric>
## 1 28.0982971194917 2.00632614550801e-11 2.80885660371121e-08
## 2  21.896790034549 2.81940312381832e-10 1.97358218667282e-07
## 3 19.5544339191185 9.27623188594239e-10 4.32890821343978e-07
## 4  17.168891891842  3.6178443913468e-09  1.0129964295771e-06
## 5 17.2866326692267 3.36909959252136e-09  1.0129964295771e-06
## 6  14.747677031564 1.74847387209903e-08 4.07977236823107e-06
##                  B
##          <numeric>
## 1 16.5325732055535
## 2 14.1125967440274
## 3 13.0760620282177
## 4 11.7112780918237
## 5  11.730411432573
## 6 10.1205460966827

# calculate number of significant detected DS cluster-marker combinations at 
# 10% false discovery rate (FDR)
threshold <- 0.1
table(topClusters(res_DS, all = TRUE)$p_adj <= threshold)

## 
## FALSE  TRUE 
##   581   819

5.1 Overview

As described in our paper introducing the diffcyt framework (Weber et al., 2018), the results from a diffcyt analysis are presented to the user in the form of a set of significant detected high-resolution clusters (for DA tests) or cluster-marker combinations (for DS tests). The detected clusters or cluster-marker combinations can then be interpreted using visualizations; for example, to interpret the marker expression profiles in order to match detected clusters to known cell populations, or to group the high-resolution clusters into larger cell populations with a consistent phenotype.

Extensive plotting functions to generate both exploratory visualizations and visualizations of results from differential testing are available in the CATALYST package (Chevrier, Crowell, Zanotelli et al., 2018). Several of these plotting functions were originally developed by Malgorzata Nowicka for the CyTOF workflow available from Bioconductor (Nowicka et al., 2017), and have been adapted by Helena Crowell for inclusion in the CATALYST package. Additional plotting functions were developed during the development of the diffcyt package. Heatmaps are generated using the ComplexHeatmap Bioconductor package (Gu et al., 2016).

Here, we generate heatmaps to illustrate the results from the differential analyses above. Note that the CATALYST plotting functions can accept diffcyt results objects in either SummarizedExperiment format (from options 1 and 3 above) or CATALYST daFrame format (option 2).

For more examples of visualizations (in particular exploratory visualizations to explore the data prior to formal differential testing, including plots of the number of cells per sample, multi-dimensional scaling plots, and t-SNE plots), see the ‘Differential analysis with CATALYST’ vignette from the CATALYST package, available from the Bioconductor website.

5.2 Heatmap: DA test results

This heatmap illustrates the phenotypes (marker expression profiles) and signals of interest (cluster abundances by sample) for the top (most highly significant) detected clusters from the DA tests. See ?plotDiffHeatmap (from the CATALYST package) for more details.

Rows represent clusters, and columns represent protein markers (left panel) or samples (right panel). The left panel displays median (arcsinh-transformed) expression values across all samples for cell type markers, i.e. cluster phenotypes. The right panel displays the signal of interest: cluster abundances by sample (for the DA tests). The right annotation bar indicates clusters detected as significantly differential at an adjusted p-value threshold of 10%.

As mentioned previously, the DA tests are not particularly meaningful for the Bodenmiller_BCR_XL dataset, since the main signals of interest in this dataset are differential expression of pS6 and other signaling markers in B cells and several other cell populations. However, we include the plot here for illustrative purposes, to show how to use the functions.

Note: using plotHeatmap from the diffcyt package for now. This will be updated to use the plotting functions from the CATALYST package instead. See ?plotHeatmap for more details.

# Heatmap for top detected DA clusters

# note: use optional argument 'sample_order' to group samples by condition
sample_order <- c(seq(2, 16, by = 2), seq(1, 16, by = 2))

plotHeatmap(out_DA, analysis_type = "DA", sample_order = sample_order)

5.3 Heatmap: DS test results

This heatmap illustrates the phenotypes (marker expression profiles) and signals of interest (median expression of cell state markers by sample) for the top (most highly significant) detected cluster-marker combinations from the DS tests. See ?plotDiffHeatmap (from the CATALYST package) for more details.

Rows represent cluster-marker combinations, and columns represent protein markers (left panel) or samples (right panel). The left panel displays median (arcsinh-transformed) expression values across all samples for cell type markers, i.e. cluster phenotypes. The right panel displays the signal of interest: median expression of cell state markers by sample (for the DS tests). The right annotation bar indicates cluster-marker combinations detected as significantly differential at an adjusted p-value threshold of 10%.

The heatmap shows that the diffcyt pipeline has successfully recovered the main differential signal of interest in this dataset. As discussed above, the Bodenmiller_BCR_XL dataset contains known strong differential expression of several signaling markers (cell state markers) in several cell populations. In particular, the strongest signal is for differential expression of pS6 in B cells.

As expected, several of the top (most highly significant) detected cluster-marker combinations represent differential expression of pS6 (labels in right annotation bar) in B cells (identified by high expression of CD20, left panel). Similarly, the other top detected cluster-marker combinations shown in the heatmap correspond to other known strong differential signals in this dataset (see Nowicka et al., 2017, Figure 29; or the description of the results for dataset BCR-XL in our paper introducing the diffcyt framework (Weber et al., 2018).

Note: using plotHeatmap from the diffcyt package for now. This will be updated to use the plotting functions from the CATALYST package instead. See ?plotHeatmap for more details.

# Heatmap for top detected DS cluster-marker combinations

# note: use optional argument 'sample_order' to group samples by condition
sample_order <- c(seq(2, 16, by = 2), seq(1, 16, by = 2))

plotHeatmap(out_DS, analysis_type = "DS", sample_order = sample_order)