2

私は I/GRanges Views オブジェクトを持っています

** データの単純化されたバージョンであり、実際のデータは膨大です

Views on a 10000000-length Rle subject

 views:
      start      end   width
 [1]      1     1000    1000 [100 100 100 100 100 100 100 100 100 100 ...]
 [2]   1001     2000    1000 [190 190 190 190 190 190 190 190 190 190 ...]
 [3]   2001     3000    1000 [280 280 280 280 280 280 280 280 280 280 ...]
 [4]   3001     4000    1000 [370 370 370 370 370 370 370 370 370 370 ...]
 [5]   4001     5000    1000 [460 460 460 460 460 460 460 460 460 460 ...]
 ...    ...      ...     ... ...
 [9996] 995001  9996000 9001000 [89650 89650 89650 89650 89650 89650 ...]
 [9997] 996001  9997000 9001000 [89740 89740 89740 89740 89740 89740 ...]
 [9998] 997001  9998000 9001000 [89830 89830 89830 89830 89830 89830 ...]
 [9999] 998001  9999000 9001000 [89920 89920 89920 89920 89920 89920 ...]
[10000] 999001 10000000 9001000 [90010 90010 90010 90010 90010 90010 ...]

各ビュー (行) の幅は 1000 で、それぞれ 100 個のデータポイントが 1000 個あることを意味します。ここで、データポイントのセットを 20 個のビン (この場合はビンごとに 50 個) に分割して平均を取りたいと思います。その結果、出力は 20 個の数値のベクトルになり、それぞれがそのビンの平均になります。

出力:

[1] 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

現在、実際の状況では、そのようなビューが 20 以上あり、各行の幅が異なり、一部の行が 5K を超えています。私のコードは正常に動作しますが、私のデータでは非常に遅く、各行で 20 個のビンのベクトルを返すのに約 1.5 秒かかり、30K 行を超えると約 12.5 時間かかります。

これらの計算を高速化する方法があると確信しています。そうでない場合は、何らかの方法でクラスターの並列ノードを使用できます。何を指示してるんですか。

データを生成するテスト コード:

library('GenomicRanges')
# generating data frame 
df=data.frame(chrom=rep('Chr1',100000),start=seq(1,1000000,by=1000),end=seq(1000,10000000,by=1000),strand=rep("+",100000))

# making GRanges object
gr=GRanges(seqnames=as.vector(df[,1]),IRanges(start=df[,2],end=df[,3]),strand=df[,4])
# obtaining coverage using function coverage in the form of RLE object
gr.cov=coverage(gr)
# generating views for specific start and end
gr.views=Views(gr.cov[[1]],start=seq(1,1000000,by=1000),end=seq(1000,10000000,by=1000))
# putting in temp variable
d=gr.views

# this following code calculates the matrix (where each line is 20 points) for 10 lines
# reduce or increase the number in the outermost sapply loop to increase/decrease the lines to be calculated

sapply(1:10,function(j)
   sapply(1:20,
   function(i)as.numeric(
     format(
       mean(
         as(d[[j]][(
           seq(0,length(d[[j]]),floor(length(d[[j]])/20))+1)[i]:
             c((seq(0,length(d[[j]]),floor(length(d[[j]])/20)))[
               -length((seq(0,length(d[[j]]),floor(length(d[[j]])/20))))
               ],length(d[[j]]))[i+1]],
            "RangedData")$score),
       digits=2)
     )
   )
)
4

1 に答える 1

1

遺伝子に基づいてビューを作成するのではなく、計算を行いたいウィンドウに基づいてビューを作成し、またはビューの統計を計算するために使用してみませんviewSumsかviewMaxs? GRanges「遺伝子」（転写物？）の開始と終了を記述したとします。

genes <- GRanges(seqnames, IRanges(geneStarts, geneEnds))

ウィンドウの始まりと終わりを計算するかもしれません

n <- 50L
starts0 <- Map(function(...) floor(seq(...)), start(genes), end(genes),
               MoreArgs=list(length.out=n + 1L))
ends <- lapply(starts0, function(x) floor(x)[-1])
starts <- lapply(starts0, function(x) x[-length(x)])

次に、ビューを作成します

v <- Views(gr.cov[[1]], start=unlist(starts), end=unlist(ends))

( ?RleViews「Views,RleList-method」を参照) 統計を計算し、遺伝子ごとに分割する

split(viewMeans(v), rep(seq_along(genes), each=n))

Bioconductor メーリング リストで質問すると、多くの巧妙な解決策が得られる可能性があります。

v「RleViews」クラスのインスタンスです。v[[1]]のインスタンスですRle。mean(v[[1]])の確認として計算するviewMeansか、さらに一歩進めv[[1]]て単純な古いベクトルに強制し、それを計算することができmean(as.vector(v[[1]])))ます。(フードの下をのぞくのではなく、適切なアクセサーを使用するrunValue(v[[1]])のと同じです) Rle の値を返します。v[[1]]@values

> (x <- Rle(c(rep(100, 10), rep(200, 10))))
numeric-Rle of length 20 with 2 runs
  Lengths:  10  10
  Values : 100 200
> runValue(x)
[1] 100 200
> runLength(x)
[1] 10 10

そして明らかにmean(runValue(x)) != mean(x)。

于 2012-11-27T02:39:21.760 に答える