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

## ---- eval = TRUE, message = FALSE--------------------------------------------
library(celda)

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

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

## -----------------------------------------------------------------------------
dim(simCounts$counts)

## -----------------------------------------------------------------------------
table(simCounts$z)

## -----------------------------------------------------------------------------
table(simCounts$y)

## -----------------------------------------------------------------------------
table(simCounts$sampleLabel)

## ---- warning = FALSE, message = FALSE----------------------------------------
celdaModel <- celda_CG(counts = simCounts$counts,
    K = 5, L = 10, verbose = FALSE)

## -----------------------------------------------------------------------------
table(clusters(celdaModel)$z, simCounts$z)
table(clusters(celdaModel)$y, simCounts$y)

## -----------------------------------------------------------------------------
factorized <- factorizeMatrix(counts = simCounts$counts, celdaMod = celdaModel)
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------------------------------
celdaHeatmap(counts = simCounts$counts, celdaMod = celdaModel, nfeatures = 10)

## -----------------------------------------------------------------------------
tsne <- celdaTsne(counts = simCounts$counts, celdaMod = celdaModel)

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
plotDimReduceCluster(dim1 = tsne[, 1],
    dim2 = tsne[, 2],
    cluster = clusters(celdaModel)$z)

plotDimReduceModule(dim1 = tsne[, 1],
    dim2 = tsne[, 2],
    celdaMod = celdaModel,
    counts = simCounts$counts,
    rescale = TRUE)

plotDimReduceFeature(dim1 = tsne[, 1],
    dim2 = tsne[, 2],
    counts = simCounts$counts,
    features = "Gene_1")

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
celdaProbabilityMap(counts = simCounts$counts, celdaMod = celdaModel)

## ---- eval = TRUE, fig.width = 7, fig.height = 7------------------------------
moduleHeatmap(counts = simCounts$counts, celdaMod = celdaModel,
    featureModule = 1, topCells = 100)

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

## ---- fig.width = 8, fig.height = 8-------------------------------------------
plotHeatmap(counts = normCounts,
    z = clusters(celdaModel)$z,
    y = clusters(celdaModel)$y,
    featureIx = geneIx,
    showNamesFeature = TRUE)

## ---- message=FALSE-----------------------------------------------------------
diffexpClust1 <- differentialExpression(counts = simCounts$counts,
    celdaMod = celdaModel,
    c1 = 1,
    c2 = NULL)

head(diffexpClust1, 5)

## ---- message=FALSE-----------------------------------------------------------
diffexpClust1vs2 <- differentialExpression(
    counts = simCounts$counts,
    celdaMod = celdaModel,
    c1 = 1,
    c2 = 2)

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

## -----------------------------------------------------------------------------
diffexpGeneIx <- which(rownames(simCounts$counts) %in% diffexpClust1vs2$Gene)

normCounts <- normalizeCounts(counts = simCounts$counts, scaleFun = scale)

## -----------------------------------------------------------------------------
plotHeatmap(counts = normCounts[, clusters(celdaModel)$z %in% c(1, 2)],
    clusterCell = TRUE,
    z = clusters(celdaModel)$z[clusters(celdaModel)$z %in% c(1, 2)],
    y = clusters(celdaModel)$y,
    featureIx = diffexpGeneIx,
    showNamesFeature = TRUE)

## ---- message = FALSE---------------------------------------------------------
moduleSplit <- recursiveSplitModule(counts = simCounts$counts,
    initialL = 2, maxL = 15)

## -----------------------------------------------------------------------------
plotGridSearchPerplexity(celdaList = moduleSplit)

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

cellSplit <- recursiveSplitCell(counts = simCounts$counts,
    initialK = 3,
    maxK = 12,
    yInit = clusters(moduleSplitSelect)$y)

## -----------------------------------------------------------------------------
plotGridSearchPerplexity(celdaList = cellSplit)

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

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

## ---- eval = TRUE-------------------------------------------------------------
cgs <- resamplePerplexity(counts = simCounts$counts,
    celdaList = cgs, resample = 5)

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

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

## ---- message=FALSE-----------------------------------------------------------
cgs <- celdaGridSearch(simCounts$counts,
    paramsTest = list(K = seq(4, 6), L = seq(9, 11)),
    cores = 4,
    model = "celda_CG",
    nchains = 2,
    maxIter = 100,
    bestOnly = FALSE)

cgs <- resamplePerplexity(counts = simCounts$counts,
    celdaList = cgs,
    resample = 2)

cgsK5L10 <- subsetCeldaList(celdaList = cgs, params = list(K = 5, L = 10))

celdaModel1 <- selectBestModel(celdaList = cgsK5L10)

## -----------------------------------------------------------------------------
featureModuleLookup(counts = simCounts$counts, celdaMod = celdaModel,
    feature = c("Gene_99"))

## -----------------------------------------------------------------------------
celdaModelZRecoded <- recodeClusterZ(celdaMod = celdaModel,
    from = c(1, 2, 3, 4, 5), to = c(2, 1, 3, 4, 5))

## -----------------------------------------------------------------------------
table(clusters(celdaModel)$z, clusters(celdaModelZRecoded)$z)

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