How to write FASTA records using "Bio.SeqIO.write()"

Question

Update: Biopython document says that "Bio.SeqIO.FastaIO.FastaWriter" class is obsolete. Now my question becomes how to I use Bio.SeqIO.write() fucntion do 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()

Answer 1

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

Answer 2

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

Answer 3

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