Hola,
recientemente me vi en la necesidad de eliminar redundancia de un fichero con millones de secuencias de nucleótidos. En mi caso se trataba de secuencias repetidas del genoma de cebada (n=4.638.834), y lo intenté con CD-HIT-EST, una herramienta que he probado muchas veces, y que tenía por muy eficiente.
Sin embargo, para esta tarea, a pesar de reservar más de 10GB de RAM, no terminá en varias horas usando 20 cores.
Por tanto, me puse a buscar opciones:
- dnaclust, usa demasiada RAM
- SigClust, muy rápido, pero usa el algoritmos K-medias y por tanto le tienes que decir cuántos clusters quieres, que en mi caso es lo que quiero averiguar
- swarm, es demasiado fino separando, solamente es eficiente para eliminar copias idénticas
Mi colega Pablo Vinuesa me habló de otras opciones que se ocupan de calcular distancias entre genomas o metagenomas, un problema relacionado:
Después de mucho buscar encontré linclust, cuyo código fuente está en https://github.com/soedinglab/MMseqs2. Este programa se describió originalmente para secuencias de péptidos (artículo) pero los autores, entre ellos J Soding, añadieron después la posibilidad de agregar nucleótidos. En mis manos este es la mejor opción ahora mismo.
A continuación os muestra un banco de pruebas con secuencias repetidas de la planta Arabidopsis thaliana (n=66.752), analizadas con el binario más rápido distribuido por los autores (AVX2).
El programa lo ejecuté de esta manera:$ command time -v mmseqs/bin/mmseqs easy-linclust \
arabidopsis_thaliana.repeats.nondeg.fasta --threads 4 \
arabidopsis_thaliana.48.repeats.nr0.95 ./ --min-seq-id 0.95
Verás que guarda datos temporales en el directorio actual (./) que luego debes eliminar.
Lo que observé en mis pruebas con cebada y A. thaliana es que a partir de un umbral de identidad del 70% el programa converge:
| umbral | RAM (kbytes) |
segundos |
secuencias |
|---|---|---|---|
| 0.99 | 311164 |
4.31 |
58285 |
| 0.95 | 303764 |
4.22 |
54353 |
| 0.90 | 307804 |
4.10 | 51359 |
| 0.80 | 306648 |
4.05 | 48722 |
| 0.70 | 304736 |
3.85 |
47770 |
| 0.50 | 307140 |
3.34 |
46433 |
Hasta pronto,
Bruno
