martes, 12 de julio de 2016

Botas, batas y bits

La verdad es que hace no tanto tiempo que si alguien decía -soy bióloga- lo primero que se pensaba era, ¡de bota! y se podía imaginar a la intrépida investigadora, amiga de los animales y las plantas, recorriendo el monte, identificando, clasificando y recogiendo muestras de todo tipo. Hacia la última década del siglo XX algunos biólogos comenzaron a matizar -eh que yo soy de bata, no de bota- y se les podía ver cuán jóvenes CSI captando la pequeña muestra en el escenario adecuado que en el laboratorio llevaría hasta el ADN y de ahí a la prueba o identificación certera... Lo cierto es que, ya en pleno siglo XXI, los biólogos seguimos, por supuesto, siendo amigos de plantas y animales, pero ya no sólo somos de bota o de bata. Este es el siglo de la tercera B biológica: de bota, de bata y... ¡de bit!



Como en tantos otros campos los ordenadores han supuesto una revolución también en la biología. Pero esta es más profunda de lo que muchos piensan. Al hablar de biología e informática, quizás el término bioinformática nos viene a la mente. Este es un concepto actualmente muy ligado a la genética y el análisis molecular del ADN. Los ordenadores permitieron en primera instancia mejorar muchísimo la eficiencia del proceso de secuenciación. En consecuencia la cantidad de datos disponibles de ADN de muchas especies, poblaciones e individuos pertenecientes a ellas es actualmente enorme, y sigue creciendo. Cuando hablamos de ADN sabemos que estamos tratando con moléculas de ácido desoxirribonucleico en las que subyace la información genética de los organismos. Una de las informaciones más interesantes que podemos extraer de estas moléculas es la secuencia de bases nitrogenadas (Adenina, Citosina, Guanina y Timina) que caracterizan los diferentes genes y que pueden variar entre especies, poblaciones e incluso, a veces, entre individuos de la misma población. Para representar estas moléculas una vez secuenciadas sólo necesitaríamos papel y lápiz, o un editor de texto en un ordenador, y anotar la secuencia de bases nitrogenadas. Por ejemplo, si hemos secuenciado un fragmento de 5 bases, esta secuencia podría ser AATCA. Pero si en vez de una secuencia de 5 bases hablamos de varias secuencias con muchos miles, cientos de miles o incluso millones de bases entonces los investigadores necesitan utilizar pequeños programas (scripts) para manejarlas de modo que las pueden copiar de unos ficheros a otros, contar el número de bases de un determinado tipo, buscar una secuencia o patrón concreto de bases, comparar unas secuencias con otras, etcétera.

Créditos de la imagen: Flickr/CC BY-SA 2.0, www.yourgenome.org, 
 www.jax.org, en.wikipedia.org, www.foodsafetynews.com
                                                 
Sin embargo, como ya dijimos, la revolución del bit en biología, es más profunda y por tanto no sólo se trata de los programas y bases de datos para manejar las secuencias. A comienzos del siglo XX, cuando la integración de la genética mendeliana y la teoría evolutiva dio lugar a la disciplina conocida como genética de poblaciones, ello incentivó importantes desarrollos estadísticos para la resolución de los problemas que la nueva disciplina planteaba. Actualmente, en vez de tratar con unos pocos genes tratamos con los genomas enteros. La integración de estos datos genómicos con los modelos poblacionales y ecológicos, plantea nuevos problemas y requiere el desarrollo de nuevos métodos y modelos.

Si estudiamos un proceso biológico que cambia a lo largo del tiempo (dinámico) podemos considerar al menos 3 tipos de modelos. En primer lugar, modelos de tipo biológico (por ejemplo, datos de apareamientos entre dos poblaciones vecinas de una especie),  que consisten en información (datos, mediciones,…) que está sometida a ruido. Este ruido proviene del muestreo, errores de medida o simplemente de eventos aleatorios de la naturaleza. El segundo tipo, el modelo matemático, sería una representación formal del modelo biológico en la cual se establecen relaciones entre las características de interés. En estos modelos los números no están sometidos a ruido. Si iteramos el modelo matemático podríamos resolver las variables fundamentales y obtener soluciones predictivas para el mismo. Si el ajuste entre el modelo matemático y el biológico es aceptablemente bueno entonces estas predicciones podrán ser muy útiles. En nuestro ejemplo, podríamos relacionar frecuencias de apareamiento con determinadas características de los individuos y de estos con su composición genética y así estaríamos en posición de conocer la evolución de los fenotipos y genotipos de esas poblaciones tanto hacia adelante como hacia atrás en el tiempo. El problema de estos modelos matemáticos es que, o bien son muy simples (no se ajustan al modelo biológico) o bien no somos capaces de resolverlos cuando el modelo biológico subyacente es muy complejo. Entra aquí en juego el tercer tipo, los modelos de simulación, los cuales resuelven los cálculos de modelos dinámicos matemáticos complejos. Es decir, los modelos de simulación nos permiten plantear y resolver modelos matemáticos complejos de modo que mejore el ajuste respecto a los modelos biológicos.



Y hemos llegado al mensaje clave de esta entrada. La necesidad cada vez más acuciante de desarrollar simulaciones de ecosistemas evolutivos a la luz de los datos genómicos y poblacionales. Si bien es cierto que este tipo de modelos llevan ya un tiempo aplicándose en biología, por ejemplo los llamados modelos basados en individuos (IBM por sus siglas en inglés) que se utilizan en ecología y genética de poblaciones, es ahora cuando emerge una extensión de los mismos combinando información genómica, ecológica y demográfica. Hablamos de modelos eco-genéticos o mejor aún eco-genómicos. Estos nos permitirán avanzar en el conocimiento teórico de aspectos muy diversos e importantes de la ecología y la evolución de las especies a la vez que realizar y contrastar predicciones. Estos modelos se aplican tanto en conservación, como en biomedicina, también para poner a punto estimadores de estadísticos de interés, en estudios de impacto de actividades humanas sobre la diversidad genética, pesquerías, etc., etc.

A pesar de que existen algunas herramientas de simulación eco-genética, pocas veces se ajustan exactamente a nuestras necesidades específicas. Hace tiempo escuché una gran verdad que personalmente he experimentado muchas veces -entendemos realmente algo cuando sabemos programarlo, en caso contrario sólo creemos entenderlo-. Surge entonces la cuestión clave. ¿Qué podemos hacer si somos biólogos y queremos programar nuestros propios bio-simuladores? Como todas las disciplinas que se mueven en la espumeante cresta de la ola del conocimiento, el campo de la simulación biológica, que enmarcaré dentro de la biología computacional, es una ciencia multidisciplinar donde se echa mano de conocimientos de genética, ecología y evolución pero también de programación, estadística y matemáticas. Respecto a los lenguajes de programación, aunque mi favorito es C++ quizás Python (https://www.python.org/) es el que mejor combina simplicidad con buenas prestaciones en términos de estructura y posibilidades del lenguaje. Dispone de  funciones predefinidas de utilidad en biología siendo relativamente sencillo realizar los primeros programas útiles. Existen además infinidad de tutoriales, comunidades y foros donde la persona interesada encuentra ayudas de todo tipo para dar los primeros pasos y avanzar en su conocimiento “pythonico”. Así que finalizo este post con una exhortación para los  nuevos biólogos de bits, os digo: ¡simulad, simulad malditos!

Para saber más

Haddock, S. H. D., and C. W. Dunn. 2011. Practical computing for biologists. Sinauer, Sunderland, Mass.
Hartmann, A. K. 2015. Big Practical Guide to Computer Simulations. World Scientific.
Hoban, S. 2014. An overview of the utility of population simulation software in molecular ecology. Molecular Ecology 23:2383-2401.
Perkel, J. M. 2015. Programming: Pick up Python. Nature 518:125-126


martes, 3 de noviembre de 2015

Generación de variables aleatorias con la librería Random de C++11


A partir de la nueva versión del estándar de C++ conocida como C++11 (o en su última actualización C++14 http://isocpp.org/blog/2014/08/we-have-cpp14)  la librería estándar aporta un buen número de clases y tipos relacionados con la generación de números aleatorios y sus distribuciones. Esta es una novedad de interés para la biología computacional puesto que gran cantidad de algoritmos y simulaciones se apoyan en el modelo de Monte Carlo cuyo pilar fundamental es la generación de números aleatorios. Antes de esta inclusión era el propio desarrollador quién debía encargarse de construir su propia librería de funciones aleatorias o importar alguna de distribución libre. Algunas librerías de este tipo eran las antiguas genlib.c y randlib.c,  o, más recientemente, ALGLIB http://www.alglib.net/  y la librería científica de la GNU (GSL) http://www.gnu.org/software/gsl/ que incluían, además de otras muchas funciones matemáticas, también las relacionadas con la generación de números aleatorios y distribuciones de probabilidad.

Sin embargo, la nueva librería de números aleatorios incluida en C++11 además de algunas de las más importantes distribuciones, aporta también algunos de los mejores y más modernos generadores de aleatoriedad. Sin entrar en detalles que nos llevarían más allá de la intención de este post, podemos aclarar que en la generación de variables aleatorias distinguimos dos aspectos fundamentales. Primero, necesitamos un mecanismo generador de aleatoriedad (random engine) es decir una fuente que genere patrones "impredecibles" de bits. En segundo lugar, necesitamos funciones que aprovechando la susodicha fuente de aleatoriedad generen números según alguna distribución de probabilidad deseada. De manera un poco confusa, a la combinación de un generador de aleatoriedad y una distribución de probabilidad se le llama generador de números aleatorios (random number generator, RNG).

Debe notarse que el generador de aleatoriedad es en realidad un generador de pseudoaleatoriedad (PRNG) es decir se genera aleatoriedad desde el punto de vista estadístico pero basándose en un modelo causal y determinista. Una simple función matemática que para un valor de entrada (semilla) genera una secuencia (muy larga) de valores de salida. Estos valores deben ser estadísticamente uniformes e independientes. La impredecibilidad radica en que estas funciones generan resultados muy distintos para pequeñas variaciones en los estados iniciales. Por tanto, si se cambia la semilla, el estado inicial del mecanismo generador cambia y la secuencia subsiguiente de números aleatorios será muy distinta a la anterior. La ventaja del uso de procesos pseudoaleatorios es que el desarrollador puede, para una misma semilla, generar exactamente la misma secuencia de números una y otra vez, lo cual es extremadamente útil a la hora de comprobar y depurar código, corregir errores, etc.

La clave para conseguir emular la deseada impredecibilidad de muchos fenómenos estocásticos  será conseguir semillas del modo menos predecible posible. Una práctica muy común es usar los milisegundos (o nanosegundos) del reloj del procesador. Esta práctica sigue siendo perfectamente válida para la mayoría de aplicaciones en biología computacional si bien ha quedado obsoleta para algunos usos como la seguridad informática (por ejemplo el hackeo del sitio web de tecnología Hacker News se pudo hacer gracias a que usaban números pseudolaeatorios incializados con una semilla basada en el reloj del sistema en milisegundos: https://blog.cloudflare.com/why-randomness-matters/).

Por otro lado, para conseguir generadores realmente aleatorios se pueden usar procesos basados en hardware como alteraciones debidas al ruido electrónico o efectos cuánticos, sin embargo, desde el punto de vista de la simulación, se pierde eficiencia, reproducibilidad y portabilidad.  Por tanto es el uso de generadores pseudoaleatorios el que está más extendido, siendo el elemento fundamental de la simulación Monte Carlo y por tanto una parte importantísima de la biología computacional.

El lector interesado puede encontrar un buen tutorial sobre la generación de números aleatorios en C++11 en el siguiente enlace: http://www.open-std.org/jtc1/sc22/wg21/docs/papers/2013/n3551.pdf.

En lo que sigue haré un breve resumen sobre los generadores aleatorios en C++11 y algunos ejemplos de su uso en simulación de biología de poblaciones.

Mecanismos generadores de aleatoriedad en la librería estándar de C++11

Según el tipo de algoritmo que se utilice para generar los patrones aleatorios de bits distinguimos entre diferentes clases de generadores. Por ejemplo, los generadores congruentes lineales (http://en.wikipedia.org/wiki/Linear_congruential_generator) y uno de  los más usados actualmente, el algoritmo "Mersenne Twister" http://en.wikipedia.org/wiki/Mersenne_Twister.

Recordad que necesitábamos dos componentes para generar variables aleatorias: un mecanismo generador y una distribución. En C++11 vienen predefinidos varios mecanismos generadores de pseudoaleatoriedad, lo mismo que diferentes distribuciones. Un resumen de todo ello se puede consultar en el siguiente enlace http://en.cppreference.com/w/cpp/numeric/random.

Personalmente suelo utilizar uno de los generadores más comunes como es el std::mt19937 o su versión de 64 bits std::mt19937_64 aunque para casos puntuales donde se van a generar pocos valores, como por ejemplo dentro de una función de muestreo sin reemplazamiento, uso el generador por defecto std::default_random_engine.

El siguiente paso es generar la semilla (seed) que inicializa la serie. Una opción es

auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();

que guarda en seed el número de tics de reloj desde el tiempo Unix o epoch (https://es.wikipedia.org/wiki/Tiempo_Unix). Como este número depende de la precisión del reloj del sistema (microsegundos en mi máquina) podemos querer definirlo con una precisión distinta. Si por ejemplo  queremos milisegundos haríamos:
auto seed =
chrono::duration_cast < chrono::milliseconds>(chrono::high_resolution_clock::now().time_since_epoch()).count()


O nanosegundos:
chrono::duration_cast< chrono::nanoseconds> (chrono::high_resolution_clock::now().time_since_epoch()).count()

Por supuesto, hay mejores (y más complicadas) formas de generar la semilla. Algunas sugerencias pueden consultarse aquí: http://stackoverflow.com/questions/24334012/best-way-to-seed-mt19937-64-for-monte-carlo-simulations.

Distribuciones de probabilidad en C++11

Existen varios tipos genéricos de distribuciones tanto discretas como continuas. Algunos ejemplos son: uniforme (incluye uniforme entera y real), Bernoulli (incluye Bernoulli, binomial, binomial negativa y geométrica), Poisson (incluye Poisson, exponencial y gamma entre otras) y normal (incluye normal, log-normal y chi-cuadrado entre otras).

Hay detalles a los que vale la pena estar atento, como por ejemplo el rango de los valores generado. Así, la uniforme entera genera valores en el intervalo cerrado [a,b] mientras que la uniforme real lo hace en el semicerrado [a,b).

Algunos ejemplos de generación y uso de números aleatorios

Lo primero que hacemos al comienzo de nuestro programa es definir el mecanismo generador deseado, la semilla y las distribuciones que necesitamos:

mt19937_64 gener; // este es el mecanismo generador que usarán las distribuciones
// la semilla:
auto seed = chrono::high_resolution_clock::now().time_since_epoch().count();
gener.seed(seed); // inicializa el generador con la semilla
// uniforme entre 0 (incluido) y 1 (no incluido):
uniform_real_distribution
< double> uniform(0.0, 1.0);
// la normal de media m y desviación típica s:
double m(0.0), s(1.0);
normal_distribution
< double> normal(m,s);
// Poisson de media avemut=5.2
double avemut(5.2);
poisson_distribution< int> poisson(avemut);


Después se utilizarán unas u otras según sea necesario.

Ensayos Bernoulli: supervivencia y  fertilidad en una población biológica

Un ensayo Bernoulli con probabilidad p es un experimento con sólo dos resultados posibles donde la probabilidad de uno de ellos es p y la del otro 1 – p. Podemos imaginar muchos procesos biológicos como ensayos Bernoulli. Por ejemplo, la probabilidad de sobrevivir dado un valor de eficacia w o la probabilidad de ser progenitor dado un valor de fertilidad f etc etc. Consideremos la eficacia de un individuo en su aspecto de viabilidad es decir, la probabilidad de sobrevivir hasta alcanzar la madurez sexual. Entonces, la supervivencia de un individuo X cualquiera con valor de eficacia w se puede modelar como una Bernoulli. La probabilidad de que ese individuo sobreviva vendrá dada por

if(X.w > uniform(gener)) ... ; // sobrevive

Nótese que si X.w = 1 el individuo siempre sobrevive (recuerda que el intervalo de la uniforme real era semicerrado) pero si X.w es 0 nunca lo hará.

Inicializar los individuos de una población con un fenotipo normalmente distribuido

Esta es una asunción muy común en biología. Una población P en la que estamos considerando un carácter cualquiera: color, tamaño, altura, etc... se distribuye de manera normal para ese carácter. Si queremos generar una población inicial de media m y desviación s hacemos:

for(auto indiv=0; indiv< P.size(); ++indiv)
    P[indiv].caracter= normal(gener);


Generar los individuos que van a mutar en una población

Dado que la ocurrencia de una mutación es un suceso de muy baja probabilidad se suele modelar mediante la distribución de Poisson. Si la probabilidad de mutación en un genoma haploide es 0.01 y tenemos N individuos diploides entonces esperamos 0.02N nuevas mutaciones en cada generación. Podemos forzar la ocurrencia de ese número exacto de nuevos mutantes en cada generación o, si queremos ser más realistas, dejar que varíe según una Poisson definida con el parámetro 0.02N. De este último modo el número de mutantes oscilará en cada generación aunque el valor medio a lo largo de las generaciones tenderá a 0.02N.

mutantes=poisson(gener);

Estos son sólo algunos ejemplos bastante sencillos de entre la multitud de posibilidades que los generadores aleatorios nos ofrecen para simular procesos estocásticos en biología.

Un último apunte respecto a la semilla es que, si estamos realizando varias réplicas de un proceso biológico relativamente largo, puede ser conveniente iniciar cada réplica con una nueva semilla de modo que actualizamos el generador:

// actualiza el generador si se quieren realizar nuevas simulaciones independientes
gen.seed(chrono::high_resolution_clock::now().time_since_epoch().count());


Contrastar la corrección de un generador de números pseudoaleatorios

Es siempre conveniente comprobar que  hemos implementado correctamente nuestros generadores. Para ello podemos comparar valores esperados de cada distribución, como su media y varianza, con los obtenidos en una secuencia dada de números pseudoaleatorios.

La media de una uniforme U(0,1) es 1/2 y la varianza 1/12. Si generamos una cantidad n de valores bajo esta distribución, la media y varianza de esos valores deberían ir aproximándose a los esperados según aumenta n.

Lo mismo podemos hacer con otras distribuciones. Por ejemplo si un individuo tiene viabilidad 0.7 esa es su probabilidad de supervivencia. Si repetimos el proceso n  veces esperamos que el número de veces que sobreviva se aproxime a 0.7n y la varianza de los valores de sobrevivir (1) o morir (0) se aproximará a 0.21n. Igualmente, si inicializamos una población de fenotipos utilizando una normal de media m y desviación s, el fenotipo medio y la desviación en la población obtenida deberían aproximarse a m y s según aumenta el tamaño de la población.

Para una información más detallada sobre estos y otros tipos de test, puede consultarse http://www.euclides.dia.uned.es/aurquia/Files/Simulacion%20Teoria%202008_09.pdf

Bibliografia


sábado, 14 de marzo de 2015

Una bibliografía C++ comentada para biólogos (computacionales)

Libros de C++ y recursos online hay infinitos. Lo que en este post traslado es sólo mi experiencia sobre una serie muy limitada de recursos bibliográficos y/o electrónicos que a mí personalmente me han sido, y siguen siendo, de gran utilidad en mi desempeño profesional. Algunas de mis recomendaciones pueden ser discutibles dado que a día de hoy cualquier persona interesada dispone de ingente cantidad de información vía blogs, manuales online, foros, etc., etc. Pero, ya digo, es mi experiencia y, en mi opinión, algunos buenos libros, cuando se hace el esfuerzo de consultarlos con detenimiento, siguen conteniendo información que a menudo es difícil, sino prácticamente imposible, encontrar online.

Antes de mencionar cualquier material o recurso vale la pena dedicar unas líneas al entorno de desarrollo a utilizar (IDE, por sus siglas en inglés). En mi caso, por diversas circunstancias, he tenido que programar bajo sistemas Windows, Mac y Linux. Por ello se que, un paso importante a la hora de disponer de un método de trabajo eficiente es encontrar un entorno que nos permita generar no sólo código portable sino proyectos de desarrollo intercambiables entre los diferentes sistemas operativos. Hoy día hay herramientas gratuitas y muy potentes, como pueden ser NetBeans y Eclipse, pero personalmente mis preferencias se inclinan, en el caso de C++, por Code::Blocks el cual puede integrarse con diferentes compiladores y que se distribuye gratuitamente bajo licencia GPL v3.0 l. Lo podemos bajar con el compilador GCC ya incluido por defecto. Esta última opción en el caso de Windows implica la instalación (la hace Code::Blocks de manera automática) de MinGW (Minimalist GNU for Windows) es decir, el sistema GNU (sistema operativo tipo-Unix pero libre) mínimo para poder trabajar bajo Windows. MinGW incorpora la colección de compiladores GNU que son más conocidos bajo el acrónimo GCC (GNU Compiler Collection) incluyendo los compiladores de C, C++, ADA y Fortran. Utilizando Code::Blocks, el que esto escribe, ha compilado y ejecutado sin problemas y sin cambiar ni un punto y coma, el mismo proyecto bajo Windows, Linux y Mac. Es cierto que para proyectos de gran complejidad y posiblemente enfocados a sistemas operativos concretos, puede haber otros IDEs más completos y potentes. Sin embargo, para el desarrollo de proyectos de mediana complejidad en biología computacional me parece que Code::Blocks supone un compromiso más que correcto entre facilidad de configuración y optimización, manejo y eficiencia.

Recursos online

Hay infinidad de manuales disponibles gratuitamente online. Para comprobarlo basta teclear en vuestro buscador favorito (no os voy a engañar, el mío es google): "c++ pdf" y obtendréis inmediatamente un listado de manuales en el formato indicado que podréis bajar sin ningún problema a vuestro ordenador. Pero aparte de todos esos pdfs, si disponemos de una conexión a internet, podemos acceder una serie de páginas que dan información y ejemplos rápidamente.

1.- Aquí se puede ver por encima la estructura de la Biblioteca Estándar de C++

http://es.wikipedia.org/wiki/Biblioteca_est%C3%A1ndar_de_C%2B%2B

2.- La siguiente NO es la mejor referencia de C++

http://www.cplusplus.com/reference/
Pues ha sido criticada por contener algunas inexactitudes conceptuales o código no eficiente en algunos de sus ejemplos pero, sin embargo, es una página útil para consulta rápida de métodos asociados a contenedores como vector o métodos de string o algoritmos. Es ligera y carga relativamente rápido. Útil si deseamos obtener información concreta sobre algunas clases de la biblioteca. En mi caso suelo hacer consultas rápidas sobre, por ejemplo, vector:
http://www.cplusplus.com/reference/vector/vector/

3.- El mítico stackoverflow

http://stackoverflow.com/questions/tagged/c%2B%2B
En realidad cuando busco una forma eficiente, o simplemente una forma, de hacer algo, lo que hago es "googlear" el término de búsqueda junto al término c++ y de la lista de resultados consulto en primer lugar los de stackoverflow. Por ejemplo si googleamos: median c++
http://stackoverflow.com/questions/2114797/compute-median-of-values-stored-in-vector-c

4.- El curso de C++ de Zator Systems

http://www.zator.com/Cpp/index.htm
Es un curso con algunos muy buenos ejemplos y explicaciones. Su resumen de operadores de C++ está en mi pestaña de favoritos.

Libros

1.- C++ básico. Autora: Carmen Fernández. Ed. Starbook (http://www.starbook.es/libros/C-BASICO/2245/978-84-936896-8-1)

En mi caso este es un descubrimiento reciente y es el libro que recomiendo actualmente para comenzar en C++ (no en C) desde 0. Tiene, a mi modo de ver, dos ventajas importantes: la primera, es que es un manual breve, sencillo y directo (de muy buen precio) y la segunda es que, a diferencia de muchos otros libros y manuales, menciona dos elementos fundamentales de la Librería Estándar de C++, como son las clases string y vector. Ambas clases facilitan inmensamente la vida del programador con respecto a las herramientas de funcionalidad similar disponibles en C. Otro aspecto interesante de este libro es que los ejemplos del mismo han sido desarrollados utilizando Code::Blocks de modo que dedica un par de páginas a presentar el IDE, explicar cómo se crea un proyecto etc.

2.- Programación en C/C++. Guias prácticas de Anaya multimedia

¿Quién no ha consultado alguna vez en sus inicios alguna de estas guías? Suelen ser manuales de bajo coste, cómodos y útiles para consultas rápidas.

3.- Enciclopedia del lenguaje C++ (C++ estándar). Autor: Fco. Javier Ceballos. Ed. Ra-Ma

Este autor es muy prolífico y tiene varios libros sobre C++ y en otros lenguajes. A mí me gusta este en particular. Es un libro muy completo (y "gordo" 1000 páginas), muy bien explicado y con buenos ejemplos.

4-. Numerical recipes in C++. The Art of Scientific Computing. Autores: Press et al. Ed. Cambridge University Press

Que yo sepa, no hay traducción al español. La edición en inglés que tengo es la del 2002 en su reimpresión de 2005. En cualquier caso, lo digo así de claro, este libro es una de mis biblias. Ello, a pesar de ser consciente de que ha sido bastante criticado, tanto por su pobre diseño desde el punto de vista de la programación orientada a objetos (es poco más que una traducción de C a C++ aprovechando las ventajas de la librería estándar), como desde el campo de la computación numérica, por presentar varios algoritmos que distan de ser los más eficientes en el estado actual del arte de algunos de los muchos campos que toca. Sin embargo, este es un riesgo a asumir, al menos en mi opinión, ante una obra enciclopédica que, en apenas 1000 páginas, introduce conceptos, explica brevemente los métodos y muestra el código completo en C++ de algoritmos que abarcan cuestiones tan dispares, a la par que útiles, como pueden ser la descomposición de Cholesky de una matriz cuadrada (no os riais que sí que sí que es muy útil, cualquiera que haya tenido que generar series de variables correlacionadas lo sabe), la generación de números aleatorios, ordenación, autocorrelación, momentos de una distribución, cálculo de probabilidades de una Chi cuadrado, la tansformada de Fourier, etc, etc. Es un libro que siempre tengo a mano, a partir del cual he programado infinidad de funciones para calcular probabilidades, medias, varianzas, covarianzas... En otras ocasiones me ha servido como una guía de información rápida sobre determinados conceptos y métodos. Una obra de consulta fundamental.

5.- The C++ Standard Library. A Tutorial and Reference. Autor: Nicolai Josuttis. Ed. Addison-Wesley 2012

Este no es un libro para aprender C++ ni algoritmos numéricos. Como su nombre indica es un tutorial y una referencia sobre la Librería Estándar de C++. Personalmente, la mayoría de mis consultas sobre la Librería Estándar las hago utilizando los recursos online, pero este libro incluye algunas explicaciones sobre determinadas características y ventajas de la librería que es difícil (aunque supongo que no imposible) encontrar online. Es un libro claro y profundo, con explicaciones muy buenas, por ejemplo, sobre las ventajas e inconvenientes de los distintos tipos de contenedores (arrays, deques, lists, vectors...), las nuevas características del nuevo estándar C++11 (incluyendo la nueva y magnífica librería de números aleatorios), la biblioteca de algoritmos, etc. Una obra magistral indicada para aquellos que ya saben programar en C++ y desean sacarle todo el partido a la librería estándar.

Y por ahora esto es todo. Hay, por supuesto, muchos más libros y recursos. Pero, como ya dije, estos me parecen útiles para biólogos computacionales o al menos lo han sido en mi caso, a la hora de programar modelos de simulación, algoritmos genéticos y métodos de contraste estadístico en C++. Otros libros de referencia, más avanzados, que quedan lejos (o muy lejos) de la intención de este post y que por tanto sólo me limito a mencionar serían:

Programming. Principles and Practice Using C++ (Stroustrup)
C++ Templates. The Complete Guide (Vandermonde & Josuttis)
Exceptional C++ (Herb Sutter)
Modern C++ Design (Alexandrescu) (muy lejos)

domingo, 22 de febrero de 2015

Paso de argumentos por valor y por referencia. Aspectos de interés en simulación biológica.

En C++ cuando pasamos los argumentos en un procedimiento o función decimos que el paso puede ser por valor o por referencia. Antes de continuar vamos a clarificar nuestra terminología recordando que en la definición de una función hablamos de los parámetros que recibe la función. Es decir si defino la función suma en C++

void suma(int x, int y, int z){
    z=x+y;
}


entonces x, y, z son los parámetros o argumentos formales de la función en los que además se ha indicado el tipo de dato al que corresponden, que en este caso es un entero (int). Cuando dentro de cualquier programa utilizo la función suma pasándole como argumentos distintos variables o valores como por ejemplo

suma(a,b,c); suma(2,3,c);

tanto a,b,c como 2 y 3 son argumentos de la función que se corresponderán dentro del cuerpo de la función a sus parámetros según el orden indicado. Por tanto a se corresponderá con xb con y, c con z.

Paso por valor

El paso por valor es el mecanismo por omisión en C++. En el paso por valor, las variables que se pasan en una llamada a la función son copiadas a nuevas variables que se corresponden a los parámetros formales y que existirán durante el tiempo que se ejecuta la función. Por ejemplo en el siguiente código

http://pastebin.com/1z8AJvz9

en la línea 6 la variable c tiene valor 7. Cuando se llama a la función suma pasándole como tercer argumento la variable c, lo que ocurre es que la función suma copia el valor de c en una nueva variable que se corresponde con el parámetro z. Todo lo que se le haga a z le ocurre a la copia de c no a c. Lo mismo podemos decir de las variables a y b. Por tanto al ejecutarse suma(a,b,c) la copia de c que llamamos z pasa a ser la suma de las copias de a y b es decir z = x+y = 0+0 = 0. Pero cuando salimos de la función la copia se destruye y la variable c sigue teniendo su valor 7. Debido a que en la llamada suma(a,b,c) no se pasan las variables a, b y c sino una copia de su valor, hablamos de paso por valor.

Paso por referencia

En este caso lo que se le pasa a la función no es una copia del valor de la variable sino una refeencia a esa misma variable. Decimos que el parámetro se convierte en un alias para el argumento, es decir, es la misma variable pero con otro nombre. Como ahora el valor que toman los parámetros de la función no es una copia sino que corresponden a las variables externas a la función los cambios que ocurran dentro de la función tendrán alcance fuera de ella (algo que no ocurría en el paso por valor). Siguiendo con el ejemplo anterior si cambiamos la manera de expresar los parámetros en la definición de suma de modo que en vez de

void suma(int x, int y, int z) 

escribimos

void suma(int x, int y, int& z) 

lo que le estamos indicando al compilador de C++ es que el parámetro z se pase por referencia. Si volvemos al ejemplo anterior, al llamar a suma con el argumento c ya no se haría una copia de c sino que sería la misma variable c que dentro dela función tendría un segundo nombre (como un apodo o alias) llamado z y por tanto al tomar dentro de la función el valor z = 0+0=0 al salir de la función c valdría 0.

Biología

Veamos ahora un detalle importante relacionado con la eficiencia computacional en relación al paso por valor o por referencia. Imaginemos que tenemos que modelar una población biológica y deseamos programar una simulación donde los individuos se reproducen de modo que una pareja de padres tienen un hijo. Si quisiéramos programar una función reproducción podríamos a priori pensar en algo como

void reproduce(Individuo padre, Individuo madre, Individuo& hijo){..código de la función..}

Donde el tipo de los parámetros es una clase llamada Individuo. Además hemos sido hábiles y nos dimos cuenta de que podemos pasar el padre y la madre por valor puesto que no necesitan cambiar al pasar por la función. Sin embargo el hijo lo vamos a calcular dentro y presumiblemente querremos conservar el valor calculado al salir de la función, por tanto hemos decidido, sabiamente, pasar el parámetro hijo por referencia. Esto hecho así, es correcto y posiblemente obtengamos el comportamiento deseado.

Hay sin embargo un posible problema relacionado con la eficiencia. Imaginemos que nuestra clase Individuo es una abstracción de un organismo diploide, es decir con dos series de cromosomas y que además de los cromosomas le queremos asignar un valor constante para el carácter que sea (por ejemplo color) y otros valores variables como altura, edad etc. Respecto a los cromosomas, vamos a considerar 10 pares de cromosomas cada uno con un millón de posiciones que pueden tener valor 0 o 1. Es decir, el Individuo padre consta de 20 vectores de longitud un millón de posiciones cada uno además de una constante y algunas variables. Si recordamos lo anteriormente dicho para el paso por valor, cuando llamemos a la función reproduce el primer argumento de la llamada se copiará en el primer parámetro (padre). Esta copia exigirá copiar cada una de las millones de posiciones de los vectores así como las constantes, variables y cualquier procedimiento que contengan los objetos de la clase Individuo. Total para usar, sin cambiarlos, unos valores que ya existían fuera de la función. Incurrimos por tanto en un gasto innecesario de memoria y tiempo de computación. En modelos de simulación con miles de individuos y procesos que se repiten durante muchos miles de generaciones el coste puede ser prohibitivo. ¿Cuál es la solución? Pues obviamente es el paso por referencia. La definición adecuada de la función sería

void reproduce(const Individuo& padre, const Individuo& madre, Individuo& hijo){..código de la función..}

Donde todos los argumentos se pasan por referencia de modo que no se realiza la copia de los millones de posiciones cromosómicas sino que simplemente se pasa una referencia al comienzo de la dirección del objeto Individuo correspondiente. El modificador const delante de padre y madre indica que estos parámetros no se modifican dentro de la función e impiden que se haga por error puesto que el compilador se quejará si intentamos cambiarlo en el código.

Otra alternativa que quizás genera un código más claro sería no pasar ningún argumento que se modifique y hacer que la función devuelva un nuevo individuo que será el hijo:

Individuo reproduce(const Individuo& padre, const Individuo& madre){ Indiv hijo; /* codigo....*/ return hijo; }

Sin embargo esta última opción no será posible si es necesario cambiar más de un objeto (si queremos por ejemplo generar dos hijos de cada pareja de padres) y tampoco es lo más razonable cuando queremos modificar vectores. En estos casos el paso por referencia no constante es más conveniente.

Punteros (paso por dirección)

Una última opción para el paso de argumentos en C++ es que el parámetro de la función sea un puntero. En este caso estamos pasando la dirección de la variable y cualquier cambio en el contenido de esa dirección afectará al valor de la variable. En el paso por dirección obtenemos por tanto un comportamiento igual al del paso por referencia. Pero ahora el nombre del parámetro ya no es un alias de la variable sino que corresponde a la dirección de la variable y por tanto hay que desreferenciarla para acceder al contenido de esa dirección (usando el operador *). Dejamos aquí está discusión pues queda más allá de la intención del post. Para saber más sobre referencias y punteros y cuando sería más conveniente el uso de unos u otros se puede consultar el libro de Stroustrup citado en la bibliografía. Personalmente suelo usar siempre referencias puesto que simplifica el código de la función ya que no es necesario desreferenciar la variable (una referencia equivale a un puntero que se desreferencia automáticamente). El paso de punteros puede generar problemas adicionales si ocurren llamadas a la función usando punteros nulos.

Bibliografía

Lenguajes de programación. Principios y práctica. K.C. Louden.
Programming: Principles and Practice Using C++. B. Stroustrup
Object-Oriented Programming in C++. R. Lafore.
Effective C++. S. Meyers.
Hola. Mi nombre es Antonio Carvajal-Rodríguez y soy profesor en la Universidad de Vigo. Me dedico a la docencia e investigación en Genética y Biología Computacional, lo cual es una suerte porque considero que ambas materias son apasionantes. Concretamente, la biología computacional es un campo de estudio que considero especial. Por un lado es un ámbito multidisciplinar donde entran en juego conocimientos de biología, estadística-matemática y programación. Por otro lado es una carrera con una gran demanda actual, y previsiblemente futura, de profesionales. Sin embargo, y aunque parezca extraño, existe en España muy poca gente formada o formándose en biología computacional, bioinformática o disciplinas afines. Las causas de esta situación no las voy a discutir aquí aunque podrían ser el tema de una futura entrada en este blog. La idea del mismo nació ya a finales de 2008 pero diversas vicisitudes personales y profesionales la fueron diluyendo. Mi intención es retomar ahora la idea original e ir compartiendo inquietudes, conocimientos e ideas que han ido surgiendo a lo largo de 15 años de programación de modelos y simulaciones en biología, especialmente en genética evolutiva. Es también una suerte y un privilegio el haberme formado en dos campos a priori bien diferentes como son la biología y la informática. Quizá por eso trato de enfocar mi trabajo siempre desde una doble perspectiva: por un lado, el problema biológico que subyace en un programa o modelo determinado y por otro, el lenguaje, la eficiencia y la técnica de programación utilizada. Pero como dice el refrán "quién mucho abarca, poco aprieta" y seguramente muchos programadores podrán plantear soluciones más eficientes a cuestiones de código que aquí se expongan. Y lo mismo ocurrirá desde el ámbito de la biología. Así que ya sabéis, toda crítica constructiva será siempre muy bien recibida.

sábado, 29 de agosto de 2009

Algoritmos genéticos y simulación en genética de poblaciones. Dinámica de sistemas multiloci

En genética de poblaciones modelar la dinámica, esto es la evolución a lo largo del tiempo, de un cierto número de genes en presencia de selección, mutación y recombinación exige usar el método de simulación de Montecarlo. Nótese que si no hay recombinación los genes se comportan como uno sólo y estaríamos hablando de un modelo de un sólo gen (o locus si nos referimos al lugar que ocupa el gen en el genoma). Este es el único caso en que la dinámica del sistema, el conocimiento de puntos de equilibrio y su estabilidad, es perfectamente conocida.

Sin embargo, en cuanto consideramos más de un gen, incluso en el caso de sólo dos, la situación se complica. Se podría pensar que a estas alturas el caso de un mísero par de genes bialélicos debería estar perfectamente resuelto para cualquier modelo de eficacia biológica o fitness. Lo cierto es que a pesar de que sí ha sido muy estudiado, sobre todo para casos especiales de la función de “fitness” (véase una interesante revisión en [1]), dista de estar completamente resuelto incluso en el caso sin mutación pudiendo aparecer dinámicas de comportamiento bastante complejo [2]. Que no decir si hablamos de tres o más genes [3, 4]. En cuanto a la dinámica de estos sistemas multiloci podemos considerar la selección a nivel haploide (sobre los gametos) o diploide (sobre los zigotos) teniendo en cuenta que la correspondencia entre unos modelos y otros dependerá del esquema de fitness asumido [5]. En cualquier caso y debido a la dificultad para analizar estos modelos multigénicos o multiloci, surge el uso de la herramienta de simulación la cual permite estudiar la evolución de los susodichos sistemas bajo distintas condiciones de mutación, selección, recombinación, tamaño poblacional etc.

El algoritmo de simulación clásico, que implementa el modelo de reproducción de Wright-Fisher y sus diferentes variantes incluyendo selección, mutación y recombinación es similar a un algoritmo genético (AG) en su versión básica, el conocido SGA (Simple Genetic Algorithm), con la peculiaridad de que los genes no son resultado de una codificación para representar un problema dado sino que podríamos decir que se representan a sí mismos o más precisamente a un fenotipo P = f(G,E) donde G es el genoma o conjunto de genes y f es la función, incluyendo la identidad, que mapea el fenotipo a partir de un genotipo y ambiente E dados. La fitness será a su vez una función sobre P. Es decir, el modelo de simulación en genética de poblaciones es formalmente un SGA pero que evita el problema de la representación asociado al caso de codificación binaria porque no es necesario inventarse una codificación genética para los puntos del espacio de soluciones de un problema dado pues los individuos bajo estudio ya incorporan su propio código genético. A parte de esto, lo que diferencia a ambos enfoques es simplemente la finalidad con la que utilizamos el algoritmo. En el caso del SGA nosotros dirigimos el proceso evolutivo mediante la función de fitness de modo que esta es máxima en la dirección en que se resuelve nuestro problema ya sea este de maximización, minimización, estimación, control, etc. Esto es, ignoramos el valor exacto de la solución pero conocemos la dirección y el sentido del vector solución y definimos la función de fitness en consecuencia. Lo que esperamos del algoritmo es que el proceso evolutivo vía selección, mutación y recombinación baraje la combinación adecuada de genes que nos lleve al punto o región deseada. Esa combinación genética al decodificarla nos dará la solución que buscamos. En genética de poblaciones el tipo de problema es bien distinto puesto que de entrada no dirigimos el proceso evolutivo si no más bien deseamos comprender ese mismo proceso y las diferentes relaciones entre los parámetros implicados y de estos con la función de fitness, el ambiente e incluso determinados factores estocásticos generados por el tamaño poblacional finito. En este sentido los algoritmos genéticos (AGs) en general son interesantes porque se caracterizan por manejar tamaños de población muy pequeños y con un componente selectivo grande. Desde el punto de vista evolutivo este es un terreno adecuado para el efecto de Hill-Robertson [6, 7] el cual predice que tenderán a aparecer combinaciones genéticas intermedias es decir gametos con alelos favorables en unos loci y desfavorables en otros. El que, de un modo u otro, cierto tipo de combinaciones se vean favorecidas genera lo que se conoce como desequilibrio de ligamiento, en este caso negativo (al favorecerse combinaciones intermedias) lo que a su vez provoca que la recombinación tenga un efecto beneficioso sobre la fitness poblacional.

Dicho lo anterior, podemos esperar que en los algoritmos genéticos el efecto de la recombinación sea especialmente importante. En cuanto a la selección, el modelo clásico en el SGA es el de selección proporcional a la fitness que garantiza que, en promedio, los individuos serán seleccionados en proporción a su fitness de modo que cuanto mayor es esta, mayor es la opción de reproducirse. El problema aquí es lo comentado sobre los tamaños poblacionales pequeños, debido a la elevada estocasticidad (deriva genética en la terminología de la genética de poblaciones) que ello genera, se pierde la garantía de que las mejores combinaciones sean seleccionadas. Por ello se utilizan a menudo AGs con otros esquemas de selección alejados del modelo biológico como la selección elitista (nos quedamos con los mejores) o por torneo (combina elitista y proporcional). El problema de estos modelos es la pérdida de variabilidad que subyace en sus estrategias; obviamente si encontrar la solución fuera tan simple como ir siguiendo las variantes óptimas de cada generación apenas necesitaríamos usar la estrategia evolutiva. En realidad la clave en el uso de la estrategia de selección proporcional a la fitness sería una adecuada combinación de los parámetros de intensidad de selección, tamaño poblacional, mutación y recombinación. Justamente, comprender la relación entre estos parámetros es el objetivo de la genética de poblaciones teórica. En el caso específico de los AGs y dado su bajo tamaño poblacional, el estudio de las condiciones del efecto de H-R y la recombinación deberían ser los aspectos principales a tener en cuenta. En este sentido, aunque desde una perspectiva puramente heurística, se han desarrollado algoritmos meta-genéticos que autoajustan los propios parámetros de búsqueda (mutación y recombinación), así como metapoblaciones de AGs realizando búsquedas independientes que se comunican vía migración esporádica de soluciones de una población a otra.

Otro aspecto que se plantea desde el punto de vista de la teoría de control y optimización, donde los AGs inicialmente surgen [8], es porqué y cómo funcionan los algoritmos genéticos. El porqué es, al menos desde el punto de vista de la genética de poblaciones, obvio. Funcionan porque las combinaciones peores se reproducen menos y las mejores más por tanto, dadas unas fuentes de variabilidad primaria, la mutación y secundaria, la recombinación, con el paso del tiempo obtenemos un conjunto de soluciones óptimas o cuasi-óptimas. Respecto al cómo, se han hecho incursiones teóricas en la dinámica de los algoritmos usando herramientas de matemática avanzada como la mecánica estadística de la física cuántica [9, 10] sin embargo, tal y como hemos comentado al principio, si ya la dinámica de unos pocos genes es compleja, difícil es obtener soluciones explícitas para el caso de muchos genes. Con todo, es interesante percatarse de que mientras la genética de poblaciones se ha centrado siempre en descripciones de la evolución de la frecuencia de las distintas combinaciones genéticas y la conexión de estas con los parámetros del modelo y con los correspondientes valores de fitness, en el ámbito de los AGs el enfoque ha sido, quizás no sorprendentemente, mucho más mecanicista. Es decir, la explicación del camino a la solución o a la fitness final se da en términos de la construcción de los bloques genéticos adecuados y no del contexto que lleva a ello. Una de las hipótesis más significativas, aunque todavía sometida a debate, sobre el funcionamiento de los AGs es el teorema de los esquemas que viene a decir que la solución se obtiene construyendo con pequeños bloques o esquemas formados por grupos de unos pocos buenos genes no muy distantes unos de otros. Lo de pocos sería para que la probabilidad de que la mutación altere los genes buenos sea baja y lo de próximos para que la recombinación no rompa esas combinaciones.

La diferencia entre ambos planteamientos, genética de poblaciones versus AGs, es fiel reflejo de las distintas finalidades en ambos campos, como ya se dijo, una se centra en el propio proceso evolutivo, el otro quiere llegar a un objetivo sin que importe el contexto en que ello ocurre. Sin embargo, la cuestión sobre cómo trabajan los AGs para encontrar la solución a un problema, es la misma que se mencionó al comienzo, esto es, conocer las pautas de comportamiento a lo largo del tiempo de un conjunto de genes, dados unos valores de mutación, recombinación y selección con la dificultad añadida de que el tamaño poblacional es apreciablemente finito lo que incorpora una importante componente estocástica al problema.

Es pues este un ejemplo de cómo un modelo que es sencillo de simular en el ordenador, se muestra tremendamente complejo a la hora de analizar su dinámica. Lo interesante es percibir cómo dos ámbitos absolutamente separados, a saber, la genética de poblaciones en biología y las áreas relacionadas con problemas de optimización, automatización, y control en ingeniería, han tratado de teorizar sobre esta dinámica realizando, curiosamente, muy poco trabajo interdisciplinar. En mi opinión, todos tienen mucho que ganar aquí: los teóricos de los AGs con el enfoque de la genética de poblaciones, pues el contexto se mostrará clave en cuanto a la resolución eficiente de problemas difíciles. Así mismo, los biólogos tienen en el campo de los algoritmos genéticos, por un lado una aplicación muy interesante de los modelos evolutivos y por otro la posibilidad de combinar sus conocimientos del contexto, las interacciones entre selección, deriva, mutación, recombinación etc con el nuevo enfoque donde se trata de ver el modo en que los diferentes bloques genotípicos son parcialmente construidos mediante esquemas.

Bibliografía

1. Li CC: Genetic equilibrium under selection. Biometrics 1967, 23(3):397-484.

2. Bürger R: The Mathematical theory of selection, recombination, and mutation. Chichester: John Wiley; 2000.

3. Robinson WP, Asmussen MA, Thomson G: Three-locus systems impose additional constraints on pairwise disequilibria. Genetics 1991, 129(3):925-930.

4. Baake M, Baake E: Erratum to: An exactly solved model for mutation, recombination and selection. Canadian Journal of Mathematics 2008, 60(2):264-265.

5. Kirzhner V, Lyubich Y: Multilocus dynamics under haploid selection. J Math Biol 1997, 35(4):391-408.

6. Hill WG, Robertson A: The effect of linkage on limits to artificial selection. Genet Res 1966, 8(3):269-294.

7. Zeng Z-B: The Hill-Robertson effect is a consequence of interplay between linkage, selection and drift: a commentary on The effect of linkage on limits to artificial selection by W. G. Hill and A. Robertson. Genet Res 2008, 89(Special Issue 5-6):309-310.

8. Goldberg DE: Genetic algorithms in search, optimization, and machine learning. Reading, Mass.: Addison-Wesley; 1989.

9. Prugel-Bennett A, Shapiro J: The dynamics of a genetic algorithm for simple random Ising systems. Phys D 1997, 104(1):75-114.

10. Shapiro J, Prugel Bennett A, Rattray M: A Statistical Mechanical Formulation of the Dynamics of Genetic Algorithms. In: Selected Papers from AISB Workshop on Evolutionary Computing. Springer-Verlag; 1994.