Я пытаюсь сделать простое геномное пересечение дорожки в R и сталкиваюсь с главными проблемами производительности, вероятно, связанными с моим использованием для циклов.
В этой ситуации я предопределил окна с промежутками в 100bp, и я пытаюсь вычислить, сколько из каждого окна охвачено аннотациями в mylist. Графически, это выглядит примерно так:
0 100 200 300 400 500 600
windows: |-----|-----|-----|-----|-----|-----|
mylist: |-| |-----------|
Таким образом, я написал некоторый код, чтобы сделать просто, что, но это довольно медленно и стало узким местом в моем коде:
##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
Естественно, это используется на наборах данных, которые намного больше, чем пример, который я обеспечиваю здесь. Посредством некоторого профилирования я вижу, что узкое место находится в для циклов, но моя неуклюжая попытка векторизовать его, использование *применяет функции, привела к коду, который выполняет порядок величины более медленно.
Я предполагаю, что мог записать что-то в C, но я хотел бы избежать этого, если это возможно. Кто-либо может предложить другой подход, который ускорит это вычисление?
"Правильно" будет использовать пакет bioconductor IRanges
, который использует структуру данных IntervalTree для представления этих диапазонов.
Имея оба объекта в своих собственных объектах IRanges
, вы должны затем использовать функцию findOverlaps
, чтобы выиграть.
Получите это здесь:
http://www.bioconductor.org/packages/release/bioc/html/IRanges.html
Кстати, внутренняя часть пакета написана на C, поэтому его сверх быстрый.
РЕДАКТИРОВАТЬ
Если подумать, это не так уж и сложно, как я предлагаю (однострочный), но вам обязательно стоит начать использовать эту библиотеку, если вы вообще работаете с геномными интервалами ( или другие типы) ... вам, вероятно, потребуется выполнить некоторые операции с наборами и прочее. К сожалению, у меня нет времени дать точный ответ.
Я просто подумал, что важно указать вам на эту библиотеку.
Ладно, я потратил на это СЛИШКОМ много времени, и все же получил только трехкратное ускорение. Может ли кто-нибудь победить это?
Код:
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
Тест:
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())
Результат:
> 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
Очень неприятно, что системное время в основном равно 0, но наблюдаемое время настолько велико. Бьюсь об заклад, если вы спуститесь до C, вы получите ускорение в 50-100 раз.
Мне кажется, я значительно усложнил задачу... System.time не помог мне в оценке производительности на таком маленьком наборе данных.
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
Модификация для устранения 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])
}
У меня нет блестящей идеи, но вы можете избавиться от внутреннего цикла и немного ускорить процесс. Обратите внимание, что если окно полностью выпадает в интервале mylist, вам просто нужно добавить 100 к соответствующему элементу windows
. Таким образом, особой обработки требуют только -ое
-ое и sp
-ое окна.
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
}
}
Я создал более крупный список:
cuts <- sort(sample(1:10000, 70)) # random interval endpoints
mylist <- split(cuts, gl(35,2))
и получил 1,08 секунды для 1000 копий этой версии по сравнению с 1,72 секунды для 1000 копий для оригинала. С реальными данными ускорение будет зависеть от того, будут ли интервалы в mylist
намного больше 100 или нет.
Кстати, можно было бы переписать внутренний цикл как отдельную функцию, а затем перекрыть
его mylist
, но это не ускоряет его работу.
Так что я не совсем уверен, почему третье и четвертое окна не 100 и 20, потому что это имело бы для меня больше смысла. Вот один лайнер для этого поведения:
Reduce('+', lapply(mylist, function(x) hist(x[1]:x[2], breaks = (0:6) * 100, plot = F)$counts))
Обратите внимание, что вам нужно указать верхнюю границу в breaks
, но не должно быть сложно сделать еще один проход, чтобы получить его, если вы этого не знаете. заранее.