从Rle向量高效构造GRanges/IRanges

Question

我有一个游程编码的向量，按顺序表示基因组上每个位置的一些值。举个简单的例子，假设我只有一个长度为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))

这是可行的，但速度非常慢。有没有更快的方法来构造相同的对象？

Answer 1

正如@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

在过去的几天里，我一直在断断续续地工作，我肯定不满意。我希望有人会接受我所做的(因为这是一种不同的方法)，并创造出更好的东西。