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.

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.

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