Exemplo de funcións recursivas

Donato publicaba sobre recursividade o mes pasado. Todos os exemplos de funcións de recursividade son sobre sucesións como a de Fibonacci, factorial, etc. e nunca se ven casos máis mundanos de recursividade.
Resulta que me tocou meter unha función de recursividade no noso servidor de modelado de GPCRs, e penso que como exemplo pode ser bastante útil.

O problema

No formulario de entrada aceptamos secuencias crudas en FASTA, ou ben identificadores Uniprot. Así que a función inicial era básicamente

def retrieve_url(url):
    import urllib2

    #url is something like
    #http://www.uniprot.org/uniprot/IDXXX.fasta

    try:
        my_response = urllib2.urlopen(url)
        return "".join(my_response.readlines())
    except urllib2.HTTPError:
        return False

Bastante simple: unha función anterior pásalle a esta o url depurado, despóis de comprobar si efectivamente é un ID válido de Uniprot. Esta función chama a Uniprot (que sigue o estándar REST), e obtén a secuencia en FASTA.
Bastantes secuencias empezaron a fallar, porque o servidor de Uniprot contestaba con un "300", a secuencia xa non se atopaba ahí, senón que a moveran a outra ID. Por outra parte, algo bastante habitual nesas bases de datos.

A solución

Fácilmente se deduce que hai que obter a nova ubicación do recurso movido, e ir outra vez a por él. Como xa temos a función que recupera FASTAS de Uniprot, ímola chamar desde sí mesma, e se o recurso volve dar un erro 300, vólvese a chamar ata que chega a un erro diferente ou ó FASTA apropiado

def retrieve_url(url):
    import urllib2

    try:
        my_response = urllib2.urlopen(url)
        return "".join(my_response.readlines())
    except urllib2.HTTPError, error:
        if error.code == 300: #Moved with multiple choices
            new_resource = error.info().getheader("Location")
            my_response = urllib2.urlopen(server + new_resource).readline()
            new_url = "http://www.uniprot.org/uniprot/uniprot/" + \
                      my_response.split("|")[1] + ".fasta"

            return retrieve_url(new_url)

    else:
        return False

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.