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

regenerando DNA degenerado

Hola,
durante la reciente visita de mi colega Pablo Vinuesa al laboratorio pasamos ratos escribiendo código en Perl con el fin de diseñar de manera automática, para grandes conjuntos de secuencias de DNA previamente alineadas, parejas de primers (cebadores de PCR) adecuadas para estudios de microbiología ambiental, siguiendo principios bien conocidos ya. En cuanto probamos el código con unos cuantos ejemplos observamos que a menudo los primers para algunas regiones de los alineamientos múltiples eran degenerados. Como se muestra en la figura, eso significa que algunas posiciones de los cebadores no podían ser nucleótidos únicos si no que tenían que ser combinaciones de 2 o más para poder hibridarse con las secuencias utlizadas para su diseño.

Pareja de primers degenerados que definen un amplicón a partir de secuencias de DNA alineadas. De acuerdo con la nomenclatura IUPAC, el de la izquierda (fwd) puede escribirse como GAYTST y el de la derecha (rev) como GHGKAG. Figura tomada de http://acgt.cs.tau.ac.il/hyden.

A la hora de evaluar parejas de primers nos encontramos con que los degenerados en realidad tendrían que ser examinados por medio de cada una de las secuencias implícitas en el código degenerado. Y por tanto, tuvimos que buscar una manera de regenerar las secuencias de nucleótidos correspondientes a un primer degenerado, que es de lo que va esta entrada. El siguiente código en Perl hace precisamente esto:

 #!/usr/bin/perl  
 use strict;  
 my $degenerate_sequence = $ARGV[0] || die "# usage: $0 <degenerate DNA sequence>\n";  
 my $regenerated_seqs = regenerate($degenerate_sequence);  
 foreach my $seq (@$regenerated_seqs)  
 {  
    print "$seq\n";  
 }  
 # returns a ref to list of all DNA sequences implied in a degenerate oligo  
 # adapted from Amplicon (Simon Neil Jarman, http://sourceforge.net/projects/amplicon)  
 sub regenerate  
 {  
    my ($primerseq) = @_;  
    my %IUPACdegen = (   
    'A'=>['A'],'C'=>['C'], 'G'=>['G'], 'T'=>['T'],   
    'R'=>['A','G'], 'Y'=>['C','T'], 'S'=>['C','G'], 'W'=>['A','T'], 'K'=>['G','T'], 'M'=>['A','C'],  
    'B'=>['C','G','T'], 'D'=>['A','G','T'], 'V'=>['A','C','G'], 'H'=>['A','C','T'],   
    'N'=>['A','C','G','T']   
    );  
    my @oligo = split(//,uc($primerseq));   
    my @grow = ('');  
    my @newgrow = ('');  
    my ($res,$degen,$x,$n,$seq);  
    foreach $res (0 .. $#oligo){  
       $degen = $IUPACdegen{$oligo[$res]};   
       if($#{$degen}>0){  
          $x = 0;  
          @newgrow = @grow;  
          while($x<$#{$degen}){  
             push(@newgrow,@grow);  
             $x++;     
          }     
          @grow = @newgrow;  
       }  
       $n=$x=0;   
       foreach $seq (0 .. $#grow){  
          $grow[$seq] .= $degen->[$x];   
          $n++;  
          if($n == (scalar(@grow)/scalar(@$degen))){  
             $n=0;  
             $x++;  
          }     
       }  
    }  
    return \@grow;     
 }            

Podemos probarlo con las secuencias de la figura:

$ perl degenprimer.pl GAYTST
GACTCT
GATTCT
GACTGT
GATTGT

 y

$ perl degenprimer.pl GHGKAG
GAGGAG
GCGGAG
GTGGAG
GAGTAG
GCGTAG
GTGTAG

Hasta otra,
Bruno

Cuando Blastn no es Blastn

BLAST cambió hace ya algún tiempo a mejor con su versión BLAST+ pero por el camino se olvidó de algún detalle que puede confundir a más de uno.

El antiguo BLAST se ejecutaba con el comando 'blastall':

Blastall
--------

Blastall may be used to perform all five flavors of blast comparison. One
may obtain the blastall options by executing 'blastall -' (note the dash). A
typical use of blastall would be to perform a blastn search (nucl. vs. nucl.) 
of a file called QUERY would be:

blastall -p blastn -d nr -i QUERY -o out.QUERY

The output is placed into the output file out.QUERY and the search is performed
against the 'nr' database.  If a protein vs. protein search is desired,
then 'blastn' should be replaced with 'blastp' etc.

De esta forma un alineamiento de proteínas comenzaría como 'blastall -p blastp' y uno de ácidos nucleicos como 'blastall -p blastn' y si queremos usar MEGABLAST tenemos el comando diferente 'megablast'.

Sin embargo en la 'nueva' versión BLAST+, se separaron el alineamiento de proteínas y el de ácidos nucleicos en dos comandos: 'blastp' y 'blastn' (ver manual). Hasta aquí todo parece lógico y normal, lo que no todo el mundo sabe es que LA OPCIÓN POR DEFECTO DE BLASTN ES MEGABLAST, si ejecutamos 'blastn -help' encontraremos lo siguiente:

 *** General search options
 -task                 'megablast' 'rmblastn' >
   Task to execute
   Default = `megablast'

Y es que no todo el mundo está interesado en la velocidad de búsqueda de alineamientos, muchos de los que todavía usamos Blastn es porque apreciamos su gran sensibilidad para detectar alineamientos. En la actualidad usamos Blast para alinear miles de secuencias en tiempos muy razonables de minutos e incluso segundos. Para ganar velocidad en alineamientos de millones de secuencias existen otras mejores alternativas como Bowtie2.

La CONCLUSIÓN de todo esto, si usamos Blastn y nos interesa la sensibilidad deberemos ejecutarlo como:

 blastn -task blastn

Si lo hacemos sin añadir esta opción estaremos ejecutando MEGABLAST y correremos el peligro de perder una gran sensibilidad y no encontrar los alineamientos que deseamos. Por ejemplo, busquemos homología entre la 2'beta microglobulina humana (NM_004048.2) y la de ratón (NM_009735.3) usando la herramienta online de Blastn con las opciones por defecto ('Highly similar sequences (megablast)'):

Sin embargo si cambiamos la opción de búsqueda a 'Somewhat similar sequences (blastn)':

La diferencia es considerable, ¿no creéis? pasamos de no encontrar similaritud a un alineamiento con E-valor de 1.5E-56!!!!

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

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.