6 de abril de 2016

Calculando experimentos de secuenciación

Buenas,
hoy necesitábamos calcular cuántos individuos (de una especie monocotiledónea) podríamos secuenciar con cierta profundidad en un secuenciador Illumina, pensando en el HiSeq2500 en concreto. Al final decidimos apostar por una profundidad promedio de 80x, para is sobre seguro, usando parejas de lecturas de 2x125b. Buscando en Internet encontré rápidamente una calculadora del propio fabricante que igual algunos no conocéis y puede ayudar a hacer esto rápidamente.

Figura tomada de http://www.danielecook.com/calculate-depth-coverage-bam-file.


Vayamos con un ejemplo con la calculadora
[ http://support.illumina.com/downloads/sequencing_coverage_calculator.html ]:

0. tipo de secuenciación: DNA             [se puede elegir RNA también]
1. protocolo: whole-genome sequencing  [otras: Nextera, Truseq, custom]
2. tamaño del genoma: 320Mbp
3. profundidad deseada: 80x
4. % de duplicados: 2%                [valor por defecto]
5. instrumento: HiSeq 1500/2500

Volumen total de secuenciación requerido: 26,1Gb   [26.122.448.980b]

En mi ejemplo, usando el protocolo v4, esto equivale a 0.42 líneas o lanes, lo que significa que podría poner hasta 2 muestras por línea.

Hasta luego,
Bruno

PD Me dicen mis colegas Dave Des Marais y Pat Edger que la longitud de un genoma (de plantas en este caso)  puede estimarse aproximadamente a partir del contenido en DNA del núclo usando la fórmula long = 1C * 980.

4 de marzo de 2016

R one-liners

Hola,
en esta entrada quería compartir unos ejemplos para aprovechar las capacidades del lenguaje R para analizar de manera rápida datos desde el terminal, espoleado por el primero de ellos, que me pasó Carlos Cantalapiedra la semana pasada. Pero vayamos al grano usando como ejemplo un fichero de salida tabular de BLAST, al que llamarmos 'all.blast'.



Receta 1. Tienes un fichero con columnas de números y quieres saber la media de una de ellas, por ejemplo la tercera. Para la la desviación estándar cambia 'mean' por 'sd':

$ cut -f 3 all.blast | Rscript -e 'data=scan(file="stdin"); mean(data)'


Receta 2. Tienes un fichero con columnas de números y quieres calcular un histograma de una de ellas, por ejemplo la tercera:

$ cut -f 3 all.blast | \
Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)' 

Se puede modificar mostrando el PDF creado inmediatamente:

$ cut -f 3 all.blast | \
Rscript -e 'data=scan(file="stdin"); pdf("hist.pdf");hist(data)'; evince hist.pdf 

Receta 3. Tienes un fichero con columnas de números y quieres calcular boxplots con dos de ellas, por ejemplo la primera y la segunda:

cut -f 1,2 all.blast | \
Rscript -e 'd=matrix(scan(file="stdin"),ncol=2); pdf("box.pdf");boxplot(d[,1],d[,2])'

Espero que os sirvan de inspiración,
un saludo,
Bruno

13 de febrero de 2016

Expresión regular de la familia de las O-fucosiltransferasas



Hola, 
el pasado 8 de febrero se publicó en la revista Nature Chemical Biology  (http://dx.doi.org/10.1038/nchembio.2019) un artículo donde se  describen las bases moleculares de la reacción de O-fucosilación. Ésta es una modificación postraduccional poco frecuente, que realizan las enzimas O-fucosiltransferasas, como nos explica la investigadora de nuestro grupo Inmaculada Yruela, una de las autoras del trabajo [reseña completa en www.eead.csic.es]:

"Esta reacción resulta esencial en algunas rutas metabólicas de los organismos eucariotas, incluidas las plantas, para mantener las funciones básicas de las células. El artículo describe el mecanismo por el cual la enzima O-fucosiltransferasa 2 (POFUT2) reconoce sin errores una secuencia de aminoácidos (TSR) de la proteína receptora y le transfiere una molécula de azúcar tipo fucosa –así resulta fucosilada–. Hasta la fecha no se conocía cómo estas enzimas reconocen y se unen a sus sustratos proteicos. La O-fucosilación es esencial para el correcto plegamiento y estabilidad del dominio TSR y el reconocimiento molecular de POFUT2."
Como se ve en el alineamiento, el dominio TSR contiene tres puentes disulfuro y una secuencia consenso en torno a los dos primeras cisteínas CX{2,3}[S|T]CX{2}G , lo que se llama un motivo, como los del repositorio Prosite:
Fragmento de un alineamiento múltiple de dominios TSR, tomado de http://dx.doi.org/10.1038/nchembio.2019

Tras alinear algunas secuencias de ejemplo y perfeccionar la expresión regular, ésta se empleó para identificar todas las secuencias con dominios TSR en los proteomas de Homo sapiens y Caenorhabditis elegans. Más generalmente, el siguiente trozo de código permite localizar dentro de un fichero FASTA todas las secuencias reconocidas por una expresión regular:
 #!/usr/bin/perl  
 use strict;  
   
 # script that takes a FASTA file(s) and scans all protein sequences   
 # looking for matches of a chosen motif expressed as a regular expression  
 # 11042015: edited following comments by JF
   
 # constant as indices to access read_FASTA_file arrays  
 use constant NAME => 0;   
 use constant SEQ => 1;  
   
 #my $motifRE = '.*C[^C]{0,21}.C[^C]{2,6}[S|T].C[^C]{4,75}C[^C]C[^C]{4,15}C.*'; 
 my $motifRE = 'C[^C]{0,21}.C[^C]{2,6}[ST].C[^C]{4,75}C[^C]C[^C]{4,15}C';
   
 if(!@ARGV){ die "# usage: $0  ... \n"; }  
 else{ print "# motif: $motifRE\n"; }  
   
 my ($n_of_matches,@matches) = (0);  
 foreach my $infile (@ARGV)  
 {  
   next if($infile !~ /.fa/);  
   my (@FASTA,$start,$end,$l);  
   my $n_of_sequences = -1;  
   open(FASTA,"<$infile") || die "# cannot read $infile $!:\n";  
   while()  
   {  
    next if(/^$/);  
    if(/^\>(.*?)[\n\r]/)  
    {  
      $n_of_sequences++;   
      $FASTA[$n_of_sequences][NAME] = $1;  
    }                
    else  
    {  
      s/[-\s\n.]//g; #s/[\s|\n|\-|\.]//g;  
      $FASTA[$n_of_sequences][SEQ] .= $_;  
    }  
   }  
   close(FASTA);  
   
   foreach my $seq ( 0 .. $n_of_sequences )  
   {  
    while($FASTA[$seq][SEQ] =~ m/($motifRE)/gi)   
    {  
      #$start=length($`);  
      #$l=length($1);  
      #$end=$start+$l;  

      $start = $-[1];
      $end   = $+[1];
      $l     = $end - $start + 1;

      print ">$FASTA[$seq][NAME] match=($start,$end,$l)\n$1\n";           
    }  
   }  
 }  

Un saludo,
Inma y Bruno

18 de enero de 2016

Word + bioinformática = CRISPR

Hola,
este será un artículo muy breve, porque tengo poco que escribir yo mismo, a parte de volver a comprobar una vez más que la ciencia básica, tan denostada a menudo, es a menudo fundamental para los avances más espectaculares de la ciencia aplicada.  En vez de contar yo aquí una historia hoy os invito a leer la increíble historia del descubrimiento de los sistemas CRISPR-Cas9, que ya están revolucionando la ciencia y pronto la medicina.

La fuente original, Cell, bajo subscripción:
http://www.cell.com/cell/fulltext/S0092-8674%2815%2901705-5

Un resumen en español en El Confidencial:
http://www.elconfidencial.com/tecnologia/2016-01-16/crispr-francis-mojica-charpentier-doudna-edicion-genomica_1136337/

Reconocimiento de una hebra sencilla de DNA por hibridación de un RNA guía de 32 bases (crRNA) por parte del complejo de proteínas Cascade. Tomado de [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC4427192/]. Como se ve en nuestro repositorio 3D-footprint, la especificidad de reconocimiento es de las más altas observadas entre proteínas naturales: http://floresta.eead.csic.es/3dfootprint/complexes/4qyz_ABCDEFGHIJ.html


 Ah, y por cierto, esta historia una cura de humildad para los que nos dedicamos a inventar y probar complicados algoritmos. A veces, la herramienta "Buscar y reemplazar" de Word es todo lo que hace falta para hacer descubrimientos en Biología, como a veces me recuerda mi colega Ana Casas, y eso también es Biología computacional,
un saludo,
Bruno

11 de enero de 2016

Bioinformática en C (y perl6 async)

Hola, y feliz año!
tras la última entrada, y antes de que escriba algo más completo sobre perl6,
parece que la gestión de tareas asíncronas tiene mucho potencial...

Pero vamos a lo que vamos. Hoy quería compartir un documento sobre cómo programar en C con los compiladores actuales, porque algunos todavía seguimos usando sintaxis anticuada y esto es relevante en nuestro campo, dado que hay muchas aplicaciones bioinformáticas escritas en C, normalmente por cuestiones de eficiencia.



El documento en cuestión es: https://matt.sh/howto-c

y ofrece un montón de trucos y recomendaciones para modernizar y simplificar la escritura de código en C, con ejemplos sencillos como éste para vectores dinámicos:

 uintmax_t arrayLength = strtoumax(argv[1], NULL, 10);
    void *array[arrayLength];

 /* Hace innecesario liberar el puntero al final.
    Ojo: debe ser un tamaño razonable para el sistema donde se ejecute. */

Ha sido traducido al español aquí y comentado/criticado aquí y en menéame,
un saludo,
Bruno