1

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

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

4

4 回答 4

4

这是一些数据

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

一个IRange概念是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（染色体）和链。

于 2015-01-31T22:41:27.150 回答
3

您正在寻找的数据结构称为区间树，它是一种排序二叉树，包含（猜猜是什么）区间，每个区间通常都有开始和结束位置。

我从来没有使用区间树来存储你需要的点，但我想你可以将你的区间定义为interval.start = interval.end.

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

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

于 2015-01-31T18:52:31.227 回答
1

我拿了那个矩阵并检查了重叠，其中只有五个有任何重叠的间隔，没有一个有 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
于 2015-01-31T23:31:56.110 回答
0

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

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

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

如果您只想将计数放在一个位置i，只需调用countChunks(i).

于 2015-01-31T18:00:40.223 回答