K-мер, считая в R с использованием идеального хеширования

Question

K-мер, считая в R с использованием идеального хеширования

Я работал в подсчете K-mers ДНК и готовлю эту формулировку, чтобы решить подсчет, используя идеальную хеш-таблицу:ссылка на сайт.
Я использовал Rcpp API (C ++) для интеграции кода в R:

#include <Rcpp.h>
using namespace Rcpp;
/*
this code can be used with c++ by replacing IntegerVector by std::vector<int>
*/
//************************************************
inline const short int V (const char x){
switch(x){
case 'A':case 'a':
return 0;
break;
case 'C':case 'c':
return 1;
break;
case 'G':case 'g':
return 2;
break;
default:
return 3;
break;
}

}

inline unsigned int  X0( const std::string A,const int k ,const int n){
unsigned int  result=0;
int j=k;
for( int i=n-1;i>n-k-1;i--) {
result+= pow(4,k-j)*V(A[i]);
j--;
}
return result;
}

// [[Rcpp::export]]
inline IntegerVector kmer4(const std::string A,const int n,const int k)
{

IntegerVector P(pow(4,k));
int x=X0(A,k,n);
P[x]++;
const int N=pow(4,k-1);
for( int i=n-k-1;i>-1;i--){
x=N*V(A[i])+x/4-x%4/4;
P[x]++;
}
return P;
}

Есть два вопроса:

Предполагая индекс х кмера, комплимент х тогда будет (4 ^ k) -x-1. Можем ли мы получить обратное, используя числовую операцию, как в предыдущей формуле?
Есть две проблемы во время выполнения: итерация по строке и создание вектора, где k больше 8.
Есть идеи для решения этих проблем?

-5

bioinformatics c++math r rcpp

Решение

Другие решения

Других решений пока нет …

Источник

Accepted Answer

Элемент 1: Я не разбираюсь в подсчете k-мер и определении x немного шатко Если бы вы могли обновить свой вопрос, то я попытаюсь обратиться к пункту 1. Но мне кажется, что вам просто нужно использовать алгоритм, который преобразует из базы 10 обратно в базу 4?

Пункт 2: Да. Выполняемая итерация является неоптимальной, так как вы постоянно конвертируете из string в int и так далее. Кроме того, почему функции объявляются inline? Кроме того, насколько большие данные передаются в A?

Чтобы решить эту проблему, мы переносим строку полностью в std::vector<int>, На самом деле, я бы сказал, просто перенести его прямо на std::vector<int> и избегайте перечисления, что это может быть переключено на код не-rcpp. Кроме того, мы решили использовать парадигму передачи по ссылке вместо просто константной переменной. Наконец, я сократил количество inline декларация.

#include <Rcpp.h>
using namespace Rcpp;

//************************************************
inline const short int V (const char x){
switch(x){
case 'A':case 'a':
return 0;
break;
case 'C':case 'c':
return 1;
break;
case 'G':case 'g':
return 2;
break;
default:
return 3;
break;
}

}

std::vector<int> conv_A2V(const std::string & A){
unsigned int obs = A.length()
std::vector<int> out(obs);
for(unsigned int i = 0; i < obs; i++){
out(i) = V(A[i]);
}
return out;
}

unsigned int X0( const std::vector<int> & V_A, const int k, const int n){
unsigned int  result=0;
int j=k;
for( int i=n-1;i>n-k-1;i--) {
result+= pow(4,k-j)*V_A[i];
j--;
}
return result;
}

// [[Rcpp::export]]
IntegerVector kmer4(const std::string A, const int n,const int k)
{
// Convert
std::vector<int> V_A = conv_A2V(A);

IntegerVector P(pow(4,k));
int x=X0(V_A,k,n);
P[x]++;
const int N=pow(4,k-1);
for( int i=n-k-1;i>-1;i--){
x=N*V_A[i]+x/4-x%4/4;
P[x]++;
}
return P;
}

0