9
$\begingroup$

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?

$\endgroup$
8
$\begingroup$

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")
$\endgroup$
10
  • $\begingroup$ 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... :) $\endgroup$ – Cleb Nov 10 '17 at 15:52
  • $\begingroup$ (I wait with acceptance for a while just to see if there is already some kind of built-in for that) $\endgroup$ – Cleb Nov 10 '17 at 16:15
  • $\begingroup$ 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. $\endgroup$ – Bioathlete Nov 10 '17 at 16:24
  • $\begingroup$ 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 $\endgroup$ – Chris_Rands Nov 10 '17 at 17:22
  • $\begingroup$ 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. $\endgroup$ – Bioathlete Nov 10 '17 at 17:34
4
$\begingroup$

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
$\endgroup$
2
  • $\begingroup$ I think you mean it creates a generator not an iterator. $\endgroup$ – Bioathlete Nov 13 '17 at 15:06
  • $\begingroup$ @Bioathlete: Of course, fixed, thanks for spotting it. $\endgroup$ – Cleb Nov 13 '17 at 15:15
2
$\begingroup$

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
$\endgroup$
1
$\begingroup$

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.

$\endgroup$
1
$\begingroup$

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

$\endgroup$
2
  • 1
    $\begingroup$ 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 :-). $\endgroup$ – Kamil S Jaron Apr 18 '19 at 13:14
  • $\begingroup$ 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! $\endgroup$ – Cleb Apr 18 '19 at 13:22
0
$\begingroup$

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")
$\endgroup$

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.