Comprime y consulta ficheros FASTA con agc

Hola,  hoy voy a escribir sobre el software AGC (Assembled Genomes Compressor), escrito en C++, que sirve para comprimir ensamblajes genómicos en múltiples ficheros FASTA (un pangenoma) y luego extraer de manera eficiente secuencias de nucleótidos arbitrarias o contigs completos. El artículo que lo describe está en https://doi.org/10.1093/bioinformatics/btad097 y el código en https://github.com/refresh-bio/agc. 

AGC comprime los ensamblajes/ficheros FASTA en varias etapas usando por defecto k=31:

Principales etapas de la compresión: (a) selección de k-meros distribuidos de manera uniforme en genoma de referencia cada 60Kb, splitters;  (b)  compresión del genoma de referencia tras partirlo en grupos y segmentos; (c y d)  compresión del resto de genomas. Los segmentos flanqueados por los mismos splitters se comprimen juntos; los de referencia con zstd y los demás con LZSS para marcar las diferencias respecto al de referencia correspondiente para comprimir solamente esas partes con zstd. Figura tomada de https://doi.org/10.1093/bioinformatics/btad097

En nuestro caso, probamos AGC con 20 genomas de cebada del pangenoma V2, cada uno de unos 4GB, obteniendo una estructura de datos comprimida de 3.4GB con el siguiente comando:

$ agc create <ref.fa> <genoma2.fa> ... <genoma20.fa> assemblies.agc

# nota: se pueden añadir nuevos genomas más tarde con agc append 

Comprobamos el contenido del archivo comprimido:

$ agc listset assemblies.agc
Akashinriki
B1K-04-12
Barke
Chiba
Du_Li_Huang
GoldenPromise
HOR_10350
HOR_13821
HOR_13942
HOR_21599
HOR_3081
HOR_3365
HOR_7552
HOR_8148
HOR_9043
Hockett
Igri
MorexV3
OUN333
Planet
 

Finalmente, lo más interesante, podemos obtener secuencias de interés por medio de sus coordenadas:

$ agc getctg assemblies.agc chr1H@Barke:1-20 

>chr1H sampleName=Barke:1-20
ATGCTATTAGTCACTAATTT

Toda la documentación y más ejemplos están en https://github.com/refresh-bio/agc, que lo disfrutéis,

Bruno

Jugando a leer y extraer secuencias de un fichero FASTQ comprimido con GZIP (.fq.gz)

Un fichero con extensión '.fq.gz' es un fichero en formato FASTQ comprimido en formato GZIP. Los ficheros FASTQ con datos de los nuevos métodos de secuenciación (NGS) suelen ocupar decenas de Gigabytes de espacio de disco (una cantidad indecente para la mayoría de los mortales) y comprimidos con GZIP se reducen.

Para empezar intentaremos saber cuanto espacio ocupa realmente un archivo 'fq.gz' comprimido, para ello usaremos el mismo comando 'gzip':

> gzip --list reads.fq.gz
         compressed        uncompressed  ratio uncompressed_name
        18827926034          1431825024 -1215.0% reads.fq

Parece que GZIP en vez de comprimir expande, pero no es verdad, simplemente que la opción '--list' de 'gzip' no funciona correctamente para archivos mayores de 2GB. Así que hay que recurrir a un método más lento y esperar unos minutos:

> zcat reads.fq.gz | wc --bytes
61561367168

Si queremos echar un vistazo al contenido del archivo podemos usar el comando 'less' o 'zless':

> less reads.fq.gz
> zless reads.fq.gz

Y para saber el número de secuencias que contiene simplemente hay que contar el número de líneas y dividir por 4 (le costará unos minutos):

> zcat reads.fq.gz | echo $((`wc -l`/4))
256678360

Podemos buscar una determinada secuencia en el archivo y contar cuántas veces aparece, por ej. ATGATGATG:

> zcat reads.fq.gz | awk '/ATGATGATG/ {nlines = nlines + 1} END {print nlines}'

398065

A veces nos interesará tomar un trozo del fichero para hacer pruebas, por ejemplo las primeras 1000 secuencias (o 4000 líneas):

> zcat reads.fq.gz | head -4000 > test_reads.fq
> zcat reads.fq.gz | head -4000 | gzip > test_reads.fq.gz

O extraer un rango de líneas (1000001-1000004):

> zcat reads.fq.gz | sed -n '1000001,1000004p;1000005q' > lines.fq 

También nos puede interesar dividir el fichero en varios más pequeños de por ejemplo 1 miĺlón de secuencias (o 4 millones de líneas):

> zcat reads.fq.gz | split -d -l 4000000 - reads.split
> gzip reads.split*
> rename 's/\.split(\d+)/.$1.fq/' reads.split*

Y posteriormente reunificarlos:

> zcat reads.*.fq.gz | gzip > reads.fq.gz


Finalmente os recomiendo leer un post similar en inglés:
http://darrenjw.wordpress.com/2010/11/28/introduction-to-the-processing-of-short-read-next-generation-sequencing-data/

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

Leyendo archivos comprimidos .gz

Hola,
tras la entrada 'Compresión de secuencias de ADN', que ya utilizaba el módulo estándar Compress::Zlib, en este ejemplo se muestra como leer línea a línea un archivo comprimido con los algoritmos de la librería zlib. En particular este ejemplo lee un archivo FASTA de gran tamaño sin necesidad de descomprimirlo entero.  Lo he probado con éxito en Linux (con el intérprete Perl que viene instalado en Ubuntu 10.04) y también en Windows (con ActiveState Perl).

 use strict;  
 use Compress::Zlib;  
 # http://perldoc.perl.org/Compress/Zlib.html  
   
 my $filename = '/path/to/swissprot.gz';  
 # ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/  
   
 my $gz = gzopen($filename,'rb') ||   
    die "# Cannot open $filename: $gzerrno\n" ;  
 while($gz->gzreadline($_) > 0)   
 {  
   # haz algo con esta linea, guardada en $_, por ejemplo   
   # imprime los identificadores GenBank encontrados  
    if($_ =~ /^>/ && /(gi\|\d+)/){ print "$1\n" }  
 }  
   
 if($gzerrno != Z_STREAM_END)  
 {  
    die "# Error reading from $filename: $gzerrno\n"  
 }  
   
 $gz->gzclose();  

Hasta pronto, Bruno

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