为什么你画的Seurat包PCA图与别人的方向不一致？

Posted on by

事情是这个样子的，老板扔给我一篇《单细胞数据挖掘》文献要我重复这个文章中的结果，然后，就然后，我发现我画出来的PCA图与作者的方向颠倒了。如下所示：

image-20201013074041221

但是我看了看《单细胞天地》的优秀学员， 他的教程：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包重复某一文章或者别人的教程还是官网的示例时，发现自己画出来的图与原有的方向呈镜像或者上下颠倒，可以试着改一下这个随机种子。

​

Comments are closed.