19 de diciembre de 2014

kseq::klib to parse FASTQ/FASTA files

Buenas,
éste tiene pinta de ser el último artículo del año, no? Desde hace tiempo venimos usando en el laboratorio la librería kseq.h , integrante de [ https://github.com/attractivechaos/klib ], para las desagradecidas tareas de manipular grandes archivos en formato FASTQ. Esta librería es un único archivo .h en lenguaje C que de manera muy eficiente permite leer archivos FASTQ, y también FASTA, incluso comprimidos con GZIP, con zlib. La eficiencia es fundamental para estas aplicaciones, porque modificar archivos .fastq de varios GB puede consumir mucho tiempo.

Fuente: http://drive5.com/usearch/manual/fastq_files.html


Como resultado de esta experiencia hace unos meses colgamos de la web del labo un programita, que llamamos split_pairs. Este software facilita trabajos tediosos como pueden ser modificar las cabeceras de una archivo FASTQ con expresiones regulares para que un programa posterior no se queje, y de paso emparejar o interlazar los reads según convenga.  Está escrito en C++ y tiene asociado un script Perl.

De lo que quería hablar hoy es de que es posible utilizar kseq.h directamente desde nuestro lenguaje preferido. Para ellos es necesario emplear herramientas que permitan incluir, y compilar, código C mezclado con el de nuestro lenguaje interpretado. En el caso de Python estaríamos hablando por ejemplo de cython , aquí os muestro cómo hacerlo en Perl con el módulo Inline::CPP, dotando a kseq de una interfaz orientada a objetos para evitar pasar tipos de datos con tipos no estándar:

 package klib::kseq;  
 use Inline CPP => config => namespace => 'klib';    
 use Inline CPP => config => libs => '-lz';  
 use Inline CPP => <<'END_CPP';  
 /*  
 The Klib library [ https://github.com/attractivechaos/klib ], is a standalone and   
 lightweight C library distributed under MIT/X11 license.   
 Klib strives for efficiency and a small memory footprint.  
 This package contains C++ wrappers for components of Klib to make them easy to use from Perl scripts.  
 */  
 using namespace std;  
 #include <iostream>  
 #include <string>  
 #include <zlib.h>  
 #include "kseq.h"   
 /*   
 kseq is generic stream buffer and a FASTA/FASTQ format parser from Klib.  
 Wrapper written by Bruno Contreras Moreira.  
 */  
 KSEQ_INIT(gzFile, gzread)  
 class kseqreader
 {  
 public:  
    kseqreader()  
    {  
       fp = NULL;  
       fromstdin=false;  
       openfile=false;  
       filename = "";  
    }  
    ~kseqreader()  
    {  
       close_stream();  
    }  
    void rewind_stream()  
    {  
       if(openfile == true)  
       {  
          if(fromstdin == false) kseq_rewind(seq);  
          else  
             cerr<<"# reader: cannot rewind STDIN stream"<<endl;  
       }  
    }           
    void close_stream()  
    {  
       if(openfile == true)  
       {  
          kseq_destroy(seq);  
          if(fromstdin == false) gzclose(fp);  
          fp = NULL;  
          fromstdin=false;  
          openfile=false;  
          filename = "";  
       }     
    }  
    // accepts both FASTA and FASTQ files, which can be GZIP-compressed  
    // returns 1 if file was successfully opened, 0 otherwise  
    int open_stream(char *infilename)  
    {  
       filename = infilename;  
       if(filename.empty())  
       {   
          cerr<<"# reader: need a valid file name"<<endl;  
          filename = "";  
          return 0;  
       }     
       else  
       {  
          FILE *ifp = NULL;  
          if(filename == "stdin")  
          {  
             ifp = stdin;  
             fromstdin = true;  
          }  
          else ifp = fopen(filename.c_str(),"r");  
          if(ifp == NULL)     
          {  
             cerr<<"# reader: cannot open input file: "<<filename<<endl;  
             filename = "";  
             return 0;  
          }  
          else // open file and initialize stream  
          {  
             fp = gzdopen(fileno(ifp),"r");   
             seq = kseq_init(fp);   
             openfile=true;  
             return 1;  
          }  
       }     
    }  
     // returns length of current sequence, which can be 0,  
    // -1 in case of EOF and -2 in case of truncated quality string  
    int next_seq()  
    {   
       return kseq_read(seq);   
    }  
     const char *get_sequence()  
    {  
       return seq->seq.s;  
    }  
    // returns FASTA/FASTQ header up to first found space  
    const char *get_name()  
    {  
       return seq->name.s;  
    }  
    const char *get_quality()  
    {  
       if(seq->qual.l) return seq->qual.s;  
       else return "";   
    }  
    const char *get_comment()  
    {  
       if(seq->comment.l) return seq->comment.s;  
       else return "";   
    }  
    const char *get_full_header()  
    {  
       _header = seq->name.s;  
       if(seq->comment.l)  
       {   
          _header += " ";  
          _header + seq->comment.s;  
       }     
      return _header.c_str();  
    }  
 private:  
    gzFile fp;  
    kseq_t *seq;  
    bool fromstdin;  
    bool openfile;  
    string filename;  
    string _header;  
 };  
 END_CPP  
 1;  

Ahora estamos en disposición de probar a leer un archivo con este módulo:

 use strict;  
 use klib;  
 my ($length,$header,$sequence,$quality);  
 my $infile = "sample.fq.gz";  
 # my $infile = "stdin";  
 my $kseq = new klib::kseq::kseqreader();  
 if($kseq->open_stream($infile))  
 {  
   while( ($length = $kseq->next_seq()) != -1 )  
   {  
     ($header,$sequence,$quality) =  
       ($kseq->get_full_header(),$kseq->get_sequence(),$kseq->get_quality());   
   }  
 }  
 $kseq->close_stream()  

Qué ganancia en tiempo de procesamiento de archivos FASTQ se obtiene? Una prueba realizada con el módulo Benchmark , con diez repeticiones y un archivo .fastq.gz con 250K secuencias sugiere que kseq::klib es el doble de rápido:

            s/iter  perlio   kseq
perlio    1.32      --     -50%
kseq    0.660  100%     --

Hasta el próximo año, un saludo,
Bruno

25 de noviembre de 2014

Una buena gráfica científica

Hola,
hoy quisiera compartir una imagen que demuestra, a pesar de toda la literatura al respecto, que no siempre es tan difícil hacer una buena gráfica con los datos de un experimento. En esta ocasión bastó con juntar todos los tomates cherry cosechados junto a la planta de la que proceden, y por supuesto una gráfica debajo para mostrar la varianza natural de estas plantas:

Fuente: http://www.nature.com/ng/journal/v46/n12/fig_tab/ng.3131_F5.html.

El artículo completo está en http://www.ncbi.nlm.nih.gov/pubmed/25362485 . Por cierto, tendrán algo de sabor esos tomates?
Hasta luego,
Bruno

12 de noviembre de 2014

Ofertas de trabajo en bioinformática

Actualización: a partir del 2 de Diciembre de 2014 las ofertas de trabajo que nos lleguen las pondremos en una página propia [http://bioinfoperl.blogspot.com.es/p/ofertas-de-trabajo.html]
separada del resto de entradas del blog.
 


1.  para desarrolladores en PerkinElmer (Granada, España)

"We have two open positions in Spain requiring bioinformatics domain expertise.

Please share the information or submit your CV at:



Thanks!"

2. staff bioinformatician  at Institute of Molecular Biology (Mainz, Alemania)

 http://www.verwaltung.personal.uni-mainz.de/Dateien/Bioinformatiker_CFBI08_docx-zi.pdf






26 de octubre de 2014

Cómo recuperar datos de un disco dañado

Primero aclarar que esta entrada será únicamente útil a los administradores de ordenadores, bioinformáticos o no, que hayan estado haciendo cambios en la distribución de su disco duro (tabla de particiones) y por algún error humano o de software hayan perdido los datos de su ordenador. Si nunca has tocado la tabla de particiones de tus discos, mejor no la toques, aunque puedes seguir leyendo...

Además el tema está de moda, porque muchos políticos y empresarios corruptos borran sus ordenadores antes de que llegue la policía... pero veremos cómo existen muchas herramientas para recuperar esos datos y meterlos en la cárcel ;)

El caso es que cuando instalamos Linux/Windows, o lo desinstalamos, o no nos arranca, o queremos crear una nueva unidad de disco dentro de nuestro propio disco, o queremos modificar las ya existentes... tenemos que modificar la tabla de particiones del disco. No explicaré lo que es la tabla de particiones, la Wikipedia lo hace mejor que yo, pero intentaré contar mi mala experiencia dar unos consejos muy útiles en caso de catástrofe.

Mi historia comienza el jueves pasado cuando me dio por ampliarme la partición de sistema (/) de mi Kubuntu, la creé de 20GB (era pequeña porque mi disco SSD es de 256GB y me gustan ligeras para poder hacer backups rápidamente), y la quería ampliar a 40GB quitando un poco de espacio a mi partición de (/home). Tenía en 10 minutos que asistir a un seminario, así que como puede tardar horas la redistribución era el momento perfecto. Y como el servidor con el que trabajo tarda mucho en arrancar, decidí hacerlo con mi portátil. Arranqué el portátil con un LiveUSB de Kubuntu 14, instalé y ejecuté GParted (ver nota 1), una aplicación que me gusta particularmente para modificar los discos con Linux y sistemas de archivos ext4 (Kubuntu también posee su propia herramienta de administración de discos). Hice los cambios de tamaño de particiones de disco que quería y pulsé la opción de 'Aplicar'. Hasta aquí todo normal, Gparted comenzó a mover datos en una partición para dejar espacio y poder amplicar la otra. Pero me di cuenta que no había conectado mi disco SSD al portátil y que estaba redimensionando mi disco interno NTFS de Windows (que también es de 256GB)... por las prisas decidí pararlo y dejarlo para después.
ERROR Nº 1: nunca detener una herramienta de administración de discos a mitad de un cambio, aunque sólo lleve 1 minuto.

Pensé que con la herramienta fdisk de Linux podía volver a asignar el tamaño original a la partición y el tipo NTFS.
ERROR Nº 2: Si hemos dañado la tabla de particiones, no hacer más cambios en ella, ni intentar formatear, si formateamos no garantizo que las herramientas que explico a continuación nos salven.

¿Qué deberíamos hacer si hemos llegado a este punto en que no podemos leer nuestro disco? Mi consejo es ir a la tienda más cercana, comprar un disco externo y hacer una copia exacta de nuestro disco 'estropeado' con un live CD de Clonezilla por ejemplo (ver nota 2). Así podemos intentar recuperar datos y estar seguros de que no puede ir la cosa a peor.

Leyendo por internet, descubrí una aplicación maravillosa y gratuita llamada TestDisk que analiza nuestro disco descubre particiones modificadas o borradas, lista los archivos que contenían y permite intentar restaurarlas. En este tutorial se describe cómo usar el programa y otras alternativas para recuperar datos. Recomiendo iniciar el ordenador con un LiveCD con TestDisk, ejecutarlo, buscar las particiones eliminadas con las opciones 'QuickSearch' y 'Advanced Search' (si la primera falla) y copiar todo a un disco externo. Posteriormente restaurar la partición con la opción 'Write'.




En mi caso encontré la partición borrada, pude copiar la mayoría de archivos del disco (algunos recientes no) y la opción de 'Write' para restaurar la partición no funcionó... El pánico había pasado, ha pero todavía no había recuperado todos mis datos. Entonces decidí probar el tan odiado chkdsk de Windows XP, aquí se explican todas sus opciones (recomiendo 'chkdsk /r').
ERROR Nº 3:  Si detenemos chkdsk a mitad de chequeo, no se arreglará el disco y al volverlo a ejecutar volverá a empezar desde el principio.

Tras 24 horas ejecutándose chkdsk y corrigiendo archivos, mi ordenador ha resucitado y parece que he recuperado casi todo, pero no todo, estoy tan feliz que he decidido escribir esta entrada :)

MORALEJA: Más vale perder unas horas haciendo copias de seguridad todos los meses, semanas o días, que tener que rezar al cielo para recuperar nuestros datos.

Tras 24 horas ejecutándose chkdsk y corrigiendo archivos, mi ordenador ha resucitado y parece que he recuperado casi todo, estoy tan feliz que he decidido escribir esta entrada :)
MORALEJA: Más vale perder unas horas haciendo copias de seguridad todos los meses, semanas o días, que tener que rezar al cielo para recuperar nuestros datos.

Otro buen consejo es usar almacenamiento redundante en la nube para nuestros datos más valiosos con servicios como Dropbox, aunque sabemos que nos los pueden espiar en cualquier momento.
 
NOTA 1: para modificar la tadesde Windows recomiendo MiniTool Partition Wizard.
 
NOTA 2: mi herramienta favorita de clonado y backup de discos es FSArchiver, pero Clonezilla es más fácil de usar con su interfaz gráfica.

10 de octubre de 2014

Apuntes de genómica microbiana

Hola,
esta semana hemos participado en un taller de 6 días organizado por Mirko Rossi (Universidad de Helsinki, Finlandia) donde hemos repasado los principales aspectos de la genómica comparativa de bacterias. Todo esto desde un punto de vista muy aplicado, dado que la mayor parte de los alumnos eran investigadores de la esfera veterinaria. Nuestra aportación fue una sesión de tres horas y media dedicada a la determinación del pangenoma de un grupo de bacterias de interés, con la ayuda de nuestro software GET_HOMOLOGUES.

Estimación de las dimensiones del core-genoma (izq) y pan-genoma (der)
de una colección de 12 plásmidos plncA/C, genomas de ejemplo en el capítulo
"Robust identification of orthologues and paralogues for microbial pan-genomics
using GET_HOMOLOGUES: a case study of pIncA/C plasmids
", parte del libro

http://www.springer.com/biomed/medical+microbiology/book/978-1-4939-1719-8 .





A parte de esto, en esta entrada quería sobre todo contar algunas de las cosas que otros profesores a los que pude escuchar contaron, por si os pueden ayudar en vuestros proyectos:
  • Con Andrey Prjibelski, del Algorithmic Biology Lab de San Petersburgo, aprendimos a ensamblar genomas bacterianos con el software SPAdes, y luego a comparar y a evaluar distintos ensamblajes con ayuda de QUAST. Nos recomendó http://nucleotid.es para comparar de una manera objetiva las virtudes y defectos de distintos ensambladores disponibles ahora mismo, aunque comentó que dada la variedad de escenarios de los usuarios y de los tipos de lecturas/reads existentes no hay ahora mismo ninguno claramente superior a los demás. Sí remarcó que algunos están mal documentados y son por tanto más difíciles de utilizar. Me quedé un poco con la idea de que Velvet y SPAdes son los mejores para genomas bacterianos. Del primero menciona que es muy estable, de calidad moderada, pero muy limitado con lecturas del tipo mate pair. El segundo es el que desarrollan ellos y fue el que usamos en la sesión. SPAdes trabaja siempre sobre una base de lecturas Illumina/IonTorrent, que por defecto se corrigen antes de ser utilizadas, pero puede hacer ensamblajes híbridos con secuencias de otras plataformas, incluso Sanger, útiles para resolver secuencias repetidas, como mate pairs convencionales o Nextera (hqmp, de alta calidad). El programa tiene además un módulo (--careful) para corregir por mapeo errores (indels,mismatches) en los contigs generados al ensamblar. Su mayor ventaja respecto a Velvet es probablemente que ensambla iterativamente con K-meros diferentes para superar las limitaciones de K pequeños (ciclos) con los grandes y viceversa (falta de solapamientos cortos, ver este blog y este otro para una comparación más amplia).  
  • Bastien Chevreux, de DSM Nutritional Products AG, filosofó acerca de distintos problemas y estrategias para secuenciar y ensamblar genomas bacterianos con su software MIRA. La diferencia fundamental de este ensamblador respecto a los anteriores es que no se basa en grafos de De Bruijn, como Velvet o SPAdes, si no en la acumulación y alineamiento de lecturas solapantes, razón por la que es menos escalable y por tanto más lento que ellos. Este software a mi me parece sin duda más complejo de usar que SPAdes, pero tiene a su favor una lista de correo muy activa y documentación extensa. MIRA se ejecuta simplemente sobre un fichero manifiesto, editado por el usuario, que puede tener una sintaxis bastante compleja, pero que para un mapeo es algo sencillo, como por ejemplo:
    project= nombre_de_proyecto
    job= mapping,genome,accurate
    
    #input data comes next
    readgroup
    is_reference
    data= /path/to/file.gbk
    strain= Ecoli_strainA
    
    readgroup
    technology= solexa
    autopairing # busca PEs por ti
    data= /path/to/reads?.fastq
    strain= Ecoli_strainK12
     
  • Rauni Kivistö, de la Universidad de Helsinki, nosexplicó el uso de los programas ACT y BRIG para visualizar datos de comparación de genomas basados en BLAST.
  •  Veronika Vonstein y Ross  Overbeek, de la Fellowship for the Integration of Genomes, presentaron su plataforma para la anotación automática de genomas RAST que se construye sobre SEED, una colección de componentes, algo parecido a operones o trozos de rutas metabólicas, que llevan curando a mano desde hace dos décadas y que les permiten proyectar en tiempo real, según se van secuenciando nuevos genomas, anotaciones de alta calidad, que tienen en cuenta el contexto genómico y no se basan en BLAST, sino en eficientes búsquedas de K-meros específicos de familias de proteínas prealineadas que llaman FIGFams. Me sorprendió que un trabajo de esta envergadura no tuviera más repercusión, de hecho yo no lo conocía, pero se explica porque durante mucho tiempo trabajaron en el sector privado. Ahora mismo tienen financiación pública y se apoyan en hardware del Argonne National Laboratory.
Creo que esto es todo, un saludo y hasta la próxima,
Bruno



29 de septiembre de 2014

La ciencia huye, el espectáculo continúa

Hace tiempo que tenía ganas de escribir una entrada de opinión en el blog y al final me he animado tras leer el artículo del periódico El País titulado "El show de la ciencia".

Resumiendo, dicho artículo más que hablar de ciencia desvela las cifras de dinero gastadas en organizar el congreso de astrofísica, show o festival Starmus en Tenerife. En la primera edición tuvo pérdidas de más de 200000 euros y este año va a costar más de 300000 euros y seguramente siga teniendo pérdidas a pesar de contar con unos 800 asistentes que pagan 300 euros para asistir a las conferencias de Stephen Hawking o algunos viejos astronautas de la NASA.

El hecho es que hace unos días asistí a un congreso de bioinformática en Sevilla con un presupuesto mucho más modesto. No importaba que pocos de los 160 científicos allí presentes fuera mundialmente conocido, la máxima gloria de muchos de nosotros es salir una vez en la vida en un periódico o televisión regional contando brevemente nuestra investigación.

Tampoco podemos compararnos con Hawking que viaja con un súbdito de 10 personas en crucero de lujo, nosotros viajamos con Ryanair y AVE tarifa promo (que si algún día perdemos el tren o el avión tenemos que pagarlo de nuestro sueldo porque la 'Administración' no entiende las condiciones de dichos billetes).

Las conversaciones en el café eran mucho más mundanas que las que pueda tener un astronauta, recuerdo frases como "ahora en mi centro de investigación no tenemos que hacer cola en el comedor", "despidieron a un científico de cada grupo", "me he pagado el viaje de mi bolsillo porque prefiero usar el poco dinero del proyecto para investigar".

No se discutió sobre la formación del universo o de la existencia de Dios. Los temas eran más mundanos, por ejemplo la medicina personalizada, muy interesante si tenemos en cuenta que la mayoría de lectores de este blog padeceremos cáncer en el futuro.

Y ¿porqué esta reflexión? porque una vez más demuestra que España es un país de bombo y pandereta, de políticos corruptos, de obras multimillonarias inservibles, de congresos de "estrellas"... pero permanecen olvidados muchos jóvenes investigadores (médicos, ingenieros...), con brillantes curriculums, con contratos temporales, en el paro o con sueldos mileuristas, a los que nunca se les dio la oportunidad de devolver a la sociedad todo lo que se invirtió en su educación.

Habrán traído a Stephen Hawking a dar un par de charlas, pero estamos olvidando que muchos de nuestros científicos en el futuro seguramente ya no las darán en nuestro país...

Saludos desde el extranjero.

23 de septiembre de 2014

apuntes de python científico

Me envía mi colega Daniel Bellomo de la UNRC, en Argentina, un enlace a una recopilación de material didáctico que están usando los alumnos de Químicas para aprender a ser autónomos en cálculo científico con software Open Source, iniciativa que apoyamos y por eso difundimos también aquí:
La idea básica es que los alumnos puedan usar python para sus trabajos
diarios y dejen de usar herramientas con licencias pagas (matlab,origin, etc).
Éste es el  temario, y los ejercicios están en el link al final de la pág:
https://github.com/pewen/cpc

Abrazo,
Daniel

9 de septiembre de 2014

contrato en mejora y genómica de cebada

El Departamento de Genética y Producción Vegetal de la Estación Experimental de Aula Dei – CSIC ofrece un contrato de 4 años para la realización de una tesis doctoral ligada al proyecto AGL2013-48756-R “Descubrimiento y utilización de la variabilidad genética que determina la adaptación de la cebada mediante herramientas genéticas y genómicas”.


El trabajo persigue el descubrimiento de caracteres y genes que ayuden a la obtención de variedades de cebada mejoradas para condiciones de sequía. Consiste en una combinación de los más modernos enfoques genéticos (diversidad intraespecífica, mapeo de QTL en poblaciones de cebada), genómicos (marcadores por métodos de secuenciación) y bioinformáticos (integración de datos de captura de exoma y otras fuentes de secuencia). El resultado de este trabajo permitirá identificar regiones cromosómicas y, potencialmente, genes de variedades tradicionales de cebada que contribuirán directamente a la mejora de variedades (el grupo también participa en un programa de mejora de variedades). Además, estos estudios permitirán  profundizar en la naturaleza de los procesos de adaptación de las plantas al clima, contribuyendo en general a la mejora de variedades para condiciones de cambio climático.

Distribución de haplotipos de HvFT1 en líneas de la Colección Nuclear de Cebadas Españolas colectas en la Península Ibérica. Tomada de http://link.springer.com/article/10.1007%2Fs00122-011-1531-x. 


Los candidatos deben tener vocación por la investigación, un buen expediente académico, excelente nivel de inglés, y deben estar admitidos en un programa de doctorado al finalizar el plazo de subsanación de la convocatoria. Se valorará tener cursado un master oficial relacionado con el tema de trabajo y/o experiencia en bioinformática.

Las solicitudes deberán ser presentadas por los candidatos del 10 de septiembre de al 26 de septiembre de 2014 a las 15:00 horas (hora peninsular española). Se ruega a las personas interesadas se pongan en contacto cuanto antes con los investigadores responsables (acasas at eead.csic.es, igartua at eead.csic.es, pgracia at eead.csic.es, bcontreras at eead.csic.es) para enviar una carta de interés y una copia del CV.

Referencias relevantes:
http://www.eead.csic.es/EEAD/barley
http://floresta.eead.csic.es/barleymap

11 de julio de 2014

blast eficiente con alias de nr

Hola,
en esta entrada quisiera mostraros una manera de buscar secuencias de proteínas dentro de la base de datos estrictamente no redundante nr del NCBI, restringiendo la búsqueda a grupos taxonómicos arbitrarios, sin necesidad de duplicar secuencias en otros archivos. De hecho, por su gran tamaño, 45 millones de secuencias a dia de hoy, desde hace tiempo nosotros ya no trabajamos con el archivo FASTA de nr, si no con la base de datos ya formateada para blast que descargamos regularmente con update_blastdb.pl desde ftp://ftp.ncbi.nlm.nih.gov/blast/db .

La receta para realizar búsquedas sobre subconjuntos de nr la hemos sacado de https://www.biostars.org/p/6528, y aquí os mostraré un ejemplo con secuencias de plantas:

1) Si no la tienes ya descárgate con el script update_blastdb.pl una copia de nr. Ojo, son una serie de archivos que en total ocupan 42G:
$ perl update_blastdb.pl nr


2) Encontrar el identificador taxonómico que mejor represente al subconjunto de nr que quieres definir, en nuestro caso Viridiplantae:
http://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=33090


3) Buscar todas las proteínas de plantas anotadas en el NCBI y exportar (Send to File => GI List) sus identificadores de GenBank (GIs) a un archivo, por ejemplo 'Viridiplantae_protein.gi.txt':
http://www.ncbi.nlm.nih.gov/protein/?term=txid33090[ORGN]



4) Hacer un alias de nr que contenga sólo proteínas de plantas:
$ ncbi-blast-2.2.29+/bin/blastdb_aliastool -gilist Viridiplantae_protein.gi.txt -db nr -out nr_plants -title nr_plants

En mi ejemplo produce el siguiente output:

Converted 3741334 GIs from Viridiplantae_protein.gi.txt to binary format in nr_plants.p.gil
Created protein BLAST (alias) database nr_plants with 2041649 sequences


Ahora veamos qué diferencias y ventajas tiene usar este subconjunto de nr en comparación con la base de datos completa:

[ NR ]
$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr

 Database: All non-redundant GenBank CDS translations+PDB+SwissProt+PIR+PRF
excluding environmental samples from WGS projects
           45,360,603 sequences; 16,206,973,370 total letters

Query= test

Length=312
                                                                      Score     E
Sequences producing significant alignments:                          (Bits)  Value

dbj|BAK07057.1|  predicted protein [Hordeum vulgare subsp. vulgare]    156    8e-47
dbj|BAJ89957.1|  predicted protein [Hordeum vulgare subsp. vulgar...  62.0    5e-10
ref|WP_014745933.1|  asparagine synthase [Tistrella mobilis] >ref  42.7    0.020
ref|WP_006720063.1|  hypothetical protein [Collinsella stercoris]...  35.0    5.0 
ref|XP_001817163.2|  hypothetical protein AOR_1_2924174 [Aspergil...  35.0    5.3 
dbj|BAE55161.1|  unnamed protein product [Aspergillus oryzae RIB40]   35.0    6.9 
ref|WP_021178774.1|  Phosphate transport system permease protein ...  34.7    7.6 
ref|WP_021807128.1|  Phosphate transport system permease protein ...  34.7    7.9 
ref|WP_026166948.1|  glucan biosynthesis protein D [Novosphingobi...  34.3    8.8

[...]

real    5m28.923s
user    5m24.564s
sys    0m3.596s

[ NR_plants ]

$ time ncbi-blast-2.2.29+/bin/blastx -query test.fna -db nr_plants

Database: nr_plants

           2,041,649 sequences; 772,729,335 total letters

Query= test

Length=312
                                                                      Score     E

Sequences producing significant alignments:                          (Bits)  Value

dbj|BAK07057.1|  predicted protein [Hordeum vulgare subsp. vulgare]    156    4e-48
dbj|BAJ89957.1|  predicted protein [Hordeum vulgare subsp. vulgar...  62.0    2e-11
gb|EMS57556.1|  hypothetical protein TRIUR3_28581 [Triticum urartu]   33.1    0.95 
ref|XP_002311179.2|  hypothetical protein POPTR_0008s05850g [Popu...  33.5    1.1  
ref|XP_004509944.1|  PREDICTED: uncharacterized protein LOC101513...  32.0    3.0  
ref|XP_004509945.1|  PREDICTED: uncharacterized protein LOC101513...  32.0    3.0  
gb|EMT15944.1|  hypothetical protein F775_17477 [Aegilops tauschii]   32.0    3.3  
dbj|BAK00277.1|  predicted protein [Hordeum vulgare subsp. vulgare]   32.0    3.7  
ref|XP_004234300.1|  PREDICTED: uncharacterized protein LOC101244...  31.2    4.7  
ref|XP_006358717.1|  PREDICTED: LRR receptor-like serine/threonin...  31.2    5.2  
ref|XP_002316273.2|  hypothetical protein POPTR_0010s20870g [Popu...  31.2    6.6  
ref|XP_003571199.1|  PREDICTED: probable tocopherol cyclase, chlo...  30.8    8.4  
ref|XP_008451998.1|  PREDICTED: respiratory burst oxidase homolog...  30.4    9.9

[...]

real    0m16.805s
user    0m16.332s
sys    0m0.372s

Como se puede ver en ambos casos las dos primeras secuencias encontradas son las mismas, con la misma puntuación en bits y distinto E-valor, al ser universos de búsqueda de tamaños  distintos.

Además,  el tiempo de búsqueda con un sólo procesador es del orden de 20x más rápido con el alias nr_plants, y a la vez el consumo de RAM en tiempo real mucho menor.

Hasta luego,
Bruno

9 de julio de 2014

regenerando DNA degenerado

Hola,
durante la reciente visita de mi colega Pablo Vinuesa al laboratorio pasamos ratos escribiendo código en Perl con el fin de diseñar de manera automática, para grandes conjuntos de secuencias de DNA previamente alineadas, parejas de primers (cebadores de PCR) adecuadas para estudios de microbiología ambiental, siguiendo principios bien conocidos ya. En cuanto probamos el código con unos cuantos ejemplos observamos que a menudo los primers para algunas regiones de los alineamientos múltiples eran degenerados. Como se muestra en la figura, eso significa que algunas posiciones de los cebadores no podían ser nucleótidos únicos si no que tenían que ser combinaciones de 2 o más para poder hibridarse con las secuencias utlizadas para su diseño.


Pareja de primers degenerados que definen un amplicón a partir de secuencias de DNA alineadas. De acuerdo con la nomenclatura IUPAC, el de la izquierda (fwd) puede escribirse como GAYTST y el de la derecha (rev) como GHGKAG. Figura tomada de http://acgt.cs.tau.ac.il/hyden.
A la hora de evaluar parejas de primers nos encontramos con que los degenerados en realidad tendrían que ser examinados por medio de cada una de las secuencias implícitas en el código degenerado. Y por tanto, tuvimos que buscar una manera de regenerar las secuencias de nucleótidos correspondientes a un primer degenerado, que es de lo que va esta entrada. El siguiente código en Perl hace precisamente esto:

 #!/usr/bin/perl  
 use strict;  
 my $degenerate_sequence = $ARGV[0] || die "# usage: $0 <degenerate DNA sequence>\n";  
 my $regenerated_seqs = regenerate($degenerate_sequence);  
 foreach my $seq (@$regenerated_seqs)  
 {  
    print "$seq\n";  
 }  
 # returns a ref to list of all DNA sequences implied in a degenerate oligo  
 # adapted from Amplicon (Simon Neil Jarman, http://sourceforge.net/projects/amplicon)  
 sub regenerate  
 {  
    my ($primerseq) = @_;  
    my %IUPACdegen = (   
    'A'=>['A'],'C'=>['C'], 'G'=>['G'], 'T'=>['T'],   
    'R'=>['A','G'], 'Y'=>['C','T'], 'S'=>['C','G'], 'W'=>['A','T'], 'K'=>['G','T'], 'M'=>['A','C'],  
    'B'=>['C','G','T'], 'D'=>['A','G','T'], 'V'=>['A','C','G'], 'H'=>['A','C','T'],   
    'N'=>['A','C','G','T']   
    );  
    my @oligo = split(//,uc($primerseq));   
    my @grow = ('');  
    my @newgrow = ('');  
    my ($res,$degen,$x,$n,$seq);  
    foreach $res (0 .. $#oligo){  
       $degen = $IUPACdegen{$oligo[$res]};   
       if($#{$degen}>0){  
          $x = 0;  
          @newgrow = @grow;  
          while($x<$#{$degen}){  
             push(@newgrow,@grow);  
             $x++;     
          }     
          @grow = @newgrow;  
       }  
       $n=$x=0;   
       foreach $seq (0 .. $#grow){  
          $grow[$seq] .= $degen->[$x];   
          $n++;  
          if($n == (scalar(@grow)/scalar(@$degen))){  
             $n=0;  
             $x++;  
          }     
       }  
    }  
    return \@grow;     
 }            

Podemos probarlo con las secuencias de la figura:

$ perl degenprimer.pl GAYTST
GACTCT
GATTCT
GACTGT
GATTGT

 y

$ perl degenprimer.pl GHGKAG
GAGGAG
GCGGAG
GTGGAG
GAGTAG
GCGTAG
GTGTAG

Hasta otra,
Bruno

17 de junio de 2014

postdoc computational systems biology Luxemburg

Postdoctoral Fellow in Computational Systems Biology at the Luxembourg Centre for Systems Biomedicine in collaboration with the Black Family Stem Cell Institute of the Mount Sinai School of Medicine

http://www.nature.com/naturejobs/science/jobs/425291-postdoctoral-fellow-in-systems-biology

The Computational Biology Group seeks a highly skilled and motivated Postdotoral Fellow to work on an exciting project on the application of network biology approaches to study the role of gene expression stochasticity in cell fate commitment. In particular, the selected candidate shall develop a computational model that integrates transcriptomics and epigenomics data to describe heterogeneity in embryonic stem cells and its implication in differentiation. Single cell and cell perturbation experiments will be performed in order to validate the predictions generated by the model. This project will be carried out in collaboration with Profs. Ihor Lemischka and Kateri Moore at the Black Family Stem Cell Institute of the Mount Sinai School of Medicine. The selected applicant will have the opportunity to visit the experimental labs of our collaborators at the Mount Sinai School of Medicine.

Requirements of the ideal candidate:

  • Ph.D. in Bioinformatics, Computer Science, Biology or a related discipline
  • Strong computational skills
  • Prior experience in mathematical modelling of biological networks, especially in network inference and analysis
  • A strong first-author publication record in the fields of Bioinformatics and Computational Biology
  • Excellent working knowledge in English.

We offer:

  • Opportunity to do applied research to medical problems within a highly dynamic research institution (LCSB) and in collaboration with internationally recognized partners
  • An exciting international environment
  • A very competitive salary

For further information, please contact:

Prof. Dr. Antonio del Sol, Luxembourg Centre for Systems Biomedicine
E-mail: antonio.delsol@uni.lu

30 de mayo de 2014

posdoc bioinformática/proteómica (Oxford)


Sir William Dunn School of Pathology, South Parks Road, Oxford
The Central Proteomics Facility at the Sir William Dunn School of Pathology, University of Oxford is looking for an enthusiastic team member experienced in computational and statistical methods to provide bioinformatics support to researchers using proteomics (http://www.proteomics.ox.ac.uk).
The successful candidate will work collaboratively with existing bioinformaticians, facility staff and academic researchers. The position will cover two main roles. The first is to advise research scientists on statistically rigorous experimental design and data analysis procedures. The second is to maintain and expand the range of computational and statistical tools the facility uses for the analysis of proteomics data. 
You will have a degree in statistics, computing, bioinformatics or equivalent with considerable demonstrable practical experience. You should have strong communication and organisational skills, experience with analysing high dimensional datasets, statistical programming and experience with at least one scripting language. While it is desirable to have previous experience with proteomics and mass spectrometry data this is not essential as training will be provided.

The post is available as a fixed-term contract for 2 years in the first instance. If you are interested in this role, please apply online. You will be required to upload a CV and supporting statement as part of your online application. 

Interviews will be held during the week commencing 23 June 2014.

http://www.jobs.ac.uk/job/AIU252/postdoctoral-bioinformatician-biostatistician

15 de mayo de 2014

Cuando Blastn no es Blastn


BLAST cambió hace ya algún tiempo a mejor con su versión BLAST+ pero por el camino se olvidó de algún detalle que puede confundir a más de uno.

El antiguo BLAST se ejecutaba con el comando 'blastall':

Blastall
--------

Blastall may be used to perform all five flavors of blast comparison. One
may obtain the blastall options by executing 'blastall -' (note the dash). A
typical use of blastall would be to perform a blastn search (nucl. vs. nucl.) 
of a file called QUERY would be:

blastall -p blastn -d nr -i QUERY -o out.QUERY

The output is placed into the output file out.QUERY and the search is performed
against the 'nr' database.  If a protein vs. protein search is desired,
then 'blastn' should be replaced with 'blastp' etc.

De esta forma un alineamiento de proteínas comenzaría como 'blastall -p blastp' y uno de ácidos nucleicos como 'blastall -p blastn' y si queremos usar MEGABLAST tenemos el comando diferente 'megablast'.

Sin embargo en la 'nueva' versión BLAST+, se separaron el alineamiento de proteínas y el de ácidos nucleicos en dos comandos: 'blastp' y 'blastn' (ver manual). Hasta aquí todo parece lógico y normal, lo que no todo el mundo sabe es que LA OPCIÓN POR DEFECTO DE BLASTN ES MEGABLAST, si ejecutamos 'blastn -help' encontraremos lo siguiente:


 *** General search options
 -task                 'megablast' 'rmblastn' >
   Task to execute
   Default = `megablast'

Y es que no todo el mundo está interesado en la velocidad de búsqueda de alineamientos, muchos de los que todavía usamos Blastn es porque apreciamos su gran sensibilidad para detectar alineamientos. En la actualidad usamos Blast para alinear miles de secuencias en tiempos muy razonables de minutos e incluso segundos. Para ganar velocidad en alineamientos de millones de secuencias existen otras mejores alternativas como Bowtie2.

La CONCLUSIÓN de todo esto, si usamos Blastn y nos interesa la sensibilidad deberemos ejecutarlo como:

 blastn -task blastn

Si lo hacemos sin añadir esta opción estaremos ejecutando MEGABLAST y correremos el peligro de perder una gran sensibilidad y no encontrar los alineamientos que deseamos. Por ejemplo, busquemos homología entre la 2'beta microglobulina humana (NM_004048.2) y la de ratón (NM_009735.3) usando la herramienta online de Blastn con las opciones por defecto ('Highly similar sequences (megablast)'):


Sin embargo si cambiamos la opción de búsqueda a 'Somewhat similar sequences (blastn)':


La diferencia es considerable, ¿no creéis? pasamos de no encontrar similaritud a un alineamiento con E-valor de 1.5E-56!!!!






9 de mayo de 2014

2nd CNB Course on Introduction to Research

Hola,
pego un anuncio de un curso para estudiantes interesados en investigar que se realizará en Madrid del  30 de Junio al 25 de Julio en el Centro Nacional de Biotecnología-CSIC, que me parece puede ser una excelente oportunidad para algunos de los seguidores del blog. Los alumnos elegidos tendrán la oportunidad de interaccionar con investigadores de todas las áreas del centro, y en concreto en el departamento de Biología de Sistemas hay varios grupos con tradición en la formación en biología computacional.




Toda la información en http://tinyurl.com/lnqtmgp.
 
Un saludo,
Bruno

25 de marzo de 2014

XII Jornadas de Bioinformática / XII Symposium on Bioinformatics

Hola,
hoy damos difusión a las pŕoximas Jornadas de Bioinformática, el mayor evento científico sobre biología computacional en España. Esta es la información que tengo de momento:


The XII Symposium on Bioinformatics (XII Jornadas de Bioinformática) will take place on 21-24 September in Sevilla, Spain, at cicCartuja (CSIC-US). 

The 21st will be the student symposium, and the main conference will start on Monday the 22nd. The URL for the meeting is: 


http://www.bioinformaticsconference2014.org/

(programme/committees are still preliminar, but will be updated periodically)
With the aim of encouraging the participation of younger bioinformaticians this year’s symposium main theme is “Bioinformatics: The New Breed”.Abstracts topics include but are not limited to:
  • Integrative Biology (NGS, -omics technologies...)
  • Structural Bioinformatics and function prediction
  • Algorithms, method, and tools development
  • Metagenomics
  • Medical Informatics 

    Abstract submission closes Thursday, July 31, 2014 

PD 25 de Agosto: nuestro laboratorio presentará dos charlas seleccionadas en las secciones de Metagenómica y Estructura y Función. Además Álvaro Sebastián, colaborador habitual del blog presentará un libro de texto sobre Bioinformática en español, en la sección The Unworkshop format.

20 de marzo de 2014

Copia de seguridad de todas nuestras bases de datos MySQL, usuarios y demás

Esta mañana me he puesto a migrar el servidor MySQL y quiero compartir un par de comandos que me han salvado horas de tiempo y muchos dolores de cabeza...

El primero permite hacer una copia de seguridad de todas nuestras bases de datos, incluso las internas de MySQL en un solo archivo:
  • mysqldump -u username -p -–all-databases > file.sql
No preocuparse si aparece el mensaje: "-- Warning: Skipping the data of table mysql.event."

Si queremos comprimir el archivo podemos ejecutar el comando anterior de la siguiente forma:
  • mysqldump -u username -p -–all-databases | gzip > file.sql.gz

Y para instalar en el nuevo servidor nuestras bases de datos y configuraciones (o para recuperar la copia de seguridad en nuestro ordenador) basta con ejecutar:
  • mysql -u username -p < file.sql

Finalmente deberemos reiniciar el servicio MySQL para recargar usuarios, permisos y demás:
  • sudo /etc/init.d/mysql restart

 Espero que estas simples indicaciones le salven la base de datos a más de uno.


27 de febrero de 2014

contrato: Anotación y diagnóstico molecular de polimorfismos en secuencias genómicas

Anotación y diagnóstico molecular de polimorfismos en secuencias genómicas

El Grupo de Biología Computacional y Estructural de la EEAD-CSIC oferta un CONTRATO de personal investigador PREDOCTORAL para la formación de doctores, renovable hasta 4 años, cofinanciado por el Gobierno de Aragón.

Plazo de solicitud finaliza el 10 de marzo de 2014.

El proyecto plantea el desarrollo de un entorno bioinformático eficiente, escalable y sencillo para el usuario final, para la anotación de secuencias genómicas obtenidas de cualquier especie y los polimorfismos observados. Por medio de algoritmos de inteligencia artificial este software deberá además aprender de las secuencias analizadas previamente para hacer predicciones de fenotipo, por ejemplo de mutaciones en un gen. Los resultados del proyecto serán directamente aplicables a los trabajos del laboratorio  en genómica de plantas y también a enfermedades humanas donde el diagnóstico molecular es una herramienta clave, como el cáncer de mama o la fibrosis quística. Para ello esta propuesta cuenta con la participación de la empresa local Blackhills Diagnostic Resources, que desarrolla este tipo de kits en Zaragoza, y que suministrará experiencia y secuencias para el adecuado desarrollo del proyecto en su vertiente clínica.

Los candidatos deben cumplir los requisitos de la convocatoria publicada en el BOA 17.02.2014 (http://tinyurl.com/nfulqbe) y estar empadronados en la Comunidad Autónoma de Aragón. Buscamos i) ingenieros o licenciados con Máster o ii) graduados con 300 créditos ECTS en Biología, Bioquímica, Biotecnología, Química, Veterinaria o Farmacia, Informática o Agronomía.

Para más información sobre el grupo consulta:
www.eead.csic.es/compbio , bioinfoperl.blogspot.com.es (este blog)

Contacto

Bruno Contreras      (bcontreras at eead.csic.es)                            
Inmaculada Yruela  (yruela at eead.csic.es)

Ubicación
Estación Experimental de Aula Dei-CSIC,
Av Montañana 1005, Zaragoza

21 de febrero de 2014

Naïve Bayes con pseudoconteos (Laplace)


Hola,
tras leer recientemente algunos capítulos del libro Doing Bayesian Data Analysis, que trae todo el código escrito en R, se me ocurrió explorar cómo construir en Perl un clasificador bayesiano ingenuo (Naïve Bayes), del que recuerdo haber leído su aplicación hace tiempo en un trabajo de Beer & Tavazoie para correlacionar secuencias reguladoras en cis y patrones de expresión de genes en levaduras.

Figura recortada de http://www.sciencedirect.com/science/article/pii/S0092867404003046



Enseguida me crucé con un módulo de CPAN (Algorithm::NaiveBayes) adecuado para entrenar un clasificador en basea  observaciones codificadas en forma de vectores de atributos y la clase a la que pertenecen:

observación = (atributo1, atributo2, atributo3 : clase)



Sin embargo,  mirando el código fuente no parece incluir la corrección de Laplace para frecuencias nulas. En esencia, esta corrección consiste en sumar a los datos empíricos, las observaciones, datos ficticios para segurarse que no sobreestimamos el modelo y para permitir clasificar muestras  con atributos no vistos en el entrenamiento, sin alejarse de lo que indican los datos reales. Si todo esto te suena a chino te redirijo a este vídeo y a este blog. Esto se ha hecho muchas veces en Bioinformática, por ejemplo para construir matrices a partir de secuencias de DNA, como en el algoritmo CONSENSUS de Hertz y Stormo.

Sigo buscando y me encuentro una versión muy compacta de clasificador Bayesiano basada en el módulo Moose, que podéis ver en http://asciirain.com/wordpress/2012/12/03/naive-bayes-in-42-lines . Tomo por tanto este código como punto de partida y lo modifico para permitir el uso de pseudoconteos de Laplace. Primero defino la clase myNaiveBayes:

 package myNaiveBayes;   
    
 # modification of http://asciirain.com/wordpress/2012/12/03/naive-bayes-in-42-lines  
 # with Laplace smoothing (K param)  
 # allows missing features, which should be encoded as 'ND'  
   
 use Moose; # turns on strict and warnings  
   
 has class_counts => (is => 'ro', isa => 'HashRef[Int]', default => sub {{}});  
 has class_feature_counts => (is => 'ro', isa => 'ArrayRef[HashRef[HashRef[Num]]]', default => sub {[]});  
 has feature_counts => (is => 'ro', isa => 'ArrayRef[HashRef[Num]]', default => sub {[]});  
 has total_observations => (is => 'rw', isa => 'Num', default => 0);  
 has K => (is => 'ro', isa => 'Num', writer => '_set_K', default => 0);  
   
 sub insert   
 {  
    # insert observation, a vector of 'features', with last element being 'class'  
    # example: ('chrome','yahoo','us','good') , where class = 'good'  
   
   my( $self, @data ) = @_;  
   my $class = pop( @data );  
   $self->class_counts->{$class}++;  
   $self->total_observations( $self->total_observations + 1 );  
       
   for( my $i = 0; $i < @data; $i++ )   
     {  
       next if($data[$i] eq 'ND');  
     $self->feature_counts->[$i]->{$data[$i]}++;  
     $self->class_feature_counts->[$i]->{$class}->{$data[$i]}++;  
   }  
     
     return $self;  
 }  
   
 sub classify   
 {  
    # takes a feature vector (a new observation) of unknown class and returns a hash reference with   
    # probabilities of being associated to all previously seen classes  
   
    my( $self, @data ) = @_;  
    my ($i,$class,%probabilities,$feature_count,$class_feature_count);  
    my ($feature_probability,$conditional_probability,$class_count,$class_probability );  
       
    printf("# classify: training data = %d , K = %d \n",$self->total_observations,$self->K);  
       
    for $class ( keys %{ $self->class_counts } ){  
       $class_count = $self->class_counts->{$class};   
       $class_probability = $class_count / $self->total_observations;   
          
       ($feature_probability, $conditional_probability) = (1.0,1.0);  
       for($i = 0; $i < @data; $i++){  
          $feature_count = $self->feature_counts->[$i]->{$data[$i]} + $self->K;  
          $class_feature_count = $self->class_feature_counts->[$i]->{$class}->{$data[$i]} + $self->K;  
         
          # if $self->K==0 (no pseudocounts) zero counts are omitted  
          next unless($feature_count && $class_feature_count);   
               
          $feature_probability *= $feature_count /   
             ($self->total_observations + (keys(%{$self->feature_counts->[$i]}) * $self->K));  
               
          $conditional_probability *= $class_feature_count /   
             ($class_count + (keys(%{$self->class_feature_counts->[$i]->{$class}}) * $self->K));   
       }  
       #                     p(class) * p(features|class)   
       # p(class|features) = ----------------------------   
       #                     p(features)  
       $probabilities{$class} = ($class_probability * $conditional_probability) / $feature_probability;  
    }  
    return \%probabilities;  
 }  
   
 no Moose;  
   
 __PACKAGE__->meta->make_immutable;  
   
 1; 

Y ahora lo probamos con los mismos ejemplos del código en que me basé, creando un programa de nombre nbayes.pl:
 #!/usr/bin/env perl  
   
 use myNaiveBayes;  
   
 use constant PSEUDOCOUNTS => 1; # Laplace smoothing  
   
 my $K = $ARGV[0] || PSEUDOCOUNTS;  
   
 my $nb = myNaiveBayes->new( K => $K);  
   
 $nb->insert('chrome' ,'yahoo'  ,'us', 'good');  
 $nb->insert('chrome' ,'slashdot','us', 'bad');  
 $nb->insert('chrome' ,'slashdot','uk', 'good');  
 $nb->insert('explorer','google' ,'us', 'good');  
 $nb->insert('explorer','slashdot','ca', 'good');  
 $nb->insert('firefox' ,'google' ,'ca', 'bad');  
 $nb->insert('firefox' ,'yahoo'  ,'uk', 'good');  
 $nb->insert('firefox' ,'slashdot','us', 'good');  
 $nb->insert('firefox' ,'slashdot','us', 'bad');  
 $nb->insert('firefox' ,'slashdot','uk', 'bad');  
 $nb->insert('opera'  ,'slashdot','ca', 'good');  
 $nb->insert('opera'  ,'yahoo'  ,'uk', 'bad');  
 $nb->insert('opera'  ,'yahoo'  ,'uk', 'bad');  
   
 my $ref_classes = $nb->classify('opera','slashdot', 'uk');  
 foreach my $class (sort { $ref_classes->{$a} <=> $ref_classes->{$b} } keys(%$ref_classes))  
 {  
   printf("%-20s : %5.5f\n", $class, $ref_classes->{$class} );  
 }  

Si lo ejecuto en el terminal obtengo:

$ perl nbayes.pl  

# classify: training data = 13 , K = 1
good                 : 0.33287
bad                  : 0.68883



Hasta luego,
Bruno


14 de febrero de 2014

trabajo en PDB Europa

The Protein Data Bank in Europe (PDBe) is seeking to recruit a 
structural analyst/programmer to join the team at the European 
Bioinformatics Institute located on the Wellcome Trust Genome Campus 
near Cambridge in the UK.

For further information and to apply, please see the full EMBL-EBI job 
listing at http://bit.ly/1jvNjIO
To discuss this post informally, please contact PDBe's Dr Sameer 
Velankar (sameer@ebi.ac.uk).

About the position:

Established in 1994, CATH and SCOP are the world?s most comprehensive 
resources classifying protein-domain structures into evolutionary 
superfamilies. They are currently being combined in a collaborative 
project - Genome3D - that aims to provide predicted 3D structures for 
sequences assigned to SCOP and CATH superfamilies. Combining SCOP and 
CATH based predictions allows us to identify more accurately regions 
that agree between the two methods. We aim to provide these 
Genome3D-predicted structures via the PDBe resource. To improve the 
assessment of the reliability of the predictions it is necessary to 
develop a mapping between SCOP and CATH and to remove any conflicts.

We are now seeking to recruit a structural analyst/programmer to assist 
with this task. The structural analyst/programmer will also build an 
automatic pipeline to generate putative domain assignments (from CATH) 
for new PDB structures, prior to classification in SCOP or CATH. The 
structural analyst/programmer will be based at PDBe, and will be jointly 
supervised by Prof. Christine Orengo (UCL, London) and Dr Alexey Muzin 
(MRC-LMB, Cambridge), together with Prof. Gerard Kleywegt and Dr Sameer 
Velankar at PDBe.

About PDBe:

The Protein Data Bank in Europe (PDBe; pdbe.org) is part of the 
Worldwide Protein Data Bank organisation (wwPDB; wwpdb.org), which 
maintains the global archive of 3D structural data on biomacromolecules. 
The PDBe team also maintains a number of databases that support 
deposition and advanced search services for structural biologists and 
the wider scientific community. The team consists of an international 
and inter-disciplinary mix of professionals (scientists and IT specialists).

Best regards,
Gary.

-- Gary Battle Protein Data Bank in Europe (PDBe) EMBL-EBI Wellcome Trust Genome Campus Hinxton, Cambridge CB10 1SD http://www.facebook.com/proteindatabank http://twitter.com/PDBeurope

11 de febrero de 2014

Jugando a leer y extraer secuencias de un fichero FASTQ comprimido con GZIP (.fq.gz)

Un fichero con extensión '.fq.gz' es un fichero en formato FASTQ comprimido en formato GZIP. Los ficheros FASTQ con datos de los nuevos métodos de secuenciación (NGS) suelen ocupar decenas de Gigabytes de espacio de disco (una cantidad indecente para la mayoría de los mortales) y comprimidos con GZIP se reducen.

Para empezar intentaremos saber cuanto espacio ocupa realmente un archivo 'fq.gz' comprimido, para ello usaremos el mismo comando 'gzip':

> gzip --list reads.fq.gz
         compressed        uncompressed  ratio uncompressed_name
        18827926034          1431825024 -1215.0% reads.fq


Parece que GZIP en vez de comprimir expande, pero no es verdad, simplemente que la opción '--list' de 'gzip' no funciona correctamente para archivos mayores de 2GB. Así que hay que recurrir a un método más lento y esperar unos minutos:


> zcat reads.fq.gz | wc --bytes
61561367168


Si queremos echar un vistazo al contenido del archivo podemos usar el comando 'less' o 'zless':

> less reads.fq.gz
> zless reads.fq.gz

Y para saber el número de secuencias que contiene simplemente hay que contar el número de líneas y dividir por 4 (le costará unos minutos):

> zcat reads.fq.gz | echo $((`wc -l`/4))
256678360

Podemos buscar una determinada secuencia en el archivo y contar cuántas veces aparece, por ej. ATGATGATG:

> zcat reads.fq.gz | awk '/ATGATGATG/ {nlines = nlines + 1} END {print nlines}'
398065

A veces nos interesará tomar un trozo del fichero para hacer pruebas, por ejemplo las primeras 1000 secuencias (o 4000 líneas):
> zcat reads.fq.gz | head -4000 > test_reads.fq
> zcat reads.fq.gz | head -4000 | gzip > test_reads.fq.gz

O extraer un rango de líneas (1000001-1000004):

> zcat reads.fq.gz | sed -n '1000001,1000004p;1000005q' > lines.fq 

También nos puede interesar dividir el fichero en varios más pequeños de por ejemplo 1 miĺlón de secuencias (o 4 millones de líneas):

> zcat reads.fq.gz | split -d -l 4000000 - reads.split
> gzip reads.split*
> rename 's/\.split(\d+)/.$1.fq/' reads.split*

Y posteriormente reunificarlos:

> zcat reads.*.fq.gz | gzip > reads.fq.gz


Finalmente os recomiendo leer un post similar en inglés:
http://darrenjw.wordpress.com/2010/11/28/introduction-to-the-processing-of-short-read-next-generation-sequencing-data/

24 de enero de 2014

breeding informatics job Dupont-Pioneer (Seville, Spain)


Buenas, 
os paso una oferta que me ha llegado por si a alguien le interesa, el contacto es  Eva Guzmán (eva dot guzman at pioneer dot com), hasta luego,
Bruno


Senior Research Associate position

Description

Provide computer and software support for researchers in crop product development programs world-wide. The senior research associate will be an integral part of the crop development team located in, and supported by, the researchers in Seville, Spain. The successful candidate will provide support for IM tools and processes and may interact with project leaders, SCRUM team(s) and members of the non-corn hybrid crop breeding staff to assist in the development and implementation of software solutions to meet research and production needs.

Provide training and client support for software applications. The successful candidate will troubleshoot and resolve technical problems primarily in support of sunflower and canola breeding in European countries. The successful applicant will work directly with breeders to become familiarized with research activities and processes relevant to the crop and region they support. They will interact directly with the app support representative in Johnston to implement software solutions, train clients worldwide on the use of software, and offer timely solutions to client IM needs. She/he will prepare and maintain training materials and provide training for new staff. 

Ensure that production bugs are routed through the recommended channels for prompt resolution and that clients’ needs are addressed in a timely and efficient manner. The successful candidate will work directly with clients and the development team(s) to promptly address urgent production issues. She/he may contribute to the establishment of a knowledge base to assist other client support representatives in addressing IM support needs, log production bugs, document bug fixes and ensure they are communicated to relevant RIM Support and Development teams.

Under the direction of the CPD RIM Support Manager, serve as the primary support contact on key initiatives and software development projects directly related to the crops and regions he/she supports. Responsibilities will include working with research clients to identify and document IM needs, working with software development staff to provide IM solutions when necessary, validating and testing newly developed software, and developing end-user documentation. Adapt and test compatibility of internal systems with locally available hardware.

Qualification

• B.S. degree in Computer Science, Bioinformatics, Information Systems, Biology, Agronomy, Molecular Genetics, Plant Breeding or related field with 2 – 4 years of industry related experience or M.Sc. degree in Computer Science, Agronomy or related field with 0 – 2 years of industry related experience or B.S.
o Science competencies;
 Experience in lab/field/greenhouse plant breeding research. Knowledge of experimental design and/or statistics. Understanding of molecular genetics and molecular breeding applications as it pertains to cultivar development.

o IM competencies;
 Excellent database, web, and/or PC skills, both functional and technical skills preferred. Experience supporting or strong familiarity with Windows Operating System. Experience with database and/or network support. Familiarity with personal computer software, client-server technology, and internet technology. Familiarity with software application to plant breeding preferred.  Demonstrated understanding of relational databases and ability to extract, summarize, and report large amounts of data.  Strong desire to contribute to the discovery and implementation of IM solutions in support of breeding activities.

o General competencies;
 Excellent written and oral communication skills and demonstrated ability to communicate complex information in a clear and concise manner.  Experience coordinating technical training and support programs to clients. Experience with customer-service preferred. Demonstrated ability to effectively manage priorities in a high-demand, time-sensitive work environment.  Demonstrated ability to work cooperatively in a diverse team environment.  Strong interpersonal and communication skills.