overlapping segments R



有一个我正在使用的数据框,它看起来像这样

这两列表示块的开始和结束。我需要知道从 0 到 23110906 的每个位置都存在多少个这样的块。有时块重叠,有时可能有一个区域根本没有块覆盖。这就像 R 中的段,但我不需要可视化,我只需要一种方法来快速找到每个位置的块数。有没有简单的方法?

这里有一些数据

m = matrix(c(10, 20, 25, 30), 2)

一个IRanges的概念是coverage()

> cvg = coverage(IRanges(start=m[,1], end=m[,2]))
> cvg
integer-Rle of length 30 with 4 runs
  Lengths:  9 10  6  5
  Values :  0  1  2  1

这是一个紧凑的运行长度编码;在第 i 个位置查询

> cvg[22]
integer-Rle of length 1 with 1 run
  Lengths: 1
  Values : 2
> runValue(cvg[22])
[1] 2

在Rle上做数学

> cvg > 1
logical-Rle of length 30 with 3 runs
  Lengths:    19     6     5
  Values : FALSE  TRUE FALSE

或强制到整数向量

> as(cvg, "integer")
 [1] 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 1 1 1 1 1

> cumsum(tabulate(m[,1], 30)) - cumsum(tabulate(m[,2], 30))
 [1] 0 0 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 1 1 1 1 1 0

也会相当快。

请注意它们之间的细微差异,与范围中是否包括端点的概念差异(IRanges:是;表格:否)。如果这些实际上是基因组坐标,那么GenomicRanges就是要去的地方,以解释seqname(染色体)和链。

您正在寻找的数据结构称为区间树,这是一种包含(猜猜看)区间的排序二叉树,每个区间通常都有开始和结束位置。

我从来没有使用间隔树来存储点,但我想你可以将间隔定义为interval.start = interval.end

构建树需要线性时间,查询数据框的间隔需要对数时间,这比 pteetor 的二次时间方法快得多。

Bioconductor 的 R 包 IRanges 可能会对您有所帮助。我会尝试findOverlaps()函数,然后table()结果。我邀请您阅读文档,看看它是否符合您的特定需求。

我拿了那个矩阵并检查了重叠,其中只有五个区间有任何重叠,没有 2 个区间,假设它们是按起始位置排序的:

> sum( mat[1:28,2] > mat[2:29,1] )
[1] 5
> sum( mat[1:27,2] > mat[3:29,1] )
[1] 0

那么他们是哪些呢?

> which( mat[1:28,2] > mat[2:29,1] )
[1] 19 21 23 25 28

因此,创建一个 2300 万个项目的向量似乎相当浪费机器资源和时间,并且简单地构建一个函数来计算任何特定位置所在的间隔数会容易得多:

 fchunk <- function(pos) {sum( mat[ , 1] <= pos & mat[,2] >= pos)}
#--------
> fchunk(16675330)
[1] 2
> fchunk(16675329)
[1] 1

这些是有 2 个的间隔:

sapply( which( mat[1:28,2] > mat[2:29,1] ) , 
       function(int1) c( mat[int1+1, 1], mat[int1, 2] ) )
#--------
       [,1]     [,2]     [,3]     [,4]     [,5]
n7 16675330 18097680 20233612 21288777 22847516
n8 16724700 18445265 20741145 22780817 22967567

如果你真的想要每个位置的计数 - 所有23,110,906个位置 - 这个代码会告诉你。

countChunks = function(i) sum(dfrm$n7 <= i & i <= dfrm$n8)
counts = sapply(1:23110906, countChunks)

但它非常慢。更快的代码需要一些巧妙的优化来消除这两行的(非常)冗余倒计时。

如果您只想在一个位置计数,i,只需调用countChunks(i)

相关内容

  • 没有找到相关文章

最新更新