为什么你画的Seurat包PCA图与别人的方向不一致？
为什么你画的Seurat包PCA图与别人的方向不一致？
下面是转录组讲师实战单细胞的投稿
事情是这个样子的，老板扔给我一篇《单细胞数据挖掘》文献要我重复这个文章中的结果，然后，就然后，我发现我画出来的PCA图与作者的方向颠倒了。如下所示：
但是我看了看《单细胞天地》的优秀学员， 他的教程：Seurat包基本分析实战—文献图表复现，并没有遇到类似的问题。
其实吧，这个发现自己画出来的图与官方中的不一致，这种情况已经不是第一次了。多遇到了几次，就非常让人抓狂，老是一根刺卡在心里特别不舒服，想要一探究竟。还跟老板吐槽过，但是老板他也没遇到过。只能自己更生了。
老板也不想
后来有我们的《单细胞转录组CNS图表复现交流群》一位同行也遇到过，他告诉我可能是随机种子的原因，一下子就找到了方向不是。如果你也想加入交流群，自己去：你要的rmarkdown文献图表复现全套代码来了（单细胞）找到我们的拉群小助手哈。
插个话题：关于随机种子
set.seed：设置R的随机数生成器的种子，这对于创建可复制的模拟或随机对象非常有用。
举个例子，创造可复制的模拟价值。
set.seed(5)
rnorm(5)
[1] -0.84085548 1.38435934 -1.25549186 0.07014277 1.71144087
set.seed(5)
rnorm(5)
[1] -0.84085548 1.38435934 -1.25549186 0.07014277 1.71144087
随机数是一样的，不管我们在序列中走了多远，它们都是一样的。
Tip：在运行模拟时使用set.seed函数，以确保所有结果、图形等都是可复制的。
经过初步探索，发现将seed设置为NULL就可以与文章中的图一致：
后面我发现只要seed大于2就会相反，小于2设置为2，比如1或者-1等都可以保持一致，这就很诡异了，作者本身的默认值42难道不是为了给大家在运行这个结果的时候保持一致的结果用的么？
#不修改，图就相反，函数默认参数是seed.use = 42
gbm <- RunPCA(gbm, features = VariableFeatures(object = gbm))
# 修改seed.use会出现与文章中一致的PCA图
gbm <- RunPCA(gbm, features = VariableFeatures(object = gbm),seed.use = NULL)
两个图很明显的对比
首先找到RunPCA的脚本查看作者代码，看一下有么有什么随机因素导致：
代码地址：https://github.com/satijalab/seurat/blob/master/R/dimensional_reduction.R，很清楚的看到作者设置的默认参数：seed.use = 42 ，全部代码如下：
RunPCA.Seurat <- function(
object,
assay = NULL,
features = NULL,
npcs = 50,
rev.pca = FALSE,
weight.by.var = TRUE,
verbose = TRUE,
ndims.print = 1:5,
nfeatures.print = 30,
reduction.name = "pca",
reduction.key = "PC_",
seed.use = 42,
...
) {
assay <- assay %||% DefaultAssay(object = object)
assay.data <- GetAssay(object = object, assay = assay)
reduction.data <- RunPCA(
object = assay.data,
assay = assay,
features = features,
npcs = npcs,
rev.pca = rev.pca,
weight.by.var = weight.by.var,
verbose = verbose,
ndims.print = ndims.print,
nfeatures.print = nfeatures.print,
reduction.key = reduction.key,
seed.use = seed.use,
...
)
object[[reduction.name]] <- reduction.data
object <- LogSeuratCommand(object = object)
return(object)
}
上面的代码基本上没有给啥信息，一层层的跟套娃一样，还是得看里面的RunPCA函数的，而不是RunPCA.Seurat ，我们要看的是RunPCA.default函数，代码如下：
再看RunPCA.default的代码：
RunPCA.default <- function(
object,
assay = NULL,
npcs = 50,
rev.pca = FALSE,
weight.by.var = TRUE,
verbose = TRUE,
ndims.print = 1:5,
nfeatures.print = 30,
reduction.key = "PC_",
seed.use = 42,
approx = TRUE,
...
) {
if (!is.null(x = seed.use)) {
set.seed(seed = seed.use)
}
if (rev.pca) {
npcs <- min(npcs, ncol(x = object) - 1)
pca.results <- irlba(A = object, nv = npcs, ...)
total.variance <- sum(RowVar(x = t(x = object)))
sdev <- pca.results$d/sqrt(max(1, nrow(x = object) - 1))
if (weight.by.var) {
feature.loadings <- pca.results$u %*% diag(pca.results$d)
} else{
feature.loadings <- pca.results$u
}
cell.embeddings <- pca.results$v
}
else {
total.variance <- sum(RowVar(x = object))
if (approx) {
npcs <- min(npcs, nrow(x = object) - 1)
pca.results <- irlba(A = t(x = object), nv = npcs, ...)
feature.loadings <- pca.results$v
sdev <- pca.results$d/sqrt(max(1, ncol(object) - 1))
if (weight.by.var) {
cell.embeddings <- pca.results$u %*% diag(pca.results$d)
} else {
cell.embeddings <- pca.results$u
}
} else {
npcs <- min(npcs, nrow(x = object))
pca.results <- prcomp(x = t(object), rank. = npcs, ...)
feature.loadings <- pca.results$rotation
sdev <- pca.results$sdev
if (weight.by.var) {
cell.embeddings <- pca.results$x
} else {
cell.embeddings <- pca.results$x / (pca.results$sdev[1:npcs] * sqrt(x = ncol(x = object) - 1))
}
}
}
rownames(x = feature.loadings) <- rownames(x = object)
colnames(x = feature.loadings) <- paste0(reduction.key, 1:npcs)
rownames(x = cell.embeddings) <- colnames(x = object)
colnames(x = cell.embeddings) <- colnames(x = feature.loadings)
reduction.data <- CreateDimReducObject(
embeddings = cell.embeddings,
loadings = feature.loadings,
assay = assay,
stdev = sdev,
key = reduction.key,
misc = list(total.variance = total.variance)
)
if (verbose) {
msg <- capture.output(print(
x = reduction.data,
dims = ndims.print,
nfeatures = nfeatures.print
))
message(paste(msg, collapse = '\n'))
}
return(reduction.data)
}
seed.use = 42，我们看到了set.seed使用的地方
但是整个函数也没看出来哪里使用了随机功能呀？
if (!is.null(x = seed.use)) {
set.seed(seed = seed.use)
}
根据默认参数，PCA结果的运行的核心代码
pca.results <- irlba(A = t(x = object), nv = npcs, ...)
查看irlba函数
irlba是来自irlba包的一个函数（哭晕了，又是一个套娃）
irlba(A, nv = 5, nu = nv, maxit = 1000, work = nv + 7, reorth = TRUE,
tol = 1e-05, v = NULL, right_only = FALSE, verbose = FALSE,
scale = NULL, center = NULL, shift = NULL, mult = NULL,
fastpath = TRUE, svtol = tol, smallest = FALSE, ...)
看到这里，终于发现使用随机的是这个函数，随机参数为maxit
maxit：maximum number of iterations
但是发现这个函数最后使用的C/C++代码…
除了RunPCA函数之外，Seurat包中使用了随机种子的还有RunTSNE函数，默认为seed.use = 1，RunUMAP，默认为seed.use = 42，这两个函数再使用RunUMAP时回遇到画出来的图不一致，RunTSNE倒是没有遇见过很明显不一样的时候。
总结
挖到最后，发现还是有点说不通，没给找到一个合理的解释。
总之，如果你发现自己在使用Seurat包重复某一文章或者别人的教程还是官网的示例时，发现自己画出来的图与原有的方向呈镜像或者上下颠倒，可以试着改一下这个随机种子。