$\,\,\,\,\,\,$Seurat 提供了一个 FindAllMarkers
方法用于在单细胞RNA测序数据中寻找差异表达基因。然而,对于大型数据集的DE分析,使用Seurat软件包的FindAllMarkers方法 在数据集的全部细胞上执行DE搜索将变得非常缓慢。即使开启了future
并行,速度也并不会显著提升,因为多线程模式只能在某些特定的计算环节中发挥作用,而在其他计算环节中仍然需要使用单线程模式。在加速搜索差异基因, FindAllMarkers
函数本身提供了对基因或细胞进行预过滤的参数,特别是设置 max.cells.per.ident
对每个类别进行降采样,使其中用于比较的细胞数量不超过设置的阈值。虽然通常会出现一定的性能损失,但速度增加会比较显,并且差异表达最高的基因可能仍会排列在最前面。
如果计算平台允许,并不推荐使用
FindAllMarkers(max.cells.per.ident = *)
通过降采样加速来执行DE 搜索,而是推荐基于 Rcpp 实现的SeuratWrappers::RunPrestoAll
方法,具有与Seurat::FindAllMarkers
一样的功能,并提供更迅速的搜索。
1、FindAllMarkers 降采样DE搜索 --性能损失评估
降采样DE搜索代码示例
pacman::p_load(Seurat,dplyr,ggplot2,patchwork,future.apply)
plan(multisession(workers = 10));options(future.globals.maxSize = 100 * 1024^4)
### load data
cur_seu <- readRDS("</USER FILE/>")
Idents(cur_seu) <- "</USER CODE/>"
### 提交异步任务
R_bg.1 <- callr::r_bg(function(obt){FindAllMarkers(obt,only.pos = T)},args = list(cur_seu),package = c("Seurat"))
### 提交并行任务
future_lapply(c(30,50,70,seq(100,1000,by = 100)), function(sm.size){print(sm.size)
degs <- FindAllMarkers(cur_seu,only.pos = T,max.cells.per.ident = sm.size,verbose = F)
degs$sm.size = sm.size
return(degs)
}) %>% bind_rows() %>% saveRDS("downsample.Rds")
### 等待异步任务完成写出数据
while(R_bg.1$is_alive()){Sys.sleep(.2)}; R_bg.1$get_result() %>% saveRDS("sourcedata.Rds")
性能损失评估 :
将每个细胞身份降采样到500 进行DE分析的时候,大部分细胞类型的 top_10,30,50,100 差异基因开始显示出比较全细胞集的DEs较高的一致性。然而,即使将采样水平提到1k,一些可能受到批次等实验技术影响的细胞类型的 top50和100 差异基因仍会表现出比较全细胞较低的一致性(最差的一致性下探到 $\sim$90%)。因此,对于需要考虑较多差异基因数量的分析中,有必要注意慎用 降采样的DE搜索方式。
#### 汇总
downsample.rslt <- readRDS("downsample.Rds")
sourcedata.rslt <- readRDS("sourcedata.Rds")
lapply(setNames(c(10,30,50,100),c("top.10","top.30","top.50","top.100")), function(top_n){
temp_dat <- data.frame()
for (batch in c(30,50,70,seq(100,1000,by = 100))) {
for (i in as.vector(unique(dat.fullcells$cluster))) {
y = sourcedata.rslt %>% filter(cluster == i) %>% slice_max(order_by = avg_log2FC,n = top_n,with_ties = F) %>% pull(gene) %>% unique()
x = downsample.rslt %>% filter(sm.size == batch) %>% filter(cluster == i) %>% slice_max(order_by = avg_log2FC,n = top_n,with_ties = F) %>% pull(gene) %>% unique()
temp_dat <- rbind(temp_dat,data.frame(cell=i, overlap_score = length(intersect(x,y))/top_n, sample.size = batch))
}
}
return(temp_dat)
}) -> summary_rslt
lapply(1:length(summary_rslt), function(i){
summary_rslt[[i]] %>% ggplot() +
geom_boxplot(aes(x = factor(sample.size), y = overlap_score,group = sample.size,fill = factor(sample.size))) +
theme_bw() + labs(title = sprintf("%s consistance",names(summary_rslt)[i]), x = "max.cells.per.ident") + NoLegend()
}) -> plt
(plt[[1]] |plt[[2]])/(plt[[3]] |plt[[4]])
----
2、RunPrestoAll 加速DE 搜索
SeuratWrappers::RunPrestoAll
是 Seurat::FindAllMarkers
基于 Presto 的实现,功能参数和输出与 Seurat::FindAllMarkers
基本一致,由于 Presto 底层是基于 Rcpp 重编译的,可以非常迅速的进行 Wilcoxon 秩和检验。
实测在【Intel Core Processor (Broadwell) Linux】 平台上,6w 细胞量23类细胞类型的数据集上,开启10线程后,
Seurat::FindAllMarkers
函数完成搜索耗时约 40min,而SeuratWrappers::RunPrestoAll
最快可以在 10min 内完成搜索。
推荐代码
pacman::p_load(Seurat,SeuratWrappers,dplyr,future);
plan(multisession(workers = 10));sprintf("FUTURE CURRENT WORKERS = %s", nbrOfWorkers())
### First Load Your Seurat Object : cur_seu
SeuratWrappers::RunPrestoAll(cur_seu,only.pos = T,verbose =0)
![RunPrestoAll输出格式]
Reference
scRNA_workshop_part3_differential_expression (compbiocore.github.io)
How can we speed up FindMarkers · satijalab/seurat · Discussion #4433 (github.com)
RunPresto: A Presto-based implementation of FindMarkers that runs... in satijalab/seurat-wrappers: Community-Provided Methods and Extensions for the Seurat Object (rdrr.io)