lapply e do.call funzionano molto lentamente?

Question

Ho un frame di dati che è di circa 35.000 righe, di 7 colonne. sembra che questo:

testa (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 è un fattore, che ha circa 3.500 livelli unici. Voglio, per ogni livello di gene_id ottenere min(start), max(end), mean(pctAT), mean(pctGC) e sum(length).

Ho provato ad usare lapply e do.call per questo, ma ci vogliono sempre +30 minuti per funzionare. il codice che sto utilizzando è:

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) 

Sono certo che sto facendo qualcosa di sbagliato per rallentare questo giù. Non ho aspettato che finisse perché sono sicuro che può essere più veloce. Qualche idea?

Answer 1

Dal momento che sono in uno stato d'animo di evangelizzazione ... ecco cosa la soluzione rapida data.table sarebbe simile:

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 può essere estremamente lento su oggetti di grandi dimensioni. Penso che questo sia dovuto al modo in cui costruisce la chiamata, ma non ne sono sicuro. Un'alternativa più veloce sarebbe il pacchetto data.table. Oppure, come suggerito da @Andrie in un commento, usa tapply per ogni calcolo e cbind i risultati.

Una nota sull'implementazione corrente: anziché eseguire la subsetting nella funzione, è possibile utilizzare la funzione split per suddividere il data.frame in un elenco di data.frames che è possibile eseguire il loopover.

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

Come altri hanno detto - do.call ha problemi con oggetti di grandi dimensioni, e di recente ho ha scoperto esattamente quanto è lento su set di dati di grandi dimensioni. Per illustrare il problema, ecco un benchamark utilizzando una semplice chiamata sintesi con un oggetto di regressione di grandi dimensioni (una regressione di Cox con il RMS-package):

> 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 

Mentre la soluzione data.table è un ottimo approccio al di sopra di essa non contiene la piena funzionalità del do.call e ho quindi pensato che condividerò la mia funzione fastDoCall - una modifica di Hadley Wickhams suggested hack sulla mailing list R. È disponibile nella versione 1.0 del pacchetto Gmisc (non ancora rilasciato su CRAN ma lo puoi trovare here). Il riferimento è:

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

Il codice completo per la funzione è la seguente:

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