## ----load-example-data, message = FALSE--------------------------------------- library(Cepo) library(SingleCellExperiment) data("cellbench", package = "Cepo") cellbench cellbench = cellbench[!duplicated(rownames(cellbench)),] ## ----visualize colData-------------------------------------------------------- colData(cellbench)[1:5,] ## ----check_cell_types--------------------------------------------------------- unique(cellbench$celltype) ## ----differential-stability-analysis------------------------------------------ ds_res = Cepo(exprsMat = logcounts(cellbench), cellType = cellbench$celltype) ## ----------------------------------------------------------------------------- ds_res ## ----filtering-lowly-expressed-genes------------------------------------------ ds_res_zprop = Cepo::Cepo(exprsMat = logcounts(cellbench), cellTypes = cellbench$celltype, exprsPct = 0.5) ## ----filtering-logfc---------------------------------------------------------- ds_res_logfc = Cepo(exprsMat = logcounts(cellbench), cellTypes = cellbench$celltype, logfc = 1) ## ----------------------------------------------------------------------------- nrow(ds_res$stats) nrow(ds_res_zprop$stats) nrow(ds_res_logfc$stats) ## ----return-p-values1--------------------------------------------------------- ds_res_pvalues = Cepo(exprsMat = logcounts(cellbench), cellType = cellbench$celltype, computePvalue = 200, prefilter_pzero = 0.4) ## ----------------------------------------------------------------------------- idx = rownames(ds_res_pvalues$stats) par(mfrow=c(1,3)) for (i in unique(cellbench$celltype)) { plot(rank(ds_res_pvalues$stats[[i]]), rank(-log10(ds_res_pvalues$pvalues[idx, i])), main = i, xlab = "rank Cepo statistics", ylab = "rank -log10 p-values") } par(mfrow=c(1,1)) ## ----return-p-values2--------------------------------------------------------- ds_res_pvalues = Cepo(exprsMat = logcounts(cellbench), cellType = cellbench$celltype, # we use a low value for demonstration purposes computePvalue = 100, computeFastPvalue = FALSE) ds_res_pvalues ## ----upset-plot--------------------------------------------------------------- library(UpSetR) res_name = topGenes(object = ds_res, n = 500) upset(fromList(res_name), nsets = 3) ## ----plot-densities----------------------------------------------------------- plotDensities(x = cellbench, cepoOutput = ds_res, nGenes = 2, assay = "logcounts", celltypeColumn = "celltype") ## ----plot-densities-genes----------------------------------------------------- plotDensities(x = cellbench, cepoOutput = ds_res, genes = c("PLTP", "CPT1C", "MEG3", "SYCE1", "MICOS10P3", "HOXB7"), assay = "logcounts", celltypeColumn = "celltype") ## ----------------------------------------------------------------------------- data("sce_pancreas", package = "Cepo") sce_pancreas ## ----fig.height=5, fig.width=5------------------------------------------------ library(scater) sce = sce_pancreas sce = scater::logNormCounts(sce) sce = scater::runPCA(sce) scater::plotPCA( sce, colour_by = "cellTypes", shape_by = "batch") ## ----------------------------------------------------------------------------- library(scMerge) data("segList", package = "scMerge") head(segList$human$human_scSEG) corrected <- scMerge( sce_combine = sce, ctl = segList$human$human_scSEG, kmeansK = c(2, 2), assay_name = "scMerge", cell_type = sce$cellTypes) ## ----fig.height=5, fig.width=5------------------------------------------------ corrected = runPCA(corrected, exprs_values = "scMerge") scater::plotPCA( corrected, colour_by = "cellTypes", shape_by = "batch") ## ----------------------------------------------------------------------------- ds_res = Cepo::Cepo(exprsMat = assay(corrected, "scMerge"), cellType = corrected$cellTypes) ## ----------------------------------------------------------------------------- ds_res_batches = Cepo::Cepo(exprsMat = logcounts(sce), cellTypes = sce$cellTypes, block = sce$batch, minCelltype = 2) ## ----------------------------------------------------------------------------- names(ds_res_batches) ## ----fig.height=10, fig.width=10---------------------------------------------- idx = Reduce(intersect, lapply(ds_res_batches, function(x) names(x$stats[, 1]))) combinedRes = as.data.frame(do.call(cbind, lapply(ds_res_batches, function(x) x$stats[idx,] ))) batches = rep(names(ds_res_batches), sapply(ds_res_batches, function(x) length(x$stats))) cty = unlist(lapply(ds_res_batches, function(x) names(x$stats)), use.name = FALSE) colnames(combinedRes) = gsub("[.]", "_", colnames(combinedRes)) annot = data.frame( batch = batches, celltype = cty ) rownames(annot) = colnames(combinedRes) pheatmap::pheatmap(cor(combinedRes), annotation = annot) ## ----marker-genes, fig.height=10, fig.width=15-------------------------------- cepo_genes = Cepo::topGenes(ds_res_batches$average, n = 3) markersPlot = lapply(cepo_genes, function(x) { pp = lapply(x, function(gene) { p = scater::plotPCA( corrected, colour_by = gene, shape_by = "cellTypes") return(p) }) pp = patchwork::wrap_plots(pp, ncol = 3) + patchwork::plot_layout(guides = "auto") return(pp) }) patchwork::wrap_plots(markersPlot, nrow = 2) ## ----gsea, message = FALSE---------------------------------------------------- library(escape) library(fgsea) hallmark <- getGeneSets(species = "Homo sapiens", library = "H") hallmarkList <- geneIds(hallmark) fgseaRes <- fgsea(pathways = hallmarkList, stats = sort(ds_res_batches$average$stats[,"beta"]), minSize = 15, maxSize = 500) enriched_beta <- -log10(fgseaRes[order(pval), "padj"][[1]]) names(enriched_beta) <- fgseaRes[order(pval), "pathway"][[1]] ## ----------------------------------------------------------------------------- enriched_beta[1:5] ## ----fig.height=5, fig.width=5------------------------------------------------ plotEnrichment(hallmarkList[["HALLMARK_PANCREAS_BETA_CELLS"]], sort(ds_res_batches$average$stats[, "beta"])) + labs(title="HALLMARK_PANCREAS_BETA_CELLS") ## ----cepo-delayed------------------------------------------------------------- library(DelayedArray) library(HDF5Array) da_matrix = DelayedArray(realize(logcounts(cellbench), "HDF5Array")) class(da_matrix) class(seed(da_matrix)) da_output = Cepo(exprsMat = da_matrix, cellType = cellbench$celltype) ## ----cepo-parallel------------------------------------------------------------ library(BiocParallel) BPPARAM = if (.Platform$OS.type == "windows") { BiocParallel::SnowParam(workers = 2) } else { BiocParallel::MulticoreParam(workers = 2) } DelayedArray::setAutoBPPARAM(BPPARAM = BPPARAM) ## Setting two cores for computation da_output_parallel = Cepo(exprsMat = da_matrix, cellTypes = cellbench$celltype) DelayedArray::setAutoBPPARAM(BPPARAM = SerialParam()) ## Revert back to only one core ## ----sessionInfo-------------------------------------------------------------- sessionInfo()