eliminar redundancia en grandes ficheros de nucleótidos (linclust)

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:

• mash
• dRep
  

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

agrupando secuencias de proteínas con Linclust

Hola,
de vez en cuando tengo que revisar un viejo script para actualizar mi copia local del Protein Data Bank (PDB). El programa descarga solamente las estructuras que han cambiado mediante rsync y otros ficheros de un servidor FTP.
Sin embargo, las rutas a las respectivas carpetas van cambiando y yo tengo que actualizarlas. En concreto, hoy habían cambiado las listas de secuencias no redundantes, que ahora se pueden encontrar en ftp://resources.rcsb.org/sequence/clusters

Leyendo descubro que en el PDB ahora agrupan sus secuencias usando MMseq2 / Linclust, dos métodos relacionados que calculan de manera muy eficiente la similitud entre secuencias a partir de su composición de K-meros con un alfabeto reducido, temas de los que ya hemos hablado por ejemplo aquí y aquí. Me centraré en Linclust.

Algoritmo de clustering de coste lineal. Fuente: https://www.nature.com/articles/s41467-018-04964-5

Según el banco de pruebas publicado por sus autores, a diferencia de otras alternativas, el algoritmo Linclust tiene un coste lineal pero un comportamiento parecido, con pérdidas controlables de sensibilidad. Consta de varias fases:

1. Transformación de las secuencias originales a una alfabeto reducido de 13 letras. Obtienen resultados óptimos haciendo las siguientes simplificaciones: (L, M), (I, V), (K, R), (E, Q), (A, S, T), (N, D), (F, Y)
2. Generación de una tabla de K-meros con K entre 10 y 14. De cada secuencia solamente guardan 20 K-meros, elegidos por su frecuencia alta con una función hash.
3. Búsqueda de secuencias con idénticos K-meros
4. Pre-clustering en varios pasos, de más a menos eficientes: distancia de Hamming con alfabeto completo, alineamientos locales sin y con gaps.
5. Clustering voraz con las secuencias ordenadas por longitud

En sus pruebas Linclust es mucho más escalable, al ser lineal, que alternativas como CD-HIT o UCLUST, y obtiene buenos resultados para cortes de identidad entre 90 y 50%. Esto es ideal para exploraciones de metagenomas por ejemplo,
hasta pronto,
Bruno

pseudoalineamientos y conteos de tránscritos

Hola,
durante el análisis de unos experimentos de RNAseq mi colega Carlos Cantalapiedra y yo nos desesperábamos al calcular los valores de expresión de cada transcrito en una serie de condiciones, dado que para cada una de ellas era necesario alinear las lecturas/reads originales (en formato FASTQ) contra los transcritos ensamblados. Al menos éste era el protocolo habitual antes de aparecer el algoritmo genérico de pseudoalineamiento, que se describe en la siguiente figura:

Pseudoalineamiento de reads a partir de su composición en k-meros (k=3 en el ejemplo) y un  vector de sufijos de un transcriptoma. Figura tomada de http://bioinformatics.oxfordjournals.org/content/32/12/i192.full.  

Este tipo de algoritmos, implementados en software como kalllisto o RapMap (éste último con licencia GPL), permiten estimar de manera muy eficiente y precisa con qué transcritos son compatibles las lecturas de un archivo FASTQ, es decir, a qué secuencias de un transcriptoma mapean, sin necesidad de alinear base a base los reads. Los alineamientos se pueden calcular, si fuera preciso, después del mapeo. 

Hasta la fecha se han propuesto al menos dos implementaciones del algoritmo genérico de pseudoalineamiento, que se diferencian fundamentalmente en las estructuras de datos utilizadas, como se explica en estos blogs (1 , 2) [en el segundo se hace un repaso a la evolución de este tipo de algoritmos y a su nomenclatura].

En los próximos párrafos muestro dos aplicaciones usando ejecutables de RapMap y como datos un fichero (tr.fna) con 67K transcritos de Arabidopsis thaliana, ensamblados de novo con trinity, y otro (1.f1) con 89M de lecturas Illumina SE en formato FASTQ.

1. Mapeo de reads (20 cores, fichero de salida SAM)

          operación   tiempo   comando      

BWAmem       índice    1m28s   bwa-0.7.12/bwa index tr.fna    

BWAmem        mapeo    9m57s   bwa-0.7.12/bwa mem -t 20 tr.fna 1.fq > 1.sam 

RapMap       índice      58s   RapMap-0.2.2/bin/rapmap quasiindex -t tr.fna -i ./ 

RapMap        mapeo    5m56s   RapMap-0.2.2/bin/rapmap quasimap -t 20 -i ./ -r 1.fq \

                                 -o 1.sam 

Por comparación con BWA queda claro que incluso generando alineamientos SAM estos algoritmos son mucho más rápidos. 


2. Conteo de transcritos (20 cores, implementación Sailfish)

          operación   tiempo   comando      

Sailfish     índice    1m06s   SailfishBeta-0.10.0/bin/sailfish index -t tr.fna -o qindex \

                                  -p 20   

Sailfish      mapeo    1m09s   SailfishBeta-0.10.0/bin/sailfish quant -i qindex/ -r 1.fq \

                                  -p 20 -o 1.bur-0.sf --auxDir tmpsf --dumpEq --libType SF
  

Como se muestra en la tabla, la estimación de valores de expresión de transcritos por pseudoalineamiento y conteo de reads es muy eficiente. En este ejemplo se genera un archivo de salida (quant.sf) con este contenido:

Name      Length EffectiveLength TPM NumReads
TR1|c0_g1_i1 517 316.623      5.3984 214.782
TR2|c0_g1_i1 390 191.501     3.36608     81
TR3|c0_g1_i1 699 498.611      10.502 658
...

Entre las opciones del program está la posibilidad de calcular los conteos con bootstraping, como hace kallisto, aunque en la versión que yo he probado es experimental.

En la siguiente entrada veremos más aplicaciones,
un saludo,
Bruno

HOWTO install Grid Engine on multi-core Linux box to run GET_HOMOLOGUES

Hi,
I post this in English so that all users of GET_HOMOLOGUES can benefit. It is inspired on a previous tutorial posted in scidom and the expertise of David Ramírez from our provider SIE.

The aim of this entry is to demonstrate how to set up a Grid Engine queue manager on your local multi-core Linux box or server so that you can get the most of that computing power during your GET_HOMOLOGUES jobs. We will install Grid Engine on its default path, /opt/, which requires superuser permissions. As I'll do this un Ubuntu, I will do 'sudo su' to temporarily get superuser privileges; this should be replaces simply with 'su' in other Linux flavours. Otherwise you can install it elsewhere if you don't have admin rights.

1) Visit http://gridscheduler.sourceforge.net , create a new user and download the latest 64bit binary to /opt:

$ cd /opt
$ sudo su 
$ useradd sgeadmin
$ wget -c http://dl.dropbox.com/u/47200624/respin/ge2011.11.tar.gz $ tar xvfz ge2011.11.tar.gz
$ chown -R sgeadmin ge2011.11/
$ chgrp -R sgeadmin ge2011.11/

$ ln -s ge2011.11 sge

2) Set relevant environment variables in /etc/bash.bashrc  [system-wide, can also be named /etc/basrhc] or alternatively in ~/.bashrc for a given user:

export arch=x86_64
export SGE_ROOT=/opt/sge
export PATH=$PATH:"$SGE_ROOT/bin/linux-x64"
export LD_LIBRARY_PATH="$LD_LIBRARY_PATH:$SGE_ROOT/lib"

And make this changes live:

$ source /etc/bash.bashrc

3) Set your host name to anything but localhost by editing /etc/hosts so that the first line is something like this (localhost or 127.0.x.x IP addresses are not valid):

172.1.1.1   yourhost

4) Install Grid Engine server with all defaults except cluster name, which usually will be 'yourhost':
$ ./install_qmaster

5) Install Grid Engine client with all defaults:
$ ./install_execd

6) Optionally configure default all.q queue:
$  qconf -mq all.q

7) Add your host to list of admitted hosts:
$ qconf -as yourhost
$ exit

You should now be done. A test will confirm this. Please open a new terminal and move to your GET_HOMOLOGUES installation folder:

$ cd get_homologues-x86-20150306
$  ./get_homologues.pl -d sample_buch_fasta -m cluster

If the jobs finishes successfully then you are indeed done. I hope this helps,
Bruno

Creando un cluster Rocks virtual

Hola,
mientras actualizaba el material didáctico de Programación en clusters Rocks a la versión actual de Linux Rocks (5.4) me pareció buena idea ir probando el código en un cluster virtual, creado como máquinas virtuales corriendo bajo VirtualBox.

(fuente: http://www.rocksclusters.org/rocks-documentation/4.1/getting-started.html)

La ventaja que tiene hacerlo virtual es que cualquiera con unos cuantos Gb de disco libres y más de 2 Gb de RAM puede probar a programar en un cluster y así aprender a manejarse en ese entorno antes de tener acceso a uno real.

Éstos son los pasos necesarios, probados en mi Ubuntu 10.04 LTS y con la versión 4.0.8 de VirtualBox:

• Instala VirtualBox para tu sistema Linux (puedes descargarlo de
  http://www.virtualbox.org/wiki/Downloads)
  
  
• Descarga e instala el 'VirtualBox Extension Pack', que puedes obtener en el mismo lugar
  
  
• Descarga ISO del DVD Jumbo de Rocks, que ya incluye los rolls esenciales
  (puedes descargarlo de http://www.rocksclusters.org,
  comprobando tras la descarga la suma MD5). En este tutorial usaremos la versión 5.4, que se llama area51+base+bio+ganglia+hpc+kernel+os+sge+web-server+xen-16.12.2009-15.16.53.i386.disk1.iso,
  que guardamos en la carpeta /home/pepe/soft/
  
  
• Elige un directorio donde ubicar los discos duros virtuales, como por ejemplo /home/pepe/
  
  
• Crea el nodo principal tecleando (copia y pega) en el terminal estos comandos:
  

   VBoxManage createvm --name "vRocks54-Frontend" --ostype RedHat --register  
     
   VBoxManage modifyvm "vRocks54-Frontend" --memory 1000 --vram 32 --nic1 intnet --nic2 nat --audio none --nictype1 82540EM --nictype2 82540EM --boot1 dvd --boot2 disk --boot3 none --boot4 none  
     
   VBoxManage createhd --filename "/home/pepe/vRocks54-Frontend.vdi" --size 50000 --variant Standard   
     
   VBoxManage storagectl "vRocks54-Frontend" --name "SATA Controller" --add sata --controller IntelAhci  
     
   VBoxManage storageattach "vRocks54-Frontend" --storagectl "SATA Controller" --port 0 --device 0 --type hdd --medium "/home/pepe/vRocks54-Frontend.vdi"  
     
   VBoxManage storagectl "vRocks54-Frontend" --name "IDEcontrol" --add ide  
     
   VBoxManage storageattach "vRocks54-Frontend" --storagectl "IDEcontrol" --port 0 --device 0 --type dvddrive --medium /home/pepe/soft/area51+base+bio+ganglia+hpc+kernel+os+sge+web-server+xen-16.12.2009-15.16.53.i386.disk1.iso  
     
   VBoxManage startvm "vRocks54-Frontend"  
  

  
  Tras estos comandos da comienzo el proceso normal de instalación del nodo maestro o frontend , documentado en la web de Rocks. El proceso es casi automático, pero debes recordar dos cosas:
  
  
  
  • escribe frontend cuando aparezca el terminal de arranque boot:
    
  • cuando llegue el momento elige como mínimo estos rolls : kernel, base, web server, os
    
  
  
• Una vez instalado el frontend puedes añadir el paquete
  VirtualBox guest additions para mayor comodidad
  
  
• En el terminal del nodo maestro teclea como superusuario $ insert-ethers
  
  
• Para cada nodo compute que quieras crear, empezando por el 0-0, hasta el 0-N teclea en el terminal de tu Linux:
  
  

   VBoxManage createvm --name "vRocks54-Compute-0-0" --ostype RedHat --register  
     
   VBoxManage modifyvm "vRocks54-Compute-0-0" --memory 1000 --vram 32 --nic1 intnet --audio none --nictype1 82540EM --boot1 net --boot2 disk --boot3 none --boot4 none  
     
   VBoxManage createhd --filename "/home/pepe/vRocks54-Compute-0-0.vdi" --size 50000 --variant Standard   
     
   VBoxManage storagectl "vRocks54-Compute-0-0" --name "SATA Controller" --add sata --controller IntelAhci  
     
   VBoxManage storageattach "vRocks54-Compute-0-0" --storagectl "SATA Controller" --port 0 --device 0 --type hdd --medium "/home/pepe/vRocks54-Compute-0-0.vdi"  
     
   VBoxManage storagectl "vRocks54-Compute-0-0" --name "IDE Controller" --add ide   
     
   VBoxManage startvm "vRocks54-Compute-0-0"  
  

  
  
• En el terminal del nodo maestro teclea como superusuario la tecla F8 para terminar de añadir nodos
  

Y ésto es todo, con esto tienes ya un cluster en miniatura para probar tus programas, que puedes borrar con dos golpes de ratón desde el VirtualBox,
un saludo,
Bruno