lapply y do.call funcionando muy lentamente?

Question

Tengo un marco de datos de aproximadamente 35,000 filas, por 7 columnas. se parece a esto:

cabeza (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 es un factor, que tiene alrededor de 3.500 niveles únicos. Quiero, para cada nivel de gene_id obtener el min(start), max(end), mean(pctAT), mean(pctGC) y sum(length).

Intenté usar lapply y do.call para esto, pero me llevaría más de 30 minutos ejecutarlo. el código que estoy usando es:

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) 

Estoy seguro de que estoy haciendo algo mal para frenar esto. No he esperado a que termine, ya que estoy seguro de que puede ser más rápido. ¿Algunas ideas?

Answer 1

Desde que estoy en un estado de ánimo evangelizador ... esto es lo que la solución rápida data.table se vería así:

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 puede ser extremadamente lento en objetos grandes. Creo que esto se debe a cómo se construye la llamada, pero no estoy seguro. Una alternativa más rápida sería el paquete data.table. O, como @Andrie sugirió en un comentario, use tapply para cada cálculo y cbind los resultados.

Una nota sobre su implementación actual: en lugar de hacer la subconjunto en su función, puede usar la función split para dividir su data.frame en una lista de marcos de datos sobre los que puede pasar el bucle.

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

Como otros han mencionado - do.call tiene problemas con objetos de gran tamaño, y yo recientemente descubrió exactamente qué tan lento es en grandes conjuntos de datos. Para ilustrar el problema, aquí hay un benchamark mediante una simple llamada de resumen con un objeto de regresión grande (una regresión de Cox utilizando el rms-paquete):

> 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 

Mientras que la solución data.table es una excelente aproximación al anterior que no contiene la funcionalidad completa de do.call y, por lo tanto, pensé que compartiría mi función fastDoCall - una modificación de Hadley Wickhams suggested hack en la lista de distribución de correo. Está disponible en la versión 1.0 del paquete Gmisc (no publicada aún en CRAN, pero puede encontrarla en here). El punto de referencia es:

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

El código completo para la función es la siguiente:

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) 
}