4

我有一个长度为 5,000,000 的数字向量

>head(coordvec)
[1] 47286545 47286546 47286547 47286548 47286549 472865

和一个 3 x 1,400,000 的数字矩阵

>head(subscores)
        V1       V2     V3
1 47286730 47286725  0.830
2 47286740 47286791  0.065
3 47286750 47286806 -0.165
4 47288371 47288427  0.760
5 47288841 47288890  0.285
6 47288896 47288945  0.225

我想要完成的是，对于 coordvec 中的每个数字，找到 V1 和 V2 包含 coordvec 中的数字的子分数中的行的 V3 平均值。为此，我采用以下方法：

results<-numeric(length(coordvec))
for(i in 1:length(coordvec)){
    select_rows <- subscores[, 1] < coordvec[i] & subscores[, 2] > coordvec[i]
scores_subset <- subscores[select_rows, 3]
results[m]<-mean(scores_subset)
}

这非常慢，需要几天才能完成。有更快的方法吗？

谢谢，

担

4

3 回答 3

6

我认为这个问题有两个具有挑战性的部分。首先是寻找重叠。我会使用BioconductorIRanges的包（在基本包中也可能有用）?findInterval

library(IRanges)

创建表示坐标向量的宽度为 1 的范围，以及表示分数的范围集；为方便起见，我对坐标向量进行排序，假设重复的坐标可以被同等对待

coord <- sort(sample(.Machine$integer.max, 5000000))
starts <- sample(.Machine$integer.max, 1200000)
scores <- runif(length(starts))

q <- IRanges(coord, width=1)
s <- IRanges(starts, starts + 100L)

在这里我们找到哪些query重叠subject

system.time({
    olaps <- findOverlaps(q, s)
})

这在我的笔记本电脑上大约需要 7 秒。有不同类型的重叠（请参阅 参考资料?findOverlaps），因此这一步可能需要一些改进。结果是一对索引查询和重叠主题的向量。

> olaps
Hits of length 281909
queryLength: 5000000
subjectLength: 1200000
       queryHits subjectHits 
        <integer>   <integer> 
 1             19      685913 
 2             35      929424 
 3             46     1130191 
 4             52       37417

我认为这是第一个复杂部分的结尾，找到 281909 重叠。（我不认为其他地方提供的 data.table 答案解决了这个问题，尽管我可能弄错了......）

下一个具有挑战性的部分是计算大量均值。内置方式类似于

olaps0 <- head(olaps, 10000)
system.time({
    res0 <- tapply(scores[subjectHits(olaps0)], queryHits(olaps0), mean)
})

这在我的计算机上大约需要 3.25 秒，并且似乎是线性缩放的，所以 280k 重叠可能需要 90 秒。但我认为我们可以使用data.table. 原坐标start(v)[queryHits(olaps)]为

require(data.table)
dt <- data.table(coord=start(q)[queryHits(olaps)],
                 score=scores[subjectHits(olaps)])
res1 <- dt[,mean(score), by=coord]$V1

所有 280k 重叠大约需要 2.5 秒。

通过识别查询命中是有序的，可以加快速度。我们想要计算每次查询命中的平均值。我们首先创建一个变量来指示每个查询命中运行的结束

idx <- c(queryHits(olaps)[-1] != queryHits(olaps)[-length(olaps)], TRUE)

然后计算每次跑步结束时的累积分数，每次跑步的长度，以及跑步结束时和开始时累积分数之间的差值

scoreHits <- cumsum(scores[subjectHits(olaps)])[idx]
n <- diff(c(0L, seq_along(idx)[idx]))
xt <- diff(c(0L, scoreHits))

最后，平均值是

res2 <- xt / n

所有数据大约需要 0.6 秒，并且与 data.table 结果相同（尽管比？）

> identical(res1, res2)
[1] TRUE

均值对应的原始坐标为

start(q)[ queryHits(olaps)[idx] ]
于 2013-01-20T01:27:20.003 回答
2

这样的事情可能会更快：

require(data.table)
subscores <- as.data.table(subscores)

subscores[, cond := V1 < coordvec & V2 > coordvec]
subscores[list(cond)[[1]], mean(V3)]

list(cond)[[1]]因为：“当 i 是单个变量名时，它不被视为列名的表达式，而是在调用范围内进行评估。” 来源：?data.table

于 2013-01-20T00:00:22.560 回答
0

由于您的答案不容易重现，即使是，您也没有一个subscores符合您的布尔条件，我不确定这是否完全符合您的要求，但您可以使用其中一个apply系列和函数。

myfun <- function(x) {
  y <- subscores[, 1] < x & subscores[, 2] > x
  mean(subscores[y, 3])
}

sapply(coordvec, myfun)

你也可以看看mclapply。如果你有足够的内存，这可能会显着加快速度。但是，您也可以查看foreach具有类似结果的包。通过分配而不是增长它，你已经得到了你的for loop“正确” results，但实际上，你正在做很多比较。很难加快速度。

于 2013-01-19T23:51:11.630 回答