How do I find identical sequences in a FASTA file?

Question

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?

Answer 1

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

Answer 2

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

Answer 3

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 4

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.

Answer 5

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

Answer 6

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