<p><strong>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.</strong></p>
<p><strong>Please let us know why this is not doing what you intend.</strong></p>
<p>I haven't checked it, but something like this would work:</p>
<pre><code>from Bio import SeqIO
import time
from collections import defaultdict

start = time.time() 

seq_counts = defaultdict(0)
records = []

 for record in SeqIO.parse(&quot;b4r2.fasta&quot;, &quot;fasta&quot;):  
   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, &quot;no_dupes_b4r2.fasta&quot;, &quot;fasta&quot;)
end = time.time()

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

print(f&quot;Run time is {(end- start)/60}&quot;) 
</code></pre>
<p>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.</p>
<p>The code that you have would probably not work as you would want the <code>records.append()</code> call to be only run in the if statement testing for non-membership.</p>
<p>There is absolutely no reason to have to write more than one script for this.</p>

CC BY-SA 4.0
2020-07-06T22:17:45.197Z
There is absolutely no need to write even this script when tools like [seqkit rmdup](https://bioinf.shenwei.me/seqkit/usage/#rmdup) already exist.
Ram RS
CC BY-SA 4.0
2020-07-06T23:50:07.967Z
Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.
M__
CC BY-SA 4.0
2020-07-07T01:52:56.900Z
Thank you! I've seen some other scripts that can do similar things but I wanted to make my own... :) Appreciate it a lot.
User
CC BY-SA 4.0
2020-07-07T09:00:26.617Z
i am indeed upset that you use `collections.defaultdict`, would prefer `collections.Counter` ;) e.g. `seq_counts = Counter(record.seq for record in SeqIO.parse(...`
Chris_Rands
CC BY-SA 4.0
2020-07-07T17:09:36.627Z
@RamRS I'll note that `seqkit rmdup` does not appear to fully satisfy the asker's question, which asked for a specific count of each sequence.
Maximilian Press
CC BY-SA 4.0
2020-07-07T17:23:25.703Z
@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...
Maximilian Press
CC BY-SA 4.0
2020-07-07T21:59:39.613Z
@MaximilianPress `-D` and `-d` do that. If OP wants to save count against each sequence (as opposed to against each sequence identifier), the file would be unreliable as it would need to be a FASTA with an annotation added in the header. With the `-D`, OP should be able to use `seqkit grep` and `seqkit duplicate` to get their custom file.
Ram RS
CC BY-SA 4.0
2020-07-07T22:35:42.470Z
@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.
Maximilian Press
CC BY-SA 4.0
2020-07-08T16:14:17.353Z
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.
Ram RS
CC BY-SA 4.0
2020-07-06T21:56:41.063Z
2020-07-07T22:37:00.690Z
2020-07-07T22:37:00.690Z
Maximilian Press
<p>Please do not reinvent the wheel - <code>seqkit rmdup</code> does what you want.</p>
<p>From the website:</p>
<hr />
<h1>rmdup</h1>
<h2>Usage</h2>
<pre><code>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
</code></pre>
<p>How does a <code>-D</code> file look like, you ask:</p>
<pre><code>$ 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
</code></pre>

CC BY-SA 4.0
2020-07-08T07:26:54.697Z
Hi @RamRS perfect thank you very much, that is really helpful
M__
CC BY-SA 4.0
2020-07-07T22:06:04.610Z
2020-07-07T22:06:04.610Z
Ram RS
<p>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...</p>
<pre><code>from Bio import SeqIO
import time

start = time.time() 

seen = []
records = []

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


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

print(f&quot;Run time is {(end- start)/60}&quot;) 
<span class="math-container">```</span>
</code></pre>

CC BY-SA 4.0
2020-07-08T07:29:02.433Z
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__
<fasta><python>
Python: How to write duplicate sequences removed from fasta file to new file
# Python: How to write duplicate sequences removed from fasta file to new file

(Asked by: username_0 on 2020-07-06)

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

**Tags:** fasta, python
> **Comments:**
> **username_1**: 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.  

---

## 2 answers

### Answer by username_2 on 2020-07-06. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_3**: There is absolutely no need to write even this script when tools like [seqkit rmdup](https://bioinf.shenwei.me/seqkit/usage/#rmdup) already exist.  
> **username_1**: Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.  
> **username_0**: Thank you! I've seen some other scripts that can do similar things but I wanted to make my own... :) Appreciate it a lot.  
> **username_4**: i am indeed upset that you use `collections.defaultdict`, would prefer `collections.Counter` ;) e.g. `seq_counts = Counter(record.seq for record in SeqIO.parse(...`  
> **username_2**: @RamRS I'll note that `seqkit rmdup` does not appear to fully satisfy the asker's question, which asked for a specific count of each sequence.  
> **username_2**: @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...  
> **username_3**: @MaximilianPress `-D` and `-d` do that. If OP wants to save count against each sequence (as opposed to against each sequence identifier), the file would be unreliable as it would need to be a FASTA with an annotation added in the header. With the `-D`, OP should be able to use `seqkit grep` and `seqkit duplicate` to get their custom file.  
> **username_2**: @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.  
> **username_3**: 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.  

---

### Answer by username_3 on 2020-07-07. Score: 4

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

> **Comments:**
> **username_1**: Hi @RamRS perfect thank you very much, that is really helpful  

---

