Como ejecutar DALI para el alineamiento de proteínas con las librerías de Fortran originales

DALI es un programa de comparación de estructuras de proteínas de finales del siglo pasado. A pesar de ser un programa ampliamente usado, el código hace décadas que no ha sido tocado por lo que sigue usando librerías de unos 20 años de antigüedad. Librerías Fortran que ya no están en la versión actual del gcc.

El problema que me daba es:

error while loading shared libraries: libg2c.so.0: cannot open shared object file: No such file or directory

Como quería usar este programa tuve que encontrar una solución. Al final lo que conseguí es instalar una versión antigua de gcc para poder usar las librerías usando el compilador antiguo de Fortran (gfortran).

Al final encontré la única solución. Consiste en instalar el compilador de GNU y después compilar y ejecutar dalilite usando la versión antigua.

Para empezar nos bajaremos la versión antigua del compilador de GNU, la versión gcc-3.4.6

wget https://ftp.gnu.org/gnu/gcc/

descomprimiremos el archivo comprimido, y crearemos el directorio dónde pondremos las librerías

tar -zxvf gcc-3.4.6.tar.gz
cd gcc-3.4.6
mkdir build
cd build

Configuraremos gcc-3.4.6

../configure --enable-shared –prefix=/cualquier/directorio --disable-multilib --disable-bootstrap # /cualquier/directorio puede ser por ejemplo /home/roc/herramientas/gcc-3.4.6

Ahora hacemos un build y lo instalamos. Esto va a tardar un rato y se instará en el directorio previamente seleccionado

make -j 4
make install

Ahora que ya tenemos nuestra copia de gcc-3.4.6 instalada en /cualquier/directorio necesitamos hacer las librerías visibles.

export PATH=/cualquier/directorio/bin:$
export LD_LIBRARY_PATH=/cualquier/directorio/

Ahora dalilite debería funcionar. Recuerda en recompilarlo con el Makefile. Además cada vez que quieras ejecutarlo deberías usar las dos lineas anteriores

Espero que os haya servido!

Introducción a TensorFlow

La web oficial de TensorFlow tiene muy buenos recursos. En esencia lo que hay en este post proviene del “get started” de la web oficial.

En el primer ejemplo importaremos TensorFlow. Crearemos dos constantes y las imprimiremos en pantalla.

import tensorflow as tf
node1 = tf.constant(3.0, tf.float32)
node2 = tf.constant(4.0) # También tf.float32 de forma implícita
sess = tf.Session()
print(sess.run([node1, node2]))
node3 = tf.add(node1, node2)
print("node3: ", node3) #Esta linea muestra las propiedades del tensor
print("sess.run(node3): ",sess.run(node3)) # Aquí se imprime el resultado

También podemos aplicar operaciones a los tensors. Las operaciones producen más tensores. En el siguiente ejemplo sumamos dos variables (tensores) y luego las sumamos. Un tensor es capaz de procesar listas también como veremos.

a = tf.placeholder(tf.float32)
b = tf.placeholder(tf.float32)
adder_node = a + b
print(sess.run(adder_node, {a: 3, b:4.5})) # 7.5
print(sess.run(adder_node, {a: [1,3], b: [2, 4]})) # Itera sobre la lista: 3 y 7

Podemos multiplicar

add_and_triple = adder_node * 3.
print(sess.run(add_and_triple, {a: 3, b:4.5}))

E incluso podemos crear modelos lineales

W = tf.Variable([.3], tf.float32)
b = tf.Variable([-.3], tf.float32)
x = tf.placeholder(tf.float32)
linear_model = W * x + b

Las constantes son inicializadas con tf.constant, y su valor no puede ser cambiado. Contrariamente las variables son creadas usando tf.Variable y no son inicializadas en el inicio. Para inicializar las variables  tenemos que usar una función especial como en el siguiente ejemplo:

init = tf.global_variables_initializer()
sess.run(init)

Lo importante es usar init. Es un “handle” que inicia el grafo, inicializando las variables que hasta entonces habían permanecido sin inicializar.

Ya que la X es la variable podemos evaluar el modelo usando distintos valores simultáneamente de la siguiente manera:

print(sess.run(linear_model, {x:[1,2,3,4]}))
# Output: [ 0, 0.30000001, 0.60000002, 0.90000004]

La función de pérdida (loss function) puede medir la distancia entre el modelo actual y los datos proporcionados. Usaremos la función de pérdida estándar para hacer una regresión lineal que sume los cuadrados de las distancias entre el modelo actual y los datos. El modelo crea un vector en el que cada posición corresponde a una distancia de error. Cuando usamos tf.square estamos elevando al cuadrado el error. Luego sumamos los errores para crear un simple escalar que abstrae el error de todos los puntos usando tf.reduce_sum:

y = tf.placeholder(tf.float32)
squared_deltas = tf.square(linear_model - y)
loss = tf.reduce_sum(squared_deltas)
print(sess.run(loss, {x:[1,2,3,4], y:[0,-1,-2,-3]}))

El valor producido es:  23.66

Podríamos mejorar el modelo manualmente asignando nuevos valores a W para obtener resultados perfectos de -1 y 1. Una variables es inicializada por el valor proporcionado por tf.Variable pero puede ser cambiada si usamos operaciones como tf.assign. Por ejemplo, W=-1 y b=1 son los parámetros óptimos de nuestro modelo. Podemos cambiar W y b acordemente:

fixW = tf.assign(W, [-1.])
fixb = tf.assign(b, [1.])
sess.run([fixW, fixb])
print(sess.run(loss, {x:[1,2,3,4], y:[0,-1,-2,-3]})
#resultado: 0

Hemos adivinado el valor perfecto para W i b aunque este no es el objetivo de  machine learning. El objetivo es encontrar los parámetros correctos para el modelo automáticamente. Pero esto lo vais a encontrar en la próxima entrega 😉

Bonus: En matemáticas y en física, un tensor es cierta clase de entidad algebraica de varias componentes, que generaliza los conceptos de escalar, vector y matriz de una manera que sea independiente de cualquier sistema de coordenadas elegido. ( wikipedia )

Script en python para convertir secuencias de proteínas de Stockholm a fasta

Aquí os dejo un pequeño python script que convierte “multiple sequence alignments” del formato Stockholm a Fasta de una forma sencilla y rápida.

import sys
from Bio import SeqIO
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord 

if(len(sys.argv) <3):
    print('two arguments needed: input path, output path')
    exit(2) 

with open(sys.argv[1],'r') as inFile:
    with open(sys.argv[2], "w") as output_handle:
        SeqIO.write(list(SeqIO.parse(inFile,'stockholm')), output_handle, "fasta")

Texmaker al estilo de Sublime text (monokai)

Sublime text, aunque es un editor de texto que me gusta mucho, he estado buscando alternativas. He encontrado textmaker. Texmaker es  un editor de latex que nos permite compilar el documento pdf además de añadir otras interesantes features. Lo que no me gustaba era el tema (theme) que tenía por defecto. Así que buscando di con la solución. Usar el theme de Sublime text, llamado monokai.

Para cambiarlo simplemente tenemos que ir a $HOME/.config/xm1 (tanto en linux, como en MacOS y editar el fichero texmaker.ini y reemplazar las siguientes lineas de la principio del fichero:

Color\Background=@Variant(\0\0\0\x43\x1\xff\xff\0\0++66\0\0)
Color\Command=@Variant(\0\0\0\x43\x1\xff\xff&&\x8b\x8b\xd2\xd2\0\0)
Color\Comment=@Variant(\0\0\0\x43\x1\xff\xffllqq\xc4\xc4\0\0)
Color\Highlight=@Variant(\0\0\0\x43\x1\xff\xff\a\a66BB\0\0)
Color\Keyword=@Variant(\0\0\0\x43\x1\xff\xff\xdc\xdc\x32\x32//\0\0)
Color\KeywordGraphic=@Variant(\0\0\0\x43\x1\xff\xff\x85\x85\x99\x99\0\0\0\0)
Color\Line=@Variant(\0\0\0\x43\x1\xff\xff\a\a66BB\0\0)
Color\Math=@Variant(\0\0\0\x43\x1\xff\xff**\xa1\xa1\x98\x98\0\0)
Color\NumberGraphic=@Variant(\0\0\0\x43\x1\xff\xff\xcb\xcbKK\x16\x16\0\0)
Color\Standard=@Variant(\0\0\0\x43\x1\xff\xff\x83\x83\x94\x94\x96\x96\0\0)
Color\Todo=@Variant(\0\0\0\x43\x1\xff\xff\xd3\xd3\x36\x36\x82\x82\0\0)
Color\Verbatim=@Variant(\0\0\0\x43\x1\xff\xff\xb5\xb5\x89\x89\0\0\0\0)

Reducción de dimensiones: Principal Component Analysis (PCA)

Principal Component Analysis o PCA en corto es un método de reducción de dimensiones bastante conocido y comúnmente usado. Este método transforma ortogonalmente las observaciones (quizás relacionadas) en un conjunto de puntos linealmente no relacionados. De esta forma se consigue que el primer componente tenga la varianza mayor. El siguiente componente será el que tendrá la varianza mayor de los restantes y es ortogonal al anterior componente, y así sucesivamente. Este método es sensible a la escala de las variables. Para evitar que los valores sean un problema tendremos que escalar las variables estandardizándolos antes de usar PCA.

Generar una barra de progreso en ipython notebooks

Manual con código para generar una barra de progreso en nuestro código para saber en que porcentaje de compleción estamos sin llenar el output con números. Muchos de los que usáis jupyter (el nuevo ipython notebooks) podéis imprimir por pantalla la iteración en la que vuestro loop reside. Eso es solo posible para una cantidad pequeña de iteraciones. Cuando llegamos a varios miles se puede generar un output bastante engorroso. Googleando un poco encontré la solución. Este código nos genera una clase que luego podemos llamar para obtener la barra de progreso junto con el porcentaje completado.

import sys, time
try:
    from IPython.core.display import clear_output
    have_ipython = True
except ImportError:
    have_ipython = False

class ProgressBar:
    def __init__(self, iterations):
        self.iterations = iterations
        self.prog_bar = '[]'
        self.fill_char = '*'
        self.width = 40
        self.__update_amount(int(0))
        if have_ipython:
            self.animate = self.animate_ipython
        else:
            self.animate = self.animate_noipython

    def animate_ipython(self, iter):
        try:
            clear_output()
        except Exception:
            # terminal IPython has no clear_output
            pass
        print '\r', self,
        sys.stdout.flush(),
        self.update_iteration(iter + 1)
        
    def update_iteration(self, elapsed_iter):
        self.__update_amount((elapsed_iter / float(self.iterations)) * 100.0)
        self.prog_bar += '  %d of %s complete' % (elapsed_iter, self.iterations)

    def __update_amount(self, new_amount):
        percent_done = int(round((new_amount / 100.0) * 100.0))
        all_full = int(self.width - 2)
        num_hashes = int(round((percent_done / 100.0) * all_full))
        self.prog_bar = '[' + self.fill_char * num_hashes + ' ' * (all_full - num_hashes) + ']'
        pct_place = int((len(self.prog_bar) / 2) - len(str(percent_done)))
        pct_string = '%d%%' % percent_done
        self.prog_bar = self.prog_bar[0:pct_place] + \
            (pct_string + self.prog_bar[pct_place + len(pct_string):])

    def __str__(self):
        return str(self.prog_bar)
        

Y el siguiente código vemos un ejemplo de como llamaremos la clase para obtener la funcionalidad deseada

c = ProgressBar(1000)
for i in range(1000):
    c.animate_ipython(i)

Para los interesados encontré el código aquí.

Usar listas en matplotlib en vez de un punto a la vez

TLTR: Cuando uses matplotlib en python pasa a la función dos listas con todos los puntos en vez de ir uno por uno. Así te ahorras el rendering cada vez que usas la función consiguiendo incrementos de velocidad de varias magnitudes.

Hasta hace poco lo que hacía mayoritariamente cada vez que tenia que hacer una gráfica con matlplotlib era poner un punto cada vez. ERROR!!! Hasta el momento asumía que el tiempo de ejecución del script era debido a los algoritmos que usaba. Pero no, era debido a una  malapraxis muy simple de resolver. Últimamente estoy trabajando con cantidades más grandes de datos y la ejecución de mi código en python me tardaba más de lo que debería hasta que descubrí que hacia un error de novato. Poner punto por punto a una gráfica en vez de pasarle la lista entera era lo que causaba una ejecución tan lenta de mi código. Iba punto por punto porque era como me llegaban los datos y no los tenia que guardar, pero resulta que si los pongo en una lista y luego paso la lista a la librería el resultado es infinitamente más veloz y eficaz. La librería se ahorra renderizar la gráfica cada vez si recibe una lista.

El siguiente código tiene una ejecución de más de 6 segundos y como veis es la forma incorrecta de hacer un scatter plot:

start_time = timeit.default_timer()
for x in range(100):
for y in range(10):
matplotlib.pyplot.scatter(x,y)

matplotlib.pyplot.show()
end_time = timeit.default_timer()
print(end_time - start_time)

El siguiente código es la manera rápida de hacer un scatter plot y tarda unos 300ms (si MILI-segundos).

xlist = []
ylist = []

start_time = timeit.default_timer()
for x in range(100):
for y in range(10):
xlist.append(x)
ylist.append(y)

matplotlib.pyplot.scatter(xlist, ylist)
matplotlib.pyplot.show()
end_time = timeit.default_timer()
print(end_time - start_time)

He subido un python notebooks (jupyter notebook) a github para que lo ejecutéis vosotros si queréis.

Reducción de dimensiones: Self-organizing feature map (SOFM)

Self-organizing map (SOM) o self-organizing feature map (SOFM) es un método que usa redes neuronales (neuronal networks) para reducir las dimensiones de un vector de datos. Para reducir las dimensiones lo que hace es usar los vecinos de un punto en concreto para moverlo al nuevo espacio dimensional manteniendo la misma topografía que en el espacio original.

Una de las ventajas de SOFM es que es un método unsupervised, esto quiere decir que no necesitamos un training dataset para entrenar nuestra red neuronal. Normalmente para la inteligencia artificial (y el machine learning) se requiere que el training dataset esté etiquetado. Un data set etiquetado son puntos del input que ya se sabe cual tiene que ser el resultado. De este modo la maquina es capaz de identificar que resultado tiene que producir.

Para explicar con más detalle este método voy a usar la foto de la wikipedia que lo hace más visual.

Somtraining
SOMF training

 

En la imagen vemos la malla que representa el input en un espacio multidimensional y la zona azul que representa un espacio 2D al que queremos extrapolar nuestros puntos. El primer paso es seleccionar nuestro punto en el espacio de origen que se acerque mas al espacio 2D. El segundo paso es mover ese punto en concreto al espacio 2D. El tercer paso, y ultimo, consiste en mantener las distancias entre puntos en el nuevo plano. Si se consigue mantener las distancias entre los puntos conseguiremos una reducción de dimensiones satisfactoria.

El paper original lo podéis encontrar aquí. Pero si no tenéis ninguna afiliación académica quizás os interese checkear este otro post.

Reducción de dimensiones: T-SNE

Como ya explicamos en el post anterior los ordenadores si que pueden procesar grandes cantidades de datos multidimensionales. Pero los humanos a veces necesitamos “ver” y entender los datos. Cuando estamos trabajando en un espacio multidimensional no podemos imaginarnos nuestro dataset. Para solventar este problema se ha desarrollado T-SNE. Éste es un algoritmo pensado especialmente en reducir dimensiones (a 2 o 3) para que podamos visualizar los datos. En el paper original se expresa de forma explícita que el algoritmo está pensado para la visualización de datos. Usar T-SNE para la reducción de dimensiones puede causar efectos desconocidos. Así que si lo usas para evitar “la maldición de la dimensionalidad”  allá tu 😉

El concepto principal consiste en que los puntos cercanos (en el espacio multidimensional) se atraen y los distantes se repelen. Para conseguir este objetivo el algoritmo tiene unos cuantos parámetros que se permiten alterar. La “perplexity” es la cantidad de vecinos que un simple punto puede afectar. Por lo que he visto hasta ahora un valor entre 30-50 suele ser el óptimo. Luego tenemos epsilon que nos sirve para determinar el tamaño de los pasos de aprendizaje. Valor pequeño el algoritmo le cuesta más a encontrar el óptimo. Con un valor grande te lo puedes pasar. Finalmente la cantidad de iteraciones o steps para conseguir la convergencia. A más iteraciones en teoría más cerca del valor óptimo, pero como es evidente como más iteraciones más tiempo de computación requieres.

Para los que queráis visualizar como funciona y un poco sus efectos cambiando los parámetros podéis verlo en vuestro navegador. Si por el contrario queréis implementarlo en Javascript aquí podéis encontrar una versión oficial del desarrollador principal (es la que estoy usando). Pero si sois tan malotes que queréis implementarlo vosotros mismos mejor que os leáis el paper original.

Bonus: En el paper para computar las distancias se usa la distancia euclidiana (la “normal”) quizás otro tipo de distancias puede ir mejor para vuestro problema.

Reducción de dimensiones: Introducción a los espacios multidimensionales

En inteligencia artificial y machine learning en la mayoría de ocasiones se usan espacios multidimensionales. Los espacios multidimensionales son espacios en los que los datos requieren más de un valor. Los espacios multidimensionales son espacios con puntos repartidos por el espacio. Un espacio 2D tiene dos dimensiones, las típicas X, Y. Un espacio 3D tiene tres dimensiones X, Y, Z. Las dimensiones las puedes llamar X, Y, Z pero también perro, gato, conejo. Como seres humanos estamos limitados a poder visualizar espacios a 3D máximo. Aunque no puedas visualizar dimensiones superiores a tres, puedes llegar a entenderlo. Imagínate que quieres ir a un restaurante a 10km máximo de tu casa. Esto en un mapa son dos dimensiones norte/sur y este/oeste. Ahora si decides que quieres ir a un restaurante italiano ya hemos añadido otra dimensión. Si no quieres pagar más de 10€ por la cena ya has añadido otra dimensión. Por lo que ya tenemos un espacio de 4D. Dos para las dimensiones del mapa norte/sur, este/oeste, otra para que tipo de restaurante y la cuarta para los precios. Se le podría sumar otra dimensión si añadiéramos en que planta esta el restaurante. Virtualmente podríamos añadir infinidad de dimensiones.

Un vector o un punto en nuestro espacio es un conjunto de valores para todas las dimensiones. Un valor (o dato) puede ser indefinido o ausente, no tiene porque tener un valor concreto. Un punto es asociado a la persona o el elemento de la acción. Tu puedes ser ese vector/punto. Con el ejemplo anterior nuestro punto sería (23, 31, italiano, 10) que hacen referencia a norte/sur, este/oeste, tipo de comida y precio máximo. Las variables pueden ser categóricas y cuantitativas. Las categóricas es “italiano” no tiene un valor númerico asociado. Otros ejemplos de variables categóricas pueden ser el sexo de una persona o su color preferido. Las variables cuantitativas son los números.

La pregunta del millón es: cómo podemos pensar en un espacio de más de tres dimensiones en el que una de ellas es “italiano”? Simplemente no podemos. Nuestro cerebro no puede visualizar más de tres dimensiones. Estos espacios sólo tienen sentido matemáticamente. Por lo que para trabajar en este tipo de espacios si queremos visualizarlos de algún modo tenemos que pensar en tres dimensiones máximo y extrapolarlo. Otra opción también puede ser usar distintas técnicas de reducción de dimensiones. En los próximos posts describiré distintos métodos para la reducción de dimensiones.