I have a genbank file and a list of addresses. I'd like to pull the exact nucleotide that corresponds to that address in biopython. How would I do that?

Question

[Tom]

I can't seem to figure out how to do this specific task. So I have a .gbk assembly file, and I have a simple excel sheet full of about 1000 numbers. The numbers correspond to SNP sites I can find in the genbank, but what I really want in the end is a list of what those nucleotides are in the assembly. What commands do I use to view said nucleotides given a list of addresses?

If biopython does not have this option what does?

[Tanvir Ahamed]

Could be informative :

1. Bioconductor : Load local genbank file
2. BioPython : Converting GenBank files to FASTA format with Biopython

Answer 1

[Peter]

Something like this:

from Bio import SeqIO
my_snp_list = [100, 1234]  # Python counting
record = SeqIO.read("single_contig.gbk")  # Assumes single contig!
for snp in my_snp_list:
    print("Position %i is nucleotide %s" % (snp, record.seq[snp]))

This is using Python's slice notation to pull out a single base from the sequence. Note you may need to convert your SNP coordinates to Python style zero-based counting by subtracting one.