## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  error    = FALSE,
  warning  = FALSE,
  eval     = TRUE,
  message  = FALSE
)

## ----loadlib------------------------------------------------------------------
library("DeeDeeExperiment")
library("ExperimentHub")
library("scater")
library("muscat")
library("limma")

## ----load_sce-----------------------------------------------------------------
# retrieve the data
eh <- ExperimentHub()
query(eh, "Kang")
sce <- eh[["EH2259"]]

## ----rm_undetected_genes------------------------------------------------------
# remove undetected genes
sce <- sce[rowSums(counts(sce) > 0) > 0, ]

## ----qc-----------------------------------------------------------------------
# calculate per-cell quality control (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]
dim(sce)

## ----rm_low_exp_genes---------------------------------------------------------
# remove lowly expressed genes
sce <- sce[rowSums(counts(sce) > 1) >= 10, ]
dim(sce)

## ----normalization------------------------------------------------------------
# compute sum-factors & normalize
sce <- computeLibraryFactors(sce)
sce <- logNormCounts(sce)

## ----prep---------------------------------------------------------------------
# data preparation
sce$id <- paste0(sce$stim, sce$ind)
(sce <- prepSCE(sce,
                kid = "cell", # subpopulation assignments
                gid = "stim",  # group IDs (ctrl/stim)
                sid = "id",   # sample IDs (ctrl/stim.1234)
                drop = TRUE))  # drop all other colData columns

## ----reddim-------------------------------------------------------------------
# compute UMAP using 1st 20 PCs
sce <- runUMAP(sce, pca = 20)

## ----aggregate_by_cell_type---------------------------------------------------
# aggregate by cell type
pb <- aggregateData(
  sce,
  assay = "counts",
  fun = "sum",
  by = c("cluster_id", "sample_id")
)

assayNames(pb)

## ----design-------------------------------------------------------------------
# construct design & contrast matrix
ei <- metadata(sce)$experiment_info
mm <- model.matrix(~ 0 + ei$group_id)
dimnames(mm) <- list(ei$sample_id, levels(ei$group_id))
contrast <- makeContrasts("stim-ctrl", levels = mm)

## ----pbDS---------------------------------------------------------------------
# run DS analysis
muscat_res <- pbDS(pb, design = mm, contrast = contrast)

names(muscat_res$table[["stim-ctrl"]])

## ----create_muscat_list-------------------------------------------------------
# preparing the results as muscat list
muscat_list <- muscat_list_for_dde(res = list(`stim-ctrl` = muscat_res),
                                   padj_col = "p_adj.loc")

## ----create_dde---------------------------------------------------------------
# create dde
dde <- DeeDeeExperiment(sce = sce,
                        de_results = muscat_list)
dde

## ----get_DEA_names------------------------------------------------------------
# check DEAs
getDEANames(dde)

## ----get_DEAs-----------------------------------------------------------------
# retrieve results
getDEA(dde,
       dea_name = "stim-ctrl_NK cells")|> head()

getDEA(dde,
       dea_name = "stim-ctrl_CD14+ Monocytes",
       format = "original") |> head()

## ----get_DEA_info-------------------------------------------------------------
# retrieving the DEA information
dea_name <- "stim-ctrl_B cells"
getDEAInfo(dde)[[dea_name]][["package"]]
getDEAInfo(dde)[[dea_name]][["package_version"]]

## ----add_scenario-------------------------------------------------------------
# adding info on the scenario under investigation
dde <- addScenarioInfo(dde,
                       dea_name = "stim-ctrl_Dendritic cells",
                       info = "This result contains the output of pseudobulk DE analysis performed on dendritic cells, comparing untreated samples to those stimulated with IFNβ")

## ----get_summary--------------------------------------------------------------
summary(dde, show_scenario_info = TRUE)

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