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...

Procesar argumentos con múltiples opciones en un script en Python

En los últimos tiempos estoy vendiendo mi alma al diablo y he empezado a usar Python. Personalmente prefiero Perl, Python me parece un lenguaje más anárquico y desorganizado. Sin embargo, Python está de moda en Bioinformática y muchos de los módulos más actualizados se están escribiendo en este lenguaje.

Hoy voy a hablar del problema que me he encontrado para encontrar un módulo que procese las opciones de línea de comandos que normalmente introducimos en los scripts, por ej.: ./script.py -i archivo1 -o archivo2 Estas opciones suelen ser archivos de entrada, salida, o simplemente parámetros que cambian la ejecución del programa. También sirve para pedir ayuda sobre como usar el script: ./script.py -h 

En Perl existe el módulo Getopt::Long para procesar dichos parámetros de línea de comando de una forma sencilla y robusta.

En Python el módulo equivalente es getopt. A continuación se muestra un ejemplo de uso. Sin embargo este módulo no contempla la posibilidad de procesar múltiples opciones para un mismo argumento, por ejemplo, varios archivos de entrada.

 #!/usr/bin/python  
 # -*- coding: utf-8 -*-   
 """  
      Script que usa el módulo 'getopt' para parsear los argumentos y opciones de línea de comandos  
 """  
 import sys, re, os, copy  
 # Importar el módulo 'getopt'  
 import getopt  
 # Imprimir en pantalla el script y sus argumentos y opciones  
 print 'ARGV:', " ".join(sys.argv[0:]),"\n"  
 # Información de ayuda de uso del script  
 def help() :  
      """Ayuda sobre opciones del script en línea de comandos"""  
      print "usage:",sys.argv[0], "[options]\n"  
      print " -h this message\n"  
      print " -i input file/s\n"  
      print " -o output file\n"  
 # Definir los argumentos obligatorios (:) y opcionales, en sus formas abreviada y completa  
 try:  
      options, args = getopt.getopt(sys.argv[1:], "hi:o:", ["help", "input", "output"])  
 except getopt.GetoptError as err:  
      # Imprimir ayuda y salir si hay algún error o argumento incorrecto  
      print str(err)  
      help()  
      sys.exit(2)  
 output = None  
 verbose = False  
 # Definir el tipo de las variables que guardarán las opciones  
 INP_input = ''  
 INP_output = ''  
 # Asignar las opciones a cada argumento  
 for _opt, _arg in options:  
      if _opt in ("-i", "--input"):  
           INP_input = _arg  
      elif _opt in ("-o", "--output"):  
           INP_output = _arg  
      elif _opt in ("-h", "--help"):  
           help()  
           sys.exit()  
      else:  
           assert False, "unhandled option"  
 # Imprimir en pantalla todo lo anotado  
 print "Script:",sys.argv[0],"\n"  
 print "Argumentos y opciones:\n",  
 for _opt,_arg in options :  
      print " -"+_opt+": "+_arg  

Para solucionar este problema de Python de procesar argumentos de línea de comandos con múltiples opciones he creado mi propio código, puesto que me parece fundamental para cualquier script bioinformático. En el siguiente ejemplo se pueden definir listas para guardar múltiples opciones para un único argumento (ej. varios archivos). El código es provisional y le faltarían muchas opciones que sí tiene el módulo getopt, sin embargo, ofrece una solución al problema y espero ir mejorándolo con el tiempo.

 #!/usr/bin/python  
 # -*- coding: utf-8 -*-   
 """  
      Script que parsea los argumentos y opciones de línea de comandos  
      permitiendo usar listas para guardar argumentos con múltiples opciones  
 """  
 import sys, re, os, copy  
 from types import *  
 # Imprimir en pantalla el script y sus argumentos y opciones  
 print 'ARGV:', " ".join(sys.argv[0:]),"\n"  
 # Información de ayuda de uso del script  
 def help() :  
      """Ayuda sobre opciones del script en línea de comandos"""  
      print "usage:",sys.argv[0], "[options]\n"  
      print " -h this message\n"  
      print " -i input file/s\n"  
      print " -o output file\n"  
      exit()  
 # Definir el tipo de las variables que guardarán las opciones  
 INP_input = []  
 INP_output = ''  
 # Asignar a cada argumento una de las anteriores variables  
 options = {'h': 'help', 'help': 'help', 'i': 'INP_input', 'input': 'INP_input', 'o': 'INP_output', 'output': 'INP_output'}  
 argv_var = str()  
 # Leer los argumentos y opciones  
 for _argv in sys.argv[1:] :  
      # Leer los argumentos que comienzan con guión  
      if re.match("^-", _argv) :  
           for _opt,_var in options.iteritems() :  
                if re.compile("^-%s$" % (_opt)).match(_argv) :  
                     argv_var = copy.deepcopy(_var)  
                     if type(vars()[argv_var]) is FunctionType :  
                          vars()[argv_var]()  
                     break  
      # Anotar las opciones de cada argumento  
      else :  
           if argv_var is not None :  
                if type(vars()[argv_var]) is ListType :  
                     vars()[argv_var].append(_argv)  
                else:  
                     vars()[argv_var] += _argv  
 # Imprimir en pantalla todo lo anotado  
 print "Script:",sys.argv[0],"\n"  
 print "Argumentos y opciones:\n",  
 for _opt,_var in options.iteritems() :  
      if vars()[_var] is not None and type(vars()[_var]) is not FunctionType:  
           #print _opt,pprint(type(vars()[_var]))  
           if type(vars()[_var]) is ListType:  
                print " -"+_opt+": "+str(", ".join(vars()[_var]))  
           else :  
                print " -"+_opt+": "+str(vars()[_var])