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
27 de febrero de 2014
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.
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:
Y ahora lo probamos con los mismos ejemplos del código en que me basé, creando un programa de nombre nbayes.pl:
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
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}'
> 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/
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/