lunes, 9 de mayo de 2011

A API do NCBI

A veces teremos que consultar de forma masiva o NCBI. O normal é baixarse un xenoma que está en miles de trozos, ou que nun artigo nos digan que subiron as secuencias AE000001-AE001000.

Imos facer un pequeniño script para facer este tipo de consultas fácilmente. Primeiro conseguimos os id de referencia. Por exemplo, si queremos baixar as secuencias de Aedes aegypti, miramos no artigo http://www.ncbi.nlm.nih.gov/pubmed/17510324, que nos dí que na entrada AAGE00000000 está o xenoma. Consultando nesa entrada vemos que o WGS está en AAGE02000001-AAGE02036206.

import time
import urllib, urllib2

root = "AAGE020"
end = 36206

efetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"

output = open("AAGE02000000.fas", "w")

for i in ["{0:05}".format(x) for x in range(1, end)]:
params = {"db": "nucleotide",
"id": root + i,
"rettype": "fasta"}
GET_params = urllib.urlencode(params)
page = urllib2.urlopen("{url}?{GET}".format(url=efetch, GET=GET_params))

output.write(page.read())

time.wait(1)

output.close()


A liña ["{0:05}".format(x) for x in range(1, end)] sólo prepara o número da entrada con "0" á esquerda ata completar unha lonxitude de 5 números ("00001", "00134", etc).

Datos con ids non consecutivas



Xogando un pouco coas APIs de NCBI, cos valores da base de datos e coas ids, podemos pedir calquer dato do NCBI. Si queremos datos non consecutivos ou que teñamos nunha lista, por exemplo unha lista de PMIDs, modificamos así o script:

import time
import urllib, urllib2

my_querys = ["19226459", "18269752", "17725574", "17507679", "17203060",
"16174338", "15797612", "15545653", "15219154", "14738738",
"12716978", "12361573", "12032249", "11732090", "11108594",
"10893864"]

efetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"

output = open("articles.txt", "w")

for i in my_querys:
params = {"db": "pubmed",
"id": i,
"rettype": "abstract",
"retmode": "html"}
GET_params = urllib.urlencode(params)
page = urllib2.urlopen("{url}?{GET}".format(url=efetch, GET=GET_params))

output.write(page.read())

time.wait(1)

output.close()


A lista que cargamos en my_querys" podémola obter de calquera fonte: si temos unha base de datos de bibliografía, é moi fácil exportar os PMID ou sacalos parseando un volcado desa base de datos a un formato de texto plano.

Be polite with NCBI



Si temos pensado baixar un número ENORME de datos, é mellor facer primeiro unha búsqueda, e cos datos da búsqueda pedir as secuencias. O NCBI prefire que lle consultemos os datos desta forma, de tal modo que pode banearnos a IP si llos pedimos masivamente como víamos no primeiro snippet.

import time
import urllib, urllib2
from xml.dom import minidom

term = "AAGE02000000"

esearch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/esearch.fcgi"
efetch = "http://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi"

search = urllib2.urlopen(esearch+ "?" + \
urllib.urlencode({"db": "nuccore",
"term": term,
"usehistory": "yes"}))

xml = minidom.parse(search)
webenv = xml.getElementsByTagName("WebEnv")[0].firstChild.data
q_key = xml.getElementsByTagName("QueryKey")[0].firstChild.data
count = xml.getElementsByTagName("Count")[0].firstChild.data

output = open("AAGE02000000.fas", "w")

ret_max = 500
for ret_start in range(0, int(count), ret_max):

params = {"db": "nuccore",
"WebEnv": webenv,
"query_key": q_key,
"retstart": ret_start,
"retmax": ret_max,
"rettype": "fasta"}
GET_params = urllib.urlencode(params)
page = urllib2.urlopen("{url}?{GET}".format(url=efetch, GET=GET_params))

output.write(page.read())

output.close()


No parámetro "term" podemos meter unha liña de texto coa que fagamos calquera búsqueda, coa sintaxis propia do NCBI. Máis en NCBI - Entrez, pero en perl.

No hay comentarios: