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 ;)

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, ...

Flower (Bioinformatics, 2011) es un programa desarrollado por Ketil Malde, también implicado en el desarrollo de FlowSim. Está escrito con Haskell, el lenguaje de programación funcional, y anteriormente se distribuía el paquete como tal (última versión flower_v0.7) pero hoy día viene en la librería biosff (v0.2). Se trata de una alternativa GPL a las SFF Tools que distribuye Roche con el paquete Data Analysis de su software propietario.

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

Correlación entre transcriptoma y proteoma

Hola,
una cuestión que ya dio qué hablar en los tiempos de los chips de RNA (microarrays) y que ha vuelto a resurgir recientemente es hasta qué punto podemos hacernos una idea de la actividad celular si sólo medimos la expresión de los genes en forma de mRNAs, cuando en realidad muchas de las funciones las desempeñan proteínas.
Por ejemplo, Tian et al.(2004) publicaron correlaciones del 40% en líneas celulares hematopoyéticas de mamíferos.

Sin embargo, ahora en muchos ámbitos los microarrays están siendo sustituidos por experimentos de RNAseq, y por tanto es pertinente volver a hacerse la pregunta, dado que se espera que esta tecnología sea superior. Dos artículos muy recientes nos ayudan precisamente a valorar esto.

El trabajo de Jiang et al.(2011), que propone el uso de controles para corregir los niveles de expresión medidos por RNAseq, incluye una figura que resume la distribución típica de valores de expresión génica que se obtienen por RNAseq:
Figura original de Jiang et al. (2011). En el A panel se muestra la equivalencia entre 'copias por célula' y 'fragmentos por kilobase mapeada' (FPKM). En el B se muestra el histograma de genes expresados, que se solapa a la izquierda (un hombro) con los valores de genes que se expresan menos de una copia por célula. Los datos se obtuvieron a partir de material extraído de una línea celular de Drosophila melanogaster.

El trabajo de Nagaraj et al. (2011), que se basa en material de la archiconocida línea celular humana HeLa, encuentra una distribución similar de mRNAs, lo que sugiere que debe ser una distribución estándar, pero además la correlaciona con valores de abundancia de péptidos medidos por espectrometría de masas de alta resolución. La figura siguiente resume un montón de datos experimentales:
Figura original de Nagaraj et al. (2011).
En el panel A efectivamente se reproduce el hombro observado en los datos de RNAseq de D.melanogaster y en el B se muestran los datos equivalentes de proteómica. El diagrama de Venn C es interesante, porque muestra que hay muchos mensajeros que no se observan como proteína, y algunas proteínas cuyo mensajero está ausente. El panel D muestra la distribución funcional de los genes expresados y, finalmente, el diagrama de dispersión E muestra la correlación observada entre ambos tipos de mediciones, con un coeficiente de correlación de Spearman de 0.6,
un saludo,
Bruno

Actualizaciones posteriores:
[1] En este trabajo se estudia cómo se heredan los niveles de expresión en mosca
[2] En este artículo se comenta que los genes ortólogos entre distintos organismos tienen mayores correlaciones en sus [proteína] que en [mRNA]

6 de marzo de 2012

ofertas de postdoc en filogenómica

Hola,
copio aquí dos ofertas de trabajo que me han llegado desde bioinformatibericos :

Postdoctoral Position in Evolutionary Genomics at IFCA (Universidad of Santander-CSIC)
A postdoctoral position is available immediately to work with Rafael Zardoya (Madrid), Julio Rozas (Barcelona), David Posada (Vigo), and Jesus Marco (Santander) on a common 
project related with gastropod phylogenomics. The position has a term of 1 year and mobility among the referred labs. 
Prerrequisites: A PhD in Biology, Chemistry, Computer Science or related fields, with proved skills in computational biology. 
Priority will be given to candidates with previous experience in NGS data analysis (assembling, gene annotation, etc), and with expertise in bioinformatics programing languages (Perl, Python, 
C, etc). 
It will be desirable some experience with Phylogenetics and/or Comparative Genomics or evolutionary biology in general, and with Relational databases (MySQL).
If interested please contact Rafael Zardoya, Department of Biodiversity and Evolutionary Biology, Museo Nacional de Ciencias Naturales, rafaz@mncn.csic.es. Application review 
will begin on March 5, 2012 and continue until the position is filled. To apply, please send the following:  
1. A curriculum vitae 
2. Names of 2 referees willing to provide a letter of recommendation upon request 
3. A brief statement of research interests and goals.
 
Postdoc on chromatin biophysics (Structural Genomics Group, National Center for Genomic Analysis)
The successful candidate will team with other members of our group in our efforts towards elucidating the 3D structure of genomic domains and genomes. Recent publication from the group in this field include:

- Umbarger, M. A., Toro, E., Wright, M. A., Porreca, G. J., Baù, D., Hong, S.-H., Fero, M. J., et al. (2011).  Molecular Cell44(2), 252–264
- Marti-Renom, M. A., & Mirny, L. A. (2011). PLoS Computational Biology7(7), e1002125
- Sanyal, A., Baù, D., Martí-Renom, M. A., & Dekker, J. (2011). Current Opinion in Cell Biology. doi:10.1016/j.ceb.2011.03.009
- Baù, D., Sanyal, A., Lajoie, B. R., Capriotti, E., Byron, M., Lawrence, J. B., Dekker, J., et al. (2011). Nature Structural & Molecular Biology18(1), 107–114.

I would appreciate if you could forward this announcement to whom may be interested,
Marc A. Marti-Renom, Group Leader 
Parc Cientific de Barcelona - Torre I - Baldiri Reixac, 4 - 2a.p - 08028 - Barcelona
Tel +34 934 033 743  Fax +34 934 037 279
email mmarti@pcb.ub.cat   web http://marciuslab.org & http://cnag.cat




28 de febrero de 2012

Leyendo el diagrama de Ramachandran

Hola,
la semana pasada les planteé a mis alumnos de la Licenciatura en Ciencias Genómicas el problema de asignar estructura secundaria a los residuos de una proteína, una vez medidos los ángulos dihedros del esqueleto peptídico psi y phi. Para esta tarea lo lógico es usar el diagrama de Ramachandran, como éste tomado de la wikipedia:


La solución naive consiste en delimitar las regiones del diagrama por medio de condiciones tales como 'si phi > X  AND psi < Y' entonces un residuo está en estado Z. Sin embargo, para tener mucha precisión sería necesario tener muchas condiciones como ésta, dada las formas irregulares del diagrama.


Una solución interesante, propuesta por F.Peñaloza, consiste en realmente leer el mapa, como lo hacemos nosotros. Lo que él hizo fue obtener un mapa de Ramachadran cuadrado con unos cuantos colores conocidos,

Diamagra de Ramachandran RGB de F.Peñaloza.


Diagrama anterior con los colores RGB empleados

para luego convertirlo a texto ASCII (con algo como text-image). Con este mapa en memoria, es sencillo obtener la estructura secundaria para unas coordenadas phi,psi dadas. 

Aprovechando el módulo GD de perl podemos leer la imagen directamente, como muestra el siguiente código:
 use strict;  
 use GD;  
   
 my $RAMAPLOT = 'ramachandran.png';  
 my $plot = GD::Image->newFromPng($RAMAPLOT);  
 my ($maxpsi,$maxphi) = checkPlot($plot);  
   
 my $inputPSI = 160;   
 my $inputPHI = -45;   
 printf("El estado de estructura secundaria para %s,%s es: %s\n",  
    $inputPSI,$inputPHI, getSimpleSS($maxpsi,$maxphi, $inputPSI,$inputPHI));  
   
 sub getSimpleSS  
 {  
   my ($maxpsi,$maxphi,$psi,$phi) = @_;  
   my %colorSS = (  
     '236,30,36','E',  
     '164,74,164','E',  
     '156,218,236','H',  
     '4,162,236','H',  
     #'180,230,28','L'   
     );  
   
   if($psi > 180 || $psi < -180 || $phi > 180 || $phi < -180)  
   {  
     die "# valid psi/phi are [-180,+180]\n";  
   }  
   $psi = (0.5 - (0.5 * $psi/180)) * $maxpsi;  
   $phi = (0.5 + (0.5 * $phi/180)) * $maxphi;  
     
   my @RGB = $plot->rgb( $plot->getPixel($phi,$psi) );  
   return $colorSS{"$RGB[0],$RGB[1],$RGB[2]"} || 'C';  
 }  
   
 sub checkPlot  
 {  
   my ($plot,$print_pixels) = @_;  
   my ($maxphi,$maxpsi) = $plot->getBounds();  
   if(!$print_pixels){ return ($maxphi,$maxpsi) }  
     
   foreach my $phi ( 0 .. $maxphi )   
   {  
     foreach my $psi ( 0 .. $maxpsi )   
    {  
       my @RGB = $plot->rgb( $plot->getPixel($phi,$psi));  
       print "$RGB[0] $RGB[1] $RGB[2]\n";  
    }  
   }  
   return ($maxphi,$maxpsi);  
 }  

Un saludo, Bruno