在两个区间数据之间找到重叠的范围

Question

在两个区间数据之间找到重叠的范围

6

我有一个包含约500000个片段的坐标(start,end)的表格和另一个包含60000个单一坐标的表格，我想将它与前一个表格中的片段匹配。也就是说，对于dtCoords表格中的每个记录，我需要在dtFrags表格中搜索具有相同chr和start<=coord<=end（并从此记录检索dtFrags的type）。我是否应该使用R处理这个问题，还是应该寻找其他编程语言？

以下是我的示例：

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)

- Vasily A

我现在没有运行R的方法，但这个问题对我来说很有趣。我认为首先应该按染色体划分数据，然后希望您在范围内没有重复的坐标，然后基于起始位置和位置创建坐标列，然后稍后与该坐标合并。如果到明天这个问题还没有解决，我会回来处理它。 - Ananta

2个回答

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

- Metrics

@Frank，输出几乎没问题（应该只是>=和<=而不是>和<），但我不能用它来处理我的大表格。 - Vasily A

很不幸，我的表格似乎太大了，我得到了这个错误信息：

Error in vecseq(f__, len__, if (allow.cartesian) NULL else as.integer(max(nrow(x),  :    Join results in 1447477452 rows; more than 519176 = max(nrow(x),nrow(i)). Check for duplicate key values in i, each of which join to the same group in x over and over again. If that's ok, try including 'j' and dropping 'by' (by-without-by) so that j runs for each group to avoid the large allocation.

- Vasily A

@VasilyA 这是大小的问题吗？我怀疑这个错误可以在一个小的可重现的例子中显示出来。我知道之前我使用四行数据表时也遇到过这个问题。 - Frank

是的，大小确实很重要，恐怕我不能在我的情况下使用这个解决方案：我的“chr”字段有24个级别，因此仅通过此字段合并会产生数十亿种组合... - Vasily A

嗯，也许你应该重新格式化你的dtFrags，例如，如果A有区间（250,350），B有区间（300,400），那么你可以创建四个记录（A在250-300，A在300-350，B在300-350，B在350-400）。然后你就会有非重叠的区间，可能可以使用roll=选项。（我对roll的使用不是很熟悉。）无论如何，我仍然相信你可以在一个更小的例子中重现你的错误。 - Frank

1

谢谢@Frank，我会尝试找出关于roll选项的相关信息。至于复现错误-我不知道如何用小例子来做到这一点，但我在帖子末尾添加了我的真实数据的链接（'UPD'）。 - Vasily A

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Arun · Accepted Answer

通常情况下，使用bioconductor软件包IRanges来处理与区间相关的问题非常合适。它通过实现interval tree来高效地解决这些问题。GenomicRanges是另一个构建在IRanges上的软件包，专门用于处理“基因组范围”。

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