一种快速的线性区间重叠查找方法

Question

我有一个线性间隔的data.frame (映射的RNA-seq读取的基因组坐标)，例如：

df <- data.frame(seqnames = c(rep("chr10",2),rep("chr5",8)),
                 start = c(12255935,12257004,12243635,12244009,12253879,12254395,12254506,12255142,12255229,12258719),
                 end = c(12257002,12258512,12243764,12244291,12254107,12254501,12254515,12255535,12255312,12258764),
                 read_id = c(rep("R9",2),rep("R10",8)),
                 stringsAsFactors = F)

对于某些读取，有包含或交叉在其他读取的间隔，我想合并它们。在上面针对read_id = "R10"的示例中，interval：chr5 12255229 12255312包含在interval chr5 12255142 12255535中。

对于单个读data.frame，我使用以下过程：

#defining helper functions
clusterHits <- function(overlap.hits)
{
  overlap.hits <- GenomicRanges::union(overlap.hits,t(overlap.hits))
  query.hits <- S4Vectors::queryHits(overlap.hits)
  search.hits <- S4Vectors::subjectHits(overlap.hits)
  cluster.ids <- seq_len(S4Vectors::queryLength(overlap.hits))
  while(TRUE){
    hit <- S4Vectors::Hits(query.hits,cluster.ids[search.hits],S4Vectors::queryLength(overlap.hits),S4Vectors::subjectLength(overlap.hits))
    tmp.cluster.ids <- pmin(cluster.ids,S4Vectors::selectHits(hit,"first"))
    if(identical(tmp.cluster.ids,cluster.ids))
      break
    cluster.ids <- tmp.cluster.ids
  }
  unname(S4Vectors::splitAsList(seq_len(S4Vectors::queryLength(overlap.hits)),cluster.ids))
}

mergeConnectedRanges <- function(x.gr,overlap.hits)
{
  cluster.ids <- clusterHits(overlap.hits)
  merged.gr <- range(IRanges::extractList(x.gr,cluster.ids))
  merged.gr <- unlist(merged.gr)
  S4Vectors::mcols(merged.gr)$merged.idx <- cluster.ids
  return(merged.gr)
}

#Now separate R10 and merge its intervals
df1 <- dplyr::filter(df, read_id == "R10")
gr <- GenomicRanges::GRanges(dplyr::select(df1,seqnames,start,end))
redundant.intervals <- GenomicRanges::findOverlaps(gr,ignore.strand=T)
query.gr <- redundant.intervals[S4Vectors::queryHits(redundant.intervals)]
subject.gr <- redundant.intervals[S4Vectors::subjectHits(redundant.intervals)]
as.data.frame(mergeConnectedRanges(x.gr=gr,overlap.hits=redundant.intervals))

这意味着：

  seqnames    start      end width strand merged.idx
1     chr5 12243635 12243764   130      *          1
2     chr5 12244009 12244291   283      *          2
3     chr5 12253879 12254107   229      *          3
4     chr5 12254395 12254501   107      *          4
5     chr5 12254506 12254515    10      *          5
6     chr5 12255142 12255535   394      *       6, 7
7     chr5 12258719 12258764    46      *          8

因此，merged.idx显示df1中的间隔6和7已经合并。

我正在寻找一种快速的方法，在成千上万的阅读中做到这一点。最明显的方法是在do.call中使用df中唯一的读取。

library(dplyr)
do.call(rbind, lapply(unique(df$read_id), function(r){
  read.df <- dplyr::filter(df, read_id == r)
  gr <- GenomicRanges::GRanges(dplyr::select(read.df,seqnames,start,end))
  redundant.intervals <- GenomicRanges::findOverlaps(gr,ignore.strand=T)
  query.gr <- redundant.intervals[S4Vectors::queryHits(redundant.intervals)]
  subject.gr <- redundant.intervals[S4Vectors::subjectHits(redundant.intervals)]
  as.data.frame(mergeConnectedRanges(x.gr=gr,overlap.hits=redundant.intervals)) %>%
    dplyr::mutate(read_id = r)
}))

但我想知道有没有更快的方法。请注意，实际具有这种相交间隔的读取的部分相对较小。

Answer 1

使用来自生物导体存储库的GenomicRanges包，可以通过几行代码来完成任务：

library(GenomicRanges)
makeGRangesListFromDataFrame(df, split.field = "read_id") |>
  reduce(with.revmap = TRUE) |>
  as.data.frame()

group_name公司名称启动端宽链改造图1 1 R10 chr5 12243635 12243764 130 2 1 R10 chr5 12244009 12244291 283 2 3 1 R10 chr5 12253879 12254107 229 *3 4 R10 chr5 12254395 12254501 107 *4 5 1 R10 chr5 12254506 12254515* R10 chr5 12255142 12255535 394 * 6，7 7 1 R10 chr5 12258719 12258764 46 *8 8 2 R9 chr10 12255935 12257002 1068 *1 9 R9 chr10 12257004 12258512 1509 *2

由于GenomeRanges包不在CRAN上，请参阅安装和管理生物导体封装或运行

install.packages("BiocManager")
BiocManager::install("GenomicRanges")

数据

df <- data.frame(seqnames = c(rep("chr10", 2), rep("chr5", 8)),
                 start = c(12255935, 12257004, 12243635, 12244009, 12253879, 12254395, 12254506, 12255142, 12255229, 12258719),
                 end   = c(12257002, 12258512, 12243764, 12244291, 12254107, 12254501, 12254515, 12255535, 12255312, 12258764),
                 read_id = c(rep("R9", 2), rep("R10", 8)), 
                 stringsAsFactors = FALSE)

作为