0
$\begingroup$

I have a multi fasta file similar to this (relatively new here so uncertain of best way to present this; I have gone for an output and the code i used to make it - belt and braces...):

# create a test dataset
with open('test_dat.txt', 'w') as test_dat:
    l1 = [">AAA a aa", ">BBB b bb", ">CCC c cc"]
    l2 = ["ATGATG", "GACGAC", "TTGTTG"]

    for i in range(len(l1)):
        test_dat.write(l1[i] + "\n")
        test_dat.write(l2[i] + "\n")

# which outputs something like this
>AAA a aa
ATGATG
>BBB b bb
GACGAC
>CCC c cc
TTGTTG

I want to modify the fasta headers in my fasta file by adding some text to the beginning and, ideally, stripping out spaces (as i was advised that white space in the headings could cause problems for me down the line). My plan was to do this using BioPython and SeqIO.

I started like this:

# my attempt at altering the fasta headers using SeqIO
from Bio import SeqIO

# open my test data and create an output file
with open('test_dat.txt', 'r') as inputs, open('test_dat_out.txt', 'w') as outputs:

    # create dictionary of the fasta sequences using SeqIO
    my_dict = SeqIO.to_dict(SeqIO.parse(inputs, "fasta"))

    # renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header and remove white space
    for v in my_dict.values(): # using the values() to loop through...
        v.description = "TEST~~~TEST~~~" + v.description.replace(" ", "_")

    # write result to file
    SeqIO.write(my_dict.values(), outputs, 'fasta')

which gave me a result with some duplication in the fasta header:

#result i want
>TEST~~~TEST~~~AAA_a_aa
ATGATG

#result i get
>AAA TEST~~~TEST~~~AAA_a_aa
ATGATG

If i modify both the .id and .description in exactly the same way, as suggested here, but without removing white space:

# renaming my fasta headers - adding "TEST~~~TEST~~~" to front of header
    for v in my_dict.values(): # using the values() to loop through...
        v.id = "TEST~~~TEST~~~" + v.id
        v.description = "TEST~~~TEST~~~" + v.description

I get almost what i want, but the header still has spaces:

>TEST~~~TEST~~~AAA a aa
ATGATG

I have tried to make the .id and .description match by adding a "_" to the end of .id, or by adding .replace(" ", "_") to both .id and .description (despite there being no spaces in .id), but I just cant seem to get exactly the output i want.

I would really appreciate any pointers on this.

$\endgroup$
1
$\begingroup$

You're very close, but unfortunately for legacy reasons this aspect of the Biopython parser is not very intuitive. The key point is the .id is anything up until the first white space character. One day I hope the .description and .id attributes are replaced by a single attribute like .header. Anyway, here is a solution for your case:

from Bio import SeqIO

to_add = "TEST~~~TEST~~~"

with open("out.fa", "w") as outputs:
    for r in SeqIO.parse("in.fa", "fasta"):
        r.id = (to_add + r.description).replace(" ", "_")
        r.description = r.id
        SeqIO.write(r, outputs, "fasta")

Output on your sample input:

>TEST~~~TEST~~~AAA_a_aa
ATGATG
>TEST~~~TEST~~~BBB_b_bb
GACGAC
>TEST~~~TEST~~~CCC_c_cc
TTGTTG
$\endgroup$
1
  • $\begingroup$ Much appreciated @Chris_Rands. So my problem was that as soon as i removed all white space characters from the .description, the description - in its entirety - is considered to be the .id. $\endgroup$
    – Mark
    May 22 '20 at 8:06

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.