如何从 fasta 文件中删除重复项，但根据标题每组至少保留一个

Question

我有一个如下所示的 multifasta 文件：（所有序列都>100bp，多于一行，且长度相同）>Lineage1_samplenameACGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA>Lineage2_samplenameBAAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG>Lineage3_samplenameCCGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA>Lineage3_samplenameDCGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA我需要删除重复项，但至少要保持每个谱系的顺序。因此，在上面的这个简单示例中（注意 samplenameA、C 和 D 是相同的），我只想删除 samplenameD 或 samplenameC，而不是同时删除它们。最后我想获得与原始文件中相同的标题信息。示例输出：>Lineage1_samplenameACGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA>Lineage2_samplenameBAAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG>Lineage3_samplenameCCGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA我找到了一种只删除重复项的方法。感谢皮埃尔林登鲍姆。sed -e '/^>/s/$/@/' -e 's/^>/#/'file.fasta  |\tr -d '\n' | tr "#" "\n" | tr "@""\t" |\sort -u -t '  ' -f -k 2,2  |\sed -e 's/^/>/' -e 's/\t/\n/'在我上面的例子中运行它会导致：>Lineage1_samplenameACGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA>Lineage2_samplenameBAAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG—> 所以失去了血统 3 序列现在我只是在寻找一种快速解决方案来删除重复项，但基于 fasta 标头每个谱系至少保留一个序列。我是脚本新手……欢迎使用 bash/python/R 中的任何想法。

Answer 1

在这种情况下，我可以看到两个相对较好的选择。A) 查看现有工具（例如 Biopython 库或 FASTX 工具包。我认为它们都有很好的命令来完成此处的大部分工作，因此学习它们可能是值得的。或者，B) 编写您自己的工具。在这种情况下，您可能想尝试（我会坚持使用 python）：

逐行遍历文件，并将谱系/序列数据添加到字典中。我建议使用序列作为键。这样，您可以很容易地知道您是否已经遇到过此密钥。

myfasta = {}

if myfasta[sequence]:

    myfasta[sequence].append(lineage_id)

else:

    myfasta[sequence] = [lineage_id]

这样你的键（序列）将保存具有相同序列的 lineage_ids 列表。请注意，此解决方案的烦人之处在于遍历文件、将 lineage-id 与序列分开、考虑可能扩展到多行的序列等。

之后，您可以遍历字典，并仅使用字典中列表中的第一个 lineage_id 将序列写入文件。