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

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)。有关此方法的后续问题应直接发送至 Bioconductor support forum。

假设我们输入了query 和subject。查找所有热门歌曲

hits <- findOverlaps(query, subject)

这会返回一个看起来像两列矩阵的东西。第一列是查询索引，第二列是主题索引。如果查询的第一个元素与主题的多个元素重叠，则在查询命中列下将有几行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 重叠。（Bioconductor 使用基于 1 的封闭区间；UCSC使用基于 0 的半开间隔；rtracklayer::import.bed() 在导入文件时做正确的事情。