Updating the GFF3 + Fasta to GeneBank code

Question

I'm trying to convert gff3 and fasta into a gbk file for usage in Mauve. I've found a solution but the code is outdated:

"""Convert a GFF and associated FASTA file into GenBank format.

Usage:
gff_to_genbank.py <GFF annotation file> <FASTA sequence file>
"""
import sys
import os

from Bio import SeqIO
from Bio.Alphabet import generic_dna
from BCBio import GFF

def main(gff_file, fasta_file):
    out_file = "%s.gb" % os.path.splitext(gff_file)[0]
    fasta_input = SeqIO.to_dict(SeqIO.parse(fasta_file, "fasta", generic_dna))
    gff_iter = GFF.parse(gff_file, fasta_input)
    SeqIO.write(gff_iter, out_file, "genbank")

if __name__ == "__main__":
    main(*sys.argv[1:])

But the Bio.Alphabet module is removed now. I tried to read the documentation on these changes and how to update the code (https://biopython.org/wiki/Alphabet) but I'm still confused about how to update the code. If I remove the Bio.Alphabet and the generic_dna, the code will return ValueError

ValueError: missing molecule_type in annotations

I know that I have to type the molecule_type = "DNA" but I'm confused as to where should I put it because the code directly inputs the fasta_file into the SeqIO.parse() function.

Anyone knows how to update code related to SeqIO.parse() with the removal of Bio.Alphabet? Or anyone knows how to convert GFF3 + Fasta into gbk files?

Note: I can't use the seqret because my files is just too big for it, same problem as this question.

Thanks in advance!

Answer 1

You can set the molecule type on the existing annotations attribute of the record objects just before they get to SeqIO.write.

There might be an even leaner way to do it, but if open the output file (so that you can do multiple writes to one file handle) and then iterate over each record from the GFF iterator, you can do it like this:

with open(out_file, "wt") as f_out:
    for record in GFF.parse(gff_file, fasta_input):
        record.annotations["molecule_type"] = "DNA"
        SeqIO.write(record, f_out, "genbank")

(The documentation is vague on that point but if you check you'll see that record.annotations is a dictionary that you can modify to set the molecule type.)