Probando deltablast

Hola,
hace tiempo que tenía pendiente escribir sobre una de las ramas más recientes de la familia de programas BLAST, que se llama deltablast, publicada el año pasado. La lectura del artículo original sugiere que este programa se desarrolló a partir de la publicación del algoritmo CS-BLAST (context specific BLAST).
Pero realmente la historia no comienza aquí.
Todo empezó con el algoritmo PSI-BLAST, una versión iterativa de BLASTP, que en vez comparar una secuencia de proteína contra una librería de secuencias, construye un perfil de la secuencia problema y ése es el que compara contra la librería, ganando en sensibilidad. Todo esto con un coste en tiempo de cálculo modesto.
Qué problema tiene esto? Pues que el perfil se construye en tiempo real con los homólogos encontrados en iteraciones previas de manera automática y en ocasiones se pueden contaminar y producir resultados no idóneos.

(tomada de http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2645910)

CS-BLAST, que en realidad es como un preprocesamiento de BLAST, logró superar estos problemas precalculando los perfiles, que capturan información del contexto, y de hecho mejoró la sensibilidad a la hora de detectar homólogos remotos, que han perdido claramente la similitud de secuencia.

Pues bien, deltablast es un programa del NCBI que se inspira en estas experiencias para superar a capacidad de búsqueda de CS-BLAST, apoyándose en la colección de dominios Conserved Domains Database (CDD):

(tomada de http://www.biology-direct.com/content/7/1/12)

El programa es muy sencillo de usar en la web, aquí os explico como probarlo localmente en vuestra máquina:
1) descarga la última versión de BLAST+ para tu arquitectura de
ftp://ftp.ncbi.nih.gov/blast/executables/LATEST
 
2) descomprime el software en una carpeta, por ejemplo soft/

3) descarga una copia de CDD de ftp://ftp.ncbi.nih.gov/blast/db/cdd_delta.tar.gz

4) descomprime la base de dominios en un carpeta, de nuevo por ejemplo en  soft/

5) elige una colección de secuencias protéicas contra la que buscar, por ejemplo soft/nr.faa y formatéala con:
 $ soft/ncbi-blast-2.2.27+/bin/makeblastdb nr.faa
 
6) Ya puedes buscar una secuencia de proteína, como ejemplo.faa, contra tu colección de secuencias:
$ soft/ncbi-blast-2.2.27+/bin/deltablast -query ejemplo.faa -db soft/nr.faa -rpsdb soft/cdd_delta

7) Puedes utilizar prácticamente los mismos parámetros que con BLASTP, que puedes consultar haciendo:
$ soft/ncbi-blast-2.2.27+/bin/deltablast -help

Suerte!
Bruno

Filogenias universales 16S

Un obituario en el último número de Science (8 Feb 2013) me ha hecho recordar el tremendo impacto que ha tenido el trabajo del desaparecido Carl Woese en nuestra visión de la historia evolutiva, ya que a partir de sus logros pudimos comparar por primera vez con precisión molecular las historias evolutivas de linajes tan distantes como las bacterias y los animales. Además de descubrir por el camino nada menos que a las arqueobacterias, el uso de las secuencias de los genes ribosomales 16S como fuente de información filogenética todavía es rutinario, a pesar de sus limitaciones y de los avances tecnológicos, y a pesar de que ahora se hable más de redes que de árboles. Si no que se lo pregunten a los que se dedican a la metagenómica...

En las siguientes figuras se resume el impacto de la introducción de los análisis de genes 16S:

Árbol filogenético de Haeckel (tomado de http://mmbr.asm.org/content/51/2/221.long).

Árbol de distancias basado en el alineamiento de 260 secuencias de rRNA 16S (tomado de http://mmbr.asm.org/content/51/2/221.long).

Al releer la revisión de Woese del año 1987 lo que más me ha sorprendido ha sido ver que no sólo se usó secuencias, sino que se molestó en analizar en detalle la estructura de los RNAs ribosomales, e incluso localizó a nivel de estructura secundaria las zonas que permitían distinguir los grandes linajes evolutivos, como se muestra en la siguiente figura. Está ya todo inventado también aquí?

(tomada de http://mmbr.asm.org/content/51/2/221.long)

Un saludo y buen finde,
Bruno

Tutorial en el congreso BIFI2013

Hola,
nuestro laboratorio (www.eead.csic.es/compbio) coordinará un tutorial el viernes 1 de Febrero a las 15.00 horas como parte del congreso BIFI 2013. El título es:

Taller de análisis bioinformático de proteínas reguladoras

Resumen: En este taller veremos cómo estudiar por medio de herramientas bioinformáticas propiedades de las proteínas reguladoras, como las regiones intrínsecamente desordenadas, su interfaz de reconocimiento de ADN y la predicción de posibles elementos cis.

Duración: 2 horas

Cómo llegar: http://goo.gl/maps/RFJmj 

Un saludo,
Bruno

PD Un artículo nuestro reciente donde estudiamos el desorden en plantas: http://www.biomedcentral.com/1471-2229/12/165

Indexando archivos FASTA

Hola y feliz año nuevo con retraso!
Recientemente me he visto en la necesidad de acceder de maneara aleatoria a secuencias de un archivo FASTA de gran tamaño, en mi caso cercano a los 2 Gb,
con las secuencias formateadas en una sola línea.

Primero probé con las herramientas de BLAST+, en concreto:
$ makeblastdb -in archivo.fasta -parse_seqid

Con la idea de posteriormente consultar las secuencias con algo como (más ejemplos aquí):

$ blastdbcmd -query -db archivo.fasta

Sin embargo, la generación del índice se hizo eterna y de hecho la aborté, por problemas que desconozco y no he seguido mirando.

Posteriormente, tras buscar en Google encontré dos caminos en Perl:

1)  con ayuda de scripts de Bioperl: bp_index.pl y bp_fetch.pl

2) con ayuda del módulo core Tie::File, como muestro en el ejemplo a continuación, mi solución preferida, un saludo:

 #!/usr/bin/perl -w  
 use strict;  
 use Tie::File;  
   
 my $fasta_file = '/path/to/file.fasta';  
   
 my (%index_fasta,@fasta_array);  
 open(FASTA,$fasta_file) ||   
    die "# cannot open $fasta_file\n";  
 while(<FASTA>)  
 {  
    if(/^> any typical marker in header, such as gi (\S+)/)  
    {  
       $index_fasta{$1} = $.;   
    }  
 }  
 close(FASTA);   
 printf("# indexed %d sequences in file %s\n\n",  
    scalar(keys(%index_fasta)),$fasta_file);  
   
 tie(@fasta_array,'Tie::File',$fasta_file) ||   
    die "# cannot tie file $fasta_file\n";  
   
 print "ejemplo: secuencia con identificador 'id12':\n";  
 print "$fasta_array[$index_fasta{'id12'}-1]\n".  
    "$fasta_array[$index_fasta{'id12'}]\n";  

Bioinformatic Merry Christmas!!!

My boss is very nervous when I prepare this blog entry at the end of December, and this 'crisis' year I wanted to prepare the best bioinformatic inspirated Christmas card ever..

If you are not as geek as us, here is written in flat words... MERRY CHRISTMAS AND HAPPY NEW YEAR 2013!!!

If you want to replicate the Christmas greeting, go to this link, or to the 'retrieve' section in UniProt and introduce the following protein identifiers: P0A183, O02437, Q40309, Q6L5Z2, E0ZPQ6.