1

重複領域を見つけるために、複数のベッド ファイルを処理したいと考えています。データ セットをデータ フレームとして読み取ります。オーバーラップ領域が発生した場所を検出するために、2 つのデータ セットを並行して効率的にスキャンするにはどうすればよいですか。私のアプローチは、データフレームオブジェクトの各セルのピーク領域をクエリとして取得し、intervaltree の別のデータフレームのすべての行のピーク領域を取得してから、重複領域を検索するたびです。R でこれを実装する方法がわかりません。バイオインフォマティクスでのベッド形式ファイルの処理について助けてください。誰かがこれを行う方法を私に指摘してくれれば感謝します...

これは私が達成したいことの簡単な例です:

  [1]     chr1 [10171, 10226]      * | MACS_peak_1      7.12
  [2]     chr1 [32698, 33079]      * | MACS_peak_2     13.92
  [3]     chr1 [34757, 34794]      * | MACS_peak_3      6.08
  [4]     chr1 [37786, 37833]      * | MACS_peak_4      2.44
  [5]     chr1 [38449, 38484]      * | MACS_peak_5      3.61
  [6]     chr1 [38584, 38838]      * | MACS_peak_6      4.12
  ..
  ..
  []     chrX [155191467, 155191508]      * | MACS_peak_77948      3.80
  []     chrX [155192786, 155192821]      * | MACS_peak_77949      3.71
  []     chrX [155206352, 155206433]      * | MACS_peak_77950      3.81
  []     chrX [155238796, 155238831]      * | MACS_peak_77951      3.81
  [n-1]     chrX [155246563, 155246616]      * | MACS_peak_77952      2.44
  [n]     chrX [155258442, 155258491]      * | MACS_peak_77953      5.08



  #step 1: read two bed files in R:

    bed_1 <- as(import.bed(bedFile_1), "GRanges")
    bed_2 <- as(import.bed(bedFile_2), "GRanges")
    bed_3 <- as(import.bed(bedFile_3), "GRanges")

  step 2: extract first row of the bed_1 (only take one specific interval as query). each row is considered as one specific genomic interval

    peak <- bed_1[1]      # only take one row once
    bed_2.intvl <- GenomicRanges::GIntervalTree(bed_2)

  #step 3: find overlapped regions:

    overlap <- GenomicRanges::findOverlaps(peak, bed_2.intvl)
  # step 4: go back to original bed_2 and look at which interval were overlapped with peak that comes from bed_1, what's the significance of each of these interval that comes from bed_2.

  # step 5: then iterate next interval from bed_1 to repeat above process
4

1 に答える 1

5

Bioconductorを使用して、 rtracklayerを使用してベッド ファイルをインポートします。

library(rtracklayer)
bed1 = import("foo.bed")
bed2 = import("bar.bed")

次に、「重複」を照会します。これがあなたにとって何を意味するのかははっきりしていません。

bed1OverlappingBed2 = bed1[bed1 %over% bed2]

より柔軟に、findOverlaps(bed1, bed2). このアプローチに関するフォローアップの質問は、Bioconductorサポート フォーラムに送信する必要があります。

aqueryと aを入力したとしsubjectます。すべてのヒットを見つける

hits <- findOverlaps(query, subject)

これは、2 列の行列のようなものを返します。最初の列はクエリへのインデックス、2 番目の列はサブジェクトへのインデックスです。クエリの最初の要素がサブジェクトの複数の要素と重複する場合1、クエリ ヒット列の下に複数回表示される行がいくつかあり、この範囲が重複するサブジェクトのインデックスとペアになっています。範囲の元のセットを取得し、ヒットに一致するようにそれらを「拡張」します。たとえば、

query[queryHits(hits)]

それらが重なっている領域の交点を見つけます

pintersect(query[queryHits(hits)], subject[subjectHits(hits)])

これにより、要素ごとのオーバーラップが得られましたが、反復を行わずにオーバーラップが行われました。

小さな例として、オブジェクトとして表された 'chr1' のいくつかの範囲を次に示しGRangesます (ベッド ファイルも GRanges として表されますがmcols()、ベッド ファイルからの追加情報が含まれます)。

query = GRanges("chr1", IRanges(c(10, 20, 30), width=5))
subject = GRanges("chr1", IRanges(c(10, 14), width=9))

彼らは次のように見えます

> query
GRanges object with 3 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1  [10, 14]      *
  [2]     chr1  [20, 24]      *
  [3]     chr1  [30, 34]      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths
> subject
GRanges object with 2 ranges and 0 metadata columns:
      seqnames    ranges strand
         <Rle> <IRanges>  <Rle>
  [1]     chr1  [10, 18]      *
  [2]     chr1  [14, 22]      *
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

ヒットは次のとおりです。

> hits = findOverlaps(query, subject)
> hits
Hits object with 3 hits and 0 metadata columns:
      queryHits subjectHits
      <integer>   <integer>
  [1]         1           1
  [2]         1           2
  [3]         2           2
  -------
  queryLength: 3
  subjectLength: 2

クエリの最初の範囲が件名の範囲 1 および 2 と重なっていることがわかります。ここに交差する範囲があります

> pintersect(query[queryHits(hits)], subject[subjectHits(hits)])
GRanges object with 3 ranges and 1 metadata column:
      seqnames    ranges strand |       hit
         <Rle> <IRanges>  <Rle> | <logical>
  [1]     chr1  [10, 14]      * |      TRUE
  [2]     chr1  [14, 14]      * |      TRUE
  [3]     chr1  [20, 22]      * |      TRUE
  -------
  seqinfo: 1 sequence from an unspecified genome; no seqlengths

したがって、クエリ 1 とサブジェクト 1 は 10 から 14 の位置で重複し、クエリ 1 とサブジェクト 2 は位置 14 で重複し、クエリ 2 とサブジェクト 2 は位置 20 から 22 で重複します。ベースの半開間隔;rtracklayer::import.bed()ファイルをインポートするときに正しいことを行います。

于 2016-01-06T20:01:28.113 に答える