## ----setup, echo=FALSE, results="hide"------------------------------------------------------------
knitr::opts_chunk$set(
  tidy = FALSE,
  cache = TRUE,
  echo = TRUE,
  dev = c("png", "pdf"),
  package.startup.message = FALSE,
  message = FALSE,
  error = FALSE,
  warning = FALSE
)

options(width = 100)

## ----preprocess.data------------------------------------------------------------------------------
library(dreamlet)
library(muscat)
library(ExperimentHub)
library(scater)

# Download data, specifying EH2259 for the Kang, et al study
eh <- ExperimentHub()
sce <- eh[["EH2259"]]

sce$ind <- as.character(sce$ind)

# only keep singlet cells with sufficient reads
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
sce <- sce[, colData(sce)$multiplets == "singlet"]

# compute QC metrics
qc <- perCellQCMetrics(sce)

# remove cells with few or many detected genes
ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]

# set variable indicating stimulated (stim) or control (ctrl)
sce$StimStatus <- sce$stim

## ----aggregate------------------------------------------------------------------------------------
# Since 'ind' is the individual and 'StimStatus' is the stimulus status,
# create unique identifier for each sample
sce$id <- paste0(sce$StimStatus, sce$ind)

# Create pseudobulk data by specifying cluster_id and sample_id for aggregating cells
pb <- aggregateToPseudoBulk(sce,
  assay = "counts",
  cluster_id = "cell",
  sample_id = "id",
  verbose = FALSE
)

## ----crumblr--------------------------------------------------------------------------------------
library(crumblr)

# use dreamlet::cellCounts() to extract data
cellCounts(pb)[1:3, 1:3]

# Apply crumblr transformation
# cobj is an EList object compatable with limma workflow
# cobj$E stores transformed values
# cobj$weights stores precision weights
cobj <- crumblr(cellCounts(pb))

## ----vp-------------------------------------------------------------------------------------------
library(variancePartition)

fit <- dream(cobj, ~ StimStatus + ind, colData(pb))
fit <- eBayes(fit)

topTable(fit, coef = "StimStatusstim", number = Inf)

## ----session, echo=FALSE--------------------------------------------------------------------------
sessionInfo()