lectura eficiente de ficheros comprimidos en el terminal

Al analizar datos genómicos a menudo tenemos que manejar ficheros de texto enormes, del orden de decenas de GB, que suelen estar comprimidos para no llenar el disco.

Imagen de https://the-chestnut.com

Por ejemplo:

-rw-r--r--  1 contrera contrera 2.2G Feb 23 12:43  mergedpairs.tsv

comprimido con gzip ocupa casi la décima parte:

-rw-r--r--  1 contrera contrera 255M Feb 23 12:43  mergedpairs.tsv.gz

A veces es necesario operar sobre este tipo de ficheros desde scripts o el terminal sin descomprimirlos en el disco. Por ejemplo, podemos contar el número de líneas con zcat:

time zcat mergedpairs.tsv.gz | wc -l
10205212

real    0m10.793s
user    0m10.689s
sys    0m1.726s

Otra opción un poco más rápida, qué descubrí en dos foros (1,2), es unpigz, que hay que instalar primero en Ubuntu:

sudo apt install pigz

time unpigz -dc mergedpairs.tsv.gz | wc -l
10205212

real    0m6.503s
user    0m10.742s
sys    0m3.391s
 

Otra opción que he encontrado es https://github.com/mxmlnkn/rapidgzip, que además tiene una interfaz python que permite indexar ficheros .gz, pero se puede usar en la línea de comando también:

python3 -m pip install rapidgzip

time rapidgzip -d -c -P 0 mergedpairs.tsv.gz | wc -l
10205212

real    0m2.833s
user    0m10.895s
sys    0m5.020s

Finalmente, he probado con https://github.com/facebook/zstd, que usa un formato de compresión menos estándar (.zst):

sudo apt install zstd
pzstd mergedpairs.tsv

time pzstd -d -c mergedpairs.tsv.zst | wc -l
mergedpairs.tsv.zst : 2273983394 bytes                                         
10205212

real    0m1.600s
user    0m5.973s
sys    0m3.348s
 

 

Hasta pronto,

Bruno

Algunos comandos útiles de linux para manejar ficheros FASTQ de NGS

Leyendo la entrada de Bruno sobre la librería kseq.h me entraron ganas de recopilar en el blog unos valiosos comandos de Linux que nos pueden ayudar a salvar mucho tiempo al trabajar con millones de secuencias de NGS.

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 (mucho más que con el clásico ZIP).

A continuación se mostrarán los comandos para manejar ficheros FASTQ comprimidos, pero con unas leves modificaciones podrán también procesar ficheros FASTQ sin comprimir (no recomendado), por ejemplo cambiando 'zcat' por 'cat'. También se pueden modificar para procesar ficheros FASTA comprimidos y sin comprimir.


Para empezar, vamos a comprimir en formato GZIP un archivo FASTQ:
   > gzip reads.fq
El fichero resultante se llamará por defecto 'reads.fq.gz'.

Nos puede interesar conocer 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
Contar el número de secuencias en un fichero FASTA:
   > grep -c "^>" reads.fa

Podemos contar cuántas veces aparece una determinada secuencia, por ej. ATGATGATG:
   > zgrep -c 'ATGATGATG' reads.fq.gz
   398065
   > zcat reads.fq.gz | awk '/ATGATGATG/ {nlines = nlines + 1} END {print nlines}'
   398065

O extraer ejemplos de dicha secuencia:
   > zcat reads.fq.gz | head -n 20000 | grep  --no-group-separator -B1 -A2 ATGATGATG
en ficheros FASTA:
   > zcat reads.fa.gz | head -n 10000 | grep  --no-group-separator -B1 ATGATGATG

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

Para extraer 1000 secuencias aleatorias de un archivo FASTQ:
    > cat reads.fq | awk '{ printf("%s",$0); n++; if(n%4==0) { printf("\n");} else { printf("X#&X");} }' | shuf | head -1000 | sed 's/X#&X/\n/g' > reads.1000.fq
O de un archivo FASTA:
    > cat reads.fa | awk '{if ((NR%2)==0)print prev"X#&X"$0;prev=$0;}' | shuf | head -1000 | sed 's/X#&X/\n/g' > reads.1000.fa

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

Algunos de estos comandos y otras ideas pueden consultarse en:
http://darrenjw.wordpress.com/2010/11/28/introduction-to-the-processing-of-short-read-next-generation-sequencing-data/

Podéis ayudarme a completar el post escribiendo vuestros comandos más útiles en comentarios...

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/