我有这行 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 因为它们很容易生长(因为我事先不知道我要保留多少列)。通过索引,我是否需要在末尾对 informativeSites vector +1(因为 R 索引从 1 开始,C++ 索引从 0 开始)?
谢谢, 本 W.
最佳答案
示例数据:
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
使用 base 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++ 会打败它,但这真的有必要吗?
关于c++ - 在 C++ 和 Rcpp 中重写慢速 R 函数,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16555752/