lapply和do.call运行得很慢？

Question

我有一个数据框,大约35,000行,7列.它看起来像这样:

头(NUC)

  chr feature    start      end   gene_id    pctAT    pctGC length
1   1     CDS 67000042 67000051 NM_032291 0.600000 0.400000     10
2   1     CDS 67091530 67091593 NM_032291 0.609375 0.390625     64
3   1     CDS 67098753 67098777 NM_032291 0.600000 0.400000     25
4   1     CDS 67101627 67101698 NM_032291 0.472222 0.527778     72
5   1     CDS 67105460 67105516 NM_032291 0.631579 0.368421     57
6   1     CDS 67108493 67108547 NM_032291 0.436364 0.563636     55

gene_id是一个因子,有大约3,500个独特的水平.我想,对于gene_id的每个级别获得min(start),max(end),mean(pctAT),mean(pctGC),和sum(length).

我尝试使用lapply和do.call,但它需要持续+30分钟才能运行.我正在使用的代码是:

nuc_prof = lapply(levels(nuc$gene_id), function(gene){
  t = nuc[nuc$gene_id==gene, ]
  return(list(gene_id=gene, start=min(t$start), end=max(t$end), pctGC =
              mean(t$pctGC), pct = mean(t$pctAT), cdslength = sum(t$length))) 
})
nuc_prof = do.call(rbind, nuc_prof)

我确定我做错了什么来放慢速度.我没有等到它完成,因为我确信它可以更快.有任何想法吗？

Answer 1

因为我正在传福音......这就是快速data.table解决方案的样子:

library(data.table)
dt <- data.table(nuc, key="gene_id")

dt[,list(A=min(start),
         B=max(end),
         C=mean(pctAT),
         D=mean(pctGC),
         E=sum(length)), by=key(dt)]
#      gene_id        A        B         C         D   E
# 1: NM_032291 67000042 67108547 0.5582567 0.4417433 283
# 2:       ZZZ 67000042 67108547 0.5582567 0.4417433 283

Answer 2

do.call在大型物体上可能会非常慢.我认为这是由于它如何构建调用,但我不确定.一个更快的替代方案是data.table包.或者,正如@Andrie在评论中建议的那样,tapply用于每个计算和cbind结果.

关于当前实现的注释:您可以使用该split函数将data.frame分解为可以循环的data.frames列表,而不是在函数中进行子集化.

g <- function(tnuc) {
  list(gene_id=tnuc$gene_id[1], start=min(tnuc$start), end=max(tnuc$end),
       pctGC=mean(tnuc$pctGC), pct=mean(tnuc$pctAT), cdslength=sum(tnuc$length))
}
nuc_prof <- lapply(split(nuc, nuc$gene_id), g)

Answer 3

正如其他人所提到的 -do.call大对象有问题，我最近发现它在大数据集上到底有多慢。为了说明这个问题，这里有一个基准测试，它使用一个简单的摘要调用和一个大型回归对象（使用 rms-package 的 cox 回归）：

> model <- cph(Surv(Time, Status == "Cardiovascular") ~  
+              Group + rcs(Age, 3) + cluster(match_group),
+              data=full_df, 
+              x=TRUE, y=TRUE)

> system.time(s_reg <- summary(object = model))
   user  system elapsed 
   0.00    0.02    0.03 
> system.time(s_dc <- do.call(summary, list(object = model)))
   user  system elapsed 
 282.27    0.08  282.43 
> nrow(full_df)
[1] 436305

虽然该data.table解决方案是上述的一种很好的方法，但它不包含 的全部功能，do.call因此我认为我将分享我的fastDoCall功能 - Hadley Wickhams 的修改建议在 R 邮件列表上的hack。它在 Gmisc-package 1.0-version 中可用（尚未在 CRAN 上发布，但您可以在此处找到）。基准是：

> system.time(s_fc <- fastDoCall(summary, list(object = model)))
   user  system elapsed 
   0.03    0.00    0.06 

该函数的完整代码如下：

fastDoCall <- function(what, args, quote = FALSE, envir = parent.frame()){
  if (quote)
    args <- lapply(args, enquote)

  if (is.null(names(args))){
    argn <- args
    args <- list()
  }else{
    # Add all the named arguments
    argn <- lapply(names(args)[names(args) != ""], as.name)
    names(argn) <- names(args)[names(args) != ""]
    # Add the unnamed arguments
    argn <- c(argn, args[names(args) == ""])
    args <- args[names(args) != ""]
  }

  if (class(what) == "character"){
    if(is.character(what)){
      fn <- strsplit(what, "[:]{2,3}")[[1]]
      what <- if(length(fn)==1) {
        get(fn[[1]], envir=envir, mode="function")
      } else {
        get(fn[[2]], envir=asNamespace(fn[[1]]), mode="function")
      }
    }
    call <- as.call(c(list(what), argn))
  }else if (class(what) == "function"){
    f_name <- deparse(substitute(what))
    call <- as.call(c(list(as.name(f_name)), argn))
    args[[f_name]] <- what
  }else if (class(what) == "name"){
    call <- as.call(c(list(what, argn)))
  }

  eval(call,
       envir = args,
       enclos = envir)
}