Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-06T22:17:45.197Z
Answers.Text: There is absolutely no need to write even this script when tools like [seqkit rmdup](https://bioinf.shenwei.me/seqkit/usage/#rmdup) already exist.
Answers.UserDisplayName: Ram RS
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-06T23:50:07.967Z
Answers.Text: Hi @RamRS any chance you could put this as an answer? It would be helpful for my own personal records if nobody else.
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T01:52:56.900Z
Answers.Text: Thank you! I've seen some other scripts that can do similar things but I wanted to make my own... :) Appreciate it a lot.
Answers.UserDisplayName: User
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T09:00:26.617Z
Answers.Text: 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(...`
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T17:09:36.627Z
Answers.Text: @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.
Answers.UserDisplayName: Maximilian Press
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T17:23:25.703Z
Answers.Text: @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...
Answers.UserDisplayName: Maximilian Press
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T21:59:39.613Z
Answers.Text: @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.
Answers.UserDisplayName: Ram RS
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T22:35:42.470Z
Answers.Text: @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.
Answers.UserDisplayName: Maximilian Press
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-08T16:14:17.353Z
Answers.Text: 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.
Answers.UserDisplayName: Ram RS
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-06T21:56:41.063Z
Answers.LastActivityDate: 2020-07-07T22:37:00.690Z
Answers.LastEditDate: 2020-07-07T22:37:00.690Z
Answers.UserDisplayName: Maximilian Press
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-08T07:26:54.697Z
Answers.Text: Hi @RamRS perfect thank you very much, that is really helpful
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-07-07T22:06:04.610Z
Answers.LastActivityDate: 2020-07-07T22:06:04.610Z
Answers.UserDisplayName: Ram RS
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-07-08T07:29:02.433Z
Comments.Text: 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.
Comments.UserDisplayName: M__
Tags: <fasta><python>
Title: Python: How to write duplicate sequences removed from fasta file to new file
ThreadText: # 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 @username_0 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 @username_3 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 @username_3 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**: @username_3 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**: @username_4 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**: @username_2 `-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**: @username_3 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 @username_3 perfect thank you very much, that is really helpful  
  
  ---
  
  