22 de abril de 2016

mapeo fino de genes por NGS

Buenas,
esta semana copio aquí una reseña de un trabajo recientemente publicado de Carlos P Cantalapiedra, autor habitual de este blog y próximo doctor del grupo, donde se explica el proceso para localizar un loci responsable de una resistencia a infección por parte de hongos, combinando genética clásica y secuenciación de nueva generación: http://www.eead.csic.es/spreading/showspreading?Id=416

Pongo aquí una de las figuras del artículo:

Genotipo de varias líneas de cebada en torno al locus que confiere resistencia. En naranja, genotipos como los del parental resistente. En verde, genotipos como los del parental susceptible. La captura de exoma permite reducir la zona de búsqueda al punto donde se unen ambos genotipos (punto 211721 dentro del recuadro). Adaptada de https://dl.sciencesocieties.org/publications/tpg/first-look/pdf/plantgenome2015.10.0101.pdf.

La referencia del artículo completo, en inglés, es:

Cantalapiedra CP, Contreras-Moreira B, Silvar C, Perovic D, Ordon F, Gracia MP, Igartua E, Casas A. (2016) A cluster of NBS-LRR genes resides in a barley powdery mildew resistance QTL on 7HL. The Plant Genome. Early access. DOI: 10.3835/plantgenome2015.10.0101. URL.

Hasta luego,
Bruno

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