Эффективно построить GRanges/IRanges из вектора Rle
У меня есть вектор длины, закодированный, представляющий некоторое значение в каждой позиции в геноме, по порядку. В качестве примера игрушки предположим, что у меня была только одна хромосома длиной 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))
который работает, но довольно медленно. Есть ли более быстрый способ построить тот же объект?
1 ответ
Как отметил @TheUnfunCat, решение ОП довольно солидно. Решение ниже только умеренно быстрее, чем оригинальное решение. Я попробовал почти каждую комбинацию 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)
}
WhichDiffZero
это Rcpp
функция, которая в значительной степени делает то же самое, что и 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
Я работал над этим в течение последних нескольких дней, и я определенно не удовлетворен. Я надеюсь, что кто-то возьмет то, что я сделал (поскольку это другой подход), и создаст что-то намного лучшее.