将个体基因组间隔连接到种群区域

Question

我想将各个基因组间隔连接到公共区域。

我的输入：

dfin <- "chr start end sample type
        1   10    20   NE1    loss
        1   5     15   NE2    gain
        1   25    30   NE1    gain
        2   40    50   NE1    loss
        2   40    60   NE2    loss
        3   20    30   NE1    gain"
dfin <- read.table(text=dfin, header=T)

我的预期输出：

dfout <- "chr start end samples type
        1   5     20   NE1-NE2  both
        1   25    30   NE1      gain
        2   40    60   NE1-NE2  loss
        3   20    30   NE1      gain"
dfout <- read.table(text=dfout, header=T)

dfin 中的间隔在同一动物中永远不会重叠，只是在动物之间（分别为 sample 和 samples 列）。 type 列在 dfin 中有两个因子（loss 和 gain），预计在 dfout 中有三个因子（loss、gain 和 both，这发生在 dfout 中的连接区域同时基于 loss 和 gain 时）。

有什么办法解决这个问题吗？

*更新@David Arenburg

Answer 1

这里尝试使用 data.table::foverlaps 对间隔进行分组，然后计算所有其余部分

library(data.table)
setkey(setDT(dfin), chr, start, end)
res <- foverlaps(dfin, dfin, which = TRUE)[, toString(xid), by = yid
                                           ][, indx := .GRP, by = V1]$indx
dfin[, .(
          chr = chr[1L],
          start = min(start), 
          end = max(end), 
          samples = paste(unique(sample), collapse = "-"),
          type = if(uniqueN(type) > 1L) "both" else as.character(type[1L])
         ),
       by = res]

#    res chr start end samples type
# 1:   1   1     5  20 NE2-NE1 both
# 2:   2   1    25  30     NE1 gain
# 3:   3   2    40  60 NE1-NE2 loss
# 4:   4   3    20  30     NE1 gain

Answer 2

（扩展评论）您可以使用 "IRanges" 包：

library(IRanges)

#build an appropriate object
dat = RangedData(space = dfin$chr, 
                 IRanges(dfin$start, dfin$end), 
                 sample = dfin$sample, 
                 type = dfin$type)
dat
#concatenate overlaps with an extra step of saving the concatenation mappings
ans = RangedData(reduce(ranges(dat), with.revmap = TRUE))
ans

无法弄清楚如何避免 reduce 丢失 "RangedData" 对象的列，但是保存映射后我们可以做类似的事情（可能有更合适的 - 根据 "IRanges"- 提取映射的方法，但我找不到它):

tmp = elementMetadata(ranges(ans)@unlistData)$revmap@partitioning
maps = rep(seq_along(start(tmp)), width(tmp))
maps
#[1] 1 1 2 3 3 4

有了区间连接的映射，我们可以聚合 "sample" 和 "type" 得到最终形式。例如：

tapply(dfin$sample, maps, function(X) paste(unique(X), collapse = "-"))
#        1         2         3         4 
#"NE1-NE2"     "NE1" "NE1-NE2"     "NE1"

将个体基因组间隔连接到种群区域

Concatenate individual genomic intervals into populational regions

r

overlap

overlapping

bioconductor