inst/doc/scMerge.R

In scMerge: scMerge: Merging multiple batches of scRNA-seq data

## ---- include = FALSE---------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

## ----loading packages, warning = FALSE, message = FALSE-----------------------
suppressPackageStartupMessages({
  library(SingleCellExperiment)
  library(scMerge)
  library(scater)
  })

## ----subsampling scMergeData, eval = FALSE, echo = FALSE----------------------
#  library(genefilter)
#  
#  load("~/Downloads/sce_mESC.rda")
#  data("segList_ensemblGeneID", package = "scMerge")
#  
#  set.seed(2019)
#  
#  example_sce = sce_mESC[, sce_mESC$batch %in% c("batch2", "batch3")]
#  example_sce$batch = droplevels(example_sce$batch)
#  batch2Sampled = sample(colnames(example_sce[,example_sce$batch == "batch2"]), 100)
#  batch3Sampled = sample(colnames(example_sce[,example_sce$batch == "batch3"]), 100)
#  
#  countsMat = SingleCellExperiment::counts(example_sce)
#  
#  batchTest = rowFtests(countsMat, fac = example_sce$batch)
#  celltypeTest = rowFtests(countsMat, fac = factor(example_sce$cellTypes))
#  
#  commonSegGenes = intersect(segList_ensemblGeneID$mouse$mouse_scSEG, rownames(sce_mESC))
#  
#  keepGenes = unique(c(commonSegGenes,
#  rownames(batchTest)[rank(batchTest$p.value) < 50],
#  rownames(celltypeTest)[rank(celltypeTest$p.value) < 250]
#  ))
#  
#  example_sce = example_sce[keepGenes, c(batch2Sampled, batch3Sampled)]
#  example_sce = example_sce[base::rowSums(counts(example_sce)) != 0, base::colSums(counts(example_sce)) != 0]
#  
#  table(example_sce$batch,
#  example_sce$cellTypes)
#  
#  dim(example_sce)
#  
#  example_sce = runPCA(example_sce, exprs_values = "logcounts")
#  scater::plotPCA(example_sce,
#  colour_by = "cellTypes",
#  shape_by = "batch")
#  
#  save(example_sce,
#  file = "data/example_sce.rda")

## ----loading data-------------------------------------------------------------
## Subsetted mouse ESC data
data("example_sce", package = "scMerge")

## ----checking raw data--------------------------------------------------------
example_sce = runPCA(example_sce, exprs_values = "logcounts")

scater::plotPCA(
  example_sce, 
  colour_by = "cellTypes", 
  shape_by = "batch")

## ----load SEG-----------------------------------------------------------------
## single-cell stably expressed gene list
data("segList_ensemblGeneID", package = "scMerge")
head(segList_ensemblGeneID$mouse$mouse_scSEG)

## ----t1, echo = FALSE---------------------------------------------------------
t1 = Sys.time()

## ----unsupervised_default, results='hide',fig.show='hide'---------------------
scMerge_unsupervised <- scMerge(
  sce_combine = example_sce, 
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3, 3),
  assay_name = "scMerge_unsupervised")

## ----t2, echo = FALSE---------------------------------------------------------
t2 = Sys.time()

## ----unsupervised_default_plotting--------------------------------------------
scMerge_unsupervised = runPCA(scMerge_unsupervised, exprs_values = "scMerge_unsupervised")
scater::plotPCA(
  scMerge_unsupervised, 
  colour_by = "cellTypes", 
  shape_by = "batch")

## ----unsupervised_prop1, results='hide',fig.show='hide'-----------------------
scMerge_unsupervised_all <- scMerge(
  sce_combine = example_sce, 
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3, 3),
  assay_name = "scMerge_unsupervised_all",
  replicate_prop = 1)

## ----unsupervised_prop1_plotting----------------------------------------------
scMerge_unsupervised_all = runPCA(scMerge_unsupervised_all,
                                  exprs_values = "scMerge_unsupervised_all")

scater::plotPCA(
  scMerge_unsupervised_all, 
  colour_by = "cellTypes", 
  shape_by = "batch")

## ----supervised, results='hide',fig.show='hide'-------------------------------
scMerge_supervised <- scMerge(
  sce_combine = example_sce,
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3, 3),
  assay_name = "scMerge_supervised",
  cell_type = example_sce$cellTypes)

## ----supervised_plotting------------------------------------------------------
scMerge_supervised = runPCA(scMerge_supervised,
                            exprs_values = "scMerge_supervised")

scater::plotPCA(
  scMerge_supervised,
  colour_by = "cellTypes",
  shape_by = "batch")

## ----semi_supervised1, results='hide',fig.show='hide'-------------------------
scMerge_semisupervised1 <- scMerge(
  sce_combine = example_sce,
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3,3),
  assay_name = "scMerge_semisupervised1",
  cell_type = example_sce$cellTypes,
  cell_type_inc = which(example_sce$cellTypes == "2i"), 
  cell_type_match = FALSE)

## ----semi_supervised1_plotting------------------------------------------------
scMerge_semisupervised1 = runPCA(scMerge_semisupervised1,
                                 exprs_values = "scMerge_semisupervised1")

scater::plotPCA(
  scMerge_semisupervised1, 
  colour_by = "cellTypes", 
  shape_by = "batch")

## ----semi_supervised2, results='hide',fig.show='hide'-------------------------
scMerge_semisupervised2 <- scMerge(
  sce_combine = example_sce,
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3, 3),
  assay_name = "scMerge_semisupervised2",
  cell_type = example_sce$cellTypes,
  cell_type_inc = NULL,
  cell_type_match = TRUE)

## ----semi_supervised2_plotting------------------------------------------------
scMerge_semisupervised2 = runPCA(scMerge_semisupervised2,
                                 exprs_values = "scMerge_semisupervised2")

scater::plotPCA(
  scMerge_semisupervised2, 
  colour_by = "cellTypes", 
  shape_by = "batch")

## ----segIndex1, eval = FALSE--------------------------------------------------
#  exprs_mat = SummarizedExperiment::assay(example_sce, 'counts')
#  result = scSEGIndex(exprs_mat = exprs_mat)

## ----segIndex2----------------------------------------------------------------
## SEG list in ensemblGene ID
data("segList_ensemblGeneID", package = "scMerge") 
## SEG list in official gene symbols
data("segList", package = "scMerge")

## SEG list for human scRNA-Seq data
head(segList$human$human_scSEG)

## SEG list for human bulk microarray data
head(segList$human$bulkMicroarrayHK)

## SEG list for human bulk RNASeq data
head(segList$human$bulkRNAseqHK)

## ----t3, echo = FALSE---------------------------------------------------------
t3 = Sys.time()

## ----computation_fast, results='hide',fig.show='hide'-------------------------
library(BiocSingular)
scMerge_fast <- scMerge(
  sce_combine = example_sce, 
  ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
  kmeansK = c(3, 3),
  assay_name = "scMerge_fast", 
  BSPARAM = IrlbaParam(), 
  svd_k = 20)

## ----t4, echo = FALSE---------------------------------------------------------
t4 = Sys.time()

## ----computation_svd_plotting-------------------------------------------------
paste("Normally, scMerge takes ", round(t2 - t1, 2), " seconds")
paste("Fast version of scMerge takes ", round(t4 - t3, 2), " seconds")

scMerge_fast = runPCA(scMerge_fast, exprs_values = "scMerge_fast")

scater::plotPCA(
  scMerge_fast, 
  colour_by = "cellTypes", 
  shape_by = "batch") +
  labs(title = "fast scMerge yields similar results to the default version")

## ----parallel1, eval = FALSE--------------------------------------------------
#  library(BiocParallel)
#  scMerge_parallel <- scMerge(
#    sce_combine = example_sce,
#    ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
#    kmeansK = c(3, 3),
#    assay_name = "scMerge_parallel",
#    BPPARAM = MulticoreParam(workers = 2)
#  )

## ---- eval = FALSE------------------------------------------------------------
#  library(Matrix)
#  library(DelayedArray)
#  
#  sparse_input = example_sce
#  
#  assay(sparse_input, "counts") = as(counts(sparse_input), "dgeMatrix")
#  assay(sparse_input, "logcounts") = as(logcounts(sparse_input), "dgeMatrix")
#  
#  scMerge_sparse = scMerge(
#    sce_combine = sparse_input,
#    ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
#    kmeansK = c(3, 3),
#    assay_name = "scMerge_sparse")

## ---- eval = FALSE------------------------------------------------------------
#  library(HDF5Array)
#  library(DelayedArray)
#  
#  DelayedArray:::set_verbose_block_processing(TRUE) ## To monitor block processing
#  
#  hdf5_input = example_sce
#  
#  assay(hdf5_input, "counts") = as(counts(hdf5_input), "HDF5Array")
#  assay(hdf5_input, "logcounts") = as(logcounts(hdf5_input), "HDF5Array")
#  
#  scMerge_hdf5 = scMerge(
#    sce_combine = sparse_input,
#    ctl = segList_ensemblGeneID$mouse$mouse_scSEG,
#    kmeansK = c(3, 3),
#    assay_name = "scMerge_hdf5")

## ----reference----------------------------------------------------------------
citation("scMerge")

## ----session info-------------------------------------------------------------
sessionInfo()