"Modern Statistics for Modern Biology" (libro)

Hola,
esta es mi primera entrada escrita desde el EMBL-EBI y en ella solamente quiero compartir un libro de libre acceso que se llama Modern Statistics for Modern Biology, escrito por Susan Holmes y Wolfgang Huber, que se puede visitar en https://www.huber.embl.de/msmb

Tiene una prosa sencilla y describe aproximaciones para enfrentarse a los problemas reales de la biología en general, incluyendo los que de manera habitual describimos en este blog. Además de explicar los fundamentos, el texto tiene muchos ejemplos y soluciones completas en lenguaje R. De hecho se puede descargar en http://web.stanford.edu/class/bios221/book/Rfiles el código fuente de todos los capítulos.

Un saludo,
Bruno

SequenceServer: nice local blast

Hi,
today I wanted to let you know about a tool that we discovered recently that has been very useful for us. Its name is SequenceServer (http://www.sequenceserver.com). It is simply a wrapper to let your collaborators run NCBI BLAST searches on your local sequence databases.
All you need is a copy of the NCBI folder with some BLAST+ release and a Linux distribution. Here we had ncbi-blast-2.7.1+ already in place but had to install the application, which is a ruby application, as recommended by the authors:

$ sudo gem install sequenceserver

Due to dependencies during installation I could not manage to install it in Centos5, but instead it was easy in CentOS release 7.5 (sudo yum install ruby-devel). Once this is done, and the appropriate port is open in the host, all that remains is to let the application know where the sequence databases are. You can do that with these commands:

$ sequenceserver -d /path/to/dbs  # add new databases

$ sequenceserver -l               # list installed dbs

$ sequenceserver &                # launch web application 

Now you are ready to go. Your users need to type the URL:port of your host in their browser and they can now run their searches. This is the way it looks in our server:

Cheers, Bruno

PS I will be moving to the EMBL-EBI so there might be a break in this blog, but please keep in touch

minimap2 vs BLASTN

Hola,
recientemente Heng Li publicó un trabajo (https://doi.org/10.1093/bioinformatics/bty191) describiendo un nuevo alineador genérico de nucleótidos que se llama minimap2, que podéis descargar en https://github.com/lh3/minimap2.

Figura tomada de https://doi.org/10.1093/bioinformatics/bty191

En el artículo se compara minimap2 en diferentes escenarios contra otros softwares alternativos, incluyendo su antecesor BWA mem y se destaca su velocidad y su versatilidad, ya que es capaz de alinear lecturas cortas, secuencias largas e incluso también puede alinear saltando intrones.

Yo lo que he hecho ha sido una prueba rápida para compararlo con BLASTN en el escenario habitual de GET_HOMOLOGUES-EST, donde se comparan por ejemplo todos los genes de una planta (Brachypodium distachyon) contra todos los genes de otra especie cercana (Oryza sativa). Esto es lo que he hecho:

# how many sequences
$ grep -c "^>" *fna
Bdistachyon.fna:36647
Osativa.fna:42189

# index and BLASTN search
$ ncbi-blast-2.6.0+/bin/makeblastdb -in Osativa.fna -dbtype nucl
$ ncbi-blast-2.6.0+/bin/blastn -query Bdistachyon.fna -db Osativa.fna \
   -out Bdistachyon.Osativa.blastn.tsv -dbsize 100000000 -evalue 1e-5 -outfmt 6
real	0m40.937s
user	0m40.280s
sys	0m0.636s

# index [assuming up 80% sequence identity] and minimap2 search
$ minimap2/minimap2 -x asm20 -d Oryza.mmi Osativa.fna
$ time minimap2/minimap2 Oryza.mmi Bdistachyon.fna > Bdistachyon.Osativa.minimap.paf
real	0m2.084s
user	0m3.360s
sys	0m0.300s

Ahora echemos un ojo a los alineamientos resultantes. Selecciono un par de secuencias, primero de BLASTN:

BdiBd21-3.2G0760100.1   LOC_Os01g70090.1        87.839  847     95      5       31      876     37      876     0.0     987
BdiBd21-3.2G0521100.1   LOC_Os01g37510.1        85.652  683     92      3       91      773     103     779     0.0     713

y ahora de minimap2, en formato PAF:

BdiBd21-3.2G0760100.1	876	155	776	+	LOC_Os01g70090.1	876	161	776	181	621	60	tp:A:P	cm:i:16	s
1:i:179	s2:i:0	dv:f:0.0980
BdiBd21-3.2G0521100.1	777	110	653	+	LOC_Os01g37510.1	783	122	659	87	543	60	tp:A:P	cm:i:10	s
1:i:85	s2:i:0	dv:f:0.1196

Al maneo para estos dos ejemplos podemos observar que:
i) el mejor hit de BLASTN y minimap coinciden
ii) los alineamiento de BLASTN son más largos


Hasta pronto, buenas vacaciones,
Bruno

conjunto diferencia entre listas con Perl

Hola de nuevo,
sirva esta entrada para compartir una receta eficiente para calcular el conjunto diferencia entre dos listas o arrays en lenguaje Perl5.

Conjunto diferencia, tomado de https://es.wikipedia.org/wiki/Diferencia_de_conjuntos.

Para ello podemos definar la siguiente subrutina, tomada del módulo Array::Utils:

my @a = 0..10000; 
my @b = 5000..10000; 

array_minus(@a, @b);

sub array_minus(\@\@) {
    my %e = map{ $_ => undef } @{$_[1]};
    return grep( ! exists( $e{$_} ), @{$_[0]} ); 
}

Podéis ver otras alternativas en reddit, incluyendo soluciones en Perl6 y python,
un saludo,
Bruno

estructuras del PDB parecidas en secuencia (REST)

Hola,
hace unos días un usuario preguntaba en la lista de usuarios del Protein Data Bank (PDB) cómo usar la interfaz de servicios REST, documentada en https://www.rcsb.org/pdb/software/rest.do, para hacer búsquedas BLAST.

Mientras otro usuario compartía un script escrito en Python3, llamado sequenceSimilarity.py, que requiere mmtf-pyspark para hacer consultas PSI-BLAST en tiempo real contra el PDB, a mi me llamó la atención la simplicidad del servicio sequenceCluster, que para cualquier cadena de una estructura depositada en el PDB permite obtener el clúster de secuencias con un porcentaje de identidad controlado por el usuario:

https://www.rcsb.org/pdb/rest/sequenceCluster?cluster=70&structureId=9ant.A

Esto devuelve una lista de cadenas de estructuras PDB en formato XML que podemos procesar por ejemplo como se explica en

https://developer.atlassian.com/server/fisheye-crucible/writing-a-rest-client-in-perl

Hasta luego, buen verano,
Bruno