sábado, 25 de agosto de 2012

Metro en man: medindo un PDB

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))

positions vale algo así como:

[[ 133.927  117.426  173.587]
 [ 133.405  117.299  174.43 ]
 [ 133.622  116.764  172.902]
 ..., 
 [ -27.042   56.598    9.173]
 [ -26.066   57.535    9.316]
 [ -27.61    56.209    8.   ]]

Calculando o tamaño do PDB

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.

lunes, 9 de enero de 2012

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

lunes, 24 de octubre de 2011

Diferentes páxinas para dispositivos móviles e fixos

Problema

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:
def mobile_browser(request):
    dict = {'mobile_browser': False}
    if hasattr(request, 'is_mobile'):
        dict['mobile_browser'] = request.is_mobile
    return dict
Podemos pasar a outro context máis simple aínda, como:
def mobile_browser(request):
    dict = {'mobile_browser': False}
    if our_server.django_mobile_utils.is_mobile(request):
        dict['mobile_browser'] = True
    return dict
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:
TEMPLATE_CONTEXT_PROCESSORS = (
    "our_server.django_mobile_utils.context_processors.mobile_browser",
    ...

Aplicación

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:
{% if mobile_browser %}
  <link href="/media/css/mobile_style.css" rel="stylesheet" type="text/css" />
{% else %}
  <link href="/media/css/style.css" rel="stylesheet" type="text/css" />
{% endif %}

miércoles, 13 de julio de 2011

Conseguir a secuencia de un PDB

Problema



Temos un archivo .pdb, e queremos extraer a secuencia de aminoácidos de unha soa letra. Atopei algunha solución por internet, pero implica importar a librería Bio.PDB.Polypeptide.PPBuilder, e resultábame un método algo escuro. Así que busquei ata atopar o diccionario que traduza código de tres letras a código de unha letra.

Librerías



from Bio.SCOP.Raf import to_one_letter_code
from Bio.PDB import PDBParser


O diccionario está algo escondido, o módulo Raf (ASTRAL RAF (Rapid Access Format) Sequence Maps) é a primeira vez que me cruzo con él. Penso que debería estar referenciado en Bio.Alphabet.IUPAC un diccionario similar de traducción tres<->unha letras.

Código



Creamos un "motor" para parsear archivos PDB (p), parseamos o archivo en s e logo extraemos a secuencia dese archivo s a unha secuencia y, que traducimos con to_one_letter_code:

p = PDBParser()

s = p.get_structure("id", "/path/to/file.pdb")
y = [x.get_resname() for x in s.get_residues()]
sequence = "".join([to_one_letter_code[z] for z in y])

¡Listo! En sequence temos un string coa secuencia do PDB. Non sei si é máis lento ou máis ineficiente que o método do PPBuilder, pero desde logo é o método que eu estaba buscando: sacar a secuencia tri-letra do PDB e traducila. Haberá outras situacións nas que teñas unha secuencia tri-letra que queiras levar a unha letra sin ter que chamar polo PPBuilder. Nese caso, importando o diccionario e facéndolle:

[to_one_letter_code[z] for z in y]

Xa o teremos.

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.

viernes, 29 de abril de 2011

Generic views con ForeignKey(User)

Con generic views podemos facer os accesos máis comúns á base de datos (engadir, borrar, editar, consultas...). A veces haberá que facer pequenos engadidos de código para conseguir axustes. Hoxe vamos ver cómo subir archivos ó servidor sin perder o trazo dos donos deses archivos.

models.py


Queremos aceptar un archivo que o usuario sube, e que o usuario lle dea unha descripción ó archivo. Comprobar que ese archivo sea seguro (que non nos suban un script, por exemplo) é outra tarea distinta que queda nas mans do programador.

from django import forms
from django.db import models
from django.forms import ModelForm
from django.contrib.auth.models import User
import uuid

def set_path(my_model, file_name):
'''Sets the file_path for the uploaded file'''
return "uploads/{0}".format(uuid.uuid1())

# Create your models here.
class OurModel(models.Model):
name = models.CharField(max_length=50)
user = models.ForeignKey(User, editable=False)
file_path = models.FileField(upload_to=set_path, blank=True)

def __unicode__(self):
return self.name

class OurModelForm(ModelForm):

class Meta:
model = OurModel
exclude = ("user",)


urls.py



Estando dentro da nosa aplicación, empezaremos necesitando a vista de list_detail.object_detail e a de create_update.create_object. Como imos a modificar aparte de crear o obxeto, a vista xenérica de create_object hai que desplazala a unha función de views.py (será add).

from django.conf.urls.defaults import *
from models import OurModel, OurModelForm

import views

urlpatterns = patterns('',
url(r'^(?P<object_id>\d+)/$',
"django.views.generic.list_detail.object_detail",
dict(queryset=OurModel.objects.all()),
name='ourmodel_detail'),
url(r'^add/$',
views.add,
name='ourmodel_add'),
)


views.py



Aquí creamos a función add que se chama desde urls.py. Nesta función trataremos as chamadas GET e POST, xa que o formulario ten o "action" cara sí mesmo. Si se chama vía GET, devolvemos a vista de create_object, si se chama vía POST, instanciamos o modelo cos datos que nos veñen do formulario, instanciamos un usuario (User) co usuario que está logueado e pasámoslle ambas instancias ó modelo, para que as grave.

O "truco" está en cargar o formulario e preparalo para grabar (commit=False), xa que a partir deste momento poderemos cargarlle tamén a instancia de "User". Despois da grabación, podemos redirixir o usuario cara calquer páxina, aquí mandámolo a un object_detail do recén grabado.

from django.views.generic.create_update import create_object
from django.views.generic.list_detail import object_detail
from django.contrib.auth.decorators import login_required
from django.contrib.auth.models import User

from models import OurModel, OurModelForm

@login_required
def add(request):
if request.method == "POST":
form = OurModelForm(data=request.POST,
files=request.FILES)
if form.is_valid():
our_model = form.save(commit=False)
our_model.user = User.objects.get(pk = request.user.id)
our_model.save()

return object_detail(request,
OurModel.objects.all(),
object_id=align.id)

return create_object(request,
model=OurModel,
form_class=OurModelForm,
login_required=True)


templates



Necesitamos só dous templates, /our_app/ourmodel_detail.html e /our_app/ourmodel_form.html. Con un pequeno cambio: os formularios que suben archivos teñen unha propiedade añadida

<form enctype="multipart/form-data" method="post" action="."></form>


Listo!

Resultado



Cando un usuario rexistrado accede á páxina http://ourserver.com/ourapp/add/ por primeira vez (vía GET), o urls.py lévao á función def add de views.py. Nesa función sáltase toda a parte "non POST", e simplemente saca o generic view para create_update.create_object, que presenta o formulario indicado.

Si o usuario enche o formulario correctamente (se non o fai, o generic view automáticamente volve para atrás), chámase a función add de novo, pero agora vía POST e cos datos apropiados. Agora entramos no proceso propio de grabar o modelo na base de datos cos datos apropiados.