Fetching sequences from NCBI

Here's the problem I want to solve:

I need to collect sequences from NCBI for a small database of 16S rRNA gene sequences from bacteria of interest. Many of these are from clonal isolates so the FASTA-formatted sequences have title lines like

>gi|29290057|gb|AY207063.1| Abiotrophia sp. oral clone P4PA_155 P1 16S ribosomal RNA gene, partial sequence

I can parse that title, or I also have the Genbank identifier in a database entry like this:

AY207063.1 Abiotrophia_clone_P4PA_155

I grab the corresponding sequence using PyCogent (as discussed in the post here):

from cogent.db.ncbi import EFetch ef = EFetch(id='AY207063') print ef.read()

> python fetch.py >gi|29290057|gb|AY207063.1| Abiotrophia sp. oral clone P4PA_155 P1 16S ribosomal RNA gene, partial sequence GACGAACGCTGGCGGTGTGCCTAATACATGCAAGTCGAACGAACCGCGACTAGGTGCTTGCACTTGGTCA ..

However, I also want to get some rRNA sequences from genome projects. In that case the FASTA-formated RefSeq record is huge, and I just want a small part of it.

From the prokaryote genome projects page at NCBI, I click on the RefSeq link for E. coli MG1655, and on the RefSeq page I click on Structural RNAs

I find what I'm looking for halfway down the page:

16S ribosomal RNA of rrnB operon

I copy the link for "DNA region in flatfile format" (the green diamond)

http://www.ncbi.nlm.nih.gov/nuccore/NC_000913?from=4164682&to=4166223

I tried feeding the range information (from and to) to EFetch but got nowhere. And I'm not sure if this is the best way, but here is what I did. I browsed through

class EFetch(UrlGetter) in /Library/Python/2.6/site-packages/cogent/db/ncbi.py class UrlGetter(object) in /Library/Python/2.6/site-packages/cogent/db/util.py

I just added a couple of keys to the dict which is used to construct the url. In cogent.utils.ncbi after line 133 I inserted:

'from', 'to',

Now this works:

from cogent.db.ncbi import EFetch D = {'id':'NC_000913','from':'4164682','to':'4166223'} ef = EFetch(**D) print str(ef) + '\n' print ef.read()

> python fetch.py http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?retmax=100&from=4164682&retmode=text&tool=PyCogent&db=nucleotide&email=Michael.Robeson%40colorado.edu&to=4166223&rettype=fasta&retstart=0&id=NC_000913 >gi|49175990:4164682-4166223 Escherichia coli str. K-12 substr. MG1655 chromosome, complete genome AAATTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAACGGT ..

Genbank is smart enough to figure out that the id we sent is not a real gi but refers to a record in the RefSeq database. Rather than using from and to, you'd think it would be OK to use something like:

NC_000913:4164682-4166223

but I couldn't get it past EFetch, even with URL encodings. Guess I should run this by the PyCogent crew for issues. Maybe they'll want to put it in the toolkit.