Add file name to FASTA file header using SeqIO

Question

I have several fasta files that represent the proteins from different orthogroups. However, the orthogroup that each of the proteins belongs to is not indicated within the files. I am leveraging BioPython and I want to add this information to the header of each sequence in the fasta files.

import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
    for file in os.listdir(dirname):
        seq = SeqIO.parse(file, 'fasta')
        ortog = Path(file).stem
        for record in seq:
            record.description = ortog
        SeqIO.write(record, output, 'fasta')

However, the output always appears with the last sequence duplicated. In the case of three fasta files, the output appears as:

>BBC_8_672 OG0015396-PROT
MKFSSVTAITLATVATVATAKKGEHDFTTTLTLSSDGSLTTTTSTHTTHKYGKFNKTSKS
KTPNHTGTHKYGKFNKTSKSKTP
>CAT_6_3128 OG0020551-PROT
MEKIKEKLNSLKLESESWQEKYEELREQLKELEQSNTEKENEIKSLSAKNEQLDSEVEKL
ESQLSDTKQLAEDSNNLRSNNENYTKKNQDLEQQLEDSEAKLKEAMDKLKEADLNSEQMG
RRIVALEEERDEWEKKCEEFQSKYEEAQKELDEIANSLENL
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP
>ZSA_0_3121 OG0025647-PROT
LLLQYQLHPALLPKALLLQYQLHQALLLKALLLQ*PALLLKALLLQYQLHQALLLKALLL
QYQLHPALLLKALLLQYQLHP

Could anyone tell me why this is happening?

Answer 1

This is because you are also parsing your output file. Remember, the first rule of debugging is to print all the things! Most bugs are because some variable doesn't actually have the value you think it has, and printing them is the best way of finding the problem. I reproduced what (I am guessing) is your setup here:

$ ls
OG0015396-PROT.fa  OG0020551-PROT.fa  OG0025647-PROT.fa

I then ran your script with an extra print statement:

#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              print("FF:",file)
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

Which printed:

$ foo.py
FF: OG0015396-PROT.fa
FF: OG0020551-PROT.fa
FF: output.fasta
FF: OG0025647-PROT.fa

As you can see, since the output.fasta file is created by the open() call, by the time you use listdir, it's already there and included. This results in duplicated sequences (not always the last one). I'm sure there are more elegant solutions, but a simple workaround is to skip the output file name:

#!/usr/bin/env python3
import os
from pathlib import Path
from Bio import SeqIO

with open('output.fasta', 'w') as output:
       for file in os.listdir("."):
              if file == "output.fasta":
                     continue
              seq = SeqIO.parse(file, 'fasta')
              ortog = Path(file).stem
              for record in seq:
                     record.description = ortog
              SeqIO.write(record, output, 'fasta')

Of course, this whole thing is a bit overkill for such a simple task and you could just use sed and the shell instead:

for file in *fa; do sed "1s|$| ${file//.fa}|" "$file"; done > output.fasta