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/