Xa escribín antes sobre os PDBs, hoxe toca algo máis metido co formato PDB. Un PDB trae unha lista de átomos con un montón de datos sobre a estructura en xeral e cada átomo en particular. Para esta entrada só nos interesan as liñas de átomos, que son as que comezan con "ATOM", "HETATM" ou "ANISOU".
Parseando lixeiramente o archivo
Buscamos estas liñas:
ATOM 10 1HG1 ILE A 1 1.867 30.218 56.564 1.00 0.00
Que según a especificación do PDB,
31 - 38 Son as coordenadas ortogonales de X en Angstroms.
39 - 46 Son as coordenadas ortogonales de Y en Angstroms.
47 - 54 Son as coordenadas ortogonales de Z en Angstroms. Así que ahí vamos a cazar esas liñas e gardalas nun array.
import numpy as np
atom = ["ATOM", "HETATM", "ANISOU"]
coords = []
n_atoms = 0
for line in pdb:
if line.split()[0] in atom:
coords.append(
[float(s) for s in [line[31:38], line[39:46], line[47:54]]])
#Reformateamos as posicions como unha matriz numpy
positions = np.reshape(coords, (n_atoms, 3))
Agora simplemente temos que tomar os valores máximos e mínimos en cada eixe, e teremos as medidas do PDB:
box_size = []
for col in range(3):
box_size.append(abs(positions[:,col].max() - positions[:,col].min()))
print box_size
box_size vale algo así como [220.609, 121.049, 178.267]. As medidas en Angstroms do PDB.
miércoles, 27 de junio de 2012
Regresión polinomial (extremely easy)
Usando a librería numpy podemos axustar unha curva de predicción a datos empíricos en catro liñas. Supoñamos que temos un programa que analiza texto. Funciona moi ben nos textos de proba, pero intuimos que textos de máis de 40 millóns de caracteres pode bloquear ou eternizar o script. ¿Cómo facer unha estima a priori de cánto tardará o script? Facemos unhas probas con textos controlados:
Millóns Tempo
de caracteres | en segundos
---------------+--------------------
1 | 0.65
2 | 1.76
3 | 3.34
4 | 5.4
4.5 | 6.8
5 | 7.9
7.5 | 16.15
10 | 27.2
20 | 101.5
---------------+--------------------
(Estes datos foron medidos realmente). Rápidamente podemos decidir un punto de corte, por exemplo textos de máis de 10 millóns de caracteres non se van procesar porque levan moito tempo. Pero mellor predecimos cánto vai tardar o script axustando un polinomio ós datos que temos. Ploteando os puntos anteriores, podemos intuir que necesitamos "máis" que unha recta, pero "chega" con unha parábola:
tempo = a + bx + cx²
import numpy as np
from numpy.polynomial.polynomial import polyfit
#Millions of chars
n_char = np.array([1, 2, 3, 4, 4.5, 5, 7.5, 10, 20])
#Seconds meassured
times = np.array([0.65, 1.76, 3.34, 5.4,
6.8, 7.9, 16.15, 27.2, 101.5])
#polyfit(X, Y, grado polinomio)
print polyfit(n_char, times, 2)
Salida:
[0.11502188 0.37922652 0.23444005]
Then, solution. Contamos o número de palabras, metémolo en x, e calculamos:
tempo estimado = 0.115 + 0.379x + 0.234x²
Para un texto de 15 millóns de palabras, estimamos 0.115 + 0.379 * 15 + 0.234 * 15², ~58 seg.
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
Temos páxinas que se leen con dificultade en dispositivos móviles, como iPhone e Android. Buscando por ahí adiante atopéi varios enfoques, e o máis extendido é analizar a petición do cliente (HTTP_USER_AGENT sobre todo) para ver si pertence a un dispositivo móvil. Unha vez determinado, aplicar a solución requerida
A base
Partín da base de Django mobile utils, que implementa unha solución robusta. Cando chega a petición (request), antes de chegar ó views.py pasa por un middleware que intenta averiguar si o USER_AGENT é móbil. A función está no __init__.py, e garda o request nunha variable local. Aquí é donde me falla, xa que non consigo recuperar as variables locais. Despóis de pasar polo middleware, o request orixinal está marcado como .is_mobile = True. Hai un código para cargar os "templates" alternativos, dentro de loaders.py; este loader mira dentro da variable local de request si ten o valor .is_mobile == True, porque non se lle pode pasar como parámetro o propio request. Así que si falla o request local, falla todo.
Por sorte, o autor tamén incluiu un context_processor, que simplemente se encarga de pasarlle un diccionario de variables ós templates. Por exemplo, si activamos o context_procesor de auth, en tódolos templates teremos dispoñibles un número de variables como user.
A solución
Utilizando o código anterior, podemos simplificalo para utilizar só o __init__.py e o context_processor.py. Os context_processors deben manterse simples según o manual oficial, así que si podemos pasar só unha variable, mellor que dúas. Dende o context_processor.py orixinal:
Vemos que o autor orixinal sigue tirando de request acumulado en local, pero nós imos testear directamente o request. Modificamos tamén o código do __init__.py. Das varias liñas de tipo:
request.mobile = True
return request
Cambiamos a:
return True
E a última liña que pon return request cambiámola a return False. Ou adornamos con variables, si nos gusta máis.
Si estamos pola simplificación extrema, todas as globáis deste archivo pódense eliminar, así como a configuración no settings.py, no que chega con engadir:
Agora temos no directorio django_mobile_utils os arquivos __init__.py, context_processor.py e un directorio data/mobile_agents.txt, e o settings.py modificado como corresponde. ¿Cómo ó usamos? Simplemente, nos templates utilizamos a variable booleana {{mobile_browser}} como condicional, por exemplo para cargar unha folla de estilo diferente:
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.
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:
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
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.
