Nothing
## ----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")
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.