R Optimierung: Wie kann ich eine For-Schleife in dieser Situation vermeiden?

Question

Ich versuche eine einfache genomische Track-Schnittmenge in R zu erstellen, und stürze mich in große Performance-Probleme, die wahrscheinlich mit der Verwendung von for-Schleifen zusammenhängen.

In dieser Situation habe ich vordefinierte Fenster in Abständen von 100 bp und ich versuche zu berechnen, wie viel von jedem Fenster durch die Anmerkungen in Mylist abgedeckt ist. Grafisch sieht es etwa so aus:

  0 100 200 300 400 500 600 
windows: |-----|-----|-----|-----|-----|-----| 

mylist: |-| |-----------| 

Also habe ich einige Code schrieb nur das zu tun, aber es ist ziemlich langsam und hat sich zu einem Engpass in meinem Code geworden:

##window for each 100-bp segment  
windows <- numeric(6) 

##second track 
mylist = vector("list") 
mylist[[1]] = c(1,20) 
mylist[[2]] = c(120,320) 


##do the intersection 
for(i in 1:length(mylist)){ 
    st <- floor(mylist[[i]][1]/100)+1 
    sp <- floor(mylist[[i]][2]/100)+1 
    for(j in st:sp){  
    b <- max((j-1)*100, mylist[[i]][1]) 
    e <- min(j*100, mylist[[i]][2]) 
    windows[j] <- windows[j] + e - b + 1 
    } 
} 

print(windows) 
[1] 20 81 101 21 0 0 

Dies ist natürlich zu sein verwendet auf Datensätzen, die viel größer sind als das Beispiel, das ich hier zur Verfügung stelle. Durch ein Profiling kann ich sehen, dass der Engpass in den for-Schleifen liegt, aber mein ungeschickter Versuch, ihn mit * apply-Funktionen zu vektorisieren, führte zu einem Code, der eine Größenordnung langsamer verläuft.

Ich nehme an, ich könnte etwas in C schreiben, aber das möchte ich wenn möglich vermeiden. Kann jemand einen anderen Ansatz vorschlagen, der diese Berechnung beschleunigt?

Answer 1

Das "Richtige" ist die Verwendung des Bioconductor IRanges Pakets, das eine IntervalTree Datenstruktur zur Darstellung dieser Bereiche verwendet.

Wenn Sie beide Objekte in ihren eigenen IRanges Objekten haben, würden Sie dann die findOverlaps Funktion verwenden, um zu gewinnen.

es hier:

http://www.bioconductor.org/packages/release/bioc/html/IRanges.html

Durch die durch die Einbauten des Pakets sind in C geschrieben, so dass seine super schnell.

EDIT

Am zweiten Gedanken, es ist nicht so viel von einem Slam-Dunk wie ich vorschlage (ein Motto), aber Sie sollten auf jeden Fall mit dieser Bibliothek starten, wenn Sie überhaupt mit genomischer Arbeits Intervalle (oder andere Arten) ... werden Sie wahrscheinlich einige Set-Operationen und so etwas tun müssen. Entschuldigung, habe keine Zeit, die genaue Antwort zu geben.

Ich dachte nur, dass es wichtig ist, Ihnen diese Bibliothek zu zeigen.

Answer 2

Ich glaube, ich habe es viel komplizierter gemacht ... System.time hat mir bei der Leistungsbewertung in so einem kleinen Datensatz nicht geholfen.

windows <- numeric(6) 

mylist = vector("list") 
mylist[[1]] = c(1,20) 
mylist[[2]] = c(120,320) 


library(plyr) 

l_ply(mylist, function(x) { 
sapply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){ 
    eval.parent(parse(text=paste("windows[",z,"] <- ", 
     min(z*100, x[2]) - max((z-1)*100, x[1]) + 1,sep="")),sys.nframe()) 
    })   
}) 

print(windows) 

EDIT

Eine Modifikation zu beseitigen eval

g <- llply(mylist, function(x) { 
ldply((floor(x[1]/100)+1) : (floor(x[2]/100)+1), function(z){ 
     t(matrix(c(z,min(z*100, x[2]) - max((z-1)*100, x[1]) + 1),nrow=2)) 
    })   
}) 

for(i in 1:length(g)){ 
    windows[unlist(g[[i]][1])] <- unlist(g[[i]][2]) 
} 

Answer 3

Ich habe keine gute Idee haben, aber Sie können von der inneren Schleife loszuwerden, und beschleunigen Dinge ein Bit. Beachten Sie, dass wenn ein Fenster innerhalb eines Mylist-Intervalls voll fällt, Sie nur 100 zum entsprechenden Element windows hinzufügen müssen. So brauchen nur die st -ten und sp -ten Fenster spezielle Handhabung.

windows <- numeric(100) 
    for(i in 1:length(mylist)){ 
    win <- mylist[[i]]   # for cleaner code 
    st <- floor(win[1]/100)+1 
    sp <- floor(win[2]/100)+1 
    # start and stop are within the same window 
    if (sp == st){ 
     windows[st] <- windows[st] + (win[2]%%100) - (win[1]%%100) +1 
    } 
    # start and stop are in separate windows - take care of edges 
    if (sp > st){ 
     windows[st] <- windows[st] + 100 - (win[1]%%100) + 1 
     windows[sp] <- windows[sp] + (win[2]%%100) 
    } 
    # windows completely inside win 
    if (sp > st+1){ 
     windows[(st+1):(sp-1)] <- windows[(st+1):(sp-1)] + 100 
    }  
    } 

I erzeugen eine größere Liste:

cuts <- sort(sample(1:10000, 70)) # random interval endpoints 
    mylist <- split(cuts, gl(35,2)) 

und bekam 1,08 sec für 1000 Wiederholungen dieser Version im Vergleich zu 1,72 sec für 1000 repliziert für das Original. Bei realen Daten hängt die Beschleunigung davon ab, ob die Intervalle in mylist dazu neigen, viel länger als 100 zu sein oder nicht.

Übrigens könnte man die innere Schleife als separate Funktion schreiben, und dann lapply es über mylist, aber das macht es nicht schneller arbeiten.

Answer 4

So bin ich mir nicht ganz sicher, warum das dritte und vierte Fenster nicht 100 und 20 sind, weil das für mich mehr Sinn machen würde. Hier ist ein Motto für dieses Verhalten:

Reduce('+', lapply(mylist, function(x) hist(x[1]:x[2], breaks = (0:6) * 100, plot = F)$counts)) 

Beachten Sie, dass Sie benötigen die obere Schranke in breaks angeben, aber es sollte nicht schwer sein, um einen anderen Pass zu machen es zu erhalten, wenn Sie es nicht im Voraus wissen, .

Answer 5

Okay, also habe ich viel zu viel Zeit damit verschwendet, und immer noch nur einen Faktor von 3 beschleunigt. Kann jemand das schlagen?

Der Code:

my <- do.call(rbind,mylist) 
myFloor <- floor(my/100) 
myRem <- my%%100 
#Add intervals, over counting interval endpoints 
counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2]))) 
windows[as.numeric(names(counts))+1] <- counts*101 

#subtract off lower and upper endpoints 
lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum) 
windows[as.numeric(names(lowerUncovered))+1] <- windows[as.numeric(names(lowerUncovered))+1] - lowerUncovered 
upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x)) 
windows[as.numeric(names(upperUncovered))+1] <- windows[as.numeric(names(upperUncovered))+1] - upperUncovered 

Der Test:

mylist = vector("list") 
for(i in 1:20000){ 
    d <- round(runif(1,,500)) 
    mylist[[i]] <- c(d,d+round(runif(1,,700))) 
} 

windows <- numeric(200) 


new_code <-function(){ 
    my <- do.call(rbind,mylist) 
    myFloor <- floor(my/100) 
    myRem <- my%%100 
    counts <- table(do.call(c,apply(myFloor,1,function(r) r[1]:r[2]))) 
    windows[as.numeric(names(counts))+1] <- counts*101 

    lowerUncovered <- tapply(myRem[,1],myFloor[,1],sum) 
    windows[as.numeric(names(lowerUncovered))+1] <- windows[as.numeric(names(lowerUncovered))+1] - lowerUncovered 

    upperUncovered <- tapply(myRem[,2],myFloor[,2],function(x) 100*length(x) - sum(x)) 
    windows[as.numeric(names(upperUncovered))+1] <- windows[as.numeric(names(upperUncovered))+1] - upperUncovered 

    #print(windows) 
} 


#old code 
old_code <- function(){ 
    for(i in 1:length(mylist)){ 
     st <- floor(mylist[[i]][1]/100)+1 
     sp <- floor(mylist[[i]][2]/100)+1 
     for(j in st:sp){  
      b <- max((j-1)*100, mylist[[i]][1]) 
      e <- min(j*100, mylist[[i]][2]) 
      windows[j] <- windows[j] + e - b + 1 
     } 
    } 
    #print(windows) 
} 

system.time(old_code()) 
system.time(new_code()) 

Das Ergebnis:

> system.time(old_code()) 
    user system elapsed 
    2.403 0.021 2.183 
> system.time(new_code()) 
    user system elapsed 
    0.739 0.033 0.588 

Sehr frustrierend, dass die Systemzeit grundsätzlich 0, aber die beobachtete Zeit ist so groß. Ich wette, wenn Sie nach C gehen würden, würden Sie eine 50-100-fache Beschleunigung bekommen.