创建一个矩阵以显示多个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，并使用GRanges识别第一个GRanges和任何其他findOverlaps之间的重叠范围;然后从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

我们已经完成了:-)如果需要，可以很方便地在方便功能中包裹上面的行。