3
$\begingroup$

Update: Biopython document says that "Bio.SeqIO.FastaIO.FastaWriter" class is obsolete. Now my question becomes how to I use Bio.SeqIO.write() fucntion to achive the same thing I asked.

I use the following code to extract specific FASTA records from a file based on a list of FASTA IDs given in another file. The script works just fine. However, I want to able to write the results using "Bio.SeqIO.FastaIO.FastaWriter" class into a file instead of Python's print or write function. How can I do that in my last if statement?

from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'
#pprint.pprint(my_dictionary)

for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:
        sequence = str(seq_record.seq)
        fasta_record = textwrap.fill(sequence, width=60)
        print(f">{seq_record.id}\n{fasta_record}")

my_fasta_file.close()
my_id_file.close()
$\endgroup$
1
  • 1
    $\begingroup$ I was considering modifying the title. I am glad you edited it @gringer. Thank you. $\endgroup$
    – Supertech
    Commented Jan 22, 2023 at 18:57

3 Answers 3

3
$\begingroup$

I figured it out:

from Bio import SeqIO
#import pprint
import textwrap

my_id_file     = open('sample_files2/fasta_id_records.txt','r')
my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
result_file    = open('sample_files2/results.fasta', "w")

my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:      
        #use this line to write standard fasta record using SeqIO.write
        SeqIO.write(seq_record, result_file, "fasta")

my_fasta_file.close()
my_id_file.close()
result_file.close()
$\endgroup$
3
  • $\begingroup$ I'm not 100% sure you can because that will be taken care within SeqIO.write(). The obvious seq_record.seq = textwrap.fill(str(seq_record.seq) , width=60) is unlikely to work. $\endgroup$
    – M__
    Commented Jan 21, 2023 at 22:13
  • $\begingroup$ The only way to do it is within a straight write to disk in you first code ... but (obviously) write rather than print. $\endgroup$
    – M__
    Commented Jan 21, 2023 at 22:14
  • 1
    $\begingroup$ it wraps to 60 by default in Biopython 1.8 at least, if saves as 'fasta-2line' uses a single line for the seq $\endgroup$
    – pippo1980
    Commented Jan 22, 2023 at 12:00
2
$\begingroup$

my attempt at changing the way Biopython formats fasta file sequence width, read into code to figure out standard behaviour:

from Bio import SeqIO, __version__
#import pprint
# import textwrap

print('Biopython version : ' , __version__)

# my_id_file     = open('sample_files2/fasta_id_records.txt','r')
# my_fasta_file  = open('sample_files2/fasta_file.fasta','r')
#result_file    = open('sample_files2/results.fasta', "w")
### CHANGED TO:
    
my_id_file     = open('fasta_id_records.txt','r')
my_fasta_file  = open('fasta_file.fasta','r')
result_file    = open('results.fasta', "w")


my_dictionary = {} # fasta IDs are keys, value can be anything.
for line in my_id_file:
    my_dictionary[line[:-1]] = 'value'

#expected lenght of sequences in my fasta
length = 10

#pprint.pprint(my_dictionary)
for seq_record in SeqIO.parse(my_fasta_file, "fasta"):
    if seq_record.id in my_dictionary:      
        
        
        ########################################################################################
        #use this line to write standard fasta record using SeqIO.write
        
        # SeqIO.write(seq_record, result_file, "fasta-2line") # writes fasta in 2 line format
        
        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with width of 60 as standard biopython format
        
        ### NOT POSSIBLE !!!!!!!!!!!!
        # seq_record.seq = textwrap.fill(seq_record.seq , width=10) ### --> AttributeError: 'Seq' object has no attribute 'expandtabs'
        
        # print(seq_record.seq, type(seq_record.seq), dir(seq_record.seq))
        
        sequence = str(seq_record.seq)
        
        result_file.write('>'+seq_record.id+'\n')
        
        while len(sequence) > 0:
            result_file.write(sequence[:length]+'\n') 
            sequence = sequence[length:]
        result_file.write('\n')
        
        
        ### if you try to add EOL inside seq_record.seq of course you get 
        ### TypeError: SeqRecord (id=seq1) has an invalid sequence.
        # SeqIO.write(seq_record, result_file, "fasta") # writes fasta with sequence width as 60   !!!! TypeError: SeqRecord (id=seq1) has an invalid sequence.

my_fasta_file.close()
my_id_file.close()
result_file.close()

input files :

fasta_id_records.txt:

seq1
seq2
seq3

fasta_file.fasta:

>seq0
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
>seq1
TTTT
TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
TTTTTTTTTTTTTTT
>seq3
CCCCCCCCCCCCCCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG

output file:

v :

>seq1
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTTTTTTTTT
TTT

>seq3
CCCCCCCCCC
CCCCCCCGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCC
CCCCCCCCCG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
GGGGGGGGGG
G


$\endgroup$
2
$\begingroup$

You could instead pass in a generator that yields SeqRecords. This avoids storing the entire list of sequences in memory and is usually the way to go assuming that you need to handle large files containing thousands of records. This example use a generator expression, which is a bit like a list comprehension, but instead of creating a list, it creates a generator:

import argparse
import contextlib
import pathlib
import sys

from Bio import SeqIO
def get_argument_parser():

    parser = argparse.ArgumentParser()

    parser.add_argument(
        "fasta_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Extract sequence records from FILE",
    )
    parser.add_argument(
        "ids_file",
        type=pathlib.Path,
        metavar='FILE',
        help="Sequence record IDs/names to extract, one per line",
    )
    parser.add_argument(
        "-o",
        "--output",
        type=pathlib.Path,
        metavar='FILE',
        default="-",
        help="Write extracted sequences to FILE (default: stdout)",
    )
    parser.add_argument(
        "-f",
        "--force",
        action='store_true',
        help="Overwrite the output file if it exists",
    )

    return parser
def main():

    parser = get_argument_parser()
    args = parser.parse_args()

    with contextlib.ExitStack() as stack:

        ids_file = stack.enter_context(open(args.ids_file))
        fasta_file = stack.enter_context(open(args.fasta_file))

        if args.output != '-':
            out = stack.enter_context(open(args.output, 'w' if args.force else 'x'))
        else:
            out = sys.stdout

        record_ids = {line.rstrip('\n') for line in ids_file}

        fasta = SeqIO.parse(fasta_file, "fasta")
        records = (rec for rec in fasta if rec.id in record_ids)

        SeqIO.write(records, out, "fasta")

if __name__ == "__main__":
    main()
$\endgroup$
5
  • 1
    $\begingroup$ cool never eard about ExitStack till now docs.python.org/3/library/contextlib.html#contextlib.ExitStack $\endgroup$
    – pippo1980
    Commented Jan 22, 2023 at 12:59
  • 1
    $\begingroup$ A better description about what it does is here rath.org/on-the-beauty-of-pythons-exitstack.html @pippo1980 $\endgroup$
    – M__
    Commented Jan 22, 2023 at 13:03
  • $\begingroup$ I think I need to ask another question on this code. It's certainly professional code. What I'm wondering is how the if-condition tuple comprehension is a generator because thats the "workhorse" of the code. $\endgroup$
    – M__
    Commented Jan 22, 2023 at 13:11
  • 1
    $\begingroup$ @M__ The parentheses just make it a generator expression. Brackets would have made it a list comprehension: peps.python.org/pep-0289 $\endgroup$
    – Steve
    Commented Jan 22, 2023 at 13:21
  • $\begingroup$ @Steve, gave me a code to study. $\endgroup$
    – Supertech
    Commented Jan 22, 2023 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.