9 de julio de 2012

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.

4 de julio de 2012

Matrices de sustitución y alineamiento de secuencias

Llevo unos días redactando un texto académico sobre alineamientos y he decidido publicar la parte de matrices de sustitución PAM y BLOSUM en este blog. Alguna vez hemos hablado sobre el tema sin entrar en profundidad, pero esta vez prometo una revisión más profunda.

La historia de las matrices de sustitución se remonta a los años 70, cuando la investigadora Margaret Oakley Dayhoff se afanaba en recopilar todas las secuencias de proteína existentes en su libro 'Atlas of Protein Sequence and Structure'  (Dayhoff and Schwartz 1978). Dayhoff y colaboradores estudiaron el modelo evolutivo de los cambios en los aminoácidos de las proteínas, para ello estudiaron 1572 cambios en 71 grupos de proteínas, dentro de cada grupo las secuencias compartían más del 85% de identidad. De esta forma anotaron el número de cambios para todas las combinaciones posibles de 2 aminoácidos, observando que 35 de las posibles mutaciones nunca ocurrían, estas se correspondían con aminoácidos poco frecuentes. También observaron que las mutaciones más frecuentes se daban entre aminoácidos con similares propiedades físico-químicas, como por ej. Asp y Glu. Muchos de los cambios de aminoácido esperados por modificación de un sólo nucleótido en los codones codificantes no se daban o eran infrecuentes, lo que demostró una mayor presión evolutiva a nivel de secuencia proteica que a nivel de DNA.
 
El cambio de un aminoácido por otro se denominó 'mutación puntual aceptada' (PAM). Normalizando los datos de las PAMs de acuerdo a la probabilidad de mutación de cada aminoácido en los datos estudiados (mutabilidad) se obtuvo la famosa matriz PAM1 en la que cada elemento de la matriz M{ij} cuantifica la probabilidad de que un aminoácido i sea remplazado por otro aminoácido j en el intervalo evolutivo de 1 PAM. 1 PAM se define como el intervalo evolutivo en que cambia un 1% de los aminoácidos en el alineamiento de 2 secuencias (1 cambio o PAM por cada 100 aminoácidos).

La matriz PAM1 sirve para simular cambios evolutivos en secuencias de proteínas. Para ello basta tomar un número aleatorio (entre 0 y 1) para cada aminoácido de una secuencia dada y asignarle un cambio si la probabilidad es menor que la anotada en la matriz para conservar el aminoácido. El proceso se puede repetir múltiples veces hasta alcanzar la distancia PAM deseada. Las matrices PAM también tienen unas propiedades my interesantes: i) la matriz PAM0 sólo posee unos en la diagonal y el resto son ceros; ii) la matriz se puede multiplicar por sí misma para calcular matrices de N PAMs; iii) si la matriz se multiplica infinitas veces por sí misma obtendremos la frecuencia del aminoácido j para todas las columnas de i.

Los intervalos evolutivos medidos en PAMs los podemos relacionar con porcentajes de residuos conservados idénticos por medio de la fórmula:


Siendo f{i} la frecuencia normalizada de aparición de un aminoácido y M{ii} el valor en la diagonal de la matriz PAM. Algunas equivalencias calculadas entre identidad y PAMs se pueden consultar en la siguiente tabla:


Toda la anterior explicación teórica de las matrices PAM está muy bien, pero volviendo al tema de alinear y comparar secuencias, ¿para qué nos sirven las matrices PAM? Las matrices PAM no nos son útiles directamente, pero sí el odd-ratio (R{ij}) calculado dividiendo un elemento de la matriz M{ij} entre la frecuencia normalizada de j (f{j}):

M{ij} nos da la probabilidad de que un aminoácido i sea sustituido por otro j en una distancia evolutiva definida por la matriz PAM y f{j} es la probabilidad de encontrar el aminoácido j en una posición de la secuencia por casualidad. El odd-ratio R{ij} cuantifica la probabilidad de que una sustitución se de en una posición dada. Un odd-ratio de valor 10 significaría que la sustitución es 10 veces más frecuente que la probabilidad de encontrar alineados ambos aminoácidos. Por el contrario, un odd-ratio de valor 0.5 significaría que la probabilidad de encontrar alineados ambos aminoácidos es el doble de probable que la mutación.

Podríamos puntuar un alineamiento de dos secuencias multiplicando los odd-ratios calculados para cada posición. Sin embargo, en informática las multiplicaciones son costosas y se prefieren las sumas, así que se calcula el log-odd multiplicado por 10 de R{ij}, estos números son más intuitivos y sencillos de sumar y serán la base de las puntuaciones de los alineamientos:


Las matrices de log-odds calculados con la anterior ecuación son las que habitualmente denominamos PAM y usamos para calcular valores de similitud en alineamiento de secuencias (puntuaciones). En la siguiente figura se puede consultar la matriz PAM250, una de las más usadas para puntuar alineamientos:
  
Si queremos encontrar un significado probabilístico de los valores log-odd de una matriz, bastaría con volver a calcular el odd-ratio (R{ij}):
Otras nuevas versiones de las matrices PAM han sido calculadas con un número mayor de grupos de secuencias homólogas alineadas, sin embargo no han conseguido mejorar sustancialmente las matrices originales de Dayhoff  (Gonnet, Cohen et al. 1992Jones, Taylor et al. 1992).

Otro tipo de matrices de sustitución que sí han conseguido mejorar a las PAM son las matrices BLOSUM (BLOcks of Amino Acid SUbstitution Matrix), creadas por Henikoff  (Henikoff and Henikoff 1992). Las matrices BLOSUM fueron creadas a partir de datos de más de 500 grupos de alineamientos de secuencias de proteínas y con el objetivo de mejorar los alineamientos de secuencias divergentes donde las matrices PAM fallaban. Para definir diferentes matrices BLOSUM se marcaron diferentes umbrales de identidad de secuencias, de forma que las secuencias con mayor o igual identidad que el umbral se agruparon para disminuir su contribución en la matriz. Por ejemplo, para calcular la matriz BLOSUM62 se agruparon las proteínas con identidad mayor o igual que 62%. Con los bloques de secuencias alineadas se calcula una tabla de frecuencias de cada pareja de aminoácidos alineados, obteniendo 210 parejas posibles con sus respectivas frecuencias de aparición que permitirán calcular los (R{ij}) entre las frecuencias observadas (q{ij}) y las frecuencias esperadas por casualidad (e{ij}):
Henikoff decidió calcular los log-odds (R{ij}) de una manera ligeramente diferente a Dayhoff, usando logaritmos en base 2:
En la siguietne figura se representa la matriz BLOSUM62, ésta es la matriz preferida para usar por defecto por algoritmos tan famosos como BLASTP.
Las matrices BLOSUM demostraron ser más sensibles a la hora de identificar alineamientos de proteínas homólogas (Henikoff and Henikoff 1992). Las principales diferencias entre ambos tipos de matrices es que las PAM son generadas por extrapolación de datos de alineamientos de secuencias muy conservadas y las BLOSUM, por contra, son derivadas de datos reales de alineamientos de secuencias menos conservadas. A continuación se muestra la equivalencia entre diferentes matrices PAM y BLOSUM, a menor distancia evolutiva PAM, mayor porcentaje de identidad BLOSUM y al contrario:
Equivalencia de matrices PAM y BLOSUM
Como norma general se prefiere el uso de matrices BLOSUM, sin embargo, cuando se realizan comparaciones de secuencias muy conservadas, las matrices PAM pueden conseguir mejores resultados.

Todo lo explicado hasta ahora sobre matrices de sustitución ha sido en el contexto de alineamientos proteicos. ¿Qué sucede en el caso de alineamientos de secuencias de DNA o RNA? Para los nucleótidos también se han calculado matrices PAM de forma similar a la explicada para proteínas (States, Gish et al. 1991), teniendo en cuenta las diferentes probabilidades de mutaciones por transición (A<->G, C<->T/U) o transversión (A/G<->C/T/U). Sin embargo, programas como BLAST emplean por defecto puntuaciones de 1 y -2 para evaluar coincidencia/no coincidencia de nucleótidos respectivamente. Aunque el uso de matrices PAM puede mejorar alineamientos de nucleótidos con identidades <70%, normalmente su mayor sensibilidad no compensa el mayor tiempo necesario para realizar los alineamientos, especialmente cuando estamos trabajando con genomas. Cuando se requiere alinear secuencias de DNA o RNA divergentes se prefiere traducirlas a secuencias proteicas antes de realizar su alineamiento.