How to query a region in a fasta file using biopython?

Question

I have a fasta file with some reference genome. I would like to obtain the reference nucleotides as a string given the chromosome, start and end indexes.

I am looking for a function which would look like this in code:

from Bio import SeqIO
p = '/path/to/refernce.fa'
seqs = SeqIO.parse(p.open(), 'fasta')
string = seqs.query(id='chr7', start=10042, end=10252)

and string should be like : 'GGCTACGAACT...'

All I have found is how to iterate over seqs, and how to pull data from NCBI, which is not what I'm looking for. What is the right way to do this in biopython?

Answer 1

AFAIK, biopython does not currently have this functionality. For random lookups using an index (please see samtools faidx), you'll probably want either pysam or pyfaidx. Here's an example using the pysam.FastaFile class which allows you to quickly 'fetch' sequences in a region:

import pysam
ref = pysam.FastaFile('/path/to/reference.fa')
seq = ref.fetch('chr7', 10042, 10252)
print(seq)

Or using pyfaidx and the 'get_seq' method:

from pyfaidx import Fasta
ref = Fasta('/path/to/reference.fa')
seq = ref.get_seq('chr7', 10042, 10252)
print(seq)