我正在使用一个数据框，它看起来像这样
两列表示块的开始和结束。我需要知道从 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 tree，它是一种包含（猜猜看）区间的排序二叉树，每个区间通常都有开始和结束位置。
我从未根据需要使用区间树来存储点，但我想您可以将区间定义为 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)
但是速度很慢。更快的代码需要一些巧妙的优化来消除这两行（非常）冗余的倒计时。
如果您只想在 one 位置计数，i，只需致电 countChunks(i)。
【讨论】：