5

Rで単純なゲノムトラック交差を実行しようとしていますが、おそらくforループの使用に関連する主要なパフォーマンスの問題が発生しています。

この状況では、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

当然、これは、ここで提供する例よりもはるかに大きいデータセットで使用されています。いくつかのプロファイリングを通じて、ボトルネックがforループにあることがわかりますが、* apply関数を使用してそれをベクトル化しようとすると、コードの実行速度が1桁遅くなります。

私はCで何かを書くことができると思いますが、可能であればそれを避けたいと思います。誰かがこの計算をスピードアップする別のアプローチを提案できますか？

4

5 に答える 5

6

「正しい」ことはIRanges、これらの範囲を表すためにIntervalTreeデータ構造を使用するbioconductorパッケージを使用することです。

両方のオブジェクトをそれぞれのIRangesオブジェクトに入れて、関数を使用しfindOverlapsて勝ちます。

ここで入手してください：

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

によって、パッケージの内部はCで書かれているので、超高速です。

編集

考え直してみると、これは私が提案しているほどのスラムダンクではありません（ワンライナー）が、ゲノム間隔（または他のタイプ）で作業している場合は、間違いなくこのライブラリの使用を開始する必要があります...おそらく、いくつかの設定操作などを行う必要があります。申し訳ありませんが、正確な答えを提供する時間がありません。

このライブラリをあなたに指摘することが重要だと思いました。

于 2010-03-25T20:25:44.450 に答える
4

したがって、3番目と4番目のウィンドウが100と20でない理由は完全にはわかりません。それは、私にとってより理にかなっているからです。これがその振る舞いのワンライナーです：

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

で上限を指定する必要があることに注意してください。ただしbreaks、事前にわからない場合は、別のパスを作成して取得するのは難しくありません。

于 2010-03-25T20:24:25.097 に答える
4

さて、私はこれにあまりにも多くの時間を無駄にしましたが、それでも3倍のスピードアップしか得られませんでした。誰もがこれを打ち負かすことができますか？

コード：

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倍のスピードアップが得られるでしょう。

于 2010-03-26T08:01:45.030 に答える
1

私はそれをもっと複雑にしたと思います...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)

編集

排除するための変更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])
}
于 2010-03-25T18:52:17.050 に答える
0

私には良い考えはありませんが、内側のループを取り除き、物事を少しスピードアップすることができます。ウィンドウがmylist間隔で完全に落ちる場合は、対応するwindows要素に100を追加するだけでよいことに注意してください。したがって、st-番目と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))

このバージョンの1000回の複製で1.08秒であるのに対し、元のバージョンでは1000回の複製で1.72秒でした。実際のデータでは、スピードアップは、の間隔mylistが100よりはるかに長くなる傾向があるかどうかによって異なります。

ちなみに、内部ループを別の関数として書き直してから、lapplyそれを上書きすることもできますが、それではmylist動作が速くなりません。

于 2010-03-25T19:17:05.393 に答える