创建一个矩阵以显示多个 GRanges 之间的重叠

Question

我正在尝试找到一种方法来在比较不同的 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)

有什么明智的导出方法吗？

Answer 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.table 的引用传递属性，避免深度复制，使其非常有吸引力（而且速度很快）。

这是一个可以准确再现您预期输出的解决方案

# 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

我们完成了:-) 如有必要，可以很容易地将上面几行换成一个方便的函数。

创建一个矩阵以显示多个 GRanges 之间的重叠

Create a matrix to show overlaps among multiple GRanges

r

bioinformatics

bioconductor

更新