为什么你画的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包重复某一文章或者别人的教程还是官网的示例时,发现自己画出来的图与原有的方向呈镜像或者上下颠倒,可以试着改一下这个随机种子。
- 区块链技术,如何提升网络安全?
- 趣店推“大白汽车”业务 启用域名dabaiqiche.com
- 糖果吃了那么多,你真的知道比特币分叉是咋回事吗?
- 静息态网络拓扑传输认知任务信息
- MYSQL字符串截取总结:LEFT、RIGHT、SUBSTRING、SUBSTRING
- 如何用Python提取中文关键词?
- 工信部:将加大网络提速降费力度加快百兆宽带普及
- 人工智能AI(5):线性代数之矩阵、线性空间
- iOS开发进阶篇——FRP与ReactiveCocoa的介绍(一)
- 英伟达修改GeForce软件使用条款:禁止在数据中心运行深度学习等应用
- 浅谈几种SLB技术的实现
- 史上最逼真人形机器人堪比健身教练,技能满满还会流汗
- 被监管前的疏忽?互联网金融大面积逾期,中介行为不容忽视
- 达尔文漏算的一步却让它填补,科学家认为人类最初认可的进化论不再适用
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法
- JUC学习之8锁现象
- TensorFlow惊现大bug?网友:这是逼着我们用PyTorch啊!
- 机器学习|用Q-Learning走迷宫
- C++设计模式笔记(03-02) - Template Method_模板方法(下)
- 从零开始学习华为路由交换 | OSPF修改计时器
- 【没落的985/211】Python爬取知乎8万字回答进行高校分析
- LinkedBlockingQueue源码学习
- 三歪吐血总结了各个中间件是如何实现持久化的
- ThreadPoolExecutor源码学习
- Docker六脉神剑(四) 使用Docker-Compose进行服务编排搭建lnmp环境
- 干的想喝水,一篇文章带你读懂硬盘工作原理!
- 微信小程序开发实战(11):滚动组件(picker)
- Docker六脉神剑 (五) Docker Swarm集群搭建及基础服务部署
- 思科模拟器GNS3将路由器变成交换机的方法
- docker安装nginx并配置https