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

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

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 (library(HDCytoData); ?Bodenmiller_BCR_XL).

suppressPackageStartupMessages(library(HDCytoData))

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

suppressPackageStartupMessages(library(flowCore))

# check data format
d_flowSet

## A flowSet with 16 experiments.
## 
##   column names:
##   Time Cell_length CD3(110:114)Dd CD45(In115)Dd BC1(La139)Dd BC2(Pr141)Dd pNFkB(Nd142)Dd pp38(Nd144)Dd CD4(Nd145)Dd BC3(Nd146)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 BC4(Tb159)Dd CD14(Gd160)Dd pSlp76(Dy164)Dd BC5(Ho165)Dd pBtk(Er166)Dd pPlcg2(Er167)Dd pErk(Er168)Dd BC6(Tm169)Dd pLat(Er170)Dd IgM(Yb171)Dd pS6(Yb172)Dd HLA-DR(Yb174)Dd BC7(Lu175)Dd CD7(Yb176)Dd DNA-1(Ir191)Dd DNA-2(Ir193)Dd 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

# dimensions
dim(exprs(d_flowSet[[1]]))

## [1] 2838   39

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

##       Time Cell_length CD3(110:114)Dd CD45(In115)Dd BC1(La139)Dd
## [1,] 33073          30     120.823265      454.6009     576.8983
## [2,] 36963          35     135.106171      624.6824     564.6299
## [3,] 37892          30      -1.664619      601.0125    3077.2668
## [4,] 41345          58     115.290245      820.7125    6088.5967
## [5,] 42475          35      14.373802      326.6405    4606.6929
## [6,] 44620          31      37.737877      557.0137    4854.1519

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 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(110:114)Dd           CD3         type
## 4    CD45(In115)Dd          CD45         type
## 5     BC1(La139)Dd           BC1         none
## 6     BC2(Pr141)Dd           BC2         none
## 7   pNFkB(Nd142)Dd         pNFkB        state
## 8    pp38(Nd144)Dd          pp38        state
## 9     CD4(Nd145)Dd           CD4         type
## 10    BC3(Nd146)Dd           BC3         none
## 11   CD20(Sm147)Dd          CD20         type
## 12   CD33(Nd148)Dd          CD33         type
## 13 pStat5(Nd150)Dd        pStat5        state
## 14  CD123(Eu151)Dd         CD123         type
## 15   pAkt(Sm152)Dd          pAkt        state
## 16 pStat1(Eu153)Dd        pStat1        state
## 17  pSHP2(Sm154)Dd         pSHP2        state
## 18 pZap70(Gd156)Dd        pZap70        state
## 19 pStat3(Gd158)Dd        pStat3        state
## 20    BC4(Tb159)Dd           BC4         none
## 21   CD14(Gd160)Dd          CD14         type
## 22 pSlp76(Dy164)Dd        pSlp76        state
## 23    BC5(Ho165)Dd           BC5         none
## 24   pBtk(Er166)Dd          pBtk        state
## 25 pPlcg2(Er167)Dd        pPlcg2        state
## 26   pErk(Er168)Dd          pErk        state
## 27    BC6(Tm169)Dd           BC6         none
## 28   pLat(Er170)Dd          pLat        state
## 29    IgM(Yb171)Dd           IgM         type
## 30    pS6(Yb172)Dd           pS6        state
## 31 HLA-DR(Yb174)Dd        HLA-DR         type
## 32    BC7(Lu175)Dd           BC7         none
## 33    CD7(Yb176)Dd           CD7         type
## 34  DNA-1(Ir191)Dd         DNA-1         none
## 35  DNA-2(Ir193)Dd         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.

Note that in the example shown here, we include terms for group_id and patient_id in the design matrix: group_id is the factor of interest for the differential tests, and patient_id is included because this dataset contains paired samples from each patient. (For an unpaired dataset, only group_id would be included.)

suppressPackageStartupMessages(library(diffcyt))

# Create design matrix
# note: selecting columns containing group IDs and patient IDs (for an 
# unpaired dataset, only group IDs would be included)
design <- createDesignMatrix(
  experiment_info, cols_design = c("group_id", "patient_id")
)

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(). This function requires a single argument: a numeric vector defining the contrast. In many cases, this will simply be a vector of zeros and a single entry equal to one, which will test whether a single parameter is equal to zero. If a design matrix has been used, the entries correspond to the columns of the design matrix; if a model formula has been used, the entries correspond to the levels of the fixed effect terms.

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
nrow(contrast) == ncol(design)

## [1] TRUE

data.frame(parameters = colnames(design), contrast)

##           parameters contrast
## 1        (Intercept)        0
## 2     group_idBCR-XL        1
## 3 patient_idpatient2        0
## 4 patient_idpatient3        0
## 5 patient_idpatient4        0
## 6 patient_idpatient5        0
## 7 patient_idpatient6        0
## 8 patient_idpatient7        0
## 9 patient_idpatient8        0

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 9.2 seconds

## calculating features...

## calculating DA tests using method 'diffcyt-DA-edgeR'...

# display table of results for top DA clusters
topTable(out_DA, format_vals = TRUE)

## DataFrame with 20 rows and 3 columns
##     cluster_id     p_val     p_adj
##       <factor> <numeric> <numeric>
## 27          27  5.91e-49  5.91e-47
## 94          94  8.15e-45  4.07e-43
## 6            6  1.23e-41   4.1e-40
## 79          79  3.42e-34  8.54e-33
## 98          98   1.9e-31  3.79e-30
## ...        ...       ...       ...
## 22          22  2.47e-22  1.54e-21
## 92          92  6.41e-21  3.77e-20
## 71          71   1.6e-20  8.91e-20
## 96          96  7.03e-20   3.7e-19
## 100        100  1.57e-15  7.87e-15

# calculate number of significant detected DA clusters at 10% false discovery 
# rate (FDR)
threshold <- 0.1
res_DA_all <- topTable(out_DA, 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 7.3 seconds

## calculating features...

## calculating DS tests using method 'diffcyt-DS-limma'...

# display table of results for top DS cluster-marker combinations
topTable(out_DS, format_vals = TRUE)

## DataFrame with 20 rows and 4 columns
##     cluster_id marker_id     p_val     p_adj
##       <factor>  <factor> <numeric> <numeric>
## 80          80       pS6  2.01e-11  2.81e-08
## 99          99       pS6  2.82e-10  1.97e-07
## 90          90       pS6  9.28e-10  4.33e-07
## 99          99    pPlcg2  3.62e-09  1.01e-06
## 89          89       pS6  3.37e-09  1.01e-06
## ...        ...       ...       ...       ...
## 33          33      pBtk  7.33e-07  6.04e-05
## 35          35      pBtk  7.12e-07  6.04e-05
## 7            7      pAkt     1e-06  7.81e-05
## 99          99      pAkt  1.13e-06  8.35e-05
## 88          88     pNFkB  1.29e-06  9.01e-05

# calculate number of significant detected DS cluster-marker combinations at 
# 10% false discovery rate (FDR)
threshold <- 0.1
res_DS_all <- topTable(out_DS, 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 7.4 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 table of results for top DA clusters
topTable(res_DA, format_vals = TRUE)

## DataFrame with 20 rows and 3 columns
##     cluster_id     p_val     p_adj
##       <factor> <numeric> <numeric>
## 27          27  5.91e-49  5.91e-47
## 94          94  8.15e-45  4.07e-43
## 6            6  1.23e-41   4.1e-40
## 79          79  3.42e-34  8.54e-33
## 98          98   1.9e-31  3.79e-30
## ...        ...       ...       ...
## 22          22  2.47e-22  1.54e-21
## 92          92  6.41e-21  3.77e-20
## 71          71   1.6e-20  8.91e-20
## 96          96  7.03e-20   3.7e-19
## 100        100  1.57e-15  7.87e-15

# calculate number of significant detected DA clusters at 10% false discovery 
# rate (FDR)
threshold <- 0.1
table(topTable(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 table of results for top DS cluster-marker combinations
topTable(res_DS, format_vals = TRUE)

## DataFrame with 20 rows and 4 columns
##     cluster_id marker_id     p_val     p_adj
##       <factor>  <factor> <numeric> <numeric>
## 80          80       pS6  2.01e-11  2.81e-08
## 99          99       pS6  2.82e-10  1.97e-07
## 90          90       pS6  9.28e-10  4.33e-07
## 99          99    pPlcg2  3.62e-09  1.01e-06
## 89          89       pS6  3.37e-09  1.01e-06
## ...        ...       ...       ...       ...
## 33          33      pBtk  7.33e-07  6.04e-05
## 35          35      pBtk  7.12e-07  6.04e-05
## 7            7      pAkt     1e-06  7.81e-05
## 99          99      pAkt  1.13e-06  8.35e-05
## 88          88     pNFkB  1.29e-06  9.01e-05

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

## 
## FALSE  TRUE 
##   581   819

6.1 Overview

As described in our paper introducing the diffcyt framework (Weber et al., 2019), the results from a diffcyt differential analysis are provided to the user in the form of adjusted p-values, allowing the identification of sets of significant detected 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). These plotting functions are used in our CyTOF workflow (Nowicka et al., 2019) available from Bioconductor. For more details, including further examples of visualizations, see the CyTOF workflow or the CATALYST Bioconductor vignette. (Heatmaps are generated using the ComplexHeatmap Bioconductor package; Gu et al., 2016.)

Here, we generate heatmaps to illustrate the results from the differential analyses performed 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).

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

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: For the example below, we use an earlier version of the heatmap plotting function included in the diffcyt package, instead of using CATALYST. This is done to reduce dependencies, in order to streamline installation and compilation of the diffcyt package and vignette. Alternative code to generate the heatmap using CATALYST is also shown; this version includes additional formatting, and will usually be preferred.)

See ?plotDiffHeatmap (CATALYST) or ?plotHeatmap (diffcyt) 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)

# Heatmap for top detected DA clusters (alternative code using 'CATALYST')

suppressPackageStartupMessages(library(CATALYST))

plotDiffHeatmap(out_DA$d_se, out_DA$res)

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

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

(Note: For the example below, we use an earlier version of the heatmap plotting function included in the diffcyt package, instead of using CATALYST. This is done to reduce dependencies, in order to streamline installation and compilation of the diffcyt package and vignette. Alternative code to generate the heatmap using CATALYST is also shown; this version includes additional formatting, and will usually be preferred.)

See ?plotDiffHeatmap (CATALYST) or ?plotHeatmap (diffcyt) 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)

# Heatmap for top detected DA clusters (alternative code using 'CATALYST')

suppressPackageStartupMessages(library(CATALYST))

plotDiffHeatmap(out_DS$d_se, out_DS$res)