## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## pre-load to avoid load messages in report
library(scRNAseq)
library(scater)
library(scran)
library(Glimma)
library(edgeR)

library(AnnotationHub)
setAnnotationHubOption("ASK", FALSE)

## -----------------------------------------------------------------------------
library(scRNAseq)
library(scater)
library(scran)
library(Glimma)
library(edgeR)

sce <- ZeiselBrainData(ensembl=TRUE)

## -----------------------------------------------------------------------------
sce <- logNormCounts(sce)

var_mod <- modelGeneVar(sce)
hvg_genes <- getTopHVGs(var_mod, n=500)
hvg_sce <- sce[hvg_genes, ]

hvg_sce <- logNormCounts(hvg_sce)

## -----------------------------------------------------------------------------
glimmaMDS(
    exprs(hvg_sce),
    groups = colData(hvg_sce)
)

## -----------------------------------------------------------------------------
colData(sce)$pb_group <-
    paste0(colData(sce)$level1class,
           "_",
           colData(sce)$level2class)

sce_counts <- counts(sce)
pb_counts <- t(rowsum(t(sce_counts), colData(sce)$pb_group))

pb_samples <- colnames(pb_counts)
pb_samples <- gsub("astrocytes_ependymal", "astrocytes-ependymal", pb_samples)
pb_split <- do.call(rbind, strsplit(pb_samples, "_"))
pb_sample_anno <- data.frame(
    sample = pb_samples,
    cell_type = pb_split[, 1],
    sample_group = pb_split[, 2]
)

## -----------------------------------------------------------------------------
pb_dge <- DGEList(
    counts = pb_counts,
    samples = pb_sample_anno,
    group = pb_sample_anno$cell_type
)

pb_dge <- calcNormFactors(pb_dge)

## -----------------------------------------------------------------------------
design <- model.matrix(~0 + cell_type, data = pb_dge$samples)
colnames(design) <- make.names(gsub("cell_type", "", colnames(design)))

pb_dge <- estimateDisp(pb_dge, design)

contr <- makeContrasts("pyramidal.SS - pyramidal.CA1", levels = design)

pb_fit <- glmFit(pb_dge, design)
pb_lrt <- glmLRT(pb_fit, contrast = contr)

## -----------------------------------------------------------------------------
glimmaMA(pb_lrt, dge = pb_dge)

## ----eval = FALSE-------------------------------------------------------------
# sc_dge <- DGEList(
#     counts = sce_counts,
#     group = colData(sce)$level1class
# )
# 
# sc_dge <- sc_dge[, colData(sce)$level1class %in% c("pyramidal CA1", "pyramidal SS")]
# 
# glimmaMA(
#     pb_lrt,
#     dge = sc_dge[, sample(1:ncol(sc_dge), 100)]
# )

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