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

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

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

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

Despedimos el 2015, viene Perl 6

Buenas,
quedan ya pocos días para terminar el año, o lo que es lo mismo, pocos días para que podamos probar la primera versión oficial de Perl 6, tras demasiados (15) años de desarrollo, de la mano del compilador Rakudo. Si no os suena nada de esto, podéis empezar por leer las preguntas frecuentes o directamente empezar tutoriales como perl6intro.

Ha habido mucho morbo con la muerte de Perl, y seguirá habiendo a lo largo de 2016 discusiones sobre la compatibilidad de perl 6 y 5, el que usamos ahora, sobre la eficiencia del nuevo lenguaje, o sobre la portabilidad de CPAN a perl 6. Todo esto habrá que ir viendo como se desarrolla muy pronto, pero desde luego es una noticia importante para los usuarios de Perl, esperemos que buena!

Un saludo y los mejores deseos para el próximo año, long live perl!
Bruno

Post-datas del 29 de Diciembre: 

1) Podéis ampliar esto en [http://perlweekly.com/archive/231.html]

2) Para instalar rakudo vía git debes tener abierto el puerto 9418; en caso contrario es posible hacerlo vía https añadiendo estas líneas a tu archivo ~/.gitconfig: 

[url "https://"]
    insteadOf = git://
[url "https://github.com"]
    insteadOf = git@github.com

Post-data del 1 de Enero: 

1) Tutorial exprés en [https://learnxinyminutes.com/docs/perl6]
2) Si ya usas perl5 [http://doc.perl6.org/language/5to6-nutshell]

Post-data de Junio de 2016: 

Mientras no sea pueda instalar perl6 directamente en Ubuntu prefiero usar rakudobrew para instalar perl6 como un ecosistema aislado, como se explica en: [http://rakudo.org/how-to-get-rakudo/#Installing-Rakudo-Star-Source-Rakudobrew]