25 de abril de 2012

StatSeq WorkShop 4 - Verona

La pasada semana, a mediadios de abril, tuvo lugar en Verona, en el norte de Italia, el cuarto workshop de StatSeq. Esta es una Action COST del campo FA (Food and Agriculture) para la coordinación de esfuerzos enfocados al análisis de datos de secuenciación de plantas. La sede del workshop era el auditorio Polo Zanotto, junto a uno de los márgenes del río Adigio, que da forma a Verona.


Tras el percance logístico de Michele Morgante, abrió la sesión Alberto Ferrarini, hablando sobre el ensamblado de novo de un cultivar de Vitis vinifera. Para la tarea obtuvieron de la planta un total de 45 muestras, representando 16 tejidos y varios estados de desarrollo. Por lo visto, obtuvieron bastantes intrones cuando mapearon los resultados de la secuenciación de RNA con la referencia existente en esa especie. Entre los posters presentados al evento se encontraban otros trabajos relacionados con el de Ferrarini y el grupo de Delledonne. Andrea Acquaviva presentaba un framework para caracterización de transcriptomas de cultivares distintos de la referencia genómica. Michele Perazzolli mostró un experimento de expresión para estudio de ISR (Induced Systemic Resistance) y Plasmopara viticola en vid.


Otros experimentos de expresión se presentaron entre los posters, como la comparación de expresión en corona y hojas, en respuesta al frío, en cebada, por Jaroslava Ovesná; y una comparativa de métodos de normalización, de Elie Maza. Entre las muy diversas herramientas y paquetes que se presentaron citamos el workflow engine Conveyor, presentado por Berkhard Linke; MotifLab, para análisis exhaustivo de secuencias reguladoras, con Finn Drablos; NarrowPeaks, un paquete R para análisis de picos de datos de ChIPseq, en póster presentado por Pedro Madrigal; y BioMark, también paquete de R, en éste caso para aplicar métodos que mejoren la selección de variables en problemas p >> n, típico de GWAS (Genome-Wide Assisted Selection), con Ron Wehrens como ponente. Al ser StatSeq una acción enfocada a los métodos estadísticos, el tema de selección de variables acogió también la excelente charla de Patrick Waldmann, en la que mostraba la capacidad del método elastic net para tener en cuenta el LD (Linkage Disequilibrium) entre marcadores cuando se trabaja en GWAS. Además, Willem Kruijer asoló al personal con una tira de ecuaciones que no facilitó la comprensión general del uso de algoritmos secuenciales de Monte Carlo para el mismo asunto.

Enlazando esto con otros asuntos bayesianos, Jimmy Vandel presentó el uso de una red bayesiana para modelizar una red de regulación génica basándose en RILs (Recombinant Inbred Lines) de Arabidopsis thaliana. También Martin-Magniette expuso su preferencia por métodos probabilísticos en la aplicación de técnicas de clustering para análisis de perfiles de expresión. Una ventaja que explicó radicaría en la posibilidad de que cada gen reporte una probabilidad de pertenecer a cada cluster, en lugar de la pertenencia de todo o nada típica de los métodos como k-means. Además, informó de que obtuvieron mejores resultados con unos algoritmos que otros, a la hora de estimar los parámetros y al determinar el número de clusters, siendo EM > CEM e ICL > BIC, respectivamente. Por otro lado, como método de validación, Micha Bayer, junto a David Marshall, presentaba un póster donde se trataban las deficiencias del uso de N50 para la calificación de ensamblajes, y proponen el mapeo de flcDNAs a los contigs para comprobar la integración del gene space, sin duda mucho más relevante biológicamente hablando. Más tímidamente, Julie Aubert presentó una comparativa de métodos de normalización para RNAseq.

Una de las charlas que más pareció gustar fue la de Jonathan Marchini, sobre estimación de haplotipos mediante SHAPEIT. Quizás lo más llamativo era la heurística del algoritmo y su buena escalabilidad.


También sobre haplotipos trató la charla de Jan de Boer, el sorprendente ponente de Wageningen. Utilizaron captura de secuencias con sondas de 120 bp con un overlap de 20 bp, para aproximadamente 800 genes de tomate tetraploide. Luego siguen un pipeline de GBS (Genotyping By Sequencing), a partir de reads 2x100 de Illumina. Si en esta ponencia los resultados resultaron nebulosos, más transparente y desafortunado pareció Thomas Odong en su análisis de poblaciones naturales de Arabis alpina, un potencial modelo para plantas perennes. Interesantes charlas sobre genotipado fueron también las de Jeff Glaubitz y Jaap Buntjer. El último proponía un nuevo método de mejora que hibrida ideas de MAS (Marked-Assisted Selection) y GS (Genomic Selection), denominada Genomic Breeding. Glaubitz, por su parte, presentó el pipeline de GBS en maíz que utilizan en la Cornell University. Usando captura, esta vez para análisis de CNV (Copy Number Variation), Guillem Rigaill propuso modelar la cobertura de la muestra como una función lineal de los controles, en lugar de clásico uso de logratios, que desprecian la cobertura total en un locus, llevando a la pérdida de información.

Exposiciones de proyectos de gran envergadura fueron la de Dan Bolser, sobre TransPLANT; Mark A. De Pristo, 1000 Genomes; y el poster informando sobre el estado actual de MELONOMICS, de Walter Sanseverino.

Sin duda, destacados en la conferencia fueron Michele Morgante y Lauren M. McIntyre. Esta contagió su ímpetu y expresividad a la hora de presentar lo que básicamente recoje su artículo "RNAseq: technical variability and sampling" de BMC Genomics. A ver quién se atreve a no utilizar réplicas técnicas y métodos de agreement delante de la amiga de Florida. En cuanto a Morgante, además del retraso, perfectamente entendible, y de la anécdota que protagonizó cuando tuvieron que ir a buscarle porque el sonido del micrófono era para él ruido de fondo desde hacía rato y no se percataba de la llamada a la mesa redonda; hizo una buena presentación de lo que él llama catálogos verticales, sobre LSV (Large Structural Variation), y horizontales (en éste caso metabolismo de lignina en chopo). En cuanto a LSV, utilizan el software BreakDancer para análisis de mapeo paired-end, a la búsqueda de PAV (Presence-Abscence Variation), y DOC (Depth Of Coverage) para CNV. En cuanto a los 5 genes de lignina en chopo, parece que el análisis de las frecuencias alélicas mediante pools de 64 individuos fue suficiente para analizqar 768 árboles con una buena correlación.


Finalmente, quizás la charla más sorpredente fue la de Maria Colomé-Tatché y sus EpiRILs (Epigenetic RILs). Mediante BSseq (BiSulphite sequencing) y MeDIP-chip (Methylated DNA InmunoPrecipitation chip) de una población obtenida cruzando parentales con DNA casi idéntico, pero perfil de metilación muy distinto, obtienen DMRs (Differentially Methylated Regions) en mapas genómicos. Esperan poder aplicar estos marcadores recombinantes robustos, basados en metilación del DNA, para análisis de QTL.


La reunión más informal tuvo como protagonistas el risotto, la pasta, la carne, los postres y los interesantes vinos de la llanura italiana.

19 de abril de 2012

perl + GNU Scientific Library

Hola,
leyendo el Linux Journal me he encontrado con este artículo que explica como usar la GNU Scientific Library en C, y me he acordado de que tenemos por el laboratorio, sin estrenar, la versión 1.12 del manual. Para qué vale esta librería? La respuesta corta es: para hacer cálculos de todo tipo con algoritmos eficientes que minimizan los errores de redondeo de los procesadores digitales.
A veces una librería es la mejor forma de emplear algoritmos extensivamente optimizados y que tal vez no sepamos cómo programar. A día de hoy la SGL incluye funciones en estas áreas:
Complex Numbers   Roots of Polynomials
Special Functions   Vectors and Matrices
Permutations   Sorting
BLAS Support   Linear Algebra
Eigensystems   Fast Fourier Transforms
Quadrature   Random Numbers
Quasi-Random Sequences   Random Distributions
Statistics   Histograms
N-Tuples   Monte Carlo Integration
Simulated Annealing   Differential Equations
Interpolation   Numerical Differentiation
Chebyshev Approximation   Series Acceleration
Discrete Hankel Transforms   Root-Finding
Minimization   Least-Squares Fitting
Physical Constants   IEEE Floating-Point
Discrete Wavelet Transforms   Basis splines
En algunos aspectos puede ser más limitada que las Numerical Recipes, pero al liberarse bajo GPL tenemos mayor libertad de uso, y sin pagar un duro. Además, gracias al trabajo de Jonathan Leto, disponemos en CPAN de una interfaz para Perl, llamada Math::GSL. Algunas de las aplicaciones de la librería se pueden conseguir con otros módulos de CPAN, pero a priori la GSL tiene a su favor: i) menor tiempo de ejecución y ii) mayor precisión numérica, una limitación de Perl ya discutida.

Pasemos a la práctica: cómo se instala esto? En mi Ubuntu 10.4 fue tan sencillo como decir $ sudo cpan -i Math::SGL y luego ir aceptando la instalación de sus dependencias.

Podemos ahora probar un ejemplo donde calcularemos una derivada numérica, mostrando algunas de las capacidades de la librería:


 #!/usr/bin/perl -w  
 # Ejemplo de derivada de ln(x) con GSL  
 # Adaptado de ejemplos de Jonathan Leto en:  
 # https://github.com/leto/math--gsl/blob/master/examples/deriv/basic  
   
 use strict;  
 use Math::GSL::Deriv qw/:all/;  
 use Math::GSL::Errno qw/:all/;  
   
 # incremento de derivada  
 my $h = 0.01;   
   
 # queremos X = pi, ejemplo constante matematica en GSL   
 my $x = $Math::GSL::Deriv::M_PI;   
    
 # derivada numerica, recuerda en perl log = ln  
 my ($status,$val,$err) = gsl_deriv_central ( sub { log($_[0]) }, $x, $h);  
   
 # derivada analitica, demostracion en  
 # http://www.math.com/tables/derivatives/more/es-ln.htm  
 sub dlndx{ return 1/$_[0] }  
   
 if($status == $GSL_SUCCESS)   
 {  
   printf("deriv(ln((%g)) = %.18g, error maximo esperado = %.18g\n", $x, $val, $err);  
   printf("dlndx(%g)   = %.18g\n" , $x, dlndx($x));  
   printf("error observado   = %.18g\n",abs($val-dlndx($x)));  
 }   
 else   
 {  
   my $gsl_error = gsl_strerror($status);  
   print "ERROR: $gsl_error (es derivable log(x) en ese punto?)\n";  
 }  

Se obtiene algo como:

deriv(ln((3.14159)) = 0.31830988618650752, error maximo esperado = 5.7425388097458009e-11
dlndx(3.14159)      = 0.31830988618379069
error observado     = 2.7168267635602206e-12

El autor del módulo tiene en la web una colección de ejemplos que sin duda amplian la descripción de las capacidad de GSL. Por ejemplo, en este ejemplo se explica como calcular productos escalares de dos maneras, encontrando que la implementación del módulo Math::GSL::BLAS es al menos el doble de rápida.

Ya que este ejemplo no representa las tareas típicas de la programación en Bioinformática os muestro un segundo mucho más terrenal, donde ordenamos vectores de números reales, comparando la implementación de mergesort de Perl con la heapsort inestable de GGL. Hay un capítulo del manual dedicado sólo a esto, pero probablemente la documentación del módulo Math::GSL::Sort sea suficiente:

 use Math::GSL::Sort qw/:all/;  
   
 my $numbers = [ map { rand(100) } (1..100000) ];  
 my $sorted = gsl_sort( $numbers,1,scalar(@$numbers) );     
 my @sortedp = sort {$a<=>$b}( @$numbers );  

Si hacéis pruebas veréis que la implementación GSL resulta ser más rápida, en mi caso aproximadamente un 20%.

Otras aplicaciones posibles incluyen el cálculo de permutaciones y combinaciones, algo ya discutido en este blog,
hasta otra,
Bruno

PD Si alguien ha probado a instalar este módulo en sistemas Windows le agradeceremos sus comentarios

13 de abril de 2012

Código genético en formato LaTeX

Seguro que si lees este blog, habrás visto alguna vez en un libro de texto la tabla de codones de RNA que codifican para los diferentes aminoácidos de alguna de las siguientes formas:


Pero seguramente pocos de vosotros os habéis atrevido a reproducir dicha tabla en código LaTeX, es un ejercicio recomendable que permite entrenarse con casi todas las particularidades de las tablas en LaTeX:

 \documentclass[12pt,spanish]{book}  
   
 \usepackage[utf8]{inputenc}  
 \usepackage{multirow}  
 \usepackage{rotating}  
   
 \begin{document}  
             
      \begin{table}[H]  
           \begin{center}  
           \begin{tabular}[t]{cccccccccccc}  
             
                & & \multicolumn{8}{c}{\textbf{\uppercase{Segunda base}}} \\  
                  
                & & \multicolumn{2}{c}{\textbf{U}} & \multicolumn{2}{c}{\textbf{C}} & \multicolumn{2}{c}{\textbf{A}} & \multicolumn{2}{c}{\textbf{G}} & \\ \cline{3-10}  
             
                \multirow{16}{*}{\rotatebox{-90}{\mbox{\textbf{\uppercase{Primera base}}}}} & \multirow{4}{*}{\textbf{U}} & \multicolumn{1}{|c}{UUU} & \multirow{2}{*}{Phe} & \multicolumn{1}{|c}{UCU} & \multirow{4}{*}{Ser} & \multicolumn{1}{|c}{UAU} & \multirow{2}{*}{Tyr} & \multicolumn{1}{|c}{UGU} & \multirow{2}{*}{Cys} & \multicolumn{1}{|c}{\textbf{U}} & \multirow{16}{*}{\rotatebox{90}{\mbox{\textbf{\uppercase{Tercera base}}}}}\\  
                & & \multicolumn{1}{|c}{UUC} & & \multicolumn{1}{|c}{UCC} & & \multicolumn{1}{|c}{UAC} & & \multicolumn{1}{|c}{UGC} & & \multicolumn{1}{|c}{\textbf{C}} \\  
                & & \multicolumn{1}{|c}{UUA} & \multirow{2}{*}{Leu} & \multicolumn{1}{|c}{UCA} & & \multicolumn{1}{|c}{UAA} & \multirow{2}{*}{FIN} & \multicolumn{1}{|c}{UGA} & FIN & \multicolumn{1}{|c}{\textbf{A}} \\  
                & & \multicolumn{1}{|c}{UUG} & & \multicolumn{1}{|c}{UGG} & & \multicolumn{1}{|c}{UAG} & & \multicolumn{1}{|c}{UGG} & Trp & \multicolumn{1}{|c}{\textbf{G}} \\ \cline{3-10}  
             
                & \multirow{4}{*}{\textbf{C}} & \multicolumn{1}{|c}{CUU} & \multirow{4}{*}{Leu} & \multicolumn{1}{|c}{CCU} & \multirow{4}{*}{Pro} & \multicolumn{1}{|c}{CAU} & \multirow{2}{*}{His} & \multicolumn{1}{|c}{CGU} & \multirow{4}{*}{Arg} & \multicolumn{1}{|c}{\textbf{U}} \\  
                & & \multicolumn{1}{|c}{CUC} & & \multicolumn{1}{|c}{CCC} & & \multicolumn{1}{|c}{CAC} & & \multicolumn{1}{|c}{CGC} & & \multicolumn{1}{|c}{\textbf{C}} \\  
                & & \multicolumn{1}{|c}{CUA} & & \multicolumn{1}{|c}{CCA} & & \multicolumn{1}{|c}{CAA} & \multirow{2}{*}{Gln} & \multicolumn{1}{|c}{CGA} & & \multicolumn{1}{|c}{\textbf{A}} \\  
                & & \multicolumn{1}{|c}{CUG} & & \multicolumn{1}{|c}{CGG} & & \multicolumn{1}{|c}{CAG} & & \multicolumn{1}{|c}{CGG} & & \multicolumn{1}{|c}{\textbf{G}} \\ \cline{3-10}  
             
                & \multirow{4}{*}{\textbf{A}} & \multicolumn{1}{|c}{AUU} & \multirow{3}{*}{Ile} & \multicolumn{1}{|c}{ACU} & \multirow{4}{*}{Thr} & \multicolumn{1}{|c}{AAU} & \multirow{2}{*}{Asn} & \multicolumn{1}{|c}{AGU} & \multirow{2}{*}{Ser} & \multicolumn{1}{|c}{\textbf{U}} \\  
                & & \multicolumn{1}{|c}{AUC} & & \multicolumn{1}{|c}{ACC} & & \multicolumn{1}{|c}{AAC} & & \multicolumn{1}{|c}{AGC} & & \multicolumn{1}{|c}{\textbf{C}} \\  
                & & \multicolumn{1}{|c}{AUA} & & \multicolumn{1}{|c}{ACA} & & \multicolumn{1}{|c}{AAA} & \multirow{2}{*}{Lys} & \multicolumn{1}{|c}{AGA} & \multirow{2}{*}{Arg} & \multicolumn{1}{|c}{\textbf{A}} \\  
                & & \multicolumn{1}{|c}{AUG} & Met & \multicolumn{1}{|c}{AGG} & & \multicolumn{1}{|c}{AAG} & & \multicolumn{1}{|c}{AGG} & & \multicolumn{1}{|c}{\textbf{G}} \\ \cline{3-10}  
             
                & \multirow{4}{*}{\textbf{G}} & \multicolumn{1}{|c}{GUU} & \multirow{4}{*}{Val} & \multicolumn{1}{|c}{GCU} & \multirow{4}{*}{Ala} & \multicolumn{1}{|c}{GAU} & \multirow{2}{*}{Asp} & \multicolumn{1}{|c}{GGU} & \multirow{4}{*}{Gly} & \multicolumn{1}{|c}{\textbf{U}} \\  
                & & \multicolumn{1}{|c}{GUC} & & \multicolumn{1}{|c}{GCC} & & \multicolumn{1}{|c}{GAC} & & \multicolumn{1}{|c}{GGC} & & \multicolumn{1}{|c}{\textbf{C}} \\  
                & & \multicolumn{1}{|c}{GUA} & & \multicolumn{1}{|c}{GCA} & & \multicolumn{1}{|c}{GAA} & \multirow{2}{*}{Glu} & \multicolumn{1}{|c}{GGA} & & \multicolumn{1}{|c}{\textbf{A}} \\  
                & & \multicolumn{1}{|c}{GUG} & & \multicolumn{1}{|c}{GGG} & & \multicolumn{1}{|c}{GAG} & & \multicolumn{1}{|c}{GGG} & & \multicolumn{1}{|c}{\textbf{G}} \\ \cline{3-10}  
             
           \end{tabular}  
           \end{center}  
           \caption{Código genético representado por tripletes de bases y los aminoácidos que codifican.}  
           \label{TablaCodigoGenetico}  
      \end{table}  
   
 \end{document}  

Lo que sería una sencilla tabla en Word, en LaTeX parece un complejo problema de ingeniería, la tabla compilada quedaría finalmente así:
Reto a Joaquín a hacer más sencillo el código de la tabla ;)