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

Question

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}") 
```

Answer 1

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.

Answer 2

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