7 de septiembre de 2017

acogemos un estudiante o postdoc del programa EMHE

Hola,
nuestro laboratorio tiene interés en acoger a un estudiante de doctorado o posdoc de Argentina, Costa Rica, Perú y Uruguay del programa EMHE “Enhancing Mobility between Latin-American and Caribean countries and Europe”. Las instrucciones para participar en la convocatoria están en http://www.csic.es/programa-emhe

Nuestra propuesta, que podréis encontrar en en las páginas web de las agencias financiadoras latinoamericanas MINCYT (ar), MICITT (cr), CONCYTEC (pe) y ANII (uy), gira en torno a la idea de que el pangenoma de una especie, la unión de todos los genomas de sus individuos, contiene variabilidad genética (codificante y no codificante) útil para estudiar el fenotipo, tanto en microorganismos como plantas. Sin embargo, construir y manejar pangenomas supone un reto. Este proyecto contribuirá a facilitar este tipo de análisis, en base a nuestra experiencia previa con https://github.com/eead-csic-compbio/get_homologues .

El candidato o candidata deberá tener conocimientos de Biología Molecular, Genómica y Biología Computacional, incluyendo experiencia real en el uso de lenguajes de programación en entornos Linux, como por ejemplo Perl, Python, R, Java o C++. Este blog es una buena muestra del trabajo en el grupo. Por favor contacta con bcontreras@eead.csic.es (https://digital.csic.es/cris/rp/rp02661).

Un saludo,
Bruno

5 de septiembre de 2017

one-liner for insert size histogram from BAM

Hola,
ayer necesitábamos obtener rápidamente un histograma con el tamaño de los insertos de una librería de lecturas/reads paired-end. Lo logramos con este one-liner que requiere haber instalado R:

$ samtools view -q 30 -F 3916 mapped_reads.bam | cut -f 9 | \
   Rscript -e 'data=abs(scan(file="stdin")); pdf("hist.pdf"); hist(data,xlab="insert size (bp)")'

$ evince hist.pdf

Lo explico por pasos:

1) Previamente habíamos alineado/mapeado las lecturas contra una referencia y convertido el alineamiento a formato BAM. Podemos hacernos una idea qué mapeos contiene el fichero mapped_reads.bam con ayuda de samtools:

$ samtools flagstat mapped_reads.bam 

497656 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 duplicates
497409 + 0 mapped (99.95%:-nan%)
497656 + 0 paired in sequencing
248828 + 0 read1
248828 + 0 read2
492706 + 0 properly paired (99.01%:-nan%)
...

2) con -q 30 le pedimos a samtools que nos devuelve solamente lecturas con calidad de mapeo (MAPQ) >= 30

3) con -F 3916 le pedimos que ignore, según podemos averiguar aquí, los siguientes tipos de lecturas: "read unmapped, mate unmapped, first in pair, not primary alignment, read fails platform/vendor quality checks, read is PCR or optical duplicate, supplementary alignment"

NOTA1: Si el fichero BAM es muy grande se puede hacer lo mismo con una muestra al azar de las lecturas, por ejemplo el 10%, con samtools view -s 0.10.

NOTA2: Si la distribución de mapeos contiene algunos insertos anormalmente grandes el histograma por defecto puede quedar demasiado ancho. En ese caso puede ser buena idea probar algo como:

hist(data[data < quantile(data,0.99)])

Hasta luego,
Bruno

28 de agosto de 2017

mean sequence length in FASTA

Hola,
para calcular la longitud promedio de las secuencias de un fichero FASTA estábamos usando el siguiente comando de Perl en el terminal, modificado de otro en https://eead-csic-compbio.github.io/perl_bioinformatica/node90.html:

$ time perl -lne 'if(/^(>.*)/){$h=$1}else{$fa{$h}.=$_} END{ foreach $h (keys(%fa)){$m+=length($fa{$h})}; printf("%1.0f\n",$m/scalar(keys(%fa))) }' ~/db/swissprot
376

real 0m4.943s
user 0m4.488s
sys 0m0.452s

Sin embargo, para ficheros de metagenomas era muy lento, y lo hemos reemplazado por este otro, que no guarda las secuencias en memoria (en un hash llamado %fa):

$ time perl -lne 'if(/^(>.*)/){$nseq++;$m+=$l;$l=0}else{$l+=length($_)} END{ $m+=$l; printf("%1.0f\n",$m/$nseq) }' ~/db/swissprot

376

real 0m2.013s
user 0m1.880s
sys 0m0.132s

Hasta luego,
Bruno

9 de agosto de 2017

semillas con licencia Open Source

Hola,

los que programamos en entornos académicos llevamos muchos años usando software distribuido con licencias de tipo Open Source (OS). De hecho, ahora mismo el canal de distribución de software más popular es GitHub, donde para alojar un nuevo proyecto debes elegir qué tipo de licencia OS le otorgas como creador.

A pesar de las controversias y del debate sobre la propiedad intelectual, este modelo ha tenido la ventaja de permitir que otros creadores aprovechen de manera eficiente los avances hechos por otros y creen soluciones mejores de las que al final todos nos beneficiamos. Una condición general de estas licencias es que el producto resultante sea también OS. Así funcionan las distribuciones Linux, por ejemplo, o gran parte del sistema operativo Mac OS X.

Variedad de tomate Sunviva del catálogo Open Source Seeds. Foto tomada de http://www.nature.com/nbt/journal/v35/n8/full/nbt0817-700.html.
Leo en Nature Biotechnology que estas licencias se están empezando a usar ahora para distribuir variedades y semillas, bajo el amparo de http://osseeds.org, como manera de acelerar y facilitar el intercambio legal de materiales vegetales. Veremos en los próximos años si también en este contexto las licencias OS son un motor de innovación, y como conviven con las empresas y cooperativas de semillas convencionales,
hasta luego,
Bruno


28 de julio de 2017

Tabla de estabilidad de sustituciones de aminoácidos

Hola,
a menudo nos preguntamos qué efecto puede tener un cambio de aminoácido sobre la estabilidad de una proteína. Una pregunta relacionada es si tendrá fenotipo una sustitución no sinónima (missense). De hecho hay toda una colección de herramientas que se han desarrollado para esto (ver por ejemplo el curso de bioinformática estructural). El trabajo reciente de Rocklin et al (2017), donde estudian la estabilidad de miles de miniproteínas (hasta 50aa) sintéticas y unas 500 naturales del PDB, nos proporciona la mejor respuesta hasta la fecha. Aunque es un artículo complejo, en esencia lo que hacen es medir la estabilidad de miles de proteínas expresándolas y exponiéndolas a proteasas en levadura.
De esa manera pueden estimar el efecto que tienen mutaciones individuales dependiendo de su contexto de estructura secundaria, ya sea una alfa-hélice, una lámina beta o un lazo o loop:
Efecto sobre la estabilidad de mutaciones en diferentes contextos estructurales. Los valores negativos, como los de la prolina en general, son desestabilizadores. Adaptada de http://science.sciencemag.org/content/357/6347/168.full.

En la siguiente figura se muestran ejemplos de mutaciones desestabilizadoras en amarillo:

Adaptada de http://science.sciencemag.org/content/357/6347/168.full.

Hasta pronto,
Bruno

4 de julio de 2017

rooting and laddering Newick trees

Buenas,
esta entrada es para compartir un script Perl para enraizar árboles, muchos árboles, de manera automática, y ordenar los nodos por distancia. Pongo el título en inglés para los buscadores.
Rubén Sancho y yo estamos probando el software GRAMPA, que requiere una colección de árboles de genes precalculados en formato Newick, pero además los quiere enraizados y ordenados de manera ascendente. Como son más de mil árboles no lo queríamos hacer uno a uno con FigTree, por poner un ejemplo; queríamos automatizarlo y para ello aprovechamos los módulos Bio::TreeIO y Bio::Phylo. El primero es de Bioperl y ya lo tenía instalado, el segundo lo instalé con: $ sudo cpanm Bio::Phylo

El árbol de ejemplo, contenido en el fichero unrooted.ph, es:

(1_Sbic:0.047,2_Osat:0.068,(((((((((((3_B422:0.007,(((4_Barb:0.000,((5_Bret:0.003,(6_BsyE:0.000,7_BsyE:0.000):0.000):0.000,8_BsyC:0.000):0.000):0.000,9_Bpho:0.002):0.000,(10_BsyC:0.000,11_BsyC:0.000):0.001):0.000):0.008,((12_B422:0.005,13_Bpho:0.005):0.002,14_Barb:0.000):0.005):0.000,((((15_Bdis:0.000,16_Bdis:0.000):0.000,17_Bmex:0.033):0.013,18_Bhyb:0.001):0.014,(19_Bboi:0.021,((20_Bmex:0.017,21_Bsta:0.034):0.012,22_Bsta:0.000):0.003):0.020):0.012):0.002,23_Barb:0.000):0.000,24_Brup:0.012):0.000,25_BsyG:0.000):0.000,(26_Bpin:0.000,27_BsyG:0.000):0.000):0.011,28_Bsta:0.000):0.053,29_BsyG:0.000):0.036,(30_Barb:0.005,(31_Bpho:0.024,(32_BsyC:0.000,33_BsyE:0.000):0.027):0.029):0.000):0.005,34_Bboi:0.030):0.035);

Invocando el script lo enraizamos con el taxón elegido, el outgroup '_Sbic':

$ perl reroot_tree.pl unrooted.ph > rooted.ph

Obtenemos el fichero rooted.ph, que podemos visualizar con FigTree:

Árbol enraizado y con orden creciente de nodos.
Éste es el código de reroot_tree.pl :

 #!/usr/bin/perl -w  
   
 # Re-roots an input Newick tree with a user-defined outgroup and   
 # prints the resulting tree in ascending or descending node order.  
 # Based on https://github.com/phac-nml/snvphyl-tools/blob/master/rearrange_snv_matrix.pl  
   
 use strict;  
 use Bio::TreeIO;  
 use Bio::Phylo::IO;  
 use Bio::Phylo::Forest::Tree;  
    
 my $OUTGROUPSTRING = '_Sbic'; # change as needed  
 my $NODEORDER = 1; # 1:increasing, 0:decreasing  
 my ($outfound,$outnode,$sorted_newick) = (0);  
   
 die "# usage: $0 <tree.newick>\n" if(!$ARGV[0]);  
   
 # read input tree  
 my $input = new Bio::TreeIO(-file=>$ARGV[0],-format=>'newick');  
 my $intree = $input->next_tree();  
   
 # find outgroup taxon  
 for my $node ( $intree->get_nodes() ) {   
  if(defined($node->id()) && $node->id() =~ m/$OUTGROUPSTRING/) {  
   $outnode = $node;   
   $outfound = 1;  
   last;  
  }  
 }  
 if($outfound == 0) {  
  die "# cannot find outgroup $OUTGROUPSTRING in input tree $ARGV[0]\n";  
 }  
   
 # root in outgroup and sort in increasing order  
 $intree->reroot_at_midpoint($outnode);  
   
 # sort nodes in defined order and print  
 my $unsorted_tree = Bio::Phylo::IO->parse(  
  '-string' => $intree->as_text('newick'),  
  '-format' => 'newick'  
 )->first();  
   
 $unsorted_tree->ladderize($NODEORDER);  
   
 $sorted_newick = $unsorted_tree->to_newick();  
 $sorted_newick =~ s/'//g;  
   
 print $sorted_newick;  


Un saludo,
Bruno

21 de junio de 2017

Leyendo FASTQ con CPAN

Buenas,
hace un par de años describía en otra entrada cómo leer de manera eficiente ficheros FASTQ con ayuda de kseq.h. Para ello definí una clase C++ que llamábamos desde Perl5 con Inline::CPP. La entrada de hoy es para contar que se puede hacer lo mismo instalando el módulo Bio::DB::HTS::Kseq desde CPAN, previa instalación en tu sistema de la librería htslib:

$ git clone https://github.com/samtools/htslib.git
$ cd htslib
$ make 
$ sudo make install
$ cd ..
$ sudo cpanm Bio::DB::HTS::Kseq

Ahora podemos escribir código como éste para probarlo:

use strict;
use Bio::DB::HTS::Kseq;

my ($length,$header,$sequence,$quality);

my $kseq = Bio::DB::HTS::Kseq->new("sample18MB.fq.gz");
my $iter = $kseq->iterator();
while(my $r = $iter->next_seq()) {

  ($header,$sequence,$quality) = ($r->name,$r->seq,$r->qual);
  print ">$header\n$sequence\n$quality\n";
}

En mis pruebas con un fichero FASTQ comprimido de 18MB, este código es un 20% más lento que la versión Inline::CPP,
hasta luego,
Bruno

16 de junio de 2017

Actualización de Perl en bioinformática

Buenas,
solamente quería comentar que el ya vetusto curso de Perl en bioinformático ha cambiado un poco de cara y ha sido revisado para traerlo al 2017, con algunas referencias a Perl6 y nuevos one-liners, que es lo que usaremos de este material en el máster de biotecnología cuantitativa de la Universidad de Zaragoza a partir del curso que viene.

Lo podéis encontrar en: https://eead-csic-compbio.github.io/perl_bioinformatica

Portada de un libro de one-liners, tomada de http://www.catonmat.net/blog/perl-one-liners-explained-part-one

Buen finde,
Bruno

11 de mayo de 2017

Jornada de agrigenómica 2017

Hola,
hoy he estado por la mañana en la 1ª jornada de agrigenómica organizada por QUAES en la UPV, en Valencia. Los 6 ponentes que hemos participado hemos hablado sobre las distintas aplicaciones de las tecnologías genómicas y de secuenciación en la mejora de cultivos. Éstas son mis notas.

Iraida Amaya nos habló de sus proyectos en torno a la mejora de caracteres en la fresa, usando marcadores SNP obtenidos con chips y por genotipado por secuenciación que clasifican por subgenoma. Estos proyectos tienen, además de genómica, metabolómica para el seguimiento de los 20 compuestos volátiles más importantes para el aroma de la fresa. Como hitos destaca la obtención de dos marcadores PCR para genotipar con un 91% de precisión fresas con aroma aceptable, así como su trabajo en curso para caracterizar por agroinfiltración FvMYB10 como regulador del color rojo de la fresa, y del color blanco de sus mutantes, que portan un cambio de marco de lectura. Para ello usaron la metodología QTLseq propuesta por Takagi, 2013.

Guillermo Marco, bioinformático de Sistemas Genómicos, nos dió un curso abreviado sobre diseño y análisis de experimentos de RNAseq, apoyándose en la revisión de Conesa, 2016. Como en un trabajo reciente nuestro (Cantalapiedra, 2017),  recomienda siempre validar las muestras por PCA para descartar réplicas y probar diferentes tuberías de análisis y estrategias de ensamblaje para elegir la más adecuada en cada caso. Su mensaje general es que el diseño del experimento RNAseq determina su éxito.

Patricia Pascual, directora de I+D de Agrícola El Bosque SL, nos presentó el diseño de su incipiente programa de mejora en mora y lo justificó como estrategia empresarial.

David Calvache, del INIA, nos habló sobre el uso de marcadores moleculares en los exámenes técnicos de nuevas variedades desde el punto de vista de la propiedad intelectual (título de obtención) y las licencias de producción y comercialización (registro de variedades). El exámen técnico permite comprobar la distinción, homogeneidad y estabilidad de una variedad candidata. Discute el uso de marcadores molecular para determinar caracteres y la idea de distancia mínima respecto a variedades existentes. Esta distancia puede ser fenotípica y/o molecular, en base a una batería de marcadores. Explica, si lo entiendo bien, que la distinción de variedades esencialmente derivadas la hace la justicia ordinaria apoyándose en distancias moleculares y pedigrís.

Esther Esteban, de la oficina española de variedades vegetales,nos recuerda los retos para la agricultura del futuro, en torno a la seguridad alimentaria, la sostenibilidad, la adaptación al cambio climático, las enfermedades emergentes, los estreses abióticos o las mejoras nutricionales. Explica las diferencias de normativa entre Europa y otros países sobre las nuevas tecnologías de mejora, incluyendo cisgénesis, transgénesis o las ZFN. Cita el trabajo de Lusser, 2011. Desde entonces han surgido nuevas tecnologías como la biología sintética y la edición mediante CRISPR/Cas9. Hace una reflexión sobre las políticas en torno a la biotecnología en Europa frente al resto del mundo.

Yo hablé sobre genómica de plantas y la aplicación de las herramientas de ultrasecuenciación en la mejora de cebada, apoyándonos en la colección nuclear de cebadas españolas, y usando ejemplos de trabajos recientes nuestros [ get_homologues-est , oídio , RNAseq en sequía , adaptación climática ] y figuras de la estupenda revisión de  Bevan, 2017.

Un saludo,
Bruno




26 de abril de 2017

Genome annotation with footprintDB

Hi,
some of you might have heard of our footprintDB collection, which is somewhat unique in that it annotates DNA motifs from different sources together with their cognate transcription factors (TF) and their interface residues. it was published in 2014 and is regularly updated and queried by users around the world, who usually perform interactive searches.

There is also a web services interface which is also quite useful, but slow if you have many sequences to scan (see examples in the manual). Things are even worse if you have a complete genome or proteome. And that's exactly what Teshome Mulugeta, who's visiting the lab from Norway, needed to do.

ACE2 DNA motif, taken from http://floresta.eead.csic.es/footprintdb/index.php?motif=cb6f6b343b895dfa1c3776c99fbedda7 .
So, we have made available FASTA files of all transcription factors in footprintDB, together with their cognate DNA motifs, at http://floresta.eead.csic.es/footprintdb/download . They come in three flavours (all, Metazoa and plants), and TF sequences look like this one:

>1:ACE2 [Saccharomyces cerevisiae] libs:JASPAR;CISBP; motif:vTGCTGGtym;mCCAGCa; url 
MDNVVDPWYINPSGFAKDTQDEEYVQHHDNVNPTIPPPDNYILNNENDDGLDNLLGMDYYNIDDLLTQELRDLDIPLVPSPKTGDGS
SDKKNIDRTWNLGDENNKVSHYSKKSMSSHKRGLSGTAIFGFLGHNKTLSISSLQQSILNMSKDPQPMELINELGNHNTVKNNNDDF
DHIRENDGENSYLSQVLLKQQEELRIALEKQKEVNEKLEKQLRDNQIQQEKLRKVLEEQEEVAQKLVSGATNSNSKPGSPVILKTPA
MQNGRMKDNAIIVTTNSANGGYQFPPPTLISPRMSNTSINGSPSRKYHRQRYPNKSPESNGLNLFSSNSGYLRDSELLSFSPQNYNL
NLDGLTYNDHNNTSDKNNNDKKNSTGDNIFRLFEKTSPGGLSISPRINGNSLRSPFLVGTDKSRDDRYAAGTFTPRTQLSPIHKKRE
SVVSTVSTISQLQDDTEPIHMRNTQNPTLRNANALASSSVLPPIPGSSNNTPIKNSLPQKHVFQHTPVKAPPKNGSNLAPLLNAPDL
TDHQLEIKTPIRNNSHCEVESYPQVPPVTHDIHKSPTLHSTSPLPDEIIPRTTPMKITKKPTTLPPGTIDQYVKELPDKLFECLYPN
CNKVFKRRYNIRSHIQTHLQDRPYSCDFPGCTKAFVRNHDLIRHKISHNAKKYICPCGKRFNREDALMVHRSRMICTGGKKLEHSIN
KKLTSPKKSLLDSPHDTSPVKETIARDKDGSVLMKMEEQLRDDMRKHGLLDPPPSTAAHEQNSNRTLSNETDAL

The header contains the internal accession number, the main TF name, the organism name, the source libraries, the DNA motifs (from JASPAR and CISBP in the example) and a URL where the full annotation and references are available,
cheers,
Bruno

25 de abril de 2017

rendimiento multihebra de BLAST+ 2.6.0

Hola,
hace unas semanas descubrí gracias a mi colega Pablo Vinuesa que BLAST+ del NCBI iba ya por la versión 2.6.0. Cuando miré el resumen de cambios me llamó la atención que ya desde la versión 2.4 soporta un algoritmo multihebra para la fase reconstrucción hacia atrás del alineamiento, que en la literatura de programación dinámica se llama traceback. Dado que nosotros usamos con mucha frecuencia BLAST quise probar en qué se traducían estas novedades en tiempo de cálculo, dado que ya habíamos observado que BLASTP no paralelizaba bien, razón por la cual desarrollamos split_blast.pl, que recientemente comparamos contra DIAMOND.

El experimento consistió en buscar alineamientos locales de 48.588 secuencias de la variedad de cebada Haruna Nijo entre los 7.927 factores de transcripción de nuestra colección http://floresta.eead.csic.es/footprintdb:

$ ncbi-blast-2.2.30+/bin/makeblastdb -in footprintdb.18042017.tf.fasta -dbtype prot

$ time ~/soft/ncbi-blast-2.2.30+/bin/blastp -query HarunaNijo_proteins.fa \
  -db footprintdb.18042017.tf.fasta -outfmt 6 -max_target_seqs 10 \
  -num_threads 8 > HarunaNijo_proteins.2.2.30.blast

real  53m47.482s
user  122m2.375s
sys 0m9.749s

$ perl split_blast.pl 8 1000 ncbi-blast-2.2.30+/bin/blastp \
  -query HarunaNijo_proteins.fa -db footprintdb.18042017.tf.fasta -outfmt 6 \
  -max_target_seqs 10 -output HarunaNijo_proteins.split.blast

# runtime: 836 wallclock secs ( 0.71 usr  0.20 sys + 6391.38 cusr  5.54 csys = 6397.83 CPU)
# this is ~14m

$ ncbi-blast-2.6.0+/bin/makeblastdb -in footprintdb.18042017.tf.fasta -dbtype prot

$ time ncbi-blast-2.6.0+/bin/blastp -query HarunaNijo_proteins.fa \
  -db footprintdb.18042017.tf.fasta -outfmt 6 -max_target_seqs 10 \
  -num_threads 8 >   HarunaNijo_proteins.2.6.0.blast
 
real  20m35.969s
user  194m1.715s
sys 2m41.827s

Como podéis ver, al menos para BLASTP esta versión de BLAST+ supone una ganancia clara en procesadores multicore (8 en esta prueba), a costa de un aumento de tamaño del binario, que pasa de 31MB a 38MB, pero sigue siendo más lento que split_blast.pl,
hasta pronto,
Bruno



18 de abril de 2017

UniPROBE heterodimers in footprintDB

Hi,
we'd like to let you know that footprintDB, our database of transcription factors (TF), cognate binding sites and interface residues has been updated. Álvaro recently added the latest version of UniPROBE which systematically annotates some TF heterodimers. Their DNA motifs look like this:
Consensus DNA motif recognized by the dimer (MXL1 , MDL1), taken from footprintDB.

The non-redundant collection of motifs has been also updated in RSAT::Plants,
have a good week,
Bruno

ciencia básica = ciencia aplicada (fármacos aprobados)

Buenas,
hoy quisiera comentar un tema recurrente cuando hablamos de ciencia: el de ciencia básica vs ciencia aplicada. En general mi impresión es que la primera se percibe como un esfuerzo romántico que da sentido a la vida de algunos locos, mientras que la segunda es la que vale, puesto que participan ingenieros, y acaba llegando a nuestro smartphone en poco tiempo. Obviamente exagero, pero por ahí van los tiros.

La excusa para sacar este tema hoy es un artículo publicado recientemente en la revista Science (http://science.sciencemag.org/content/356/6333/78.full) donde se analiza cómo se citan los proyectos de Biomedicina financiados por los NIH norteamericanos en patentes. El estudio cubre los años entre 1980 y 2007. Las conclusiones de este trabajo son:

1) que si tenemos en cuenta las citas indirectas, es decir, patentes que citan artículos que a su vez citan proyectos NIH, finannciados con dinero público, hasta un 31% de proyectos en ese periodo son citados en patentes.

2) los proyectos citados en patentes que protegen fármacos aprobados por la FDA son en igual proporción "básicos" y "aplicados", según las definiciones de los autores, que son conscientes de lo resbaladizo de estos términos.

Figura tomada de http://science.sciencemag.org/content/356/6333/78.full


Hasta luego,
Bruno

24 de marzo de 2017

Apuntes sobre ensamblaje de genomas de plantas

Buenas, ayer asistimos Ernesto Igartua y yo al 6th CNAG Symposium on Genome Research: Agrigenomics, organizado por el Centro Nacional de Análisis Genómico en Barcelona, donde a menudo contratamos servicios de secuenciación.


Allí presentamos nuestro trabajo con cebada, junto a otros colegas que trabajan en ganadería, piscicultura y agricultura y utilizan herramientas de la genómica contemporánea.

Como curiosidades me apunté que André Eggen, de Illumina, mencionó que comparando razas bovinas habían imputado SNPs mezclando genotipos de baja densidad (chips de ~10K SNPs), con genomas completos, alcanzando millones de SNPs. Por cierto, habían usado el software propietario DeNovoMAGIC para ensamblar genomas bovinos.

Otra cosa fue que los peces que estudian Franscesc Piferrer y su grupo tienen un mecanismo de metilación en función de la temperatura para controlar la producción de hormonas sexuales, algo que me recordó mucho a la memoria de vernalización en las plantas.

Pero además de estas charlas, y de visitar las salas de secuenciación y de servidores del CNAG, tuvimos dos sesiones casi seguidas donde repasamos los últimos métodos de ensamblaje y validación de genomas de plantas de la mano de Tyler Alioto y Gareth Linsmith. Éstas son mis notas:

Detección de contaminantes en las lecturas/reads
kraken : https://ccb.jhu.edu/software/kraken

Ensamblajes híbridos y diploides, combinando lecturas cortas y largas y estrategias más complejas para genomas de individuos heterocigotos.
  • reads cortos, generalmente Illumina, de entre 100 y 300b, para alcanzar profunidades de al menos 30X en cada tipo de librería: 
    • paired-end (PE) con insertos de por ejemplo 400 y 730pb 
    • mate-pair (MP) con insertos de 4 y 8Kb para superar la longitud de la mayoría de secuencias repetidas
  • reads largos, generalmente PacBio o de Oxford Nanopore. EN CNAG usan secuenciadores minIon para producir lecturas de 11.5Kb de media, alcanzando longitudes máximas > 100kb. Gareth comentó que en manzano necesitaron 60x, y eso que era material doble haploide. Este tipo de reads requieren consensos calculados con software como Sparc, Racon o Nanopolish.
En cuanto a ensambladores, Tyler destacó DISCOVAR de novo y Platanus, más adecuado para individuos con moderadas tasas de sitios heterocigotos. Pero advirtió del efecto negativo que tiene la heterocigosis sobre N50. En cambio, Gareth mencionó que primero ensambla las lecturas cortas con SOAPdenovo sin resolver las burbujas de Bruijn para luego luego combinar los reads largos con DBG2OLC y CANU.

Estrategias complementarias de ensamblaje
Datos de RNAseq para scaffolding con AGOUTI y Rascaf.

Pools de fósmidos como los empleados en el genoma de la ostra.
Mapas ópticos con enzimas nickasas que cortan cada 10Kb, con Bionano.
Dovetail genomics, aproximación basada en Hi-C.

Herramientas para corregir y finalizar genomas
PILON : https://github.com/broadinstitute/pilon/wiki
BESST : https://github.com/ksahlin/BESST

Estrategias para evaluar y validar genomas
Aparte del criterio clásico de sintenia respecto a especies cercanas, ambos mencionaron los problemas de evaluar un ensamblaje solamente por su N50 sin mirar por ejemplo los genes core anotados, por ejemplo con BUSCO, el sucesor de CEGMA. Gareth mencionó ALE para calcular la verosimilitud de un ensamblaje dadas las librerías de secuencias y KAT para comparar los k-meros originales de los reads con los del ensamblaje, que deberían coincidir, o para determinar la fracción de sitios heterocigotos:

Frecuencias de k-meros de los genotipos B73 y Mo17 de maíz, tomada de http://www.nature.com/articles/srep42444.

Casi se me olvida mencionar la comparación entre el mapa físico y el genético como criterio de calidad, muy útil en el genoma de manzano o en el de la cebada:

Comparación entre las posiciones de marcadores en una población de mapeo en cebada y sus posiciones en los mapas físico IBSC y POPSEQ de cebada, tomada de http://link.springer.com/article/10.1007%2Fs11032-015-0253-1.


Hasta  pronto,
Bruno















9 de marzo de 2017

Tutorial: pan-genome analysis with GET_HOMOLOGUES

Hi,
a new tutorial on the analysis pan-genomes using GET_HOMOLOGUES and GET_HOMOLOGUES-EST is now available. After a short introduction, where the main concepts are illustrated, the remaining sections cover the installation and typical operations required to analyze and annotate genomes and transcriptomes from a pan-genome perspective, in which individuals or species contribute genetic material to a pool.

The examples include both bacterial sequences in GenBank format and plant transcripts. This tutorial has been created for a two-day workshop to be held at BIOS (Manizales, Colombia) next week, with title "From genomes to pangenomes: understanding variation among individuals and species":



The tutorial can be found at: http://digital.csic.es/handle/10261/146411 

Code, sample datasets and documentation are available at:
https://github.com/eead-csic-compbio/get_homologues

Suggestions and error reports are welcome,
Bruno



13 de febrero de 2017

Actualización de Algoritmos3D

Hola,
sirva esta entrada para publicar la versión actualizada del material sobre "Algoritmos en Bioinformática estructural", que podéis encontrar en:

http://eead-csic-compbio.github.io/bioinformatica_estructural

Este material lo usamos anualmente en la licenciatura en ciencias genómicas de la UNAM en el campus de Cuernavaca, México.

Predicción de estructuras protéicas a partir de sus patrones de mutaciones correlacionadas en secuencias homólogas, tomada de http://science.sciencemag.org/content/355/6322/248.full. Ests tipo de protocolos se presentan en la sección http://eead-csic-compbio.github.io/bioinformatica_estructural/node35.html.

Un saludo,
Bruno

5 de enero de 2017

HS-BLASTN as replacement of MEGABLAST

Hi,
this is hopefully the last of a series of posts where I explored software choices that might replace BLAST+ programs in some scenarios. Today I'll write about HS-BLASTN, a parallel nucleotide local aligner reported to accelerate the default BLASTN algorithm (megablast), while producing the same results:


HS-BLASTN algorithm overview, taken from http://nar.oxfordjournals.org/content/early/2015/08/06/nar.gkv784.long
Note that megablast is a fast choice for intra-species comparisons and typically retrieves sequence matches with nucleotide identities greater than 70% (see Figure 3 in http://eead-csic-compbio.github.io/get_homologues/manual-est).

We will benchmark HS-BLASTN (version hs-blastn-0.0.5+) using the same Hordeum vulgare and Arabidopsis thaliana sequences used in a previous post.

$ hs-blastn index bur-0.fasta
#[IndexBuilder] Time elapsed: 89.002952 secs.
#-rw-rw-r-- 1 contrera contrera 442M Jan  5 09:59 bur-0.fasta.bwt
#-rw-rw-r-- 1 contrera contrera  13M Jan  5 09:58 bur-0.fasta.header
#-rw-rw-r-- 1 contrera contrera 116M Jan  5 09:59 bur-0.fasta.sa
#-rw-rw-r-- 1 contrera contrera  58M Jan  5 09:58 bur-0.fasta.sequence

$ ncbi-blast-2.2.30+/bin/makeblastdb -in bur-0.fasta -dbtype nucl
#Adding sequences from FASTA; added 67259 sequences in 5.51543 seconds.
#-rw-rw-r-- 1 contrera contrera  12M Jan  5 10:06 bur-0.fasta.nhr
#-rw-rw-r-- 1 contrera contrera 789K Jan  5 10:06 bur-0.fasta.nin
#-rw-rw-r-- 1 contrera contrera  15M Jan  5 10:06 bur-0.fasta.nsq

It can be seen that indexing a FASTA file with 67K sequences is about 16x slower with HS-BLASTN than with standard NCBI mableblastdb, and produces much larger index files (629Mb vs 28Mb). Now let's review search performance:
 
hs-blastn align -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.hsblastn
#[HS-BLASTN] done. Elpased time: 9.0055 secs.

hs-blastn align -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.hsblastn -num_threads 20
#[HS-BLASTN] done. Elpased time: 1.6599 secs.

time ncbi-blast-2.2.30+/bin/blastn -task megablast -query SBCC073_fLF.fasta \
  -db bur-0.fasta -evalue 0.00001 -soft_masking true -outfmt 6 -max_target_seqs 5 \
  -out SBCC073_fLF.bur-0.blastn
#real 0m33.943s

perl _split_blast.pl 20 2000 ncbi-blast-2.2.30+/bin/blastn \
  -task megablast -query SBCC073_fLF.fasta -db bur-0.fasta -evalue 0.00001 -soft_masking true \
  -outfmt 6 -max_target_seqs 5 -out SBCC073_fLF.bur-0.blastn 
# runtime:  7 wallclock secs ( 2.81 usr  0.50 sys + 60.47 cusr  7.28 csys = 71.06 CPU)

It can be seen that HS-BLASTN is ~4x faster using a single thread. About the same speedup is obtained when 20 threads are used and BLASTN is parallelized with help from an external script (_split_blast.pl). There are small nuisances, such as the fact that no soft-masking is available, or the ocasional non-stable output order of hits with same score, but it seems worth it nevertheless,
have a good weekend,
Bruno