Mostrando entradas con la etiqueta estructura de datos. Mostrar todas las entradas
Mostrando entradas con la etiqueta estructura de datos. Mostrar todas las entradas

3 de septiembre de 2019

cómo hacer filogenias de miles de genomas

Hola,
la acumulación de genomas completos humanos, actualmente del orden de decenas de miles, plantea problemas a la hora de calcular filogenias con las estructuras de datos y los algoritmos tradicionales. Por esa razón hay grupos desarrollando nuevas estrategias que beneficiarán también a los que, como nosotros, trabajamos en plantas, cuando lleguemos a esos números.

Hoy comento muy brevemente dos métodos que acabo de ver publicados en Nature Genetics. El primero se llama tsinfer y usa un árbol comprimido para almacenar las variantes genómicas en mucho menos espacio que una matriz VCF:

Tamaño de las estructuras de datos probadas por los autores de tsinfer, tomado de https://www.nature.com/articles/s41588-019-0480-1.

El segundo método se llama relate y se basa en reconstruir los eventos de recombinación de cromosomas ancestrales que explican los haplotipos observados. Este método calcula longitudes de ramas:


Resumen del algoritmo relate, tomado de https://www.nature.com/articles/s41588-019-0484-x.

Un saludo,
Bruno

15 de noviembre de 2011

grafos de De Bruijn

Hola,
un artículo educativo publicado recientemente en Nature Biotechnology me ha vuelto a recordar los grafos de De Bruijn, que ya había mencionado de pasada en una entrada anterior sobre vectores de sufijos. Hoy los veremos un poco más de cerca, por su importancia para la reconstrucción de secuencias genómicas a partir de fragmentos más pequeños que han sido previamente secuenciados.

Los modernos ensambladores de secuencias emplean grafos de De Bruijn para guardar en memoria los prefijos (nodos) y sufijos (aristas) de las lecturas o reads obtenidas en experimentos de secuenciación. Este tipo de representación permite reducir en parte la redundancia natural de este tipo de datos a la vez que facilita su ensamblaje posterior. La figura  compara los métodos de ensamblaje tradicionales (a) con los actuales (d). Los tradicionales se usan para pegar entre si secuencias obtenidas por el método de Sanger, que suelen ser largas, de buena calidad y en números pequeños. Los métodos que acompañan a los secuenciadores actuales se adaptan al nuevo escenario de secuencias cortas y en números varios órdenes de magnitud más elevados.

fuente: http://www.nature.com/nbt/journal/v29/n11/full/nbt.2023.html

Como se muestra en el panel d de la figura, un grafo de De Bruijn facilita la reconstrucción del genoma de partida, que se obtiene buscando un ciclo euleriano que recorra todas las aristas (sufijos) una sola vez. Esta estrategia es computacionalmente mucho más barata que la alternativa de buscar ciclos hamiltonianos, que en la práctica es imposible al ser un problema de tipo NP completo. Sin embargo, como pasa a menudo, la realidad se resiste a ser capturada por una sola teoría, e incluso los grafos de De Bruijn son todavía incapaces de resolver el problema de las secuencias repetidas, que aparecen en un genoma o región genómica no una sinó muchas veces.

En el Russell's blog encontraréis código fuente en Python para aprender cómo se construyen estos grafos. Espero haberos despertado la curiosidad, un saludo,
Bruno

21 de marzo de 2011

Vectores de sufijos para buscar patrones (suffix arrays)

Buenas,
hace poco tuve la suerte de cursar un taller práctico sobre problemas en el análisis, mapeo y ensamblado de secuencias genómicas, impartido por Paolo Ribeca (autor del software GEM). A parte de volver un poco abrumado por las dificultades que encierran los datos producidos por estas tecnologías de secuenciación de gran capacidad, sobre todo para bichos que no tienen todavía un genoma de referencia de buena calidad, he aprendido qué estructuras de datos soportan el manejo y análisis de todas estas secuencias. Además de tablas hash y grafos de De Bruijn, y la magia de Burrows-Wheeler, que se merecerá su propia entrada, hemos aprendido las ventajas de usar vectores de sufijos (suffix arrays) para el problema de buscar patrones o subcadenas exactas en textos muy largos, como son las secuencias genómicas. Pero, qué es un vector de sufijos? Trataré de explicarlo con el siguiente código en Perl:

 # calcula un vector de sufijos de manera que las diferentes subcadenas   
 # de una secuencia queden preordenadas lexicograficamente  
 sub make_suffix_array  
 {  
    my ($seq,$verbose) = @_;  
    my @suffix = (0 .. length($seq)-1); # en base 0  
   
    # ordena los length($seq) posibles sufijos lexicogr´aficamente 
    @suffix = sort {substr($seq,$a) cmp substr($seq,$b)} (@suffix);  
      
    if($verbose)  
    {  
       print "# suffix array for $seq :\n";  
       foreach my $suf (@suffix)  
       {   
          printf("%3d %s\n",$suf,substr($seq,$suf));   
       }  
       print "\n";  
    }  
    return @suffix;  
 }  

Si llamamos a la subrutina
my @suffix = make_suffix_array('TTTTAGATCGATCGACTAGACTACGACTCGA',1);
obtendremos:

30 A
22 ACGACTCGA
19 ACTACGACTCGA
14 ACTAGACTACGACTCGA
25 ACTCGA
17 AGACTACGACTCGA
4  AGATCGATCGACTAGACTACGACTCGA
10 ATCGACTAGACTACGACTCGA
6  ATCGATCGACTAGACTACGACTCGA
28 CGA
12 CGACTAGACTACGACTCGA
23 CGACTCGA
8  CGATCGACTAGACTACGACTCGA
20 CTACGACTCGA
15 CTAGACTACGACTCGA
26 CTCGA
29 GA
18 GACTACGACTCGA
13 GACTAGACTACGACTCGA
24 GACTCGA
9  GATCGACTAGACTACGACTCGA
5  GATCGATCGACTAGACTACGACTCGA
21 TACGACTCGA
16 TAGACTACGACTCGA
3  TAGATCGATCGACTAGACTACGACTCGA
27 TCGA
11 TCGACTAGACTACGACTCGA
7  TCGATCGACTAGACTACGACTCGA
2  TTAGATCGATCGACTAGACTACGACTCGA
1  TTTAGATCGATCGACTAGACTACGACTCGA
0  TTTTAGATCGATCGACTAGACTACGACTCGA

que es una lista de los 31 posibles sufijos de la cadena original, en orden lexicográfico. De hecho, si te fijas en el código, el vector realmente contiene sólo las posiciones (en base 0) de los sufijos ordenados, no su secuencia. Obviamente la construcción de este vector es costosa al necesitar de una ordenación (en Perl es por defecto un mergesort con un caso peor O(NlogN) ), pero luego permite consultas  O(logN), es decir, mucho más rápidas que simplemente recorrer la secuencia de principio a fin, al usar implícitamente un árbol binario.

Puedes probarlo mediante el siguiente ejemplo:

 match_pattern('TAG','TTTTAGATCGATCGACTAGACTACGACTCGA'); 
 
 # aprovecha el orden del vector de sufijos para buscar cualquier subcadena   
 # o patro´n (pattern) por medio de un árbol binario   
 sub match_pattern   
 {  
    my ($pattern,$seq) = @_;  
    print "# looking for pattern $pattern in sequence $seq (base 0)\n";  
    my @suffix = make_suffix_array($seq,1);  
    my $low = 0;  
    my $high = $#suffix;  
    my $patl = length($pattern);  
    my ($med,$submed);  
    while ($low <= $high)   
    {  
       my $med = int (($low+$high)/2); # punto medio de la búsqueda        
       # comparacion lexicográfica en punto medio, mira 'perldoc perlop'  
       my $comp = $pattern cmp substr($seq,$suffix[$med],length($pattern));  
       if($comp < 0){ $high = $med-1 }  # retrocedemos en @suffix   
       elsif($comp > 0){ $low = $med+1 } # avanzamos  
       else   
       {  
          my $submed = $med - 1; # sufijo inmediatamente anterior al punto medio  
          while($submed > 0 && $pattern eq substr($seq,$suffix[$submed],$patl))  
          {   
             print "# match at position $suffix[$submed]\n";  
             $submed--;   
          }  
          while ($med < $#suffix-1 && $pattern eq substr ($seq,$suffix[$med],$patl))  
          {  
             print "# match at position $suffix[$med]\n";  
             $med++;   
          }  
          last;  
       }  
    }  
 }  

Hasta otra,
Bruno

25 de febrero de 2011

Rendimiento de estructuras de datos en Perl

Buenas,
recientemente, tratando de resolver las dudas de un principiante en Perl, he podido volver a constatar la dificultad que suponen para los novatos los vectores (@arrays) y las tablas asociativas (%hashes). Sin embargo, estas estructuras de datos son muy empleadas por usuarios experimentados, puesto que permiten guardar y modelar datos complejos de manera relativamente sencilla. No obstante, después de tanto tiempo usando estas estructuras, creo que nunca me he puesto a medir el coste de un vector y de un hash para guardar los mismos datos en un programita Perl. Este coste es especialmente relevante para manejar volúmenes de datos grandes, como averiguaremos con ayuda de los módulos Benchmark y Devel::Size.

Si además invocamos el módulo GraphViz::Data::Grapher, podremos dibujar ambos tipos de estructuras de datos, lo cual puede ayudar a los más principiantes:




 #!/usr/bin/perl -w  
   
 use strict;  
 use Benchmark;    
 use Devel::Size;   
 use Scalar::Util;  
 #use GraphViz::Data::Grapher;   
   
   
 my $MUESTRAS = 100;   
   
 # 1) mide tiempo de creacion de estructura de datos  
 print "# construyendo estructuras de datos:\n";  
 print "> hash :\n";  
 timethis( $MUESTRAS, "crea_estructura_datos('hash')" );  
 print "> array:\n";  
 timethis( $MUESTRAS, "crea_estructura_datos('array')" );  
   
 # 2) mide memoria que necesita cada estructura  
 print "\n# memoria RAM usada (KB):\n";  
 my $ref_hash = crea_estructura_datos('hash');  
 printf("\n> hash: %1.1f\n",  
    Devel::Size::total_size($ref_hash)/1024);  
   
 my $ref_array = crea_estructura_datos('array');  
 printf("> array: %1.1f\n",  
    Devel::Size::total_size($ref_array)/1024);  
      
 # 3) mide tiempo de consulta de estructura de datos  
 print "\n# consultando estructuras de datos:\n";  
 print "> hash:\n";  
 timethis( $MUESTRAS, sub{ consulta_estructura_datos($ref_hash) } );  
 print "> array:\n";  
 timethis( $MUESTRAS, sub{ consulta_estructura_datos($ref_array) } );     
   
 ############################################  
   
 sub crea_estructura_datos  
 {  
    my ($hash_o_array) = @_;   
   
     my $referencia;  
      
    if($hash_o_array eq 'hash')  
    {  
        foreach my $n (1..100_000)   
        {  
          # las llaves ocupan menos como cadenas de caracteres  
          # http://codenode.com/perl-memory-usage  
          $referencia->{"$n"} = $n * 10;         #7165.6KB en CPU 64bits  
         #$referencia->{sprintf("%01.3f",$n)} = $n * 10; #7556.2KB "  
          #$referencia->{$n/0.3} = $n * 10;        #7914.3KB "  
        }  
         
       #descomenta para generar grafo de estructura como el del blog  
       # OJO: mejor que sea una estructura no muy grande  
       #my $grafo = GraphViz::Data::Grapher->new($referencia);  
       #print $grafo->as_png("hash.png");  
     }  
    else  
    {  
       foreach my $n (1..100_000)   
        {  
         push(@{$referencia}, $n * 10);  # 3367.9KB           
          #$referencia->[$n-1] = $n * 10; # lo mismo  
        }  
         
       #my $grafo = GraphViz::Data::Grapher->new($referencia);  
       #print $grafo->as_png("array.png");  
    }  
   
    return $referencia;  
 }   
      
 sub consulta_estructura_datos  
 {  
    my ($referencia) = @_;   
      
    my $index;  
      
    if(Scalar::Util::reftype($referencia) eq "HASH")  
    {  
        foreach my $n (1..100_000)   
        {  
          $index = int(rand(100_000));  
          $referencia->{$index} += 1;   
        }  
     }  
    else  
    {  
       foreach my $n (1..100_000)   
        {  
         $index = int(rand(100_000));  
          $referencia->[$index] += 1;   
        }  
    }  
   
    return $referencia;  
 }   
   

Mediante este código, ejecutado en mi máquina de 64 bits, podemos ver que un hash sencillo ocupa un poco más del doble que un array para guardar una lista de enteros, y que su tamaño varía según el tipo de llave que usemos. Además, tardamos el triple de tiempo en llenar el hash que el array equivalente. Finalmente, este pequeño experimento muestra que el acceso y modificación de datos dentro un hash es al menos dos veces más lento que en un array:

# construyendo estructuras de datos:
> hash :
timethis 100: 10 wallclock secs (10.45 usr +  0.01 sys = 10.46 CPU) @  9.56/s (n=100)
> array:
timethis 100:  3 wallclock secs ( 2.17 usr +  0.00 sys =  2.17 CPU) @ 46.08/s (n=100)

# memoria RAM usada (KB):

> hash: 7165.6
> array: 3367.9

# consultando estructuras de datos:
> hash:
timethis 100:  6 wallclock secs ( 6.08 usr +  0.01 sys =  6.09 CPU) @ 16.42/s (n=100)
> array:
timethis 100:  3 wallclock secs ( 2.61 usr +  0.00 sys =  2.61 CPU) @ 38.31/s (n=100)

Un saludo,
Bruno