PHMMER como alternativa a BLASTP

Hola,
hoy quería hablar de herramientas de búsqueda de secuencias de proteína, uno de los caballos de batalla en nuestro trabajo. Es un tema ya tocado en este blog, por ejemplo cuando hablamos de deltablast y cs-blast, pero que sigue siendo de actualidad.

En esta ocasión ha sido un artículo de Saripella et al el que me lo ha vuelto a recordar, comparando algoritmos que usan perfiles de secuencia (CS-BLAST, HHSEARCH and PHMMER) con algoritmos convencionales (BLASTP, USEARCH, UBLAST and FASTA). Es uno de esos artículos aburridos pero necesarios, donde se compara por vía independiente la validez de los mejores algoritmos para esta tarea, y se guardan los tiempos de cálculo empleados por cada uno de ellos (usando un core de CPU), para anotar dominios conocidos de secuencias de SwissProt extraídos de 3 repositorios complementarios (Pfam, Superfamily y CATH):

Figura original de http://bioinformatics.oxfordjournals.org/content/32/17/2636.full.

En el artículo se calculan áreas bajo curvas ROC para caracterizar a los diferentes algoritmos. De todos ellos destacaría dos:

1) BLASTP por ser muy rápido y muy preciso, con áreas de 0.908, 0.857 y 0.878 para las 3 colecciones de dominios y el tiempo que se muestra en la gráfica para buscar 100 secuencias al azar.

2) PHMMER por ser marginalmente más lento que BLASTP pero con una ganancia en área de aproximadamente el 3% (0.922,0.903,0.903).

Además, PHMMER es muy fácil de usar. Lo descargas de http://eddylab.org/software/hmmer3/CURRENT/ y solamente necesitas un archivo FASTA con la(s) secuencia(s) problema y otro con un repositorio de secuencias como SwissProt:

$ hmmer-3.1b1-linux-intel-x86_64/src/phmmer problema.faa uniprot_sprot.fasta

Hasta pronto,
Bruno

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

Alineamiento de perfiles con HHalign

Hola,
hoy discutiré un problema habitual cuando trabajamos con familias de proteínas, el de alinear entre si dos grupos de secuencias previamente alineadas (MSAs=multiple sequence alignments) , como se ilustra en el esquema siguiente:

Dadas las secuencias del grupo 1, que hemos alineado entre si con ayuda de programas como muscle, clustal o MAFFT, y las secuencias del grupo 2, que hemos procesado de manera similar, ahora queremos ver cómo se alinean todas juntas, por ejemplo para inferir diferentes funciones o estructuras.

En el caso trivial sería tan sencillo como poner todas las secuencias en el mismo archivo y hacer un alineamiento múltiple de una vez. Sin embargo, el caso que discutimos aquí es más complicado, el que se da cuando no es sencillo alinear una secuencia del grupo 1 con otra del 2, porque han divergido bastante. Entonces podemos probar a convertir el alineamiento MSA1 a un perfil y a alinearlo contra el perfil correspondiente a MSA2, como si fuera un alineamiento pareado. Podemos incorporar, además de la información evolutiva capturada en cada perfil, información de estructura secundaria para ayudar a guiar el alineamiento.

Todo esto es sencillo de realizar con el paquete de programas HHsearch, que podemos descargar de ftp://toolkit.lmb.uni-muenchen.de/HHsearch.

Los pasos a seguir son 3:

1) Obtener archivos FASTA de MSA1.faa (30 secuencias) y MSA2.faa (35 secuencias). Para mayor precisión en el alineamiento es recomendable añadir en ambos archivos una predicción de estructura secundaria, en el formato de PSIPRED. Por ejemplo, el contenido de MSA1.faa comenzaría con las siguientes líneas:

>aa_pred
MVVSIGVFDGVHI--GHQKVL...
>ss_pred
CEEEEECCCCEEH--HHHHHH...
>ss_conf
8888617775007--889999...
>seq1 (y todas las 29 siguintes a continuación)
MVVSIGVFDGVHI--GHQKVL...

2) Convertir MSA1.faa y MSA2.faa en perfiles de Markov (MSA1.hhm y MSA2.hhm), con el ejecutable hhmake:

$ ~/soft/HHsearch/hhmake -i MSA1.faa -seq 30
$ ~/soft/HHsearch/hhmake -i MSA2.faa -seq 35

3) Alinear MSA1 contra MSA2, ya sea de manera global o local:

$ ~/soft/HHsearch/hhalign -i MSA1.hhm -t MSA2.hhm -glob -seq 65 -ofas globalMSA.faa

La salida obtenida es algo como:

REMARK: in -mac -global mode -mact is forced to 0
Query file is in HHM format
Read in HMM seq1 with 467 match states and effective number of sequences = 3.1
Query file is in HHM format
Read in HMM seq31 with 576 match states and effective number of sequences = 1.8
Using maximum accuracy (MAC) alignment algorithm ...
Printing alignments in FASTA format to globalMSA.faa 

Aligned seq1 with seq31: Score = 218.04   P-value = 1     
Done

Un ejemplo de aplicación real de esta técnica de alineamiento lo tenéis en el artículo http://www.biomedcentral.com/1471-2148/10/311 .