AGRUPAMIENTO (CLUSTERING) DE
SECUENCIAS, CON BLAST Y SCRIPTS DE PERL
Manuel J. Gómez,
Protein Design Group, Centro Nacional de Biotecnología, CSIC, Cantoblanco
Ramón
Alonso-Allende, Protein Design Group, Centro Nacional de Biotecnología,
CSIC, Cantoblanco
Software que debiera estar instalado:
Objetivos de la práctica
-
Hacer un repaso sobre el uso de BLAST por linea de comandos;
-
Usar otras aplicaciones incluidas en el paquete BLAST, como blastclust,
formatdb y fastacmd;
-
Usar las aplicaciones mencionadas para identificar grupos de secuencias
parecidas a nivel de secuencia (clusteres), en el caso concreto de un proteoma
completo.
Parte 1:
Conseguir un conjunto de secuencias y crear de
una base de datos, en formato apropiado para
BLAST.
-
Crea, en tu "home directory", un directorio
llamado BLASTDATA.
-
Salva, en ese directorio, el fichero en
formato fasta con las secuencias de TODAS las proteínas predichas
para Mycoplasma pneumoniae, disponible en:
Después,
haz
una copia de esos ficheros en tu "home directory". Esto
es fundamental para evitar complicaciones después.
-
FORMATEA LA BASE DE DATOS
-
Para poder usar las colecciones de secuencias como una base de datos
usable por BLAST, hay que usar la aplicación formatdb, que
está incluida en el paquete BLAST.
-
Ejecuta, desde tu "home", el comando siguiente.
Debieran ser generados siete ficheros de indices y un fichero de incidencias
con nombre formatdb.log.
formatdb -o -i BLASTDATA/Mycoplasma_Pep.faa
Puede que, desde tu "home directory", haya que especificar el path
a
formatdb (por ejemplo: /blast/formatdb).
-
Intenta averiguar para que sirve la opción
-o (buscando en los ficheros README de BLAST, que debieran estar en el
directorio que contiene el programa BLAST, o usando formatdb ?).
-
PRUEBA A RECUPERAR SECUENCIAS DE LA BASE DE DATOS
-
La aplicación fastacmd, que también es parte de
BLAST, permite recuperar secuencias de una base de datos formateada con
formatdb.
-
Abre el fichero Mycoplasma_Pep.faa, selecciona una proteína
cualquiera y apunta su identificador (por
ejemplo: 13507741).
-
Ejecuta, desde tu "home directory", el
primero de los comandos siguientes para recuperar la secuencia del peptido
asociado al identificador que has seleccionado.
fastacmd -d BLASTDATA/Mycoplasma_Pep.faa -s identificador
Como antes, puede que haya que especificar el path completo a fastacmd.
-
Importante: es posible traerse un conjunto
de secuencias, incluyendo sus identificadores en un fichero y usando el
comando:
fastacmd -d BLASTDATA/Mycoplasma_Pep.faa -i fichero
Parte 2:
Identificacion de clusteres con blastclust.
-
IDENTIFICACION DE CLUSTERES CON blastclust
-
La aplicación blastclust hace una comparación entre
los miembros de un conjunto de secuencias, todos contra todos, haciendo
busquedas con BLAST.
-
blastclust usa como input un fichero con secuencias en formato
fasta. Los pares de secuencias entre las que existe una similitud caracterizada
por e-values menores de 1e-6 son aceptados (este es un valor por defecto,
que puede ser cambiado editando un fichero de configuracion). Al invocar
blastclust
se pueden pasar, como argumentos, otras restricciones que permiten especificar
el grado de parecido que debe de existir entre dos secuencias para que
sean consideradas como relacionadas (los parametros S, T y L, ver más
abajo). La información de pares de secuencias conectadas es usada
para construir clusteres según el principio de single linkage.
-
Dicho comando debiera generar una lista de clusteres de proteínas
de Mycoplasma pneumoniae. Los parametros -S 30 -L 0.0 -b T especifican
que para que un par de secuencias sean consideradas como vecinas y se consideren
"conectadas", el segmento alineado tiene que estar caracterizado por al
menos un 30% de identidad, y tiene que corresponder a un solapamiento de
las secuencias de al menos 0% (el mínimo posible).
-
El fichero que se DEBE usar como input
de blastclust es el fichero Mycoplasma_Pep.faa que fue copiado en
el "home directory", y no el que hay en el directorio BLASTDATA. La razon
de esto es que blastclust genera su propia base de datos y luego
la borra. Si usárais el fichero que hay en BLASTDATA, blastclust
borraría los índices generados por formatdb.
-
Abre el fichero Myco_blastclust.txt para
hacerte una idea de cuantos clusteres han sido definidos y cuantas secuencias
tiene cada uno. Cada línea (horizontal) es un cluster, y los identificadores
están separados por espacios en blanco.
-
Es importante darse cuenta de que lo que
se obtiene es una lista de clusteres, y NO una lista de adyacencia. De
hecho, la información sobre como están conectadas las proteínas
dentro de cada cluster no es devuelta por blastclust.
-
ALINEAMIENTO DE LAS SECUENCIAS DE UN CLUSTER CON ClustalW
-
Los clusteres generados por blastclust debieran de estar formados
por proteínas parecidas. Para visualizar ese parecido, y evaluar
si la similitud es global o se limita a un fragmento de las secuencias,
lo más apropiado sería hacer un alineamiento múltiple.
-
Copia los identificadores de todas las
proteínas de uno de los clusteres, y pégalas en un fichero
de texto, en forma de columna.
-
Utiliza entonces el fichero para recuperar todas las secuencias de dicho
cluster, en formato fasta, con la aplicacion fastacmd.
-
Compara las secuencias del cluster que
has escogido con ClustalW. Si observas que las secuencias no se alinean
demasiado bien, puedes probar a ejecutar de nuevo blastclust, con
unos parametros más restrictivos. Por ejemplo, aumenta la identidad
mínima de 30% a 60%.
blastclust -i Mycoplasma_Pep.faa -S 60 -L 0.0 -b T >
Myco_blastclust_60.txt
Debieras obtener MAS clusteres de temaño MAS pequeño,
formados por secuencias que, tras alinearlas con ClustalW, se vería
que son MAS parecidas.
-
A partir del mismo conjunto de secuencias en formato fasta, se pueden
recuperar todas las lineas de cabecera, simplemente usando:
para hacernos una idea de si todas las secuencias que pertenecen
a un cluster han sido anotadas de forma parecida.
Parte 3:
BLAST por linea de comandos.
-
blastclust devuelve información sobre clusteres, pero
no sobre cómo estan conectadas cada par de secuencias en los clusteres.
Para obtener esa informacion y generar clusteres de secuencias de una forma
más manual, podemos empezar por hacer directamente la comparación
de secuencias con BLAST.
-
PRUEBA A EJECUTAR BLAST POR LINEA DE COMANDOS
-
Primero, recupera una secuencia cualquiera
de la base de datos, por ejemplo 13507741 (la misma de antes), y sálvala
en un fichero. El comando sería:
fastacmd -d BLASTDATA/Mycoplasma_Pep.faa -s 13507741 >
BLASTDATA/Seq2.faa
-
Ejecuta los siguientes comandos, para comparar
la secuencia salvada anteriormente con todas las secuencias en la base
de datos:
blastall -p blastp -d BLASTDATA/Mycoplasma_Pep.faa -i
BLASTDATA/Seq2.faa -m 0
blastall -p blastp -d BLASTDATA/Mycoplasma_Pep.faa -i BLASTDATA/Seq2.faa
-m 8
Como siempre, puede que haya que especificar el path completo a
blastall.
-
la opción -p especifica el programa.
-
la opcion -d especifica la base de datos en la que se hará la
búsqueda.
-
la opción -i especifica la secuencia con que se buscará
(QUERY).
-
la opcion -m especifica el tipo de salida (0 para formato estandar,
con alineamientos; 8 para formato tabular).
blastall
para conseguir información sobre todas las opciones.
-
COMPARACION DE TODAS LAS PROTEINAS DE UN PROTEOMA ENTRE SI
-
En vez de usar una única secuencia como QUERY, se puede usar
el fichero que contiene todas las secuencias del proteoma de Mycoplasma
pneumoniae.
-
Ejecuta el siguiente comando:
blastall -p blastp -d BLASTDATA/Mycoplasma_Pep.faa -i BLASTDATA/Mycoplasma_Pep.faa
-m 8 > Myco_vs_Myco.blas
-
En el fichero Myco_vs_Myco.blast se habrán salvado
los resultados de la comparación mutua de todas las proteínas
de Mycoplasma pneumoniae en formato tabular (ejemplo).
Parte 4:
Identificación de clusteres a partir de
la salida de BLAST.
-
Escribe un script de Perl que extraiga de la tabla los identificadores
de cada par de proteínas y el e-value, y construya una LISTA DE
ADYACENCIA.
-
El programa necesitará dos argumentos: nombre del fichero que
contiene la tabla de BLAST y e-value máximo (UMBRAL).
-
Habrá que leer el fichero línea a línea.
-
Usar una expresion regular para extraer los identificadores de las secuencias
y el e-value.
-
Si los dos identificadores son iguales no usar la informacion de esa
linea.
-
Si son diferentes chequear el valor del e-value. Si es menor que el
UMBRAL, almacenar en un hash de hashes la
información del par de identificadores y el e-value:
$hash{$prot1}{$prot2} = e-value
Guardándolo de forma cruzada se obliga a que los datos sean
simétricos:
$hash{$prot1}{$prot2} = e-value
$hash{$prot2}{$prot1} = e-value
-
Antes de guardar un valor en el hash de hashes, tambien es conveniente
incluir una condición if que chequee si ya se ha guardado
previamente, de manera que el valor solo se actualice si es menor que el
ya existente.
-
Al guardar los valores en el hash de hashes se esta generando, de hecho,
una matriz de adyacencia. Dicha matriz de
adyacencia puede ahora ser imprimida en una variedad de formatos, por ejemplo,
como una lista de adyacencia (ejemplo),
o como la matriz en si misma (ejemplo).
-
Para generar la lista de adyacencia, recorrer
la tabla para recuperar los pares conectados. Para ello habrá que
hacer dos bucles anidados foreach, para extraer las claves de los
hashes.
-
Si necesitas inspiración puedes bajarte el fichero:
blast2adjLists_ed1.pl
blast2adjLists-ed1.pl se ejecuta como sigue:
blast2adjLists-ed1.pl blast_table e-value
-
Escribe otro script de Perl que identifique las componentes conexas
(o clusteres) a partir de la lista de adyacencia. El siguiente te
puede servir de base:
adja2cluster_v2.pl
-
adja2cluster.pl identifica componentes conexas a partir de una lista
de adyacencia como la generada por blast2adjLists_ed1.pl. Se usa
como sigue:
adja2cluster-v2.pl adjacency_list
-
La salida de adja2cluster_v2.pl (ejemplo)
es un fichero parecido al generado por blastclust y, de hecho, ajustando
los parámetros de forma apropiada debieran de definir los mismos
clusteres.
Parte 5:
Ultimas consideraciones.
-
Algunos de los clusteres identificados con cualquiera
de las dos estrategias revisadas en estas prácticas podrían
coincidir con familias ya caracterizadas de proteínas.
-
Por ejemplo, el cluster mayoritario identificado
a partir de la comparación con BLAST y el uso de blast2adjLists-ed1.pl
con un umbral de 1e-100, contiene
13 proteinas de Mycoplasma pneumoniae (el fichero con la lista de clusteres
está aquí). Si recuperais
alguna de las secuencias del cluster, usando fastacmd, y haceis
una búsqueda en Pfam, descubrireis que las 13 secuencias forman
parte de una familia exclusiva del genero Mycoplasma, del que se conocen
20 representantes: 13 son de M. pneumoniae y 7 en M. genitalium (PF03257)
. El dominio catalogado en Pfam es de unos 100 amino ácidos, que
es la longitud del bloque conservado que podríais haber detectado
al hacer alineamientos con ClustalW.
-
BioLayout
es un programa escrito en Java, diseñado para visualizar grafos.
-
Requiere Java 1.3, que podeis descargar de la pagina http://java.sun.com/
.
-
El input de BioLayout es muy simple: un fichero de texto plano, con
dos columnas al menos, en las que se especifican los nodos involucrados
en arcos o conexiones. Si hay una tercera columna, esta debe incluir los
pesos de las conexiones.
-
BioLayout acepta, también, tablas generadas con BLAST (deben
de estar guardadas con la extension .blast).