9 de diciembre de 2010

Compresión de secuencias de ADN

Un problema que me ha surgido recientemente es el de manejar en memoria RAM grandes volúmenes de secuencias de ADN, en mi caso genes, pero que podrían ser también lecturas de secuenciación masiva. Cuando el tamaño total de las secuencias supera la RAM disponible el sistema operativo se verá forzado a hacer swapping/intercambio con el disco y el proceso se ralentiza considerablemente. Una de las posibles soluciones, la que nos ocupa hoy, es la compresión de las secuencias, de manera que ocupen menos espacio físico en RAM.
El siguiente programa muestra cómo podemos realizar esto en Perl, y nos sugiere que esta estrategia vale realmente la pena a medida que el tamaño de las secuencias se aproxima a las 200 bases:




 use strict;  
 use Compress::Zlib;  
   
 # http://perldoc.perl.org/Compress/Zlib.html  
   
 my (@S,$source,$compressed,$uncompressed);  
   
 print "size(source)\tsize(compressed)\tratio\n";  
 for(my $size=25;$size<=500;$size += 25)  
 {  
    # secuencia aleatoria de longitud $size  
    @S = map { ('A','C','G','T')[rand(4)] } 1 .. $size;  
    $source = join(/ /,@S); #print "$source\n";  
   
    $compressed = compress($source);  
    $uncompressed = uncompress($compressed) ;  
      
    printf("%d\t%d\t%g\n",  
       $size,length($compressed),length($compressed)/$size);  
 }  
   

Un saludo,
Bruno

1 comentario:

  1. Hoy en día, los ordenadores pueden tener mucha memoria RAM... varios Gb.

    Y aún en los casos en los que el fichero a tratar exceda el tamaño disponible en la RAM, tampoco suele ser necesario realizar técnicas de compresión u otros trucos para ir leyendo el fichero por partes: el sistema operativo lo hará por nosotros.

    Los sistemas operativos actuales disponen de una función llamada mmap(2). Esta función permite que los programas de usuario accedan al contenido de los ficheros abiertos en el disco, sin tener que preocuparnos del tamaño que tenga el fichero o de las operaciones de lectura/escritura. Para el programa, será como si, de repente, dispusiera de una gigantesca variable que almacenase todo el contenido del fichero.

    En Perl hay varios módulos que permiten usar esa técnica. Por ejemplo, File::Map, que funciona en Linux, Windows y VMS.

    Un ejemplo de uso: tenemos la secuencia de tres mil millones de bases de un ser humano en un fichero de un poco menos de 3Gb, llamado "human.gen". Para acceder a su contenido y buscar por la secuencia "ATTCG", hacemos

    use File::Map 'map_file';
    map_file my $humano, 'human.gen';
    $número++ while $humano =~ /ATTCG/g;

    (Bueno, sí que hay límites: los de la arquitectura del sistema, y los "bugs" y "pitfalls" del módulo ;)

    Más información:
    http://en.wikipedia.org/wiki/Mmap
    http://es.wikipedia.org/wiki/Extensi%C3%B3n_de_direcci%C3%B3n_f%C3%ADsica
    http://perlenespanol.com/foro/post23785.html#p23785
    http://search.cpan.org/perldoc?File::Map

    ResponderEliminar