查看原文
其他


学徒和学员已经陆续出师,是时候把生信技能树的舞台交给后辈了!

下面是转录组讲师实战单细胞的投稿

事情是这个样子的,老板扔给我一篇《单细胞数据挖掘》文献要我重复这个文章中的结果,然后,就然后,我发现我画出来的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包重复某一文章或者别人的教程还是官网的示例时,发现自己画出来的图与原有的方向呈镜像或者上下颠倒,可以试着改一下这个随机种子

您可能也对以下帖子感兴趣

文章有问题?点此查看未经处理的缓存