<p>Just as a small variation to <a href="https://bioinformatics.stackexchange.com/a/2818/351">@Bioathlete's answer</a> in case you want to write the fasta using Biopython (e.g. to add names and description):</p>

<pre><code>from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)

# this creates a generator; if you need a physical list, replace the outer "(", ")" by "[" and "]", respectively
final_seq = (SeqRecord(Seq(seqi, IUPAC.protein), id="|".join(gi), name='', description='') for seqi, gi in dedup_records.items())

# write file
SeqIO.write(final_seq, 'test_output_final.fasta', 'fasta')
</code></pre>

<p>This will give the desired output:</p>

<pre><code>&gt;prot_id_1|prot_id_5|prot_id_7
MAWIGLDISKLFVENSRDAAA
&gt;prot_id_2|prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
&gt;prot_id_4
MEMVLSSANPSTTFADSYVV
&gt;prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
</code></pre>

CC BY-SA 3.0
2017-11-13T15:06:24.630Z
I think you mean it creates a generator not an iterator.
Bioathlete
CC BY-SA 3.0
2017-11-13T15:15:22.143Z
@Bioathlete: Of course, fixed, thanks for spotting it.
Cleb
CC BY-SA 3.0
2017-11-13T08:14:38.017Z
2017-11-13T15:14:57.760Z
2017-11-13T15:14:57.760Z
Cleb
<p>using <strong>datamash</strong> ans some awk</p>

<pre><code>$ cat input.fa | paste - -  |\
cut -c2- | sort -t $'\t' -k2,2 |\
datamash  collapse 1 -g 2  | tr "," "_" |\
awk '{printf("&gt;%s\n%s\n",$2,$1);}'

&gt;prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
&gt;prot_id_4
MEMVLSSANPSTTFADSYVV
&gt;prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
&gt;prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
</code></pre>

CC BY-SA 3.0
2017-12-05T13:44:01.903Z
2017-12-05T13:44:01.903Z
Pierre
<p>Here is my free program on Github  <strong>Sequence database curator</strong>
(<a href="https://github.com/Eslam-Samir-Ragab/Sequence-database-curator" rel="nofollow noreferrer">https://github.com/Eslam-Samir-Ragab/Sequence-database-curator</a>)</p>

<p>It is a very fast program and it can deal with:</p>

<ol>
<li>Nucleotide sequences</li>
<li>Protein sequences</li>
</ol>

<p>It can work under Operating systems:</p>

<ol>
<li>Windows</li>
<li>Mac</li>
<li>Linux</li>
</ol>

<p>It also works for:</p>

<ol>
<li>Fasta format</li>
<li>Fastq format</li>
</ol>

<p>Best Regards</p>

CC BY-SA 4.0
2019-04-18T13:14:10.157Z
Hello Eslam, thanks for your contribution. It would be nice if you could also include how to use your program to achieve what the original poster has asked for... Here we try to be as concrete as possible both when asking and answering :-).
Kamil S Jaron
CC BY-SA 4.0
2019-04-18T13:22:22.497Z
I agree with @KamilSJaron: Please show how to use your package to solve the issue from above (ideally with the same input). That will help others to get started. Thanks!
Cleb
CC BY-SA 4.0
2019-04-18T13:06:15.737Z
2019-04-18T13:06:15.737Z
Eslam Ibrahim
<p>If you have a very big fasta file(60 GB) why use SeqIO.parse and write the output file using biopython?</p>

<p>Why not use SimpleFastaParser instead (you get a tuple(name, seq)) and you dont need to wait so long as using the parser and using Seqrecords.</p>

<p>What do you think about it guys?</p>

<p>Thanks</p>

<p>Paulo</p>

<p>Kind of:</p>

<pre><code>from Bio.SeqIO.FastaIO import SimpleFastaParser
from collections import defaultdict

dedup_records = defaultdict(list)
for (name, seq) in SimpleFastaParser(filename):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[seq)].append(name)
with open(filename, 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write("&gt;{}\n".format('|'.join(ids)))
        output.write(seq + "\n")
</code></pre>

CC BY-SA 4.0
2020-01-19T23:17:45.590Z
2020-01-19T23:17:45.590Z
Paulo Sergio Schlogl
<p>This can also be done quite easily with awk:</p>

<p><code>awk '{if(NR%2){printf("%s\t",$1)}else{print}}' filename | sort -k2,2| awk '{if($2==SEQ){gsub("&gt;","",$1);ID=ID"|"$1} else{if(SEQ!=""){printf("%s\n%s\n", ID,SEQ);}SEQ=$2;ID=$1;} }END{printf("%s\n%s\n", ID,SEQ)}'</code></p>

<p>The first <code>awk</code> converts the fasta file to a tab separated file with format <code>ID\tSequence</code>, which is then sorted by sequence by <code>sort</code>. The last <code>awk</code> goes through the sorted file looking at the sequences: if the sequence in the current line is the same as that in the previous line, it appends its id to the ID variable. Otherwise (i.e. we have encountered a new sequence), it prints the previous concatenated ids and the corresponding sequence. </p>

CC BY-SA 3.0
2017-12-05T13:14:45.593Z
2017-12-05T13:14:45.593Z
Tom
<p>The trick would be to swap the key in the dictionary to be the sequence itself.</p>

<p>Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.</p>

<pre><code>from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write("&gt;{}\n".format('|'.join(ids)))
        output.write(seq + "\n")
</code></pre>

CC BY-SA 3.0
2017-11-10T15:52:27.497Z
Yes, that works fine and is more condensed than my solution (your `dedup_records` is the same as `seq_id_dict`). But true, all the additional filtering is not needed; Friday afternoon... :)
Cleb
CC BY-SA 3.0
2017-11-10T16:15:47.757Z
(I wait with acceptance for a while just to see if there is already some kind of built-in for that)
Cleb
CC BY-SA 3.0
2017-11-10T16:24:35.510Z
There is no built in function that I know of that does the id concatenation.  If it did exist it would likely be doing the same logic.
Bioathlete
CC BY-SA 3.0
2017-11-10T17:22:28.837Z
Seems like a good approach if fitting the file into memory is ok, but why not use biopython to write also? This asymmetry seems odd, plus that way you won’t dump all the sequence on a single line
Chris_Rands
CC BY-SA 3.0
2017-11-10T17:34:18.217Z
Mainly because I find the BioPython interface cumbersome to write out sequences in a fasta format when the data is not a Seq object already.  @cleb can use Biopython if he wants the main changes are in the use of the defaultdict with the swapped key.
Bioathlete
CC BY-SA 3.0
2017-11-10T18:41:08.840Z
I see, well still consider using record.description not record.id to ensure you get the whole headers. Also you could save memory by storing the sequences as a hash rather than the full sequences
Chris_Rands
CC BY-SA 3.0
2017-11-10T18:58:14.577Z
OP referred to the protein ID so thats what I used.  I agree that hashes could be useful to handle large input amounts.
Bioathlete
CC BY-SA 3.0
2017-11-10T19:59:44.263Z
@Chris_Rands: Mind adding the hash solution!? In my case I only have the IDs, no further description but if you could add a more detailed answer I am more than happy to upvote it.
Cleb
CC BY-SA 3.0
2017-11-27T20:43:16.103Z
@Cleb I meant to say, the reason I didn't add a hash solution is because it would be tricky to resolve hash collisions (although they'd be very rare)
Chris_Rands
CC BY-SA 3.0
2017-11-27T21:14:45.860Z
@Chris_Rands: No worries, the solutions posted are currently sufficient; would have just been a 'nice-to-see' type of thing...
Cleb
CC BY-SA 3.0
2017-11-10T15:45:20.237Z
2017-12-02T12:53:15.567Z
2017-12-02T12:53:15.567Z
Bioathlete
<p>I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in <a href="https://en.wikipedia.org/wiki/Biopython" rel="nofollow noreferrer">Biopython</a> for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).</p>

<p>Using the following example input</p>

<pre><code>&gt;prot_id_1
MAWIGLDISKLFVENSRDAAA
&gt;prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
&gt;prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
&gt;prot_id_4
MEMVLSSANPSTTFADSYVV
&gt;prot_id_5
MAWIGLDISKLFVENSRDAAA
&gt;prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
&gt;prot_id_7
MAWIGLDISKLFVENSRDAAA
</code></pre>

<p>I currently do it as follows:</p>

<pre><code>from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) &gt; 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')
</code></pre>

<p>which gives me the desired output:</p>

<pre><code>&gt;prot_id_4
MEMVLSSANPSTTFADSYVV
&gt;prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
&gt;prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
&gt;prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
</code></pre>

<p>Is there a smarter way of accomplishing this?</p>

<fasta><python><filtering>
How do I find identical sequences in a FASTA file?
# How do I find identical sequences in a FASTA file?

(Asked by: username_0 on 2017-11-10)

I want to create a database for a proteomics study. Therefore, the mapping from a given sequence to a protein ID has to be unique. I am wondering whether there is already a built-in function in Biopython for that, but I could not find any. The duplicated sequences should be merged into a single entry and there respective IDs should be just concatenated (see expected outcome below).
Using the following example input
>prot_id_1
MAWIGLDISKLFVENSRDAAA
>prot_id_2
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_5
MAWIGLDISKLFVENSRDAAA
>prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_7
MAWIGLDISKLFVENSRDAAA

I currently do it as follows:
from Bio import SeqIO
import numpy as np
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC

# Read entries into dictionary, gi: sequence
record_dict = SeqIO.to_dict(SeqIO.parse("test.fasta", "fasta"))

# Extract sequence as string
record_dict = {gi: str(record_dict[gi].seq) for gi in record_dict}

# Check for duplicated sequences
seq_id_dict = {}
for k, v in record_dict.items():
    seq_id_dict.setdefault(v, []).append(k)

# Filter sequences that appear more than once
filter_dict = {seqi: li for seqi, li in seq_id_dict.items() if len(li) > 1}

# Merge identifiers to a single one
merged_ids = {"_".join(sorted(li)): seqi for seqi, li in filter_dict.items()}

# Remove all entries that are associated with duplicate IDs
for gi in set(np.concatenate(filter_dict.values())):
    record_dict.pop(gi, None)

# Add the merged IDs again
record_dict.update(merged_ids)

# Create list of Seq objects
final_seq = [SeqRecord(Seq(seqi, IUPAC.protein), id=gi, name='', description='') for gi, seqi in record_dict.items()]

# Write file
SeqIO.write(final_seq, 'test_output.fasta', 'fasta')

which gives me the desired output:
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS

Is there a smarter way of accomplishing this?

**Tags:** fasta, python, filtering

---

## 6 answers

### Answer by username_1 on 2017-11-10. Score: 8 | Accepted answer

The trick would be to swap the key in the dictionary to be the sequence itself.
Also I would recommend using a different separator that "_" since that is what the current ids have so that you can easily separate the individual ids from the concatenated id. I used a pipe "|" in this example. Also I just manually wrote the FASTA output instead of using Biopython to do it.
from Bio import SeqIO
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)
with open("Output.fasta", 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")

> **Comments:**
> **username_0**: Yes, that works fine and is more condensed than my solution (your `dedup_records` is the same as `seq_id_dict`). But true, all the additional filtering is not needed; Friday afternoon... :)  
> **username_0**: (I wait with acceptance for a while just to see if there is already some kind of built-in for that)  
> **username_1**: There is no built in function that I know of that does the id concatenation.  If it did exist it would likely be doing the same logic.  
> **username_7**: Seems like a good approach if fitting the file into memory is ok, but why not use biopython to write also? This asymmetry seems odd, plus that way you won’t dump all the sequence on a single line  
> **username_1**: Mainly because I find the BioPython interface cumbersome to write out sequences in a fasta format when the data is not a Seq object already.  @cleb can use Biopython if he wants the main changes are in the use of the defaultdict with the swapped key.  
> **username_7**: I see, well still consider using record.description not record.id to ensure you get the whole headers. Also you could save memory by storing the sequences as a hash rather than the full sequences  
> **username_1**: OP referred to the protein ID so thats what I used.  I agree that hashes could be useful to handle large input amounts.  
> **username_0**: @Chris_Rands: Mind adding the hash solution!? In my case I only have the IDs, no further description but if you could add a more detailed answer I am more than happy to upvote it.  
> **username_7**: @Cleb I meant to say, the reason I didn't add a hash solution is because it would be tricky to resolve hash collisions (although they'd be very rare)  
> **username_0**: @Chris_Rands: No worries, the solutions posted are currently sufficient; would have just been a 'nice-to-see' type of thing...  

---

### Answer by username_0 on 2017-11-13. Score: 4

Just as a small variation to @Bioathlete's answer in case you want to write the fasta using Biopython (e.g. to add names and description):
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.Alphabet import IUPAC
from collections import defaultdict

dedup_records = defaultdict(list)
for record in SeqIO.parse("test.fasta", "fasta"):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[str(record.seq)].append(record.id)

# this creates a generator; if you need a physical list, replace the outer "(", ")" by "[" and "]", respectively
final_seq = (SeqRecord(Seq(seqi, IUPAC.protein), id="|".join(gi), name='', description='') for seqi, gi in dedup_records.items())

# write file
SeqIO.write(final_seq, 'test_output_final.fasta', 'fasta')

This will give the desired output:
>prot_id_1|prot_id_5|prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_2|prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH

> **Comments:**
> **username_1**: I think you mean it creates a generator not an iterator.  
> **username_0**: @Bioathlete: Of course, fixed, thanks for spotting it.  

---

### Answer by username_2 on 2017-12-05. Score: 2

using datamash ans some awk
$ cat input.fa | paste - -  |\
cut -c2- | sort -t $'\t' -k2,2 |\
datamash  collapse 1 -g 2  | tr "," "_" |\
awk '{printf(">%s\n%s\n",$2,$1);}'

>prot_id_1_prot_id_5_prot_id_7
MAWIGLDISKLFVENSRDAAA
>prot_id_4
MEMVLSSANPSTTFADSYVV
>prot_id_2_prot_id_6
MENVSRIWRYQQRRVLSRAFTHWYLMGLTKHNHPS
>prot_id_3
MRTRPSRDFHYIVGRITRESPEEKGIFYVFH


---

### Answer by username_3 on 2019-04-18. Score: 1

Here is my free program on Github **Sequence database curator**
(<https://github.com/Eslam-Samir-Ragab/Sequence-database-curator>)

It is a very fast program and it can deal with:

1. Nucleotide sequences
2. Protein sequences

It can work under Operating systems:

1. Windows
2. Mac
3. Linux

It also works for:

1. Fasta format
2. Fastq format

Best Regards

> **Comments:**
> **username_4**: Hello Eslam, thanks for your contribution. It would be nice if you could also include how to use your program to achieve what the original poster has asked for... Here we try to be as concrete as possible both when asking and answering :-).  
> **username_0**: I agree with @KamilSJaron: Please show how to use your package to solve the issue from above (ideally with the same input). That will help others to get started. Thanks!  

---

### Answer by username_5 on 2020-01-19. Score: 1

If you have a very big fasta file(60 GB) why use SeqIO.parse and write the output file using biopython?
Why not use SimpleFastaParser instead (you get a tuple(name, seq)) and you dont need to wait so long as using the parser and using Seqrecords.
What do you think about it guys?
Thanks
Paulo
Kind of:
from Bio.SeqIO.FastaIO import SimpleFastaParser
from collections import defaultdict

dedup_records = defaultdict(list)
for (name, seq) in SimpleFastaParser(filename):
    # Use the sequence as the key and then have a list of id's as the value
    dedup_records[seq)].append(name)
with open(filename, 'w') as output:
    for seq, ids in dedup_records.items():
        # Join the ids and write them out as the fasta
        output.write(">{}\n".format('|'.join(ids)))
        output.write(seq + "\n")


---

### Answer by username_6 on 2017-12-05. Score: 1

This can also be done quite easily with awk:

`awk '{if(NR%2){printf("%s\t",$1)}else{print}}' filename | sort -k2,2| awk '{if($2==SEQ){gsub(">","",$1);ID=ID"|"$1} else{if(SEQ!=""){printf("%s\n%s\n", ID,SEQ);}SEQ=$2;ID=$1;} }END{printf("%s\n%s\n", ID,SEQ)}'`

The first `awk` converts the fasta file to a tab separated file with format `ID\tSequence`, which is then sorted by sequence by `sort`. The last `awk` goes through the sorted file looking at the sequences: if the sequence in the current line is the same as that in the previous line, it appends its id to the ID variable. Otherwise (i.e. we have encountered a new sequence), it prints the previous concatenated ids and the corresponding sequence.


---

