## ----echo=FALSE, results="hide", message=FALSE--------------------------------
require(knitr)
opts_chunk$set(error=FALSE, message=FALSE, warning=FALSE)

library(BiocStyle)
self <- Biocpkg("scrapper")

## -----------------------------------------------------------------------------
library(scRNAseq)
sce.z <- ZeiselBrainData()
sce.z

## -----------------------------------------------------------------------------
library(scrapper)
is.mito <- grep("^mt-", rownames(sce.z))
results <- analyze.se(
    sce.z,
    rna.qc.subsets=list(Mito=is.mito),
    num.threads=2 # limit number of threads to keep R CMD check happy.
)
results$x # most outputs are stored in a copy of the input SingleCellExperiment. 

## -----------------------------------------------------------------------------
library(scater)
plotReducedDim(results$x, "TSNE", colour_by="graph.cluster")

## -----------------------------------------------------------------------------
# Ordering by the median AUC across pairwise comparisons involving the cluster.
previewMarkers(results$markers$rna[[1]], order.by="auc.median")

## -----------------------------------------------------------------------------
test <- sce.z
test <- quickRnaQc.se(test, subsets=list(mt=is.mito), num.threads=2)
test <- test[,test$keep]
test <- normalizeRnaCounts.se(test, size.factors=test$sum)
test <- chooseRnaHvgs.se(test, num.threads=2)
test <- runPca.se(test, features=rowData(test)$hvg, num.threads=2)
test <- runAllNeighborSteps.se(test, num.threads=2)
markers <- scoreMarkers.se(test, groups=test$clusters, num.threads=2)

## -----------------------------------------------------------------------------
sce.t <- TasicBrainData()

# Only using the common genes in both datasets.
common <- intersect(rownames(sce.z), rownames(sce.t))
combined <- combineCols(sce.t[common,], sce.z[common,])
block <- rep(c("tasic", "zeisel"), c(ncol(sce.t), ncol(sce.z)))

## -----------------------------------------------------------------------------
# No mitochondrial genes in the combined set, so we'll just skip it.
blocked_res <- analyze.se(combined, block=block, num.threads=2)

## -----------------------------------------------------------------------------
plotReducedDim(blocked_res$x, "TSNE", colour_by="block")

## -----------------------------------------------------------------------------
table(blocked_res$x$graph.cluster, blocked_res$x$block)

## -----------------------------------------------------------------------------
previewMarkers(blocked_res$markers$rna[[1]])

## -----------------------------------------------------------------------------
aggregates <- aggregateAcrossCells.se(
    blocked_res$x,
    list(cluster=blocked_res$x$graph.cluster, dataset=blocked_res$x$block)
)
aggregates

## -----------------------------------------------------------------------------
sce.s <- StoeckiusHashingData(mode=c("human", "adt1", "adt2"))
sce.s <- sce.s[,1:5000]

# Combining the various ADT-related alternative experiments into a single object.
altExp(sce.s, "all.adts") <- rbind(altExp(sce.s, "adt1"), altExp(sce.s, "adt2"))
altExp(sce.s, "adt1") <- altExp(sce.s, "adt2") <- NULL

sce.s

## -----------------------------------------------------------------------------
is.mito <- grepl("^MT-", rownames(sce.s))
is.igg <- grepl("^IgG", rownames(altExp(sce.s, 'all.adts')))
multi_res <- analyze.se(
    sce.s,
    adt.altexp="all.adts",
    rna.qc.subsets=list(mito=is.mito),
    adt.qc.subsets=list(igg=is.igg),
    num.threads=2
)

## -----------------------------------------------------------------------------
plotReducedDim(multi_res$x, "UMAP", colour_by="graph.cluster")

## -----------------------------------------------------------------------------
# By default, markers are sorted by the mean Cohen's d.
# We just replace the column name when printing the dataframe here.
previewMarkers(multi_res$markers$adt[[5]], c(effect="cohens.d.mean"))

## -----------------------------------------------------------------------------
# Using the delta-detected to focus on silent-to-expressed genes.
previewMarkers(multi_res$markers$rna[[5]], order.by="delta.detected.mean")

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