celda.R
In celda: CEllular Latent Dirichlet Allocation

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE, dev = "png")

## ---- eval = TRUE, message = FALSE--------------------------------------------
library(celda)
library(SingleCellExperiment)
data(celdaCGSim, celdaCGMod, celdaCGGridSearchRes)
scegrid <- celdatosce(celdaCGGridSearchRes, celdaCGSim$counts)
scecg <- celdatosce(celdaCGMod, celdaCGSim$counts)

## ---- eval= FALSE-------------------------------------------------------------
#  if (!requireNamespace("BiocManager", quietly = TRUE)) {
#      install.packages("BiocManager")
#  }
#  BiocManager::install("celda")

## ---- eval = FALSE------------------------------------------------------------
#  help(package = celda)

## -----------------------------------------------------------------------------
simsce <- simulateCells("celda_CG",
    S = 5, K = 5, L = 10, G = 200, CRange = c(30, 50))

## -----------------------------------------------------------------------------
dim(assay(simsce, i = "counts"))

## -----------------------------------------------------------------------------
table(colData(simsce)$celda_cell_cluster)
table(colData(simsce)$celda_sample_label)

## -----------------------------------------------------------------------------
table(rowData(simsce)$celda_feature_module)

## ---- warning = FALSE, message = FALSE----------------------------------------
simsce <- selectFeatures(simsce)

## ---- warning = FALSE, message = FALSE----------------------------------------
sce <- celda_CG(x = simsce, K = 5, L = 10, verbose = FALSE, nchains = 1)

## -----------------------------------------------------------------------------
table(celdaClusters(sce), celdaClusters(simsce))
table(celdaModules(sce), celdaModules(simsce))

## -----------------------------------------------------------------------------
factorized <- factorizeMatrix(x = sce)
names(factorized)

## -----------------------------------------------------------------------------
dim(factorized$proportions$cell)
head(factorized$proportions$cell[, seq(3)], 5)

## -----------------------------------------------------------------------------
cellPop <- factorized$proportions$cellPopulation
dim(cellPop)
head(cellPop, 5)

## -----------------------------------------------------------------------------
dim(factorized$proportions$module)
head(factorized$proportions$module, 5)

## -----------------------------------------------------------------------------
topGenes <- topRank(matrix = factorized$proportions$module,
    n = 10, threshold = NULL)

## -----------------------------------------------------------------------------
topGenes$names$L1

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
plot(celdaHeatmap(sce = sce, nfeatures = 10))

## -----------------------------------------------------------------------------
sce <- celdaTsne(sce)

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
plotDimReduceCluster(x = sce, reducedDimName = "celda_tSNE")

plotDimReduceModule(x = sce, reducedDimName = "celda_tSNE", rescale = TRUE)

plotDimReduceFeature(x = sce, reducedDimName = "celda_tSNE",
    normalize = TRUE, features = "Gene_1")

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
celdaProbabilityMap(sce)

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
moduleHeatmap(sce, featureModule = seq(10), topCells = 100)

## -----------------------------------------------------------------------------
genes <- c(topGenes$names$L1, topGenes$names$L10)
geneIx <- which(rownames(altExp(sce)) %in% genes)
normCounts <- normalizeCounts(counts = counts(altExp(sce)), scaleFun = scale)

## ---- fig.width = 8, fig.height = 8-------------------------------------------
plot(plotHeatmap(counts = normCounts,
    z = celdaClusters(sce),
    y = celdaModules(sce),
    featureIx = geneIx,
    showNamesFeature = TRUE))

## ---- message=FALSE-----------------------------------------------------------
diffexpClust1 <- differentialExpression(sce, c1 = 1, c2 = NULL)

head(diffexpClust1, 5)

## ---- message=FALSE-----------------------------------------------------------
diffexpClust1vs2 <- differentialExpression(sce, c1 = 1, c2 = 2)

diffexpClust1vs2 <- diffexpClust1vs2[diffexpClust1vs2$FDR < 0.05 &
    abs(diffexpClust1vs2$Log2_FC) > 2, ]
head(diffexpClust1vs2, 5)

## -----------------------------------------------------------------------------
diffexpGeneIx <- which(rownames(altExp(sce)) %in% diffexpClust1vs2$Gene)
normCounts <- normalizeCounts(counts = counts(altExp(sce)), scaleFun = scale)

## -----------------------------------------------------------------------------
plot(plotHeatmap(counts = normCounts[, celdaClusters(sce) %in% c(1, 2)],
    clusterCell = TRUE,
    z = celdaClusters(sce)[celdaClusters(sce) %in% c(1, 2)],
    y = celdaModules(sce),
    featureIx = diffexpGeneIx,
    showNamesFeature = TRUE))

## ---- message = FALSE---------------------------------------------------------
moduleSplit <- recursiveSplitModule(simsce, initialL = 2, maxL = 15)

## -----------------------------------------------------------------------------
plotGridSearchPerplexity(moduleSplit)

## ---- message = FALSE---------------------------------------------------------
moduleSplitSelect <- subsetCeldaList(moduleSplit, params = list(L = 10))

cellSplit <- recursiveSplitCell(moduleSplitSelect,
    initialK = 3,
    maxK = 12,
    yInit = celdaModules(moduleSplitSelect))

## -----------------------------------------------------------------------------
plotGridSearchPerplexity(cellSplit)

## ---- eval = TRUE-------------------------------------------------------------
sce <- subsetCeldaList(cellSplit, params = list(K = 5, L = 10))

## ---- eval = TRUE, message = FALSE--------------------------------------------
cgs <- celdaGridSearch(simsce,
    paramsTest = list(K = seq(4, 6), L = seq(9, 11)),
    cores = 1,
    model = "celda_CG",
    nchains = 2,
    maxIter = 100,
    verbose = FALSE,
    bestOnly = TRUE)

## ---- eval = TRUE-------------------------------------------------------------
cgs <- resamplePerplexity(cgs, resample = 5)

## ---- eval = TRUE, fig.width = 8, fig.height = 8, warning = FALSE, message = FALSE----
plotGridSearchPerplexity(cgs)

## ---- eval = TRUE-------------------------------------------------------------
sce <- subsetCeldaList(cgs, params = list(K = 5, L = 10))

## ---- eval = FALSE, message=FALSE---------------------------------------------
#  cgs <- celdaGridSearch(simsce,
#      paramsTest = list(K = seq(4, 6), L = seq(9, 11)),
#      cores = 1,
#      model = "celda_CG",
#      nchains = 2,
#      maxIter = 100,
#      verbose = FALSE,
#      bestOnly = FALSE)
#  
#  cgs <- resamplePerplexity(cgs, celdaList = cgs, resample = 2)
#  
#  cgsK5L10 <- subsetCeldaList(cgs, params = list(K = 5, L = 10))
#  
#  sce <- selectBestModel(cgsK5L10)

## -----------------------------------------------------------------------------
featureModuleLookup(sce, feature = c("Gene_99"))

## -----------------------------------------------------------------------------
sceZRecoded <- recodeClusterZ(sce,
    from = c(1, 2, 3, 4, 5), to = c(2, 1, 3, 4, 5))

## -----------------------------------------------------------------------------
table(celdaClusters(sce), celdaClusters(sceZRecoded))

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