☰

如何通过命令行工具显示两个 DNA 序列之间的差异

tag-icon • tag-icon shell-script text-processing bioinformatics
如何通过命令行工具显示两个 DNA 序列之间的差异

我有以下问题：

我的数据表看起来像这样

AAAGGGTTT AAAGGG
AAAGGGCCC GGGCCC

我想像这样在第三行显示两个序列之间的差异

AAAGGGTTT AAAGGG TTT
AAAGGGCCC GGGCCC AAA

我尝试过使用 diff。我提取了文件（f1.txt 和 f2.txt）中的各个序列并对其进行了格式化，以便可以将它们与 diff 逐行进行比较，这产生了一个问题，即它仅在序列的开头相似时才起作用（数据表的第 1 行） ）。

awk '{gsub(".","&\n");printf "%s",$0}' < f1.txt >f1a.txt
awk '{gsub(".","&\n");printf "%s",$0}' < f2.txt >f2a.txt
 
diff -y f1a.txt f2a.txt

有谁知道如何实现所需的输出？

答案1

这就是你所追求的吗？

awk '{$3=$1;sub($2,"",$3)}1' file

  • $3=$1将第一个字段复制到第三个字段并

  • sub($2,"",$3)在第三个字段中查找第二个字段。如果存在匹配，sub则用第三个字段中的空字符串替换匹配的字符串。

  • 最后1打印结果。它相当于一个{print}语句，因此您可以将其重写为{$3=$1;sub($2,"",$3);print}.

结果：

AAAGGGTTT AAAGGG TTT
AAAGGGCCC GGGCCC AAA

答案2

如果您需要显示序列的成对比对，请使用适当的生物信息学工具。序列比对将显示典型用户期望的格式差异。在这里，对于核苷酸序列的成对比对，您可以使用例如blastnBLAST 工具（例如，可以使用 来安装conda），请从线程中查看以下命令：https://www.biostars.org/p/18087/#18095:

blastn -query querySeqSet.fasta -subject targetSeqSet.fasta

相关内容