RunPrestoAll 替代FindAllMarkers 加速DE 搜索

Alex • 2023-06-12 07:48 • 杂文

$/,/,/,/,/,/,$ 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("")
Idents(cur_seu) <- ""

### 提交异步任务 
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)