子集在大矩阵中慢

Question

我有一个长度为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中的每个数字,找到子行中行的V3的平均值,其中V1和V2包含coordvec中的数字.为此,我采取以下方法:

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)
}

这很慢,需要几天才能完成.有更快的方法吗？

谢谢,

担

Answer 1

我认为这个问题有两个具有挑战性的部分.首先是找到重叠.我使用Bioconductor的IRanges包(在基础包中也可能有用)?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] ]