13 de junio de 2012

Alineamiento con transformadas de Fourier

En 1992 el grupo de Ilya Vakser publicó en PNAS el método fundamental para poder hacer simulaciones de docking de biomoléculas de manera eficiente, disminuyendo considerablemente el tiempo de cálculo aplicando transformadas rápidas de Fourier. Diez años más tarde se publicó la primera versión del exitoso programa de alineamiento múltiple de secuencias MAFFT, del que ya hemos hablado en otras entradas, que aplica ideas similares para el problema de encajar secuencias similares, como se hace al alinear secuencias.

Figura tomada de http://www.ncbi.nlm.nih.gov/pubmed/12136088. donde se muestran dos picos del análisis de Fourier que se corresponden con dos subalineamientos locales.

En esta entrada mi intención es explicar el fundamento de esta técnica recurriendo al lenguaje Octave, la versión open source de MatLab. El código es el siguiente:

 % alineamiento (sin gaps) de secuencias proteicas usando la FFT  
 % escrito en Octave, que es compatible con Matlab   
 % requiere la libreria 'bioinfo'   
   
 clear all  
   
 %% Conversor binario de secuencias   
 %% Parametros:   
 %% 1) longitud de la secuencia S  
 %% 2) longitud del fragmento F (de S)  
 function [S , F] = aa2bits( sequence , fragment )  
   
 % cada residuo se representa por un numero del 1 al 20  
 % http://www.mathworks.com/help/toolbox/bioinfo/ref/aa2int.html  
 [ seq ]= aa2int( sequence )  
 [ frag ] = aa2int( fragment )   
   
 % codificamos residuos como cadenas binarias de ancho fijo (20 columnas)  
 S=[];  
 for i=1:length(seq)  
   % cada residuo ocp  
   S=[S ones(1,seq(i)) zeros(1,20-seq(i)) ];  
 end  
   
 F=[];  
 for i=1:length(frag)  
   F=[F ones(1,frag(i)) zeros(1,20-frag(i)) ];  
 end  
   
 % completa F hasta igualar la longitud de S  
 F=[F zeros(1,length(S)-length(F))];  
   
 endfunction % no se usa en Matlab  
   
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
   
 % define la secuencia S y el fragmento F que queremos alinear  
 sequence = 'SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY';  
 fragment = 'KMLNSVCRSWWWWC';  
   
 % convertimos ambas secuencias en dos segnales binarias,  
 % de manera que cada residuo se representa por un natural del 1 al 20  
 [S , F] = aa2bits(sequence,fragment);  
   
 % calculamos el alineamiento de F con S optimizando la funcion objetivo Fobj,  
 % que en este ejemplo es simplemente el producto (como un AND binario)  
 % obtenido al desplazar F de izq a der sobre S:   
 % producto(n) = sum { F(i) * S(i-n) }  
 %        i   
 % aprovechamos que Fobj se puede expresar por medio de transformadas  
 % de Fourier para reducir las operaciones necesarias (de N^2 a 4NlogN)  
 % https://sites.google.com/site/cartografiaygeodesia/prac7.pdf  
 % http://www.pnas.org/content/89/6/2195.abstract  
   
 FTsec = fft(S);  
 FTfrag = fft(F);  
   
 FTproducto = conj(FTfrag) .* FTsec; % producto termino a termino  
   
 % Deshacemos transformacion y localizamos valor maximo,   
 % la posicion que optimiza el alineamiento   
 producto = ifft(FTproducto);   
 [valor maxpos] = max(producto);  
   
 % imprime alineamiento sobre secuencicas binarias   
 tit = sprintf("Optimal alignment, optimal position = %d",maxpos);  
 plot(S*0.75)  
 title(tit)  
 xlabel('sequence position')  
 hold on  
 plot([zeros(1,maxpos) F*0.9] , 'r');  
 axis([0 length(S) 0 1]);  
 legend('sequence','fragment');  
 print -dpng figure_align.png;  
   
 % alineamiento sobre secuencias originales, cambiando maxpos de escala  
 align = printf("Alignment:\nS %s\nF %s%s\n",sequence,blanks((maxpos-1)/20),fragment);  

Como resultado obtendremos un alineamiento como éste, que :
S SDEVRKNLMDMFRDRQAFSEHTWKMLLSVCRSWAAWCKLNNRKWFPAEPEDVRDYLLY
F                        KMLNSVCRSWWWWC
 

Que se corresponde con éste producto de Fourier:



En Perl podríamos usar el módulo Math::FFT para este fin, o la GNU Scientific Library que ya revisamos en otra entrada. Un saludo,
Bruno

5 de junio de 2012

WINTER SCHOOL: ALGORITHMS IN STRUCTURAL BIOINFORMATICS

Hola,
ha llegado a mis manos un anuncio de un curso con contenidos muy cercanos a algunas entradas de este blog e incluso a nuestro material de Algoritmos en Bioinformática Estructural. El curso tendrá lugar en Diciembre en la riviera francesa y acogerá a 20 estudiantes. Os pego los detalles en inglés, el enlace al curso es www-sop.inria.fr/manifestations/algoSB :
 
Dear Colleague,

please find below the announcement for a Winter School on Algorithms
in Structural Bioinformatics: we would appreciate your forwarding it
to PhD students and postdocs who could be interested.

With best regards,
 The organizers :
    Frédéric Cazals, INRIA Sophia Antipolis - Méditerranée
    Juan Cortés, LAAS-CNRS, Toulouse 


%%i%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ALGORITHMS IN STRUCTURAL BIOINFORMATICS : WINTER SCHOOL
   2-7 DECEMBER 2012, INRIA SOPHIA ANTIPOLIS, FRANCE  
   http://www-sop.inria.fr/manifestations/algoSB/
%%i%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We are pleased to announce a one-week school on Algorithms in
Structural Bioinformatics. The aim is to introduce advanced methods in
this domain, giving special attention to interdisciplinary
approaches. The main focus will be on methodological developments
meant to analyze and predict macromolecular assemblies, as well as on
the corresponding software. The following topics will be taught:

I.   Obtaining and organizing structural information for modeling studies (J. Janin, C. Robert)
II.  Modeling protein complexes and assemblies with Voronoi diagrams (F. Cazals)
III. Molecules as robots: mining the flexibility of proteins (J. Cortés)
IV.  Structural comparisons (R. Andonov, N. Malod-Dognin)
V.   Docking algorithms (C. Prévost, M. Zacharias)

The program will span 5.5 days. The first afternoon will consist of a
mini research symposium during which each participant will present
his/her research interests and achievements.  Each of the five
subsequent days will consist of a lecture (morning session) and a
hands-on computer lab (afternoon session).

Additional information can be found at
http://www-sop.inria.fr/manifestations/algoSB/

INTENDED AUDIENCE. The school is primarily intended for twenty
European PhD students and postdoctoral researchers, but
applications from other parts of the world will also be
considered if space is available.

DEADLINE FOR APPLICATIONS. A one-page application (pdf format) should
be sent by email to algosb-coordinators@inria.fr with Subject "AlgoSB
application", by September 7th, 2012.

The application should provide a mini-vitae or a pointer to the full
vitae, a brief presentation of the applicant's research interests, and
a few lines commenting on expected synergy between one's research and
the courses that will be taught.

Notification of acceptance will be provided by September 30th. It is
expected that all accepted students will participate in the entire
school period.

VENUE AND PARTICIPATION FEES. The Algorithmics in Structural
Bioinformatics school will take place on the French Riviera, at INRIA
Sophia Antipolis - Méditerranée.

The participation fees will be of 390 euros, which include lunches,
coffee breaks, and social events.
 
This school is sponsored by INRIA (http://www.inria.fr/), the CNRS/GDR
Bioinformatique Moleculaire (http://www.gdr-bim.u-psud.fr/), and by
SANOFI (http://en.sanofi.com/).

The organizers :
    Frédéric Cazals, INRIA Sophia Antipolis - Méditerranée
    Juan Cortés, LAAS-CNRS, Toulouse 

18 de mayo de 2012

beca FPU en genómica comparada de cereales

El laboratorio de Biología Computacional de la Estación Experimental de Aula Dei/CSIC (en Zaragoza, España) busca un candidato/a para desarrollar un proyecto de investigación sobre genómica comparativa de cereales y la planta modelo Brachypodium distachyon.

El grupo tiene amplia experiencia en diferentes áreas de la Bioinformática y actualmente participa en diversos proyectos de genómica y secuenciación de plantas, incluyendo Arabidopsis thaliana, arroz y cebada, con un interés especial en el estudio funcional y evolutivo de las redes de regulación genética y en el descubrimiento de las raíces genéticas de la biodiversidad vegetal y agrícola en relación con la adaptación al medio ambiente. Para una lista completa de las publicaciones del grupo consultar http://www.eead.csic.es/compbio y http://www.eead.csic.es/index.php?id=99 .

Las personas interesadas deberán cumplir los requisitos para solicitar una beca FPU (2 años beca + 2 años contrato, ver convocatoria en http://tinyurl.com/863o6cr y requisitos en http://tinyurl.com/7klg24g). Se requiere expediente académico superior a 1.6, preferiblemente superior a 2, y se valorará positivamente la experiencia previa en Bioinformática y el haber terminado un Máster oficial. La convocatoria especifica una fecha límite de terminación de los estudios de licenciatura/grado/ingeniería posterior al 1 de enero de 2009.

Interesados contactar con el Dr. Bruno Contreras Moreira (bcontreras@eead.csic.es) o el Dr. Ernesto Igartua (igartua@eead.csic.es) antes del 31 de Mayo de 2012.

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