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
Muchas gracias Bruno por compartir esta información. No conocía linclust del paquete MMseqs2. Está interesante.Saludos.
ResponderEliminarLo que comentas del umbral del 70% asumo se debe a saturación mutacional creciente y consiguiente caída en la eficiencia del filtrado, ¿es correcto?
ResponderEliminarEs saturación pero también consecuencia del algoritmo; mira lo que hay en el fichero de revisión por pares, que puedes ver en https://static-content.springer.com/esm/art%3A10.1038%2Fs41467-018-04964-5/MediaObjects/41467_2018_4964_MOESM2_ESM.pdf:
EliminarThe only remaining point I have is that, the authors should acknowledge in the paper that (in the authors's own words, from their point-by-point reply "neither Linclust, MMseqs2, nor UCLUST nor any other tool we know of is suitable to cluster down to 70% or 50% sequence identity, since they all use fast prefilters that can miss some similar pairs".