Hola,
en
éste artículo quería compartir algunas ideas que han surgido
durante la eliminación de secuencias redundantes de un conjunto de
ESTs. Estas cadenas de cDNA, procedentes del RNA de distintos
cultivares de cebada, se obtuvieron de experimentos de clonación
seguidos de secuenciación basada en metodología Sanger. Por tanto,
es de esperar que muchas de estas secuencias sean fragmentos
incompletos de los RNA originales, y que algunas sean incluso
subsecuencias de otras. Eliminando la redundancia se debería trabajar más cómodamente que con todas ellas. Otra implicación de su origen es que no todas las bases se habrán determinado con claridad en el proceso
de basecalling, por lo que encontraremos N (bases sin
determinar) y, posiblemente, otros caracteres como S (C o G) o W (A o
T).
Hoy
día existe gran variedad de programas de alineamiento y mapeo. Desde
el archiconocido BLAST (http://blast.ncbi.nlm.nih.gov/)
hasta programas consolidados en análisis NGS,
como BWA (http://bio-bwa.sourceforge.net/); pasando por software de creciente popularidad como
Exonerate (http://www.ebi.ac.uk/~guy/exonerate/)
. Sin embargo, para la obtención de conjuntos no
redundantes de secuencias biológicas, CD-HIT (http://weizhong-lab.ucsd.edu/cd-hit/)
, del que ya hablamos en
el blog, es sin duda la elección más habitual. Es un programa muy
versátil de clustering
(o agrupación) de secuencias, que principalmente se ha venido
utilizando para analizar conjuntos de polipéptidos. De su página oficial se extrae la afirmación “CD-HIT
helps to significantly reduce the computational and manual efforts in
many sequence analysis tasks and aids in understanding the data
structure and correct the bias within a dataset.”.
Últimamente, el paquete CD-HIT se ha diversificado enormemente. La
aparición de CD-HIT-EST demanda la misma popularidad que ya tiene
CD-HIT para agrupar proteínas, en éste caso para los ácidos
nucleicos.
Todo lo detallado en éste texto se refiere a la versión 4.6 de CD-HIT, que es la más reciente a fecha 09 de Julio de 2012.
También cuenta con una serie de herramientas auxiliares para
tareas muy concretas. Una de ellas es CD-HIT-LAP, con la que se
pueden agrupar las secuencias en base a la existencia de
solapamientos entre ellas. Una de las grandes bazas del paquete
CD-HIT es que mantiene la misma filosofía en todos sus programas, y
el mismo formato de salida, por lo que es muy fácil trabajar con las
distintas utilidades una vez que ya se ha familiarizado uno con
alguna de ellas. Por otro lado, no todos los programas cuentan con el
mismo repertorio de opciones, proporcionando CD-HIT-EST un conjunto
general y amplio de opciones que le dotan de gran versatilidad; lo
cual no ocurre con, el por otro lado muy útil, CD-HIT-LAP.
Como
ha ocurrido en muchos entornos competitivos, CD-HIT parece tener de
su lado su extrema popularidad y el trasfondo que conlleva su
continuo uso por una comunidad activa. Sin embargo, su más cercano
competidor, USEARCH (http://www.drive5.com/usearch/)
(quizás más conocido como UCLUST), asegura ser superior tanto en
velocidad como sensibilidad a CD-HIT, en cuanto a clustering; y a BLAST y Mega BLAST, en cuanto a búsqueda en bases de datos
(http://drive5.com/usearch/perf/).
Una ventaja de USEARCH frente a CD-HIT es sin duda su menos ambigua y
más trabajada documentación. También proporciona una salida más
versátil y detallada que CD-HIT, pudiendo visualizar incluso los
alineamientos en formato BLAST (sin duda, una de las grandes
flaquezas de CD-HIT-EST, que pierde de esta forma mucha
transparencia). En lugar de distribuirse en distintos paquetes,
USEARCH permite ejecutar distintos algoritmos mediante distintas
opciones, como –derep_fullseq, para agrupar secuencias idénticas,
o –derep_subseq, para hacer lo mismo con subsecuencias.
En éste artículo nos estamos siempre refiriendo a la versión actual (09 de Julio, 2012), USEARCH 5.2.32.
En
la práctica, la primera diferencia que uno encuentra entre los dos
programas es que CD-HIT-EST ordena internamente las secuencias,
mientras que con USEARCH hay que hacerlo de forma explícita con el
comando –sort. Además, USEARCH cuenta con la opción –minlen,
que nos permite descartar las secuencias que no superan un tamaño
determinado, y a la que daremos específicamente el valor 0, si queremos
reproducir el comportamiento de CD-HIT-EST, en éste sentido.
Una
vez ordenadas las secuencias nos pusimos a trabajar con los comandos
de clustering.
La primera sorpresa vino con el uso de CD-HIT-EST y CD-HIT-LAP.
Ejecutando CD-HIT-EST con la opción -c 1 (identidad del 100%)
esperaríamos capturar todas las secuencias idénticas. Sin embargo,
corriendo CD-HIT-LAP sobre los resultados reveló que aún existían
secuencias iguales que no habían sido agrupadas previamente.
Poniéndonos en comunicación con el soporte de CD-HIT nos
informaron, rápidamente, de que efectivamente existe un bug
referente al alineamiento de secuencias con Ns. Curiosamente, esto no
parece ocurrir con CD-HIT-LAP.
Otro
indicio del buen hacer de CD-HIT-LAP vino con el análisis de los
resultados de USEARCH –derep_fullseq. Con éste comando se pretende
capturar todas las secuencias idénticas. De los 469.188 ESTs,
462.948 resultaban no redundantes. Todas las secuencias idénticas
así detectadas las había localizado también CD-HIT-LAP, que además
con la opción -p 1 había capturado los solapamientos en los que se
cubría completamente una de las secuencias en comparación,
rindiendo 441.098 secuencias no redundantes.
También
fue sorprendente el resultado del uso de USEARCH –derep_subseq.
Éste lleva a cabo un algoritmo heurístico, por lo que no garantiza
localizar todas las subsecuencias. Sin embargo, de la documentación
de USEARCH podemos citar: “Since
the results are equivalent, there is no reason to use --derep_fullseq
as a preprocessing step
unless you find that using --derep_subseq
directly on your input set is computationally expensive
and you will
be repeating these processing steps repeatedly with new data.”
Esto ha resultado no ser correcto, como ha confirmado el propio
Robert Edgar, autor de USEARCH, ya que algunas secuencias agrupadas
por –derep_fullseq, aparecían como clusters independientes en
–derep_subseq.
Tanto
USEARCH como CD-HIT planean subsanar estos bugs
en próximas versiones. USEARCH 6 debería aparecer en las próximas
semanas, como mínimo con la documentación corregida, mientras que
CD-HIT nos comunicaron que están trabajando en el problema de
CD-HIT-EST con las Ns y darán noticias cuando esté resuelto.
Mientras
tanto, hemos podido sacar varias conclusiones de todo esto. Por un
lado, hay que andar con pies de plomo cuando se afronta el uso de un
nuevo programa, y conviene hacer una inspección previa de algunos de
los resultados, comparandolo con lo que podríamos esperar. Además,
no siempre la documentación dice la verdad. La verdad está en el
output. Por otro lado, más concretamente, cuando se utiliza
CD-HIT-EST conviene utilizar previamente CD-HIT-LAP para eliminar
muchas de las redundancias. Análogamente, con USEARCH conviene
utilizar en primer lugar –derep_fullseq y después, en su caso,
--derep_subseq. Yo me he decantado, para alcanzar mi objetivo de
obtener un conjunto no redundante de secuencias, con criterios de
identidad bastante estrictos, por ejecutar en primer lugar el sort de
USEARCH, a continuación el CD-HIT-LAP (pues, como hemos visto,
incluye los resultados de USEARCH --derep_fullseq), y ya finalmente
una combinación de los resultados de CD-HIT-EST y USEARCH
–derep_subseq.
Esta
entrada ha servido para documentar algunos casos prácticos en el uso
de CD-HIT-EST / CD-HIT-LAP y USEARCH para eliminar redundancia de un
conjunto de ESTs. No pretende, por tanto, ser una comparativa ni una
revisión de las opciones de estos dos programas. Ambos poseen una
buena y extensa documentación, así como una cantidad ingente de
opciones y posibilidades (a saber: creación de índices, comparación
de bases de datos, generación de salidas custom, especificar la
cobertura de las secuencias por tramos, etc.) que no nos interesa
cubrir aquí.
Sin embargo, como epílogo, me gustaría comentar
algunos comandos útiles para manejar la salida estándar de estos
dos programas. En
cuanto a CD-HIT-EST, el formato .clstr está bastante establecido, e
incluso el propio USEARCH permite generar salida de éste tipo, por
motivos de compatibilidad con otras herramientas. De alguna forma, un
fichero .clstr imita al formato FASTA, con una cabecera que comienza
con “>” y el identificador de cada cluster.
De esta forma, para seleccionar las líneas con información de cada
cluster,
se puede hacer simplemente con un “grep -v '^>'”. Sin embargo,
es algo más complicado obtener las líneas solamente de aquellos
clusters
con más de un miembro. Mi solución, estando muy lejos yo de ser un
gurú de linux, consiste en buscar en primer lugar los segundos
miembros de cada cluster,
y retroceder dos líneas: “grep -B 2 '^1[^0-9]' | grep “^>”.
Más práctico me ha resultado, a la hora de interpretar los
resultados, obtener los identificadores de las secuencias que han
sido alineadas. Para ello, de forma análoga, obtengo en primer lugar
las secuencias representativas de cada cluster
con: “grep -B 1 '^1[^0-9]' | grep '^0'”, y las paso a un fichero.
Luego capturo el resto de líneas de los clusters: “grep '^[1-9]'”
y las adjunto al mismo fichero con “>>”. Ahora, con un “cat
fichero | awk '{print $3}' | sort | uniq | sed -e 's#^>##' -e
's#...$##' > tmp && mv tmp fichero” obtengo todos los
identificadores únicos de las secuencias que forman parte de
clusters con más de un miembro.
En
cuanto al formato de USEARCH, los ficheros .uc se pueden manejar de
forma análoga. Los clusters vienen identificados empezando con la
letra “C”, mientras que las secuencias representativas de cada
uno con “S”. Para obtener los alineamientos hay que seleccionar
las líneas de las secuencias redundantes, en nuestro caso se puede
hacer mediante “grep '^H'”. Como cada línea ya incluye a qué
secuencia representativa ha sido alineada, se pueden capturar todos
los identificadores pertenecientes a clusters
con más de un miembro mediante “grep '^H' | awk '{print
$9”\n”$10}' | sort | uniq”.
Una
vez tenemos los identificadores de secuencias relevantes podemos
comparar distintos resultados mediante “cat repeats_usearch | while
read line; do grep $line repeats_cdhitlap; done | wc -l”, por
ejemplo. Podemos obtener las secuencias que no coincidan por medio de
“cat repeats_usearch | while read line; do if [ $(grep -c $line
repeats_cdhitlap) -eq 0 ]; then echo $line; fi; done”.
Eso
es todo lo que os quería contar, por supuesto a la espera de
comentarios, sugerencias, consejos y correcciones.
Saludos.
