## ----warning=FALSE, message=FALSE---------------------------------------------
# load required packages
library(CATALYST)
library(cowplot)
library(flowCore)
library(diffcyt)
library(scater)
library(SingleCellExperiment)

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

## -----------------------------------------------------------------------------
(sce <- prepData(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 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")))

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

## ----eval=FALSE---------------------------------------------------------------
#  plotMDS(sce, color_by = "condition")

## ----echo=FALSE, fig.width=5, fig.height=4.5----------------------------------
CATALYST::plotMDS(sce, color_by = "condition")

## ----fig.width=10, fig.height=6-----------------------------------------------
plotExprHeatmap(sce, 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")
sce <- cluster(sce, features = lineage_markers, 
    xdim = 10, ydim = 10, maxK = 20, verbose = FALSE, seed = 1)       

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

## ----fig.width = 12, fig.height = 6-------------------------------------------
p <- plotMedExprs(sce, 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(sce, facet = "antigen", 
    group_by = "condition", shape_by = "patient_id")
p$facet$params$ncol <- 6
p

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

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

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

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

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

## ----run-umap-----------------------------------------------------------------
set.seed(1601)
sce <- runUMAP(sce, exprs_values = "exprs")

## -----------------------------------------------------------------------------
# view & access DRs
reducedDimNames(sce)
head(reducedDim(sce, "UMAP"))

## ----plot-umap, fig.width=8---------------------------------------------------
# color by marker expression & split by condition
plotDR(sce, color_by = "pS6") + facet_wrap("condition")

## ----fig.width=8--------------------------------------------------------------
# color by 8 metaclusters & split by sample ID
plotDR(sce, color_by = "meta8") + facet_wrap("sample_id", ncol = 4)

## ----fig.width=8--------------------------------------------------------------
u <- filterSCE(sce, patient_id == "Patient1")
table(u$sample_id)

u <- filterSCE(sce, k = "meta8",
    cluster_id %in% c(1, 3, 8))
plot_grid(
    plotDR(sce, color_by = "meta8"),
    plotDR(u, color_by = "meta8"))

## ----message=FALSE, warning=FALSE, fig.show='hide'----------------------------
# create design & constrast matrix
design <- createDesignMatrix(ei(sce), cols_design = "condition")
contrast <- createContrast(c(0, 1))

# test for
# - differential abundance (DA) of clusters
# - differential states (DS) within clusters
res_DA <- diffcyt(sce, clustering_to_use = "meta10",
    analysis_type = "DA", method_DA = "diffcyt-DA-edgeR",
    design = design, contrast = contrast)
res_DS <- diffcyt(sce, clustering_to_use = "meta10",
    analysis_type = "DS", method_DS = "diffcyt-DS-limma",
    design = design, contrast = contrast)

## ----fig.width=10, fig.height=5-----------------------------------------------
plotDiffHeatmap(sce, res_DA, all = TRUE, th = 0.05)

## ----fig.width=6, fig.height=6------------------------------------------------
plotDiffHeatmap(sce, res_DS, hm1 = FALSE, top_n = 20)

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