求字符串序列之间的汉明距离答案

【问题标题】：Find the Hamming distance between string sequences求字符串序列之间的汉明距离
【发布时间】：2020-01-28 05:04:04
【问题描述】：

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

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

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

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

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]] 包含ith DNA 序列。

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

非常感谢任何帮助！

【问题讨论】：

一个快速的谷歌想出了这个：johanndejong.wordpress.com/2015/09/23/…
您可以考虑使用专门的函数，例如phangorn 中的dist.hamming，或者Bioconductor 上的Biostrings 包。
感谢您的建议。 @GordonShumway，不幸的是，我的数据集对于矩阵乘法来说太大了......
感谢您的建议。 @Axeman，我现在就试一试……

标签： r arrays string vectorization hamming-distance

【解决方案1】：

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

这是 C++ 代码（我在这里使用了Sugar 运算符）...

#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

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

【讨论】：