Search code examples
pythonbiopythongenbank

How do I edit AND SAVE the sequence of a genbank file to a NEW genbank file using biopython?

I have a .gbk file that's wrong, and I have the list of corrections that follows the format of

"Address of Nuclotide: correct nucleotide"

1:T
2:C
4:A
63:A
324:G
etc...

I know how to open and parse exact original sequence with

list(SeqIO.parse(sys.argv[1], "genbank"))[0].seq

I just need to know how to replace it with my own nucleotide corrections. I've tried

seq_records[0].seq = "".join(dna_refseq)

Where the dna_refseq is a was just a list that constitutes the entire genome

I literally cannot find this specific action anywhere in the documentation or online, and intuitively, this is something that biopython should be capable of.

Solution

  • You are assigning a string where a Bio.Seq object is expected. For me, this works:

    from Bio import Seq
from Bio import SeqIO

my_entries = list(SeqIO.parse('my_file.gb', 'genbank'))
my_entry = my_entries[0]

# Make a new Seq object and assing to my_entry.seq. 'TTT' is just an example sequence
my_entry.seq = Seq.Seq('TTT', my_entry.seq.alphabet) 

# Write back to file
SeqIO.write(my_entries, 'my_updated_file.gb', 'genbank')

    If your Genbank file has only one entry, you might consider using SeqIO.read:

    my_entry = SeqIO.read('my_file.gb', 'genbank')

my_entry.seq = Seq.Seq('TTT', my_entry.seq.alphabet)
SeqIO.write(my_entry, 'my_updated_file.gb', 'genbank')

    Alternatively, you can directly convert the sequence into a mutable sequence and manipulate it directly:

    from Bio import SeqIO

my_entry = list(SeqIO.parse('my_file.gb', 'genbank'))[0]
my_entry.seq = my_entry.seq.tomutable()

my_entry.seq[0] = 'T'  # Remember that Genbank position 1 is 0 in seq
my_entry.seq[1] = 'C'
....
SeqIO.write(my_entry, 'my_updated_file.gb', 'genbank')