从Rle载体高效构建GRanges / IRanges
我有一个运行长度编码的向量，代表基因组上每个位置的某个值。 作为一个玩具的例子，假设我只有一个长度为10的染色体，那么我会看到一个向量
library(GenomicRanges)
set.seed(1)
toyData = Rle(sample(1:3,10,replace=TRUE))
我想将它强制为一个GRanges对象。 我能拿出最好的是
gr = GRanges('toyChr',IRanges(cumsum(c(0,runLength(toyData)[-nrun(toyData)])),
width=runLength(toyData)),
toyData = runValue(toyData))
哪个起作用，但相当缓慢。 有没有更快的方法来构建相同的对象？
正如@ TheUnfunCat指出的那样，OP的解决方案非常稳固。 下面的解决方案只比原来的解决方案稍微快一些。 我尝试了base R几乎所有组合，并且无法从S4Vectors软件包中击败效率Rle ，因此我使用了Rcpp 。 这是主要功能：
GenomeRcpp <- function(v) {
x <- WhichDiffZero(v)
m <- v[c(1L,x+1L)]
s <- c(0L,x)
e <- c(x,length(v))-1L
GRanges('toyChr',IRanges(start = s, end = e), toyData = m)
}
WhichDiffZero是Rcpp函数，几乎与base R which(diff(v) != 0)完全相同。 很多功劳归于@ G.Grothendieck。
#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
IntegerVector WhichDiffZero(IntegerVector x) {
int nx = x.size()-1;
std::vector<int> y;
y.reserve(nx);
for(int i = 0; i < nx; i++) {
if (x[i] != x[i+1]) y.push_back(i+1);
}
return wrap(y);
}
以下是一些基准：
set.seed(437)
testData <- do.call(c,lapply(1:10^5, function(x) rep(sample(1:50, 1), sample(1:30, 1))))
microbenchmark(GenomeRcpp(testData), GenomeOrig(testData))
Unit: milliseconds
expr min lq mean median uq max neval cld
GenomeRcpp(testData) 20.30118 22.45121 26.59644 24.62041 27.28459 198.9773 100 a
GenomeOrig(testData) 25.11047 27.12811 31.73180 28.96914 32.16538 225.1727 100 a
identical(GenomeRcpp(testData), GenomeOrig(testData))
[1] TRUE
过去几天我一直在努力解决这个问题，我绝对不满意。 我希望有人会采取我所做的（因为它是一种不同的方法）并创造更好的东西。
链接地址: http://www.djcxy.com/p/36041.html