k-mer使用完美散列法計入R

Question

我從事DNA K-mers計數工作，並準備使用完美散列表來計算這個公式：link。 我用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; 
} 

有兩個問題：

1. 假設k聚體的索引x，x的稱讚然後（4^K）-x-1。我們可以使用像上述公式那樣的數字操作來獲得反向嗎？

2. 運行時間有兩個問題：對字符串和矢量創建進行迭代，其中k超過8。 有沒有想法來解決這些問題？

Answer 1

項目1：我不是在k鏈節計數精通和x的定義是有點不穩。如果你可以更新你的問題，那麼我會嘗試解決第一點。但是，對我來說，好像你只需要使用一個算法，從基數10轉換回基數4？

項目2：是的。正在執行的迭代並不理想，因爲您不斷從string轉換爲int等等。此外，爲什麼這些功能被宣佈爲inline？另外，數據傳輸到A有多大？

爲了解決這個問題，我們將字符串全部移植到std::vector<int>。事實上，我會說只是將其直接移植到std::vector<int>，並避免列出它可以切換到非rcpp代碼。另外，我們選擇使用通過引用範例而不是const變量。最後，我減少了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; 
}