r语言 - 识别两个序列之间的突变

Question

给定长度相等的两个序列(例如DNA序列)，我希望能够找到它们之间的突变—序列中的类型和索引。作为一个例子，如果我输入序列AGGCTAC和AGCCTTC，我想得到一个像G3C, A6T这样的列表。

我可以得到数的差异很好:

seqs <- Biostrings::readAAStringSet("PSE-1_Round20.fas")
nmut <- adist(seqs[1], seqs)

但是我想不出比循环更优雅的方法来获取位置，这看起来很笨拙——我想把这作为一个学习的机会。

我正在使用R中的Biostrings包，但是对于我想要做的任何事情，该包中似乎没有任何特殊的工具，我认为任何适用于通用字符串的解决方案也应该适用于我。事实上，如果有更优雅的python或bash脚本解决方案，我也会接受。

Answer 1

似乎有多个包应该这样做。一个是adegenet包中的findMutations函数。

关于字符串比较问题，参见这个问题。如果字符串长度相同，下面的函数将工作:

mutations <- function(str1, str2) {
str1vec <- unlist(strsplit(str1, ""))
str2vec <- unlist(strsplit(str2, ""))
iMut <- (1:length(str1vec))[str1vec != str2vec]
return(paste0(str1vec[iMut], iMut, str2vec[iMut]))
}
> mutations("AGGCTAC", "AGCCTTC")
[1] "G3C" "A6T"