r语言 - 求字符串序列之间的汉明距离

Question

我有一个包含3156个DNA序列的数据集，每个序列有98290个字符(SNP(，包括(通常(5个符号：A，C，G，T，N(间隙(。

找到这些序列之间的成对汉明距离的最佳方法是什么？

请注意，对于每个序列，我实际上想找到序列数(包括其自身(的倒数，其中每个站点的汉明距离小于某个阈值(在本例中为 0.1(。

到目前为止，我已经尝试了以下方法：

library(doParallel)
registerDoParallel(cores=8)
result <- foreach(i = 1:3156) %dopar% {
temp <- 1/sum(sapply(snpdat, function(x) sum(x != snpdat[[i]])/98290 < 0.1))
}

snpdat是一个list变量，其中snpdat[[i]]包含第i个 DNA 序列。

在具有 36GB 内存的核心 i7 - 4790 上运行大约需要 16 分钟。 我还尝试使用stringdist包，这需要更多时间来生成相同的结果。

任何帮助都非常感谢！

Answer 1

我不确定这是否是最佳解决方案，但我能够使用 Rcpp 将运行时间缩短到 15 分钟左右。我会在这里编写代码，以防有人有一天会发现它很有用......

这是C++代码(我在这里使用了糖运算符(...

#include <Rcpp.h>
using namespace Rcpp;
double test5(const List& C, const int& x){
double HD;
for(int i = 0; i < 3156; i++) if(sum(CharacterVector(C[x])!=CharacterVector(C[i])) < 9829) HD++;
return HD;
}

编译后：

library(Rcpp)
sourceCpp("hd_code.cpp")

我只是从 R 调用这个函数：

library(foreach)
library(doParallel)
registerDoParallel(cores = 8)
t =Sys.time()
bla = foreach(i = 1:3156, .combine = "c") %dopar% test5(snpdat,i-1)
Sys.time() - t

谁能想到更快的方法来做到这一点？