## ---- 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()