6

ca の座標 ( start、 )を持つテーブルが 1 つあります。end500000 個のフラグメントと、前のフラグメントと一致させたい 60000 個の単一座標を持つ別のテーブル。つまり、テーブルの各レコードについて、同じおよび<= <=を持つテーブルdtCoords内のレコードを検索する必要があります (そして、このレコードからを取得します)。これにRを使用するのはまったく良い考えですか、それとも他の言語に目を向けるべきですか?dtFragschrstartcoordendtypedtFrags

これが私の例です：

require(data.table)

dtFrags <- fread(
  "id,chr,start,end,type
 1,1,100,200,exon
 2,2,300,500,intron
 3,X,400,600,intron
 4,2,250,600,exon
")

dtCoords <- fread(
"id,chr,coord
 10,1,150
 20,2,300
 30,Y,500
")

最後に、私はこのようなものを持ちたいと思います:

"idC,chr,coord,idF,type
 10,  1,  150,  1, exon
 20,  2,  300,  2, intron
 20,  2,  300,  4, exon
 30,  Y,  500, NA, NA
"

でテーブルをサブテーブルに分割することで、タスクを少し簡素化できるchrので、座標だけに集中します

setkey(dtCoords, 'chr')
setkey(dtFrags,  'chr')

for (chr in unique(dtCoords$chr)) {
  dtCoordsSub <- dtCoords[chr];
  dtFragsSub  <-  dtFrags[chr];
  dtCoordsSub[, {
    # ????  
  }, by=id]  
}

しかし、内部でどのように作業すればよいかはまだ明確ではありません...ヒントをいただければ幸いです。

アップデート。念のため、実際のテーブルをここのアーカイブに入れました。作業ディレクトリに解凍した後、次のコードでテーブルをロードできます。

dtCoords <- fread("dtCoords.txt", sep="\t", header=TRUE)
dtFrags  <- fread("dtFrags.txt",  sep="\t", header=TRUE)
4

2 に答える 2

7

一般に、間隔に関連する問題に対処するには、 bioconductorパッケージを使用するのが非常に適切です。インターバルツリーIRangesを実装することで効率的にこれを行います。は、特に「ゲノム範囲」を処理するために、上に構築される別のパッケージです。GenomicRangesIRanges

require(GenomicRanges)
gr1 = with(dtFrags, GRanges(Rle(factor(chr, 
          levels=c("1", "2", "X", "Y"))), IRanges(start, end)))
gr2 = with(dtCoords, GRanges(Rle(factor(chr, 
          levels=c("1", "2", "X", "Y"))), IRanges(coord, coord)))
olaps = findOverlaps(gr2, gr1)
dtCoords[, grp := seq_len(nrow(dtCoords))]
dtFrags[subjectHits(olaps), grp := queryHits(olaps)]
setkey(dtCoords, grp)
setkey(dtFrags, grp)
dtFrags[, list(grp, id, type)][dtCoords]

   grp id   type id.1 chr coord
1:   1  1   exon   10   1   150
2:   2  2 intron   20   2   300
3:   2  4   exon   20   2   300
4:   3 NA     NA   30   Y   500
于 2013-11-03T12:01:39.510 に答える
3

これは機能しますか？merge最初に使用してから使用できますsubset

   kk<-merge(dtFrags,dtCoords,by="chr",all.x=TRUE)
> kk
   chr id.x start end   type id.y coord
1:   1    1   100 200   exon   10   150
2:   2    2   300 500 intron   20   300
3:   2    4   250 600   exon   20   300
4:   X    3   400 600 intron   NA    NA


 kk[coord>=start & coord<=end]
   chr id.x start end type id.y coord
1:   1    1   100 200 exon   10   150
2:   2    4   250 600 exon   20   300
于 2013-11-03T00:50:17.480 に答える