# Python: How to write duplicate sequences removed from fasta file to new file

I currently am using this code to remove duplicate sequences from the fasta file. However, I would also like to write a new file with only the removed duplicates as well as a count for how many times they appear. Please let me know if you have any suggestions, I am hoping to keep it written in one py file, but I am not opposed to running a second file to do so. I attended to write the appended strings to a new file but it looked slightly off...

from Bio import SeqIO
import time

start = time.time()

seen = []
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
if str(record.seq) not in seen:
seen.append(str(record.seq))
records.append(record)

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

print(f"Run time is {(end- start)/60}")
$$$$

• Hi @User this has been nicely answered from different perspectives. Could you please send off one of the answers to this question as "accepted". The site has certainly helped here, so to reciprocate that help, this certainly helps.
– M__
Jul 8 '20 at 7:29

Edit: user @RamRS has a more robust answer for actually using this. Please use that. Leaving answer up as OP appears to be interested in implementing as a learning activity.

Please let us know why this is not doing what you intend.

I haven't checked it, but something like this would work:

from Bio import SeqIO
import time
from collections import defaultdict

start = time.time()

seq_counts = defaultdict(0)
records = []

for record in SeqIO.parse("b4r2.fasta", "fasta"):
seq = str(record.seq)
if seq in seq_counts:  # python is smart about dict keys
records.append(record)

seq_counts[seq] += 1

#writing to a fasta file
SeqIO.write(records, "no_dupes_b4r2.fasta", "fasta")
end = time.time()

# print the counts
for seq in seq_counts:
print(seq_counts[seq], seq)

print(f"Run time is {(end- start)/60}")


It will be a bit slow for large fasta files and I am sure that some people will be very upset that I'm using a defaultdict for this.

The code that you have would probably not work as you would want the records.append() call to be only run in the if statement testing for non-membership.

There is absolutely no reason to have to write more than one script for this.

• There is absolutely no need to write even this script when tools like seqkit rmdup already exist. Jul 6 '20 at 22:17
• Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.
– M__
Jul 6 '20 at 23:50
• @Chris_Rands yeah totally, i just know the defaultdict syntax by heart and i'd have to go look up how to use a counter again... Jul 7 '20 at 17:23
• @RamRS all right then, that is what they should do if they're using this in earnest. On the other hand, it does seem like they want to implement it themselves as an exercise, so I'll leave my answer up for now with an edit. Jul 7 '20 at 22:35
• Agree with you 100%. This is great as a learning exercise (who among us did not write our own fasta deduplicator?), but these robust tools, as you rightly call them, are better tested so they catch edge cases that our code might not look out for. Jul 8 '20 at 16:14

Please do not reinvent the wheel - seqkit rmdup does what you want.

From the website:

# rmdup

## Usage

remove duplicated sequences by id/name/sequence

Usage:
seqkit rmdup [flags]

Flags:
-n, --by-name                by full name instead of just id
-s, --by-seq                 by seq
-D, --dup-num-file string    file to save number and list of duplicated seqs
-d, --dup-seqs-file string   file to save duplicated seqs
-h, --help                   help for rmdup
-i, --ignore-case            ignore case


How does a -D file look like, you ask:

$zcat hairpin.fa.gz \ | seqkit rmdup -s -i -o clean.fa.gz -d duplicated.fa.gz -D duplicated.detail.txt$ cat duplicated.detail.txt   # here is not the entire list
3   hsa-mir-424, mml-mir-424, ppy-mir-424
3   hsa-mir-342, mml-mir-342, ppy-mir-342
2   ngi-mir-932, nlo-mir-932
2   ssc-mir-9784-1, ssc-mir-9784-2
`
• Hi @RamRS perfect thank you very much, that is really helpful
– M__
Jul 8 '20 at 7:26