Эффективное построение GRanges/IRanges из Rle-вектора

Question

У меня есть кодированный вектор длины пробега, представляющий некоторое значение в каждой позиции на геноме, по порядку. В качестве примера предположим, что я игрушка была только одна хромосома длиной 10, то я бы вектор похожий

library(GenomicRanges) 

set.seed(1) 
toyData = Rle(sample(1:3,10,replace=TRUE)) 

Я хотел бы, чтобы заставить это в объект Усадьбы. Лучшее, что я могу придумать, это

gr = GRanges('toyChr',IRanges(cumsum(c(0,runLength(toyData)[-nrun(toyData)])), 
           width=runLength(toyData)), 
      toyData = runValue(toyData)) 

который работает, но довольно медленный. Есть ли более быстрый способ построения одного и того же объекта?

Answer 1

Как отметил @TheUnfunCat, решение OP довольно прочное. Решение ниже только умеренно быстрее, чем исходное решение. Я пробовал почти каждую комбинацию base R и не смог побить эффективность Rle из пакета S4Vectors, таким образом, я обратился к 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) 
} 

WhichDiffZeroRcpp является функцией, которая в значительной степени делает точно такую ​​же вещь, как which(diff(v) != 0) в base R. Большой кредит - @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 

я работаю над этим и выключения в течение последних нескольких дней, и я определенно не удовлетворены. Я надеюсь, что кто-то возьмет то, что я сделал (поскольку это другой подход) и создать что-то намного лучше.