4

長さ 5,000,000 の数値ベクトルがあります

>head(coordvec)
[1] 47286545 47286546 47286547 47286548 47286549 472865

および 3 x 1,400,000 数値行列

>head(subscores)
        V1       V2     V3
1 47286730 47286725  0.830
2 47286740 47286791  0.065
3 47286750 47286806 -0.165
4 47288371 47288427  0.760
5 47288841 47288890  0.285
6 47288896 47288945  0.225

私が達成しようとしているのは、coordvec の各数値について、V1 と V2 が coordvec の数値を含むサブスコアの行の V3 の平均を見つけることです。そのために、私は次のアプローチを取っています。

results<-numeric(length(coordvec))
for(i in 1:length(coordvec)){
    select_rows <- subscores[, 1] < coordvec[i] & subscores[, 2] > coordvec[i]
scores_subset <- subscores[select_rows, 3]
results[m]<-mean(scores_subset)
}

これは非常に遅く、完了するまでに数日かかります。もっと速い方法はありますか？

ありがとう、

ダン

4

3 に答える 3

6

この質問には 2 つの難しい部分があると思います。1 つ目は、重なりを見つけることです。BioconductorIRangesのパッケージを使用します(基本パッケージも役立つ場合があります)。?findInterval

library(IRanges)

座標ベクトルを表す幅 1 の範囲と、スコアを表す範囲のセットを作成します。便宜上、重複した座標を同じように扱うことができると仮定して、座標ベクトルを並べ替えます

coord <- sort(sample(.Machine$integer.max, 5000000))
starts <- sample(.Machine$integer.max, 1200000)
scores <- runif(length(starts))

q <- IRanges(coord, width=1)
s <- IRanges(starts, starts + 100L)

ここで、どれがどれqueryと重複しているかを見つけますsubject

system.time({
    olaps <- findOverlaps(q, s)
})

これには、私のラップトップで約 7 秒かかります。オーバーラップにはさまざまなタイプがある?findOverlapsため (「 」を参照)、この手順には少し調整が必要になる場合があります。結果は、クエリと重複するサブジェクトにインデックスを付けるベクトルのペアです。

> olaps
Hits of length 281909
queryLength: 5000000
subjectLength: 1200000
       queryHits subjectHits 
        <integer>   <integer> 
 1             19      685913 
 2             35      929424 
 3             46     1130191 
 4             52       37417

281909 のオーバーラップを見つけて、最初の複雑な部分はこれで終わりだと思います。（私は間違っている可能性がありますが、他の場所で提供されているdata.tableの回答はこれに対処しているとは思いません...）

次の難しい部分は、多数の平均を計算することです。組み込みの方法は次のようになります

olaps0 <- head(olaps, 10000)
system.time({
    res0 <- tapply(scores[subjectHits(olaps0)], queryHits(olaps0), mean)
})

私のコンピューターでは約 3.25 秒かかり、線形にスケーリングするように見えるため、280k のオーバーラップの場合は 90 秒になる可能性があります。しかし、この集計は で効率的に実行できると思いますdata.table。元の座標はstart(v)[queryHits(olaps)]ですので、

require(data.table)
dt <- data.table(coord=start(q)[queryHits(olaps)],
                 score=scores[subjectHits(olaps)])
res1 <- dt[,mean(score), by=coord]$V1

280k のオーバーラップすべてに約 2.5 秒かかります。

クエリのヒットが順序付けられていることを認識することで、速度をさらに向上させることができます。クエリヒットの実行ごとに平均を計算したいと考えています。まず、各クエリ ヒット実行の終了を示す変数を作成します。

idx <- c(queryHits(olaps)[-1] != queryHits(olaps)[-length(olaps)], TRUE)

次に、各実行の終了時の累積スコア、各実行の長さ、および実行の終了時と開始時の累積スコアの差を計算します。

scoreHits <- cumsum(scores[subjectHits(olaps)])[idx]
n <- diff(c(0L, seq_along(idx)[idx]))
xt <- diff(c(0L, scoreHits))

そして最後に、平均は

res2 <- xt / n

これは、すべてのデータに対して約 0.6 秒かかり、data.table の結果と同じです (より不可解ではありますが?)。

> identical(res1, res2)
[1] TRUE

平均に対応する元の座標は

start(q)[ queryHits(olaps)[idx] ]
于 2013-01-20T01:27:20.003 に答える
2

このようなものの方が速いかもしれません:

require(data.table)
subscores <- as.data.table(subscores)

subscores[, cond := V1 < coordvec & V2 > coordvec]
subscores[list(cond)[[1]], mean(V3)]

list(cond)[[1]]理由:「i が単一の変数名の場合、列名の式とは見なされず、呼び出しスコープで評価されます。」ソース：?data.table

于 2013-01-20T00:00:22.560 に答える
0

あなたの答えは簡単に再現することはできず、たとえそうであったとしても、subscoresあなたのブール条件を満たしていないので、これがあなたが探しているものを正確に実行するかどうかはわかりませんが、applyファミリの1つと関数を使用できます。

myfun <- function(x) {
  y <- subscores[, 1] < x & subscores[, 2] > x
  mean(subscores[y, 3])
}

sapply(coordvec, myfun)

もご覧くださいmclapply。あなたが十分なメモリを持っているなら、これはおそらく物事をかなりスピードアップするでしょう。ただし、foreach同様の結果でパッケージを確認することもできます。あなたはそれを成長させるのではなくfor loopに割り当てることによってあなたの「正しい」を手に入れました、しかし実際には、あなたはたくさんの比較をしています。これを大幅にスピードアップするのは難しいでしょう。results

于 2013-01-19T23:51:11.630 に答える