如何从fasta文件中删除重复项,但每个组基于头至少保留一个

2024-05-06 18:04:58 发布

您现在位置:Python中文网/ 问答频道 /正文

我有一个multifasta文件,它看起来像这样:

(所有序列均大于100bp,多行,长度相同)

>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage3_samplenameD
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA

我需要删除重复项,但至少保留每个谱系的序列。因此,在上面这个简单的示例中(注意samplenameA、C和D是相同的),我只想删除samplenameD或samplenameC,但不想同时删除它们。最后,我希望获得与原始文件中相同的头信息

示例输出:

>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA

我找到了一种只删除重复项的方法。感谢皮埃尔·林登鲍姆

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_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG

->;所以失去了血统3序列

现在,我正在寻找一种快速解决方案,以删除重复项,但基于fasta头,每个谱系至少保留一个序列

我不熟悉脚本编写。。。欢迎使用bash/python/R中的任何想法

谢谢


Tags: 文件示例序列sedtr谱系samplenameasamplenameb
1条回答
网友
1楼 · 发布于 2024-05-06 18:04:58

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

逐行循环文件,并将沿袭/序列数据添加到字典中。我建议使用序列作为键。这样,您就可以很容易地知道是否已经遇到了此密钥

myfasta = {}
if myfasta[sequence]:
    myfasta[sequence].append(lineage_id)
else:
    myfasta[sequence] = [lineage_id]

这样,您的密钥(序列)将保存具有相同序列的沿袭ID列表。请注意,这个解决方案令人恼火的地方是在文件上循环,将沿袭id与序列分开,解释可能扩展到多行的序列,等等

之后,您可以循环字典,并仅使用字典中列表中的第一个沿袭id将序列写入文件

相关问题 更多 >