Эффективно построить GRange/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))

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


r bioinformatics bioconductor run-length-encoding iranges
person user1356855    schedule 30.08.2016    source источник
comment
вы можете использовать start(toyData)-1, чтобы получить начало интервала, но это не улучшит скорость.   -  person NicE    schedule 13.09.2016
comment
@NicE Спасибо за совет, даже если он не быстрее, он намного понятнее и лаконичнее.   -  person user1356855    schedule 15.09.2016
comment
‹Всю cumsum можно заменить на start(toyData)-1›   -  person The Unfun Cat    schedule 16.09.2016
comment
@user1356855 user1356855, с какой типичной длиной хромосом вы столкнетесь? Кроме того, будет ли достаточно 3 вариаций в вашем реальном приложении (например, вы могли бы иметь sample(1:15,10^8,replace=TRUE))?   -  person Joseph Wood    schedule 19.09.2016
comment
@JosephWood Да, значения, хранящиеся в геномных данных, часто являются реальными числами, поэтому 3 будет недостаточно ... Но я бы принял почти любой ответ. Самый длинный геном: например, 247 249 719? Хр1 люди...   -  person The Unfun Cat    schedule 19.09.2016
comment
Может быть, проблема не в преобразовании объекта Rle. Это функция IRanges/GRanges, которая вообще является медленной частью?   -  person Roman    schedule 19.09.2016
comment
@Джимбоу, я именно об этом и думаю. Мне интересно, как с этим справятся ребята из https://bioconductor.org/.   -  person Joseph Wood    schedule 19.09.2016

Ответы (1)


arrow_upward
3
arrow_downward

Как отметил @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)
}

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

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

person Joseph Wood    schedule 19.09.2016
comment
Это может означать, что метаданные OP содержат не векторизованные данные? Наличие объектов в векторах возможно в пандах, не знаю насчет R... - person The Unfun Cat; 19.09.2016
comment
Должен признаться, я тоже не совсем доволен. Похоже, что объект GRanges достаточно похож на вектор Rle (для одной хромосомы), поэтому этап построения должен быть практически мгновенным. Вместо этого это самая медленная часть моего кода. Очевидно, я недостаточно хорошо понимаю внутренности, чтобы понять, почему это неправильно/как сделать это быстрее. Однако альтернатива Rcpp хороша и дает некоторую дополнительную скорость. Спасибо! - person user1356855; 19.09.2016