识别两个序列之间的突变

Posted 2023-03-23

【中文标题】识别两个序列之间的突变

【英文标题】：Identifying mutations between two sequences

【发布时间】：2021-12-07 12:13:03

【问题描述】：

给定两个长度相等的序列（例如 DNA 序列），我希望能够找到它们之间的突变 - 序列内的类型和索引。例如，如果我输入序列AGGCTAC 和AGCCTTC，我想得到一个类似G3C, A6T 的列表。

我可以得到数量的差异就好了：

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

但我想不出比循环更优雅的方式来获得位置，这看起来很笨拙 - 我想借此机会学习。

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

【参考方案1】：

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

关于字符串比较问题，见this question。如果字符串长度相同，下面是一个可以工作的函数：

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"