simExperiment.R
In bcTSNE: Projected t-SNE for Batch Correction

## ----echo=FALSE,results='hide'------------------------------------------------
library(knitr)
knitr::opts_knit$set(progress = TRUE, verbose = TRUE)
thm <- knit_theme$get("bclear")
knitr::knit_theme$set(thm)
knitr::opts_chunk$set(background = "#FCFCFC",
                      fig.width = 5,
                      fig.height = 3)

## ----warning=FALSE------------------------------------------------------------
pkgLoaded <- suppressPackageStartupMessages({
  c(require(bcTSNE),
    require(data.table),
    require(batchelor),
    require(kBET),
    require(splatter),
    require(scater),
    require(Rtsne),
    require(lisi),
    require(harmony),
    require(dlfUtils),
    require(xtable))
})
pkgLoaded <- all(pkgLoaded)
## Uncomment to install kBet & lisi
## Note: the packages will require compilation
# if (!require(drat)) {install.packages("drat"); library(drat)}
# drat::addRepo("daynefiler")
# install.packages(c("lisi", "kBET", "harmony", "dlfUtils"))

## ----eval=TRUE----------------------------------------------------------------
if (pkgLoaded) {
  p <- newSplatParams(seed = 1234,
                      batchCells = rep(200, 4), 
                      batch.facLoc = 0.2, 
                      batch.facScale = 0.1,
                      group.prob = c(0.1, 0.2, 0.3, 0.4),
                      de.facLoc = 0.1,
                      de.facScale = 0.4)
  sim <- splatSimulate(p, method = "groups", verbose = FALSE)
  sizeFactors(sim) <- librarySizeFactors(sim)
  sim <- logNormCounts(sim)
  sim <- logNormCounts(sim, log = FALSE)
  assay(sim, "centered") <- t(scale(t(normcounts(sim)), 
                                    center = TRUE, 
                                    scale = FALSE))
  Z <- model.matrix( ~ -1 + factor(colData(sim)$Batch))
  grp <- factor(sim$Group)
  bch <- as.integer(factor(sim$Batch))
}

## ----calcMet,echo=FALSE-------------------------------------------------------
if (pkgLoaded) {
  calcMetrics <- function(Y, bchLst) {
  calcSil <- function(x) {
    s <- batch_sil(pca.data = list(x = Y), batch = x, nPCs = 2)
    1 - abs(s)
  }
  calcKBET <- function(x) {
    kBET(Y, batch = x, do.pca = FALSE, plot = FALSE)$average
  }
  calcPCA <- function(x) {
    pcRegression(pca.data = prcomp(Y), batch = x, n_top = 2)$pcReg
  }
  sil  <- sapply(bchLst, calcSil)
  kbet <- sapply(bchLst, calcKBET)
  lisi <- compute_lisi(Y, 
                       meta_data = as.data.frame(bchLst),
                       label_colnames = names(bchLst))
  lisi <- colMeans(lisi)
  sizes <- sapply(bchLst, function(x) length(unique(x)))
  lisi <- (lisi - 1)/(sizes - 1)
  pca  <- sapply(bchLst, calcPCA)
  res <- list(sil = sil, kbet = kbet, lisi = lisi, pca = pca)
  do.call(cbind, res)
}
}

## ----pltFunc,echo=FALSE-------------------------------------------------------
if (pkgLoaded) {
  pltSimRes <- function(Y, lbl) {
  par(mar = c(4, 1, 2, 1) + 0.1)
  plot(Y, 
       ann = FALSE, 
       axes = FALSE, 
       bty = "n", 
       col = grp, 
       pch = bch)
  title(main = lbl)
  legend(x = grconvertX(0.2, "nfc"),
         y = line2user(2, 1),
         legend = rep(" ", 4),
         pch = 1:4,
         xpd = NA,
         bty = "n",
         horiz = TRUE,
         xjust = 0,
         yjust = 0.5)
  legend(x = grconvertX(0.2, "nfc"),
         y = line2user(3, 1),
         legend = rep(" ", 4),
         pch = 16,
         col = 1:4,
         xpd = NA,
         bty = "n",
         horiz = TRUE,
         xjust = 0,
         yjust = 0.5)
  text(x = grconvertX(0.2, "nfc"), 
       y = line2user(2:3, 1), 
       labels = c("Batch:", "Cell-type:"),
       xpd = NA,
       adj = c(1, 0.5))
}
}

## -----------------------------------------------------------------------------
res <- vector(mode = "list", length = 6)
names(res) <- c("btcc", "btlc", "hmlc", "hmcc", "mnn", "tsne")

## ----tsne---------------------------------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  res$tsne <- Rtsne(t(assay(sim, "centered")), inital_dims = 50)$Y
  pltSimRes(res$tsne, "t-SNE")
}

## ----bctsne-cc----------------------------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  res$btcc <- bctsne(t(assay(sim, "centered")), Z, k = 50)$Y
  pltSimRes(res$btcc, "bcTSNE-centered")
}

## ----bctsne-lc----------------------------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  res$btlc <- bctsne(t(logcounts(sim)), Z, k = 50)$Y
  pltSimRes(res$btlc, "bcTSNE-logcounts")
}

## ----harmony,message=FALSE----------------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  sim <- runPCA(sim, 50, exprs_values = "logcounts")
  sim <- RunHarmony(sim, group.by.vars = "Batch")
  res$hmlc = Rtsne(reducedDim(sim, "HARMONY"), pca = FALSE)$Y
  pltSimRes(res$hmlc, "Harmony-logcounts")
}

## ----harmony-centered,message=FALSE-------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  sim <- runPCA(sim, 50, exprs_values = "centered")
  sim <- RunHarmony(sim, group.by.vars = "Batch")
  res$hmcc = Rtsne(reducedDim(sim, "HARMONY"), pca = FALSE)$Y
  pltSimRes(res$hmcc, "Harmony-centered")
}

## ----mnn----------------------------------------------------------------------
if (pkgLoaded) {
  set.seed(1234)
  tmp <- mnnCorrect(sim, batch = factor(sim$Batch))
  res$mnn <- Rtsne(t(assay(tmp, "corrected")), initial_dims = 50)$Y
  rm(tmp)
  pltSimRes(res$mnn, "mnnRes")
}

## ----comp---------------------------------------------------------------------
if (pkgLoaded) {
  batchList <- list(batch = factor(sim$Batch), 
                    cellType = factor(sim$Group))
  met <- lapply(res, calcMetrics, bchLst = batchList)
}

## ----figs,echo=FALSE,results='hide'-------------------------------------------
if (pkgLoaded) {
  met <- lapply(met, round, 4)
  met <- lapply(met, as.data.table, keep.rownames = TRUE)
  for (i in seq_along(met)) met[[i]][ , alg := names(met)[i]]
  met <- rbindlist(met)
  met <- met[ , .(rn, alg, sil, kbet, lisi, pca)][order(rn, alg)]
  met[alg == "btcc", alg := "bcTSNE-centered"]
  met[alg == "btlc", alg := "bcTSNE-logcounts"]
  met[alg == "hmlc", alg := "Harmony-logcounts"]
  met[alg == "hmcc", alg := "Harmony-centered"]
  met[alg == "mnn", alg  := "MNN"]
  met[alg == "tsne", alg := "\\textit{t}-SNE"]
  met[rn == "batch", rn := "Batch"]
  met[rn == "cellType", rn := "Cell type"]
  met[ , tmp := seq_along(alg), by = rn]
  met[tmp != 1, rn := ""]
  met[ , tmp := NULL]
  setnames(met, c("", "", "SIL", "kBET", "iLSIS", "PcR"))
}

## ----echo=FALSE,results='asis'------------------------------------------------
if (pkgLoaded) {
  print(xtable(met, digits = 4),hline.after = c(-1, 0, 6, nrow(met)), 
        include.rownames = FALSE, 
        sanitize.text.function = I)
}

## -----------------------------------------------------------------------------
if (pkgLoaded) {
  calcMetrics <- function(Y, bchLst) {
  calcSil <- function(x) {
    s <- batch_sil(pca.data = list(x = Y), batch = x, nPCs = 2)
    1 - abs(s)
  }
  calcKBET <- function(x) {
    kBET(Y, batch = x, do.pca = FALSE, plot = FALSE)$average
  }
  calcPCA <- function(x) {
    pcRegression(pca.data = prcomp(Y), batch = x, n_top = 2)$pcReg
  }
  sil  <- sapply(bchLst, calcSil)
  kbet <- sapply(bchLst, calcKBET)
  lisi <- compute_lisi(Y, 
                       meta_data = as.data.frame(bchLst),
                       label_colnames = names(bchLst))
  lisi <- colMeans(lisi)
  sizes <- sapply(bchLst, function(x) length(unique(x)))
  lisi <- (lisi - 1)/(sizes - 1)
  pca  <- sapply(bchLst, calcPCA)
  res <- list(sil = sil, kbet = kbet, lisi = lisi, pca = pca)
  do.call(cbind, res)
}
}

## ----eval=FALSE,echo=FALSE,resuls='hide'--------------------------------------
#  ## For internal use, save the sim results to the untracked inst/notrack folder
#  # saveRDS(res, "../inst/notrack/sim.res")
#  # saveRDS(met, "../inst/notrack/sim.metrics")
#  # saveRDS(sim, "../inst/notrack/sim.sce")
#  # saveRDS(as.data.table(colData(sim)), "../inst/notrack/sim.cellMeta")