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.
Ideas y código para problemas de genómica de plantas, biología computacional y estructural
18 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.
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.
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.
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:
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:
Se obtiene algo como:
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:
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
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:
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.
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
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:
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 ;)
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 ;)
23 de marzo de 2012
Es primavera: lidiando con Flower
Hola,
vamos a hablar un poco sobre el trabajo con secuencias obtenidas de experimentos de NGS (Next Generation Sequencing; para despistados) con el secuenciador 454. En concreto, nos vamos a centrar en una herramienta que puede servir para tener un primer vistazo de los datos obtenidos, antes de realizar el ensamblaje, mapeo, ...
Para instalarlo lo más sencillo es utilizar cabal, gestor de paquetes de Haskell, y bajar además ghc, entorno para desarrollo y compilación en ese lenguaje.
sudo apt-get install ghc cabal-install
cabal install biosff
Vamos a desarrollar algunos ejemplos con los datos de la entrada SRR000001 del NCBI SRA, run que se llevó a cabo en un equipo Roche 454 GS FLX en 2007. Quizás la opción más simple es la que nos ofrece un vistazo rápido del experimento, dándonos información de índice generado (creo que es un índice a la manera de un suffix array), el número de reads, número de flows y la key usada para el run: secuencia de bases que se añade a toda la muestra durante la preparación de las librerías.
flower -i SRR000001.sff
Index: (782054672,9419708)
Num_reads: 470985
Num_flows: 400
Key: TCAG
En éste caso vemos que hay un index (se regeneró mediante sfffile -o, ya que las secuencias de SRA no suelen llevar el índice). Con aquel run obtuvieron 470K reads en 400 flows (100 ciclos: 100 nts/read de media esperada). Por último nos indica que la key es TCAG (aunque al parecer siempre ha sido la misma).
Flower (algo así como FLOW ExtracteR), puede generar FASTA y FASTQ con scores en Phred+33 o basados en la codificación de Illumina. La carencia de un comando para obtener un FASTQ es una de las cosas que más parece echar de menos la comunidad que utiliza SFF Tools. Para obtener el FASTQ en Phred+33 con Flower:
flower -Q SRR000001.sff > SRR000001.fastq
Con la opción -s se obtienen datos sobre cada uno de los reads del experimento. La salida del programa da los campos separados convenientemente por tabuladores.
flower -s SRR000001.sff > SRR000001.reads
# name........ date...... time.... reg trim_l trim_r x_loc y_loc len K2 trimK2 ncount avgQ travgQ
EM7LVYS01C1LWG 2 7-04-10 15:13:00 01 1 235 1131 1422 255 0.74 0.76 0 26.38 26.33
K2 es una métrica de calidad "K-square", que es la suma de cuadrados de las distancias desde el entero más cercano para cada flow value. Como el sistema de scoring de 454 se basa en estimar la longitud de un homopolímero en base a la señal registrada en un flow, un entero supondría la definición de una longitud de homopolímero con exactitud. ncount es el número de flow cycles sin valor determinado o basecalls ambiguos, que también llevan asociado un error si el flow value no es exactamente 0. avgQ es la calidad promedio del read. Todos los campos que comienzan con tr son equivalentes a los anteriores, pero tras realizar el trimming indicado en el SFF.
Además, para el análisis del experimento también se pueden generar datos para representación gráfica de los flow values. Con la opción -h se obtiene la distribución de flow values por nucleótido.
flower -h SRR000001.sff > SRR000001histogram.txt
Score A C G T
0.00 361230 586168 620918 409514
0.01 136705 284992 305021 168203
0.02 182110 407799 437550 228386
0.03 239742 571509 614918 306284
0.04 311296 779612 839794 401568
0.05 397364 1029955 1110715 513386
0.06 503853 1318643 1412014 648700
....
así hasta 99.99
Con la opción -H se obtiene la distribución de flow values para los flow cycle. Si representamos el resultado de esta salida, tenemos un gráfico como el siguiente, donde se muestra claramente que el sistema está optimizado para obtener la menor tasa de error posible. Como ya se ha comentado, la tasa de error está asociada a los puntos entre cada dos enteros. Para poder dar el mayor quality score posible, la mayor confianza posible de que la longitud del homopolímero es la predicha por el basecall, se debe maximizar el número de flow values lo más cercanos que sea posible a un entero.
Imagen tomada de Bioinformatics (2010) 26 (18): i420-i425 |
En un principio, Flower incluia una herramienta para seleccionar reads del SFF (flowselect). Sin embargo, en la versión de biosff esta utilidad ha desaparecido. Podría venir a suplir las opciones -i y -e de sfffile, que generalmente se usan ya durante el proceso de análisis, para generar subconjuntos de los datos originales una vez se tienen datos de mapeo o ensamblaje.
Como conclusión, Flower es más que una alternativa a SFF Tools para el pre-análisis de las secuencias de 454, en vista de la buena representación de los datos, el buen desempeño que tiene y las opciones extra que nos aporta. Por otro lado, sigue sin cubrir algunas de las utilidades que presenta SFF Tools, sfffile en especial, pero es código libre así que todos estamos invitados a mejorarlo.
saludos!
Carlos
Suscribirse a:
Entradas (Atom)