0
$\begingroup$

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}") 
```
$\endgroup$
1
  • $\begingroup$ 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. $\endgroup$ – M__ Jul 8 '20 at 7:29
3
$\begingroup$

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.

$\endgroup$
9
  • 1
    $\begingroup$ There is absolutely no need to write even this script when tools like seqkit rmdup already exist. $\endgroup$ – Ram RS Jul 6 '20 at 22:17
  • $\begingroup$ Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else. $\endgroup$ – M__ Jul 6 '20 at 23:50
  • 1
    $\begingroup$ @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... $\endgroup$ – Maximilian Press Jul 7 '20 at 17:23
  • 1
    $\begingroup$ @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. $\endgroup$ – Maximilian Press Jul 7 '20 at 22:35
  • 1
    $\begingroup$ 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. $\endgroup$ – Ram RS Jul 8 '20 at 16:14
4
$\begingroup$

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
$\endgroup$
1
  • 1
    $\begingroup$ Hi @RamRS perfect thank you very much, that is really helpful $\endgroup$ – M__ Jul 8 '20 at 7:26

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.