单细胞转录组数据分析||Seurat并行策略

分享是一种态度

作者 |  周运来

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

随着单细胞技术的成熟，单细胞数据分析往往不再是单个组织样本，这有时候在计算（资源与时间）上是一个挑战。为此，Seurat也提供了可以探索的并行策略。鉴于入门单细胞数据分析的同事大多不是计算机出身，我们借助知乎的回答来解释一下什么是并行：

 1你吃饭吃到一半，电话来了，你一直到吃完了以后才去接，这就说明你不支持并发也不支持并行。
 2你吃饭吃到一半，电话来了，你停了下来接了电话，接完后继续吃饭，这说明你支持并发。
 3你吃饭吃到一半，电话来了，你一边打电话一边吃饭，这说明你支持并行。
 4
 5并发的关键是你有处理多个任务的能力，不一定要同时。
 6并行的关键是你有同时处理多个任务的能力。
 7所以我认为它们最关键的点就是：是否是『同时』。
 8
 9作者：知乎用户
10链接：https://www.zhihu.com/question/33515481/answer/58849148
11来源：知乎
12著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

在数据分析过程中，比如我们计算差异基因，其实是单个基因的计算，一般是算完一个再算下一个，并行的意思就是同时计算，以节约时间。

在Seurat中，我们选择使用future的并行化框架。在本文中，我们将演示如何从用户的角度利用某些Seurat函数的future实现。如果您有兴趣了解更多关于future框架的内容，请参阅这个R包文档以获得全面和详细的描述。

如何使用呢？

要访问Seurat中的函数的并行版本，需要加载future的包并设置plan。该plan将指定如何执行该函数。默认行为是以非并行的方式(顺序地)计算的。为了实现并行(异步)行为，我们通常推荐“多进程”策略。默认情况下，这将使用所有可用的内核，但是您可以设置workers参数来限制并发的数量。

首先应该检查你的计算机系统是否支持R的并行。

1library(future)
2# check the current active plan
3plan()
4
5sequential:
6- args: function (expr, envir = parent.frame(), substitute = TRUE, lazy = FALSE, seed = NULL, globals = TRUE, local = TRUE, earlySignal = FALSE, label = NULL, ...)
7- tweaked: FALSE
8- call: NULL

1# change the current plan to access parallelization
2plan("multiprocess", workers = 4)
3plan()
4
5## multiprocess:
6## - args: function (expr, envir = parent.frame(), substitute = TRUE, lazy = FALSE, seed = NULL, globals = TRUE, workers = 4, gc = FALSE, earlySignal = FALSE, label = NULL, ...)
7## - tweaked: TRUE
8## - call: plan("multiprocess", workers = 4)

‘Futurized’ functions in Seurat

再看看Seurat哪些函数可以使用并行。

编写以下函数是为了利用future的框架，如果当前plan设置正确，这些函数将被并行化。重要的是，调用函数的方式不变。

1NormalizeData
2ScaleData
3JackStraw
4FindMarkers
5FindIntegrationAnchors
6FindClusters - if clustering over multiple resolutions

例如，要运行findmarker的并行版本，只需设置计划并像往常一样调用该函数。

1library(Seurat)
2pbmc <- readRDS("../data/pbmc3k_final.rds")
3
4# Enable parallelization
5plan("multiprocess", workers = 4)
6markers <- FindMarkers(pbmc, ident.1 = "NK", verbose = FALSE)

并行的比较

在这里，我们将对具有并行化和不具有并行化的相同函数调用的运行时进行一个简短的比较。请注意，虽然我们预期使用并行化策略将减少上面列出的函数的运行时，但是这种减少的幅度将取决于许多因素(例如数据集的大小、工作人员的数量、系统的规格、未来的策略等)。下面的基准测试是在运行Ubuntu 16.04.5 LTS和Intel(R) Core(TM) i7-6800K CPU @ 3.40GHz和96 GB RAM的桌面计算机上进行的。

 1timing.comparisons <- data.frame(fxn = character(), time = numeric(), strategy = character())
 2plan("sequential")
 3start <- Sys.time()
 4pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)
 5end <- Sys.time()
 6timing.comparisons <- rbind(timing.comparisons, data.frame(fxn = "ScaleData", time = as.numeric(end - 
 7    start, units = "secs"), strategy = "sequential"))
 8
 9start <- Sys.time()
10markers <- FindMarkers(pbmc, ident.1 = "NK", verbose = FALSE)
11end <- Sys.time()
12timing.comparisons <- rbind(timing.comparisons, data.frame(fxn = "FindMarkers", time = as.numeric(end - 
13    start, units = "secs"), strategy = "sequential"))
14
15plan("multiprocess", workers = 4)
16start <- Sys.time()
17pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt", verbose = FALSE)
18end <- Sys.time()
19timing.comparisons <- rbind(timing.comparisons, data.frame(fxn = "ScaleData", time = as.numeric(end - 
20    start, units = "secs"), strategy = "multiprocess"))
21
22start <- Sys.time()
23markers <- FindMarkers(pbmc, ident.1 = "NK", verbose = FALSE)
24end <- Sys.time()
25timing.comparisons <- rbind(timing.comparisons, data.frame(fxn = "FindMarkers", time = as.numeric(end - 
26    start, units = "secs"), strategy = "multiprocess"))

来看看效果

1library(ggplot2)
2library(cowplot)
3ggplot(timing.comparisons, aes(fxn, time)) + geom_bar(aes(fill = strategy), stat = "identity", position = "dodge") + 
4    ylab("Time(s)") + xlab("Function") + theme_cowplot()

FAQ

• 我的进度条去哪了?
  不幸的是，当以任何并行计划模式运行这些函数时，您将丢失进度条。这是由于一些技术限制在未来的框架和R一般。如果您想监视函数的进度，就需要放弃并行化，而使用plan(“sequential”)。

• 如果我一直看到以下错误，我应该怎么做?

1Error in getGlobalsAndPackages(expr, envir = envir, globals = TRUE) : 
2  The total size of the X globals that need to be exported for the future expression ('FUN()') is X GiB. This exceeds the maximum allowed size of 500.00 MiB (option 'future.globals.maxSize'). The X largest globals are ... 

或者某些函数，每个worker需要访问某些全局变量。如果这些值大于默认限制，您将看到此错误。要解决这个问题，可以设置选项(future.globals。maxSize = X)，其中X是允许的最大大小(以字节为单位)。因此，要将其设置为1GB，需要运行options(future.globals。maxSize = 1000 * 1024^2)。请注意，这将增加您的RAM使用量，因此请谨慎地设置这个数字。

说到底Seurat不过是个R包，想要并行计算是要懂一些R里面并行原理的：由内而外释放R的力量||摘自《R大数据分析实用指南》。特别地，当我们在R中计算的中途突然发现某个任务报错说超出内存了，怎么办呢？

1使用更少的线程进行并行；
2如果你的电脑内存非常小，有一个简单的方法确定你的最大使用线程：max cores = memory.limit() / memory.size()；
3将大量的并行分小部分进行；
4在代码中多使用rm()删除没用的变量，使用gc()回收内存空间；

References

[1] 由内而外释放R的力量||摘自《R大数据分析实用指南》: https://www.jianshu.com/p/db95385b5340
[2] https://satijalab.org/seurat/v3.1/future_vignette.html

往期回顾

---

---

---

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程

• 生信爆款入门-第6期（线上直播4周，马拉松式陪伴，带你入门） 你的生物信息入门课

• 数据挖掘学习班第4期（线上直播3周，马拉松式陪伴，带你入门） 医学生/医生首选技能提高课

• 生信技能树的2019年终总结 你的生物信息成长宝藏

看完记得顺手点个“在看”哦！

生物 | 单细胞 | 转录组丨资料每天都精彩

长按扫码可关注