创建一个矩阵以显示多个 GRange 之间的重叠答案

【问题标题】：Create a matrix to show overlaps among multiple GRanges创建一个矩阵以显示多个 GRange 之间的重叠
【发布时间】：2019-02-28 08:58:41
【问题描述】：

在比较不同的GRange 对象时，我试图找到一种方法来有效地提取显示“0”或“1”的矩阵。在我的例子中：

df <- data.frame(chr = c("chr1", "chr10"), start = c(1,4), end=c(2, 4))
gr.1 <- makeGRangesFromDataFrame(df)

df <- data.frame(chr = c("chr1", "chr10"), start = c(2,3), end=c(2, 4))
gr.2 <- makeGRangesFromDataFrame(df)

df <- data.frame(chr = c("chr1"), start = c(1), end=c(1))
gr.3 <- makeGRangesFromDataFrame(df)

我尝试findOverlaps 来评估这些区域之间的重叠，但显然它不能处理两个以上的GRanges：

> GenomicRanges::findOverlaps(gr.1, gr.2, gr.3)
> Error in IRanges:::NCList_find_overlaps_in_groups(ranges(query),
> q_space,  :    'maxgap' must be a single integer

此外，我需要的输出类似于此示例数据框：

out <- "gr.1 gr.2 gr.3
chr1-1 1  0  1
chr1-2 1  1  0
chr10-3 0 1  0
chr10-4 1 1  0"

out <- read.table(text=out, header=TRUE)

有什么明智地导出它的想法吗？

【问题讨论】：

你能提供更多的最小样本数据吗（真的有必要拥有 3 个 GRanges 每个包含 1000 个特征）吗？此外，regionR::createRandomRegions 具有进一步的依赖关系，这使得生成样本数据变得不必要地尴尬。如果您要提供较小的样本数据（包括匹配的预期输出），这将有很大帮助。为了解决您的问题，这本质上是一个多数据集合并，其解决方案通常非常简单，例如使用Reduce(findOverlaps, list_of_gr) 并处理生成的对象。
事实上，regionR::createRandomRegions 是我发现的代表一组随机基因组范围的最简单方法。这是一个相当新的功能，适合我的例子。此外，我不确定生成随机基因组区域是否直接（参见biostars.org/p/225520）。最后，我正在模拟 1000 个基因组区域以增加变化，以使一些复杂的重叠模式作为我的真实数据集。另外，Reduce(findOverlaps, glg) 只是给了我一个错误。
我强烈建议手动生成GRanges 的最小集合；我几乎可以保证，对于像你提供的那样笨重的样本数据，你将得到非常少的（积极的）回应。如果您手动生成数据，您还可以确保范围重叠；无需生成那么多范围！正如我所解释的，如果无法访问最少的样本数据，就很难（不可能）提供适当和具体的帮助；使用Reduce 绝对是一种选择，但要使其发挥作用，我们需要数据。
[更新] 请记住，您要求这里的其他人利用他们的空闲时间来帮助您解决您的问题。第一条规则是：让我们尽可能轻松地为您提供帮助。
我明白你的意思，并且无意冒犯那些在 stackoverflow 上花时间帮助他人的人。希望我的更新现在更容易消化。

标签： r bioinformatics bioconductor

【解决方案1】：

首先，这是一个部分解决方案，它仅显示第一个和任何其他 GRanges 之间的重叠区域（这应该产生类似于 bedtools intersect 的结果，它允许“识别单个查询 (-a) 文件和多个数据库文件 (-b) 一次");这应该是进一步完善的良好起点。

我们可以定义一个函数，它接受任意数量的GRanges，并使用findOverlaps 识别第一个GRanges 和任何其他GRanges 之间的重叠范围；然后从pintersect 获得相交区域。

请注意，我使用了常见的tidyverse 语法；虽然这不是绝对必要的（对于每个 purrr::map/purrr::map2 函数，都可以使用它们的基本 R lapply/mapply 等效项），但我更喜欢 tidyverse 方法以提高代码可读性。

multiOverlap <- function(...) {
    require(GenomicRanges)
    require(tidyverse)

    # Store GRanges in list
    lst <- list(...)
    names(lst) <- paste0("gr", 1:length(lst))

    # Calculate mutual overlaps
    lst.matches <- map(lst[-1L], ~ findOverlaps(lst[[1L]], .x))

    # List of intersecting regions
    lst.gr <- map2(
        lst[-1L], lst.matches,
        ~pintersect(lst[[1]][queryHits(.y)], .x[subjectHits(.y)]))
    names(lst.gr) <- paste0("gr1-gr", 2:length(lst))

    # Convert GRanges to data.frame and reshape data
    map(lst.gr, ~.x %>%
        as.data.frame() %>%
        unite(locus, seqnames, start, sep = "-") %>%
        select(locus)) %>%
        bind_rows(.id = "id") %>%
        separate(id, into = c("grx", "gry")) %>%
        gather(gr, no, -locus) %>%
        transmute(
            locus,
            no,
            val = 1) %>%
            spread(no, val, fill = 0)
}

当我们将此函数应用于三个样本GRanges 时，我们得到以下结果

multiOverlap(gr.1, gr.2, gr.3)
#    locus gr1 gr2 gr3
#1  chr1-1   1   0   1
#2  chr1-2   1   1   0
#3 chr10-4   1   1   0

更新

另一个（快速）选项可能是使用data.table；尤其是在处理基因组数据时，data.tables 的传递引用属性，避免深度复制，使其非常有吸引力（而且速度很快）。

这是一个完全重现您预期输出的解决方案

# Load the library
library(data.table)

# Convert GRanges to data.table and row-bind entries
dt <- rbindlist(
    lapply(list(gr.1 = gr.1, gr.2 = gr.2, gr.3 = gr.3), as.data.table),
    idcol = "id")

# Remove width and strand
dt[, c("width", "strand") := NULL]

# Expand rows by range using start and end
dt <- dt[, .(pos = seq(start, end, by = 1L)), by = .(id, seqnames, grp = 1:nrow(dt))]

# Remove helper group label
dt[, grp := NULL]

# Unite seqnames and pos into one column
dt <- dt[, .(locus = do.call(paste, c(.SD, sep = "-")), id, pos), .SDcols = seqnames:pos]

# Add count variable
dt[, ct := 1]

# Convert from long to wide
dcast(dt, locus ~ id, value.var = "ct", fill = 0)
#     locus gr.1 gr.2 gr.3
#1:  chr1-1    1    0    1
#2:  chr1-2    1    1    0
#3: chr10-3    0    1    0
#4: chr10-4    1    1    0

我们就完成了

【讨论】：