## ----warning=FALSE---------------------------------------------------------
# load required packages
suppressPackageStartupMessages({
    library(CATALYST)
    library(flowCore)
    library(diffcyt)
    library(SummarizedExperiment)
})

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

## ----eval=FALSE------------------------------------------------------------
#  # 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)

## --------------------------------------------------------------------------
# construct daFrame
(daf <- daFrame(PBMC_fs, PBMC_panel, PBMC_md))

## ----eval=FALSE------------------------------------------------------------
#  # 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 daFrame
#  daFrame(PBMC_fs, panel2, md2,
#      panel_cols = list(channel = "channel_name", antigen = "marker"),
#      md_cols = list(file = "file", id = "sampleID",
#          factors = c("cond1", "cond2", "patientID")))

## ----fig.width=6, fig.height=3.5-------------------------------------------
n_cells(daf)
plotCounts(daf, color_by = "condition")

## ----fig.width=5, fig.height=4.5-------------------------------------------
plotMDS(daf, color_by = "condition")

## ----fig.width=10, fig.height=6--------------------------------------------
plotExprHeatmap(daf, bin_anno = TRUE, row_anno = TRUE)

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

## ----fig.width=5, fig.height=3---------------------------------------------
# access & render delta area plot
metadata(daf)$delta_area

## ----fig.width = 12, fig.height = 6----------------------------------------
p <- plotMedExprs(daf, k = "meta8", facet = "cluster_id", 
  group_by = "condition", shape_by = "patient_id")
p$facet$params$ncol <- 4
p

## ----fig.width = 12, fig.height = 8----------------------------------------
p <- plotMedExprs(daf, facet = "antigen", 
  group_by = "condition", shape_by = "patient_id")
p$facet$params$ncol <- 6
p

## ----message = FALSE, fig.width = 12, fig.height = 8-----------------------
plotClusterExprs(daf, k = "meta8", markers = "type")

## --------------------------------------------------------------------------
data(merging_table)
head(merging_table)
daf <- mergeClusters(daf, k = "meta20", table = merging_table, id = "merging1")
head(cluster_codes(daf))

## ----fig.width=8, fig.height=6---------------------------------------------
# median pS6 expression by sample as 2nd heatmap
plotClusterHeatmap(daf, hm2 = "pS6", k = "meta12", m = "meta6")

## ----fig.width=10, fig.height=6--------------------------------------------
# population frequencies by sample as 2nd heatmap
plotClusterHeatmap(daf, hm2 = "abundances", 
    draw_freqs = TRUE, cluster_anno = FALSE)

## ----fig.width=6, fig.height=4---------------------------------------------
plotAbundances(daf, k = "meta12", by = "sample_id", group_by = "condition")
plotAbundances(daf, k = "merging1", by = "cluster_id", 
  group_by = "condition", shape_by = "patient_id")

## --------------------------------------------------------------------------
# run PCA on all cells
set.seed(1)
daf <- runDR(daf, "PCA")

# run UMAP on 200 cells per sample
set.seed(2)
daf <- runDR(daf, "UMAP", rows_to_use = 200)

## --------------------------------------------------------------------------
# view & access DRs
reducedDimNames(daf)
reducedDims(daf)
head(reducedDim(daf, "PCA"))

# all cells
nrow(reducedDim(daf, "PCA")) == nrow(daf) 

# 200 per sample
nrow(reducedDim(daf, "UMAP")) == 200 * nlevels(sample_ids(daf))

## ----fig.width=7, fig.height=3---------------------------------------------
# color by marker expression & split by condition
plotDR(daf, "UMAP", color_by = "pS6", facet = "condition")

## ----fig.width=8, fig.height=4---------------------------------------------
# color by 8 metaclusters & split by sample ID
plotDR(daf, "UMAP", color_by = "meta8", facet = "sample_id")

## ----message=FALSE, warning=FALSE, fig.show='hide'-------------------------
# create design & constrast matrix
formula <- createFormula(
    experiment_info = PBMC_md, 
    cols_fixed = "condition", 
    cols_random = "patient_id")
contrast <- createContrast(c(0, 1))
# test for
# - differential abundance (DA) of clusters
# - differential states (DS) within clusters
res_DA <- diffcyt(daf, contrast = contrast, formula = formula, 
    analysis_type = "DA", method_DA = "diffcyt-DA-GLMM", clustering_to_use = "meta10")
res_DS <- diffcyt(daf, contrast = contrast, formula = formula, 
    analysis_type = "DS", method_DS = "diffcyt-DS-LMM", clustering_to_use = "meta10")

## ----fig.width=10, fig.height=6--------------------------------------------
plotDiffHeatmap(daf, res_DA, all = TRUE, th = 0.05)
plotDiffHeatmap(daf, res_DS, top_n = 20)

## --------------------------------------------------------------------------
sessionInfo()