如何从 fasta 文件中删除重复项,但基于 header 每组至少保留一个
How to remove duplicates from fasta file but keep at least one per group based on header
我有一个 multifasta 文件,如下所示:
(所有序列>100bp,多于一行,长度相同)
>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage3_samplenameD
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
我需要删除重复项,但至少要保持每个谱系的顺序。因此,在上面的这个简单示例中(注意 samplenameA、C 和 D 是相同的),我只想删除 samplenameD 或 samplenameC,而不是同时删除它们。最后我想获得与原始文件中相同的 header 信息。
示例输出:
>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 header.
每个谱系至少保留一个序列
我是脚本新手...欢迎 bash/python/R 中的任何想法。
谢谢!!!
在这种情况下,我可以看到两个比较好的替代方案。 A) 查看现有工具(例如 Biopython 库或 FASTX 工具包。我 认为 它们都有很好的命令来完成这里的大部分工作,所以它可能值得学习它们。或者,B) 编写你自己的。在这种情况下,您可能想尝试(我会坚持 python):
逐行遍历文件,并将 lineage/sequence 数据添加到字典中。我建议使用序列作为键。这样,您就可以很容易地知道您是否已经遇到过这个密钥。
myfasta = {}
if myfasta[sequence]:
myfasta[sequence].append(lineage_id)
else:
myfasta[sequence] = [lineage_id]
这样你的键(序列)将保存具有相同序列的 lineage_ids 列表。请注意,此解决方案的烦人之处在于遍历文件、将 lineage-id 与序列分开、考虑可能扩展到多行的序列等。
之后,您可以遍历字典,并仅使用字典中列表中的第一个 lineage_id 将序列写入文件。
我有一个 multifasta 文件,如下所示:
(所有序列>100bp,多于一行,长度相同)
>Lineage1_samplenameA
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage2_samplenameB
AAATTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAG
>Lineage3_samplenameC
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
>Lineage3_samplenameD
CGCTTCAACGGAATGGATCTACGTTACAGCCTGCATAAAGAAAACGGAGTTGCCGAGGACGAAAGCGACTTTAGGTTCTGTCCGTTGTCTTTGGCGGAAA
我需要删除重复项,但至少要保持每个谱系的顺序。因此,在上面的这个简单示例中(注意 samplenameA、C 和 D 是相同的),我只想删除 samplenameD 或 samplenameC,而不是同时删除它们。最后我想获得与原始文件中相同的 header 信息。
示例输出:
>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 header.
每个谱系至少保留一个序列我是脚本新手...欢迎 bash/python/R 中的任何想法。
谢谢!!!
在这种情况下,我可以看到两个比较好的替代方案。 A) 查看现有工具(例如 Biopython 库或 FASTX 工具包。我 认为 它们都有很好的命令来完成这里的大部分工作,所以它可能值得学习它们。或者,B) 编写你自己的。在这种情况下,您可能想尝试(我会坚持 python):
逐行遍历文件,并将 lineage/sequence 数据添加到字典中。我建议使用序列作为键。这样,您就可以很容易地知道您是否已经遇到过这个密钥。
myfasta = {}
if myfasta[sequence]:
myfasta[sequence].append(lineage_id)
else:
myfasta[sequence] = [lineage_id]
这样你的键(序列)将保存具有相同序列的 lineage_ids 列表。请注意,此解决方案的烦人之处在于遍历文件、将 lineage-id 与序列分开、考虑可能扩展到多行的序列等。
之后,您可以遍历字典,并仅使用字典中列表中的第一个 lineage_id 将序列写入文件。