21 de diciembre de 2012

Bioinformatic Merry Christmas!!!

My boss is very nervous when I prepare this blog entry at the end of December, and this 'crisis' year I wanted to prepare the best bioinformatic inspirated Christmas card ever..
Merry Bioinformatic Christmas!

If you are not as geek as us, here is written in flat words... MERRY CHRISTMAS AND HAPPY NEW YEAR 2013!!!

If you want to replicate the Christmas greeting, go to this link, or to the 'retrieve' section in UniProt and introduce the following protein identifiers: P0A183, O02437, Q40309, Q6L5Z2, E0ZPQ6.
 

30 de noviembre de 2012

comparando secuencias de igual longitud

Hola,
hoy quería compartir recetas en Perl para la comparación eficiente de secuencias de igual longitud. Las dos funciones que vamos a ver devuelven como resultado listas con las posiciones (contando desde 0) en las que difieren dos secuencias.

Para qué sirve esto? En general, para calcular distancias de edición (edit distances); en bioinformática es una manera, por ejemplo, de comparar oligos de igual longitud o secuencias extraídas de un alineamiento múltilple, incluyendo gaps. El siguiente código fuente incluye una subrutina en Perl (con el operador binario ^ y la función pos) y otra en C embebido para hacer la misma operación por dos caminos:

 use strict;  
 use warnings;  
 my $seq1 = 'ACTGGA';  
 my $seq2 = 'AGTG-A';  
   
 my @Pdiffs = find_diffs($seq1,$seq2);  
 my @Cdiffs = find_diffs_C($seq1,$seq2);  
   
 printf("# version Perl\ndiferencias = %d\n%s\n",scalar(@Pdiffs),join(',',@Pdiffs));  
 printf("# version C\ndiferencias = %d\n%s\n",scalar(@Cdiffs),join(',',@Cdiffs));  
   
 sub find_diffs  
 {  
    my ($seqA,$seqB) = @_;  
    my $XOR = $seqA ^ $seqB; 
    my @diff;   
    while($XOR =~ /([^\0])/g)  
    {   
       push(@diff,pos($XOR)-1);  
    }  
    return @diff;  
 }  
   
 use Inline C => <<'END_C';  
 void find_diffs_C(char* x, char* y)   
 {                        
   int i;   
    Inline_Stack_Vars;                               
   Inline_Stack_Reset;                                                                
   for(i=0; x[i] && y[i]; ++i) {                          
    if(x[i] != y[i]) {                              
     Inline_Stack_Push(sv_2mortal(newSViv(i)));                 
    }                                       
   }   
    Inline_Stack_Done;                                                                  
 }     
 END_C  
Podéis ver en stackoverflow otras versiones y su comparación en base a su tiempo de cálculo,
un saludo,
Bruno




21 de noviembre de 2012

Comprimiendo BLAST

Hola,
ya hemos hablado antes en este blog de las crecientes aplicaciones de los algoritmos de compresión en la bioinformática. Hoy precisamente quería enlazar a dos artículos recientes que los explotan para acelerar algo que a priori parece imposible, la herramienta más universal de la biología computacional, BLAST.

En un paper en Nature Biotech, Loh, Maym y Berger nos presentan el prototipo Compression-accelerated BLAST (fuente C++ aquí), que en pruebas empíricas acelera varios órdenes de magnitud las búsquedas de BLAST simplemente haciendo las operaciones en un espacio comprimido. Como ventaja añadida se generan archivos de resultados de tamaños significativamente menores. Obviamente pagamos una pequeña pérdida de sensibilidad, pero puede ser perfectamente asumible para tareas de mapeo de lecturas (reads) de secuenciación de última generación:

Original de http://www.nature.com/nbt/journal/v30/n7/full/nbt.2241.html.

En otro paper recientemente publicado en Bioinformatics, Koskinen y Holm, aplican el marco de los vectores de sufijos (ya discutidos aquí) a la búsqueda de proteínas homólogas con una identidad superior al 50%, acelerando de nuevo varias órdenes de magnitud por encima de BLAST.

Familias de algoritmos para la búsqueda inexacta de secuencias, según Koskinen y Hol (http://bioinformatics.oxfordjournals.org/content/28/18/i438.full).
El programa que implementa estas ideas se llama SANS, escrito en FORTRAN 90, está disponible en este enlace,
un saludo,
Bruno

31 de octubre de 2012

Tutorial de modelado comparativo

Hola,
en colaboración con mis colegas Daniela Medeot, Romina Rivero y Edgardo Jofré de la Universidad Nacional de Río Cuarto (Argentina) hemos organizado recientemente un curso de posgrado titulado "Herramientas de bioinformática aplicadas al análisis de ácidos nucleicos y proteínas", al que asistieron cerca de 20 alumnos.

Además de hacer una introducción práctica al uso de Perl en diferentes tareas cotidianas de la bioinformática, mi participación se centró en un taller práctico para el que preparé el tutorial "Modelado comparativo de proteínas"  que espero pueda servir de ayuda a personas interesadas en aprender a modelar proteínas. El material puede utilzarse desde el navegador web o descargarse como un documento PDF,
un saludo, Bruno

18 de septiembre de 2012

TFcompare - a tool for structural alignment of DNA binding protein complexes

I want to introduce you the new bioinformatic contribution of our lab to the science world: TFcompare (http://floresta.eead.csic.es/tfcompare/)

TFcompare is a tool for structural alignment of DNA motifs and protein domains from DNA binding protein complexes in Protein Data Bank

The TFcompare algorithm calculates structural alignments between three dimensional structures of two DNA-protein complexes. The most interesting feature of TFcompare when compared with other methods is that it extracts individual protein domains and their recognized DNA sequences, aligning them separately and returning not only the structure superposition but the DNA sequence superposition too. In this way we can compare single domain affinity for different DNA sequences in DNA-protein complexes, especially transcription factors and their recognized cis elements.

The working schema of TFcompare is the following:




TFcompare takes as input two PDB identifiers. Structures from PDB are retrieved automatically and Pfam domains contacting DNA are calculated and trimmed from the original structure. Then all the domains from the first structure are aligned to all the domains from the second in several steps:
  1. The program MAMMOTH performs the structural alignment.
  2. The produced transformation matrices are applied to the coordinates of the DNA binding sites in order to derive the equivalent cis element superpositions.
  3. Root-mean-squared deviations of superposed coordinates are calculated with beta-carbon atoms (proteins) and with N9 (purines) and N1 (pyrimidines) atoms (DNA).
  4. Structural alignments are scored in terms of i) the number of identical superposed nucleotides (DNA Score 1-0) and ii) the sum of N9 and N1 atom pairs within 3.5 Å (DNA Score).
We can take as example the alignment of 1D5Y and 1BL0 structures, both are bacterial proteins with Helix-turn-helix (HTH) protein domains binding DNA.
We obtain the following results:

Results are ordered by structural similarity (RMSD), from both protein domain and DNA. In green colour are showed the alignment of similar structures (protein RMSD <=5.0 Å and DNA RMSD <= 3.5 Å) and in red colour dissimilar ones. 1D5Y contains two protein chains with HTH domains contacting DNA (trimmed domain structures are 1d5y_A1 and 1d5y_C1). 1BL0 have two HTH domains in its unique chain A (1bl0_A1 and 1bl0_A2). Structural alignment results show how 1bl0_A1 superposes very well with 1d5y_A1 and 1d5y_C1 (green colour). When DNA sequences recognized by these domains are aligned, they show a DNA motif conserved with three common nucleotides ‘CAC’. However, 1bl0_A2 superpositions (red colour) are not as good as previous ones, and DNA motif ‘CAC’ is not preserved when checking the resulting DNA alignment.

Each row contains an alignment of a pair of DNA binding domains, showing a picture of their structures before and after superposition. DNA alignment is also shown.

PDB files with aligned structures can be downloaded by left-clicking on the domain names and DNA sequences. Opening them with PDB viewer software (Pymol for ex.) is possible to visualize the resulting superposition after structural alignment.

Results column headers and their meaning:
Pair: Pair number
Domain_Query: PDB name, chain and domain number of the Query
Domain_Sbjct: PDB name, chain and domain number of the Sbjct
DNA_Query: DNA site recognized by the Query domain
DNA_Sbjct: DNA site recognized by the Query domain
Similar: 1 if both protein domains and DNA sites are below RMSD thresholds, 5.0 A and 3.5 A
respectively
DNA_Alignment: DNA sites structurally aligned
DNA_Aligned: Number of aligned nucleotides
DNA_Score_1-0: Number of identical nucleotides
DNA_Score: Structural alignment score
DNA_RMSD: RMSD of the structurally aligned DNA sites
PROT_RMSD: RMSD of the structurally aligned protein domains
3D_Alignment: 3D Visualization of aligned structures

9 de septiembre de 2012

Un SSD con Kubuntu 12 para gobernarlos a todos

ACTUALIZACIÓN: en la actualidad tengo instalado Kubuntu 12 de 64 bits en un disco de 256GB Samsung.

Hace ya unos meses cuando se liberó la versión beta de Ubuntu Pagolin 12.04 LTS me planteé instalarla en mi ordenador de trabajo. La principal novedad de Ubuntu 12.04 es que es una versión estable que será soportada durante 5 años. Uno de los inconvenientes que a mi parecer tiene Ubuntu es su falta de estabilidad en la migración de una versión a otra (yo siempre hago copia de seguridad antes). Siempre tengo problemas al actualizar una determinada versión y me tengo que quedar con la antigua si no quiero reconfigurar todo el ordenador de cero, cosa que con mi Windows XP nunca me ha pasado. Así que con Ubuntu 12.04 puedo estar tranquilo al menos 5 años :)



A su vez, me planteé la posibilidad del teletrabajo, aunque más que teletrabajo suele ser llevarse trabajo extra a casa. Y debido a las restricciones de acceso que tengo a mi ordenador de trabajo desde el exterior de mi centro y a la necesidad de portabilidad cuando viajo al extranjero, me surgió una vieja idea que tenía en mi cabeza... llevar mi ordenador en un pendrive. Haciendo cuentas de espacio, necesitaba un mínimo de 50GB para montar un sistema linux con todos los programas y bases de datos que suelo utilizar.
Los actuales pendrives a precios razonables no llegan a dicha capacidad, además son muy lentos comparados con los discos duros tradicionales y los modernos SSDs. Encontré por casa un disco duro viejo de portátil que tenía en una carcasa externa para hacer copias de seguridad y decidí instalar la versión de 64 bits de Kubuntu 12.04 para probar como funcionaba en un disco externo conectado por USB.

Disco duro IDE de portátil (2,5' 100GB) con carcasa USB.

Primero copié la imagen ISO de Kubuntu en un pendrive (así se evita desperdiciar 1 CD) con el programa LinuxLive USB. Entonces instalé Kubuntu tal y como se explica en los numerosos tutoriales que hay en internet, teniendo cuidado de elegir el disco correcto (en mi caso '/dev/sdc') y particionando el disco en una partición root ('/') de 20GB y una para el home ('/home') de 80GB. El espacio Swap no lo asigné a una partición sino a un archivo tras finalizar la instalación, ver más abajo.

Particionado de mi disco duro externo.
Al probar la instalación de Kubuntu en el disco externo USB quedé sorprendido que funcionaba con la misma fluided que mi anterior instalación en un disco fijo SATA. Para tareas ordinarias parece que la velocidad USB es suficientemente alta como para no notarlo en el rendimiento. Ahora llegaba la prueba de fuego... conectar el disco en otro ordenador por USB y ver si también funcionaba. La instalación de Kubuntu no arrancaba en mi portátil (un poco antiguo por cierto), la explicación es que mi portátil no soporta sistemas operativos de 64 bits.
NOTA: aunque sigo siendo partidario de los 32 bits, recomiendo instalar 64 bits a todo el que trabaje en NGS, si no programas como Bowtie o Trinity nunca funcionarán.

Volví a instalar Kubuntu, pero esta vez la versión de 32 bits en el disco y ahora sí que funcionaba, tanto en mi ordenador de sobremesa de 64 bits, en mi portátil de 32 e incluso en mi netbook Acer Aspire One, reconociendo tanto tarjetas de video, sonido y wifi en todos ellos. Con la sorpresa que en el netbook aparecía un escritorio diferente, el KDE Netbook Plasma. Dicho escritorio es algo raro al principio, pero según te acostumbras resulta más cómodo para trabajar con pantallas pequeñas. Haciendo una prueba sencilla de rendimiento, con la versión de 32 bits perdía un 10% de velocidad al ejecutar scripts de perl, una pérdida que no me importaba debido a su gran compatibilidad. Seguramente en programas optimizados de compresión de vídeo o retoque fotográfico el rendimiento será peor, pero en mi caso los programas científicos que requieren alto rendimiento están en servidores externos, y en mi ordenador simplemente pruebo scripts o tareas que más tarde puedo mandar a los servidores.

Aspecto del KDE Netbook Plasma.


Ahora tenía lo que quería... UN DISCO CON LINUX QUE PUEDO CONECTAR A CUALQUIER ORDENADOR Y TRABAJAR CON ÉL.

Tras probar "mi nuevo ordenador" durante varios meses en el trabajo, en casa, viajando... he tenido tan buenas experiencias que he decidido "ampliarlo". La ampliación ha consistido en migrar el sistema a un disco externo SSD. En teoría la mejora en rendimiento de un disco duro clásico a uno SSD es sustancial. El disco elegido ha sido un Samsung 830 de 128GB, tuve dudas entre éste y un Crucial M4, pero me decanté por el Samsung debido a la mínima diferencia de precio y la fiabilidad de la marca

Disco SSD y carcasa externa.
El primer paso de la migración fue usar un pendrive con Clonezilla (también creado con LinuxLive USB). Tras arrancar el ordenador con el pendrive y los 2 discos externos USB conectados, seleccionando la opción 'device-device' pude clonar mi antiguo disco externo en el nuevo SSD. Sin embargo, debido a las diferencias entre los discos, el nuevo disco no arrancaba. Para arreglarlo, retomé mi pendrive de Kubuntu y lo reinstalé en mi partición root ('/') copiada en el disco SSD. El instalador de Kubuntu es bastante inteligente y restaura el registro de arranque (MBR) y alguna cosa más sin alterar la partición root y los programas que ya teníamos instalados.

Pantalla ejemplo de Clonezilla.
Ya sólo quedaba comprobar que funcionaba en mis varios ordenadores...
(sólo se resiste a reconocerlo una placa base Gygabyte GA-965P-DS4 rev. 1.0)


Respecto a la mejora en el rendimiento, hasta el momento sólo puedo decir que el tiempo de arranque es menor y que LibreOffice se inicia en la mitad de tiempo usando el disco SSD.

Disco en la carcasa comparado en tamaño con un bolígrafo.
Y aún se puede optimizar más el discoy alargar su vida útil, básicamente minimizando las operaciones de lectura/escritura en disco, para ello podemos hacer las siguientes modificaciones en nuestro linux:
  • Añadir a '/etc/fstab' las opciones noatime,nodiratime y discard en todos los puntos de montaje del disco.
  • Crear y montar el archivo swap indicando que sólo se use cuando sea estrictamente necesario:
  1. Crear el archivo:
    • sudo dd if=/dev/zero of=/media/swapfile bs=1M count=4000 
  2. Activar el archivo swap:
    • sudo mkswap /media/swapfile
    • sudo swapon /media/swapfile
  3. Añadir en '/etc/fstab' la siguiente línea para que se active al arrancar:
    • /media/swapfile swap    swap    defaults        0       0  
  4. Añadir la siguiente línea en '/etc/sysctl.conf' para dar preferencia de uso a la memoria RAM:
    •  vm.swappiness = 0
  5.  Comprobar el funcionamiento:
    •  cat /proc/swaps
  • Montar en RAM el cache del sistema, para ello añadir a '/etc/fstab' la línea:
    • tmpfs /tmp tmpfs defaults,noatime,mode=1777 0 0
Más detalles sobre la optimización en:
http://askubuntu.com/questions/1400/how-do-i-optimize-the-os-for-ssds
http://www.howtogeek.com/62761/how-to-tweak-your-ssd-in-ubuntu-for-better-performance

Al instalar Kubuntu 12.04 de 64 bits el ordenador se volvió anormalmente lento y no funcionaba mi conexión a monitor externo, lo cual se solucionó actualizando los drivers de NVIDIA:
http://askubuntu.com/questions/133181/ubuntu-12-04-x64-very-slow-response-and-sluggishness

30 de agosto de 2012

Cuando no queda nada para secuenciar...

Normalmente pensamos que una DNA polimerasa puede amplificar cualquier resto de materia orgánica, sin embargo sucesos tan desafortunados como el error en la investigación policial de 2 niños "supuestamente" asesinados que ha conmocionado a la sociedad española en los últimos días hacen replantearnos los límites de la ciencia.

Cualquier juez, abogado o persona de a pie pensaría que la prueba de ADN es infalible para identificar restos humanos. Sin embargo la ciencia tiene límites, y en casos extremos como el mencionado, hay que recurrir a técnicas clásicas de reconocimiento forense de restos humanos, por ejemplo la observación al microscopio.

Leyendo la literatura científica podemos encontrar un interesantes artículo científico del laboratorio de la Dra. Nicole von Wurmb-Schwark
titulado "Reliable genetic identification of burnt human remains".  En el mismo se estudia la posibilidad de éxito en la identificación genética de cadáveres mediante restos de ADN en huesos quemados o calcinados.

Se estudian 71 fragmentos de huesos de 13 cadáveres. Los fragmentos se dividen en 5 categorías, según su grado de combustión (de menor a mayor):
  • A: Bien conservado
  • B: Semi-qiemado (200–300 °C)
  • C: Quemado negro (300–350 °C)
  • D: Quemado azul-gris (550–600 °C)
  • E: Quemado azul-gris-blanco (>650 °C)

La conclusión obtenida es que el análisis genético de los restos de ADN en los huesos quemados es muy complicada por tres motivos:
  1. La excasez y degradación del ADN en las muestras.
  2. La posible contaminación por ADN foráneo.
  3. La existencia de inhibidores de la polimerasa: colágeno, gasolinas, plásticos, componentes textiles...
Finalmente se concluyó que la identificación de forma reproducible de los cadáveres era posible únicamente en huesos bien conservados (A) o semi-quemados (B). La mayoría de los casos de huesos quemados negros (C) se podían identificar, aunque con perfiles incompletos, y también en algunos casos de los huesos quemados azul-gris (D). Sin embargo la identificación con huesos
quemados azul-gris-blanco (E) era infructuosa en la mayoría de los casos.

El estudio sugiere la amplificación de un fragmento de 220 pares de bases de la región HVI mitocondrial como último recurso para la identificación de huesos
quemados azul-gris-blanco (E), con éxito en el 30% de estos casos extremos.


21 de agosto de 2012

beca en el laboratorio de Biología Computacional

El laboratorio de Biología Computacional de la Estación Experimental de Aula Dei/CSIC (en Zaragoza, España) busca un candidato/a para desarrollar un proyecto de investigación en el marco de un programa de doctorado de la Universidad de Zaragoza.
El grupo tiene amplia experiencia en diferentes áreas de la Bioinformática y actualmente participa en diversos proyectos de genómica y secuenciación de plantas, incluyendo Arabidopsis thaliana, arroz y cebada, con un interés especial en el estudio funcional y evolutivo de las redes de regulación genética y en el descubrimiento de las raíces genéticas de la biodiversidad vegetal y agrícola en relación con la adaptación al medio ambiente. Para una lista completa de las publicaciones del grupo consultar http://www.eead.csic.es/compbio y http://www.eead.csic.es/index.php?id=99 .
Los candidatos interesados deben contactar con el Dr. Bruno Contreras Moreira (bcontreras@eead.csic.es)  antes del 15 de Septiembre de 2012 y elegir un tema de investigación de entre estos:
1) estudio de la especificidad y el desorden de factores de transcripción por medio de técnicas de dinámica molecular
2) genómica comparativa de cereales y la planta modelo Brachypodium distachyon (en colaboración con el grupo Bioflora de Pilar Catalán)

3)  estudio de las redes transcripcionales que controlan la pluripotencia y la diferenciación de células madre, usando una combinación de métodos bioquímicos y bioinformáticos (en colaboración con el grupo de Jon Schoorlemmer del  Instituto Aragonés de Ciencias de la Salud)

Las personas interesadas deberán cumplir los requisitos para solicitar una beca predoctoral del Gobierno de Aragón (2 años beca + 2 años contrato, ver convocatoria y requisitos en http://tinyurl.com/cocap6g). Se requiere un expediente académico superior a 6 y se valorará positivamente la experiencia previa en Bioinformática y/o haber terminado un Máster oficial. La convocatoria especifica una fecha límite de terminación de los estudios de licenciatura/grado/ingeniería posterior a Junio de 2011.

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.