5

私は2つのファイルを持っています

エンコード

X.pattern.name       chr       start    stop    strand  score   p.value q.value matched.sequence
1   V_CETS1P54_01   chr1    98769545    98769554    +   11.42280    8.89e-05    NA  TCAGGATGTA
2   V_CETS1P54_01   chr1    152013037   152013046   +   11.98020    4.74e-05    NA  ACAGGAAGTT
3   V_CETS1P54_01   chr1    168932563   168932572   +   11.60860    7.59e-05    NA  ACCGGATGCT

エンコード合計

    chr     start           stop
1   chr1    58708485    58708713
2   chr1    58709084    58710538
3   chr1    98766295    98766639
4   chr1    98766902    98770338
5   chr1    107885506   107889414
6   chr1    138589531   138590856
7   chr1    138591180   138593378
8   chr1    152011746   152013185
9   chr1    152014263   152014695
10  chr1    168930561   168933076
11  chr1    181808064   181808906
12  chr1    184609002   184611519
13  chr1    193288453   193289567
14  chr1    193290105   193290490
15  chr1    193290744   193291092
16  chr1    196801920   196804954

2 つのファイルを比較したいのですが、各エントリはchr、start、およびstopです。最初のファイルの開始値と終了値が同じ染色体の 2 番目のファイルの開始値と終了値の間にある場合、最初のファイルの開始値と終了値は 2 番目のファイルの開始値と終了値に置き換えられる必要があります。そのために for ループを作成しましたが、時間がかかりすぎます。代替手段は何ですか？

コード：

for(i in 1:nrow(encode))
{
     for(j in 1:nrow(encode.total))
     {
           if(encode[i,2]==encode.total[j,1])
           {
               if((encode[i,3]>=encode.total[j,2]) & (encode[i,4]<=encode.total[j,3]))
               {
                   encode[i,3]=encode.total[j,2]
                   encode[i,4]=encode.total[j,3]
               }
           }
     }
}

同じ目的で、以下のように GenomicRanges パッケージも試しました。私のデータフレームのサイズは巨大で、それらのマージ機能は非常に大きなデータフレーム（許可されていない> 20億行）を作成しますが、最終的にはデータフレームをより小さなデータフレームにサブセット化します。しかし、merge() は多くのメモリを消費し、R を終了させて​​います。

gr1<-GRanges(seqnames=encode$chr,IRanges(start=encode$start,end=encode$end))
gr2<-GRanges(seqnames=encode.total$chr, IRanges(start=encode.total$start,end=encode.total$end))

ranges <- merge(as.data.frame(gr1),as.data.frame(gr2),by="seqnames",suffixes=c("A","B"))
ranges <- ranges[with(ranges, startB <= startA & endB >= endA),]
4

1 に答える 1

12

Bioconductor GenomicRangesパッケージを使用します。

source("http://bioconductor.org/biocLite.R")
biocLite("GenomicRanges")

元の例のとのように、x0との2 つのデータ フレームがあるとします。これらを GRanges インスタンスにしたいと思います。私はこれをやったx1encodeencode.total

library(GenomicRanges)
gr0 = with(x0, GRanges(chr, IRanges(start=start, end=stop)))
gr1 = with(x1, GRanges(chr, IRanges(start=start, end=stop)))

多くの場合、単純に と言うmakeGRangesFromDataFrame(x0)か、標準の R コマンドを使用して GRanges インスタンスを「手動で」作成することができます。ここでは を使用するwith()ので、GRanges(chr, IRanges(start=start, end=stop))の代わりに書くことができますGRanges(x0$chr, IRanges(start=x0$start, end=x0$stop))。

次のステップは、クエリ (gr0) とサブジェクト (gr1) の間の重複を見つけることです。

hits = findOverlaps(gr0, gr1)

につながる

> hits
Hits of length 3
queryLength: 3
subjectLength: 16
  queryHits subjectHits 
   <integer>   <integer> 
 1         1           4 
 2         2           8 
 3         3          10

次に、関連する開始/終了座標を更新しました

ranges(gr0)[queryHits(hits)] = ranges(gr1)[subjectHits(hits)]

与える

> gr0
GRanges with 3 ranges and 0 metadata columns:
      seqnames                 ranges strand
         <Rle>              <IRanges>  <Rle>
  [1]     chr1 [ 98766902,  98770338]      *
  [2]     chr1 [152011746, 152013185]      *
  [3]     chr1 [168930561, 168933076]      *
  ---
  seqlengths:
   chr1
     NA

これは、数百万の範囲に高速化されます。

于 2013-10-01T01:29:22.297 に答える