如何处理两个床层文件，并行查找重叠区域？

Question

我想要处理多个床文件，以查找重叠的区域。我将数据集读取为数据帧，如何有效地并行扫描两个数据集，以检测重叠区域发生在何处。我的方法是每次以数据帧对象的每个单元格的峰值区域作为查询，在中间树中取另一个数据帧的所有行的峰值区域，然后搜索重叠区域。我很困惑如何在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

Answer 1

使用Bioconductor，使用rtracklayer导入床文件

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

然后查询“重叠”，这对你来说意味着什么，也许不太清楚。

bed1OverlappingBed2 = bed1[bed1 %over% bed2]

更灵活点，findOverlaps(bed1, bed2)。这个方法的后续问题应该指向生物导体support forum。

假设我们输入了一个query和一个subject。找到所有的命中

hits <- findOverlaps(query, subject)

这会返回类似于两列矩阵的东西。第一列是查询索引，第二列是主题索引。如果查询的第一个元素与主题的多个元素重叠，那么就会有几行1在query列下出现几次，并与此范围重叠的主题的索引配对。获取原始的范围集并“扩展”它们以匹配命中，例如，

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位置。(生物导体使用基于1的封闭间隔；UCSC使用基于0的半开间隔；rtracklayer::import.bed()在导入文件时做正确的事情。