将个体基因组间隔串联成群体区

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

(扩展注释)您可以使用"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“可能有一个更合适的-according --提取映射的方法，但我找不到)：

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"

Answer 2

下面是使用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