6

我有这行 R 代码：

croppedDNA <- completeDNA[,apply(completeDNA,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))]

它的作用是识别 DNA 序列矩阵（1 行 = 一个序列）中不通用（信息丰富）的位点（列），并将它们从矩阵中子集化以形成一个新的“裁剪矩阵”，即去掉所有值相同的列。对于大型数据集，这大约需要 6 秒。我不知道我是否可以在 C++ 中更快地做到这一点（仍然是 C++ 的初学者），但尝试一下对我有好处。我的想法是使用 Rcpp，遍历 CharacterMatrix 的列，拉出列（站点）作为 CharacterVector 检查它们是否相同。如果它们相同，请记录该列号/索引，继续所有列。然后最后制作一个仅包含这些列的新 CharacterMatrix。重要的是，我保留行名和列名，因为它们在矩阵的“R 版本”中，即如果列出现，

我已经写了大约两分钟，到目前为止我所拥有的是（未完成）：

#include <Rcpp.h>
#include <vector>
using namespace Rcpp;
// [[Rcpp::export]]
CharacterMatrix reduce_sequences(CharacterMatrix completeDNA)
{
  std::vector<bool> informativeSites; 
  for(int i = 0; i < completeDNA.ncol(); i++)
  {
    CharacterVector bpsite = completeDNA(,i);
    if(all(bpsite == bpsite[1])
    {
      informativeSites.push_back(i);
    }
  }
CharacterMatrix cutDNA = completeDNA(,informativeSites);
return cutDNA;
}

我在这件事上走对了吗？有没有更简单的方法。我的理解是我需要 std::vector 因为它们很容易生长（因为我事先不知道我要保留多少列）。通过索引，我是否需要在最后对 informationativeSites 向量 +1（因为 R 索引从 1 开始，C++ 从 0 开始）？

谢谢，本·W。

4

1 回答 1

14

样本数据：

set.seed(123)
z <- matrix(sample(c("a", "t", "c", "g", "N", "-"), 3*398508, TRUE), 3, 398508)

OP的解决方案：

system.time(y1 <- z[,apply(z,2,function(x) any(c(FALSE,x[-length(x)]!=x[-1])))])
#    user  system elapsed 
#   4.929   0.043   4.976

使用基础 R 的更快版本：

system.time(y2 <- (z[, colSums(z[-1,] != z[-nrow(z), ]) > 0]))
#    user  system elapsed 
#   0.087   0.011   0.098

结果是相同的：

identical(y1, y2)
# [1] TRUE

c++很有可能会打败它，但真的有必要吗？

于 2013-05-15T02:52:13.943 回答