使用完美散列将 k-mer 计入 R

Question

我从事 DNA K-mers 计数工作，我准备了这个公式来解决使用完美哈希的计数 table :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;
}

有两个问题：

假设一个kmer的索引为x，则x的补码为(4^k)-x-1。我们可以像前面的公式一样使用数字运算得到相反的结果吗？
在运行时间内有两个问题：迭代字符串和向量创建，其中 k 大于 8。有解决这些问题的思路吗？

Answer 1

项目 1: 我不精通 k-mer 计数，x 的定义有点不稳定。如果您可以更新您的问题，那么我将尝试解决第 1 点。但是，对我来说，您似乎只需要使用一种从 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;
}

使用完美散列将 k-mer 计入 R

k-mer counting into R using perfect hashing

c++

math

r

bioinformatics

rcpp