将单个基因组区间连接到种群区域答案

【问题标题】：Concatenate individual genomic intervals into populational regions将单个基因组区间连接到种群区域
【发布时间】：2015-11-16 15:08:53
【问题描述】：

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

我的意见：

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 更新

【问题讨论】：

我不认为有一个简单的解决方法可以一次性创建所有内容。我将首先按chr 和start 排序，然后创建合并的区间。之后检查 dfin 中的 each 是否包含在每个 dfout 间隔中，并相应地更新 samples 和 type。这应该涉及一些逻辑和编程，但我认为没有快速解决方法。
你会发现“IRanges”包很有帮助； reduce(RangedData(space = dfin$chr, IRanges(dfin$start, dfin$end)))
@alexis_laz，您的代码在我的真实数据集上运行良好。请写一个正式的答案，我会把它标记为正确的。

标签： r overlap overlapping bioconductor

【解决方案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

【讨论】：

您的回答非常适合我的示例，但过于简单化。我在同一个chr 中有几个间隔。因此，只要它们有重叠，就需要将它们变成一个公共区域。请检查我的新示例！

【解决方案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"

【讨论】：