Efektywnie skonstruuj GRanges/IRanges z Rle wektorowego

Question

Mam wektor zakodowany w długości Run reprezentujący pewną wartość w każdej pozycji w genomie, w kolejności. Jako przykład zabawki przypuszczać miałem tylko jeden chromosom o długości 10, wtedy mam wektor patrząc jak

library(GenomicRanges) 

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

chciałbym zmusić to do obiektu folwarków. Najlepsze, co mogę wymyślić, to:

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

który działa, ale działa dość wolno. Czy istnieje szybszy sposób na skonstruowanie tego samego obiektu?

Answer 1

Jak zaznaczył @TheUnfunCat, rozwiązanie OP jest całkiem solidne. Poniższe rozwiązanie jest tylko umiarkowanie szybsze niż oryginalne rozwiązanie. Próbowałem prawie każdej kombinacji base R i nie mogłem pokonać wydajności Rle z pakietu S4Vectors, dlatego uciekłem się do Rcpp. Tutaj jest główną funkcją:

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 jest funkcja Rcpp że dość dużo robi dokładnie to samo jak which(diff(v) != 0) w base R. Dużą zasługą jest @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); 
} 

Poniżej przedstawiamy kilka punktów odniesienia:

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 

Pracuję na to wyłączyć i przez ostatnie kilka dni, a ja na pewno nie jestem zadowolony. Mam nadzieję, że ktoś weźmie to, co zrobiłem (ponieważ jest to inne podejście) i stworzy coś znacznie lepszego.