Algunas ideas del uso de CD-HIT-EST 4.6 y USEARCH 5.2.32

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.