Curso de Python para biólogos - Lección 1. Ejecutando Python en Windows

Traicionando al blog, pero espero que beneficiando a sus lectores, vamos a publicar un curso de Python por entregas (cada 1 ó 2 semanas).

El curso está diseñado para introducir el mundo de la programación a biólogos, bioquímicos, médicos... en fin a cualquier investigador o experimentalista que sabe manejar con soltura un ordenador y le interesa aprender a programar.

Durante el curso se aprenderá a solucionar de forma más sencilla problemas comunes que son complicados de resolver con otras herramientas como hojas de cálculo o editores de texto. Por ejemplo, cambiar de formato un fichero con miles de secuencias, extraer unos pocos datos de interes de un fichero de cientos de megas u extraer y ordenar datos a nuestro antojo de una tabla con miles de resultados.

El curso original está preparado en inglés, por ello el código y ejemplos son en inglés, no obstante creo que cualquier lector de este blog puede entenderlos. Además el código serán imágenes para forzar al lector a escribir los ejemplos y no limitarse a copiar y pegar.

¿Qué es Python?

En Wikipedia encontraremos una buena definición:

Python es un lenguaje de alto nivel y de propósito general cuya filosofía de diseño se centra en la facilidad de lectura y escritura de código, permitiendo a los programadores expresar sus conceptos en pocas líneas de código. Esta filosofía contrasta con la complejidad de lenguages como C++ o Java.

¿Para qué sirve Python? 

Saber programar en biología es sumamente útil, de ahí la importancia de la bioinformática y la biología computacional. A continuación mostraremos algunos buenos ejemplos de su utilidad.

Ejemplo 1. Supongamos tenemos un fichero FASTA con cientos de secuencias de mRNA y queremos traducirlas a secuencias de proteínas. Conociendo el código genético del organimo en cuestión podemos hacer un sencillo script de Python que nos genere automáticamente otro fichero FASTA con las secuencias de proteínas.

Ejemplo 2. En este caso queremos realizar un experimento de expresión diferencial de genes en dos condiciones diferentes, por ej. en un paciente sano y en otro enfermo. Para ello usamos microarrays (o chips de DNA) y obtenemos un fichero con las intensidades de las sondas que hibridan con los genes. El valor de las miles de intensidades no nos permite sacar ninguna conclusión valiosa a mano, pero si agrupamos los genes por familias o función, podremos ver qué mecanismos celulares están afectados, por ejemplo una vía apoptótica en un caso de cáncer. Para realizar dichos análisis podemos ayudarnos de Python.

Ejemplo 3. Ahora que están de moda las nuevas tecnologías de secuenciación, éstas nos proporcionan millones de secuencias cortas (menos de 300bp) que debemos ensamblar o alinear con un genoma de referencia. Aunque existen muchas herramientas y muy eficientes para estos fines (escritas principalmente en C++, un lenguaje mucho más rápido que Python), podemos escribir scripts en Python que nos ayuden a procesar estas secuencias antes de ser alineadas o a analizar los resultados del alineamiento.

¿Perl o Python? 

En mi opinión ambos lenguajes de programación son igualmente válidos. Python tiene la ventaja de ser rápido y sencillo de escribir, y también de leer y modificar posteriormente el código. Aunque personalmente me encuentro más cómodo con Perl, mucha gente actualmente usa Python. La elección de uno u otro va a depender principalmente de la existencia previa de utilidades que nos interesen para nuestra investigación en uno u otro lenguaje. Por ejemplo, si trabajamos en biología estructural elegiremos Python porque existen muy buenos módulos para leer y manejar ficheros PDB. Y si trabajamos en filogenética posiblemente preferiremos Perl porque las utilidades existentes para este lenguaje son más completas. En la siguiente imagen podemos ver como un código sencillo de Python es más sencillo que uno de Perl, pero no nos engañemos, un código de Python mal indentado puede ser un infierno.

¿Cómo instalar Python? 

Para empezar este curso instalaremos la versión estable más reciente de Python desde su página oficial:

 

Durante la instalación será muy importante activar la opción "Añadir python.exe al PATH":

Comprobaremos la existencia del nuevo grupo de programas en nuestro escritorio:

¿Cómo ejecutar Python?

Entre las utilidades instaladas por Python encontraremos: "Python GUI" y "Python Command line", al ejecutarlas su aspecto es el siguiente:

Para empezar a programar y realizar los siguientes ejemplos recomiendo "Python GUI", pero "Python Command line" es igualmente válida. Podemos probar a escribir unas simples operaciones aritméticas y pulsar "Enter" como se muestra en la imagen:

Usando Python como una calculadora

En el anterior ejemplo vimos como Python se puede usar como calculadora. A continuación realizaremos unas cuantas operaciones aritméticas y asignaremos valores numéricos a las palabras 'width' y 'heigth' para posteriormente calcular su producto. Notar también el uso del cáracter almohadilla ´#´, todo el texto que escribamos a continuación de éste carácter no será leído ni ejecutado en el programa de Python, sirve para insertar comentarios en el código.

Palabras y frases en Python (cadenas)

En Python llamamos 'cadenas' (strings) cuando en vez de números trabajamos con palabras, frases o textos. Este tipo de datos se escribe siempres entre comillas simples ('texto') o dobles ("texto"). Podemos incluir un tipo de comillas dentro del otro, pero si queremos incluir unas comillas del mismo tipo que las que delimitan el texto, deberemos usar el carácter barra invertida '\' y a continuación la comilla (ej. "\""), la barra invertida se denomina carácter de escape y tendrá más usos como se muestra a continuación.

Si en un texto incluimos la barra invertida seguida de 'n' ('\n') insertaremos un salto de línea. Pero para que el salto de línea sea efectivo, necesitaremos utilizar nuestro primer comando de Python: 'print'. Un comando ejecuta una acción sobre unos argumentos que se incluyen entre paréntesis, en este caso 'print' ejecuta la acción de imprimir el texto que se incluye entre paréntesis. Escribir texto entre comillas en el ejemplo anterior nos ha servido para practicar, pero a partir de ahora para imprimirlo lo escribiremos dentro de los paréntesis del comando 'print' y conservando las comillas dobles o sencillas.

Variables en Python

Una variable es un valor/datos almacenados en la memoria del ordenador (ej. números o cadenas), la localización del valor se asocia a un nombre o identificador. El nombre de la variable sirve para referirnos a los datos que guarda, y estos datos pueden ser modificados o reemplazados por otros nuevos.

Los nombres de variables sólo pueden contener letras en su comienzo y letras, número y '_' (barra baja) en el resto del nombre (ver ejemplos en la siguiente imagen).

Juguemos un poco con variables, veamos como guardar textos en ellas y concatenarlos. También como guardar números y sumarlos, restarlos... además como imprimir los resultados con el comando 'print'.

Para finalizar la primera lección, guardaremos la palabra "Python" en la variable 'word'. Veamos cómo acceder a cada una de sus letras según su posición (empezando siempre desde la posición 0 para la primera letra 'P') y veamos qué pasa al intentar acceder a una posición que no existe o intentar modificar una letra de la variable.

Próxima lección

Esto es todo por hoy, en la próxima lección más y mejor...

Ofertas de trabajo en bioinformática

[He movido el contenido de esta entrada a la página de ofertas laborales]

La PCR más grande del mundo

En la entrada de hoy veremos la utilidad del comando de linux AWK, que fue utilizado en el anterior post para extraer secuencias aleatorias de un archivo FASTA o FASTQ.

AWK es más que un comando, es un poderoso lenguaje usado para manipular y procesar grandes ficheros de texto. Si queréis más información, en internet hay muy buenos tutoriales.

GAWK es la versión GNU de AWK que usan todas las distribuciones linux, se puede ver que 'awk' en Ubuntu redirige a 'gawk':

 $ which awk  
 /usr/bin/awk  
 $ ls -l /usr/bin/awk  
 /usr/bin/awk -> /etc/alternatives/awk  
 $ ls -l /etc/alternatives/awk  
 /etc/alternatives/awk -> /usr/bin/gawk  

A lo que íbamos, el título del post es 'La PCR más grande del mundo' porque la idea es mostrar cómo hacer una PCR in silico usando AWK.

Para ello prepararemos 3 ficheros, uno con los primers, otro con las secuencias con las que vamos a simular la PCR (ambos en formato FASTA) y otro con el código AWK.

Ejemplo de fichero 'primers.fa':

 >Marker_1  
 GTTGTGTCTTTAGCTCCCCTG(.+)TA[G|T]ATCAGGTGAGCCGAT  
 >Marker_2  
 ACTACAACCAGAGCGAGG(.+)AG[A|T]TC[C|G]CCCAAAGGCACA  

Observar que hay que poner el primer directo en sentido 5'->3' y el reverso en sentido 3'->5', tal como los veríamos en la secuencia 5'->3' del archivo 'secuencias.fa' donde los vamos a buscar.

Ejemplo de fichero 'secuencias.fa':

 >QK4PW:00814:02853  
 AAGACAGTTGTGTCTTTAGCTCCCCTGAGCTGAACGACATCAGTACATCAGGTCCAACTTCTACAACAAGCTGGAGCTTTTCAGGTTTGACAGCAACCTGGGGAGTTTGTTGGATACACAGAATATGGAGTGAAACAAGCTGAATACAGGAACAACAACCCGTCATATATCGCATCACTGAGAGCTCGAGAGGGCAGACCGCCTGCCTGCACAACATTGGTATTGACTACCAGAACATTCTGACTAGATCAGGTGAGCCGATTCGGTT  
 >QK4PW:03633:02817  
 TCACTCGTTGTGTCTTTAGCTCCCCTGAGCTGAAGGACATTCAGTACATCAACTCCTATTATTACAACAAGCTGGAATGGGCCAGGTTTGACAGCAACGTGGGTAGATATGTTGGATACACGAAGTTCGGAGTTTATAATGCAGAACGATGGAACAAAGACCCGTCAGAGATCGCTATGAGGATGGCTCAGAGGGAGACCTACTGCCTGCACAACATTGGTAACTGGTACCCAAACATGCTGACTAGATCAGGTGAGCCGATCGCGTT  
 >QK4PW:01993:02583  
 CACAGTGTTGTGTCTTTAGCTCCCCTGAGCTCAAAGACATCCAGTACATCAGGTCCTATTATTACAACAAGCTGGAGTTCATCAGGTTTGACAGCAACGTGGGGGAGTTTGTTGGATACACGGAGCTGGGAGTGAAGAACGCTGAGCGGCTCAACAACGACCCGTCACAGATCGCTGGGTTGAAAGCTCAGAGGGAGGCCTACTGCATGAACCATGTTACTGCTTTCTACCCAAACGCTCTGGATATATCAGGTGAGCCGATCAGCGG  
 >QK4PW:01336:01442  
 CTCATGACTACAACCAGAGCGAGGGCGGCTCTCACACCATCCAGGAGATGTATGGCTGTGACGTGGGGTCGGACGGGCGCCTCCTCCGCGGGTACAGTCAGTTCGCCTACGATGGCCGCGATTACATCGCCCTGAACGAGGACCTGACGACATGGACGGCGGCGGACACGGCGGCGCAGATCTCCCAGCGCAAGTTGGAGCAGGCTGGTGCTGCTGAGATATACAAGGCCTACCTGGAGGACGAGTGCGTGCAGTGGCTGCGCAGACACCTGGAGAACGGGAAGGGACGAGCTGCTGCGCACAGTTCGCCCAAAGGCACAATGTTA  
 >QK4PW:00383:00875  
 TCGACGACTACAACCAGAGCGAGGGCGGCTCTCACACCATCCAGGAGATGTATGGCTGTGACGTGGGGTCGGACGGGCGCCTCCTCCACGGGTATTATCAGTTCGCCTACGACGGCCGCGATTACATCGCCCCTGAACGAGGACCTGAAGACGTGGACGGCAGCGGACACGGCGGCGCAGATCACCCAGCGCAAGTGGGAGCAGGCTGGTGATGCAGAGAGTCTGAAGGCCTTCCTGGAGGGCGAGTACACGCAGTGGCTGCGCAGATTCCTGGAGATCGGGAAGGACGCGCTGCTGCGCACAGTTCCCCCAAAGGCACACACATA  

Fichero con el código AWK 'pcr.awk':

 # Primero leemos el archivo FASTA con los primers  
 # Y guardamos las secuencias de los primers en un array  
 FNR==NR{  
      if ($0~/>/){  
           primer_name=$0;  
      } else if ($0) {  
                primers[primer_name]=$0;  
      }  
      next;  
 }  
 # Después procesamos el archivo FASTA de secuencias  
 # Buscando los primers en cada secuencia  
 {  
      if ($0~/>/) {  
           query_name=$0;  
      } else if ($0) {  
           for (primer_name in primers) {  
                if (p=match($0, primers[primer_name])) {  
                     print query_name"\t"primer_name"\t"RSTART"\t"RLENGTH"\t"substr($0, RSTART, RLENGTH);  
                }  
           }  
      }  
 }  

Y ejecutaremos el siguiente comando:

 $ awk -f pcr.awk primers.fa secuencias.fa  

O en un sola línea podemos poner todo el código AWK:

$ awk 'FNR==NR{ if ($0~/>/){primer_name=$0;} else if ($0) {primers[primer_name]=$0;} next;} { if ($0~/>/) {query_name=$0;} else if ($0) { for (primer_name in primers) { if (p=match($0, primers[primer_name])) { print query_name"\t"primer_name"\t"RSTART"\t"RLENGTH"\t"substr($0, RSTART, RLENGTH); } } } }' primers.fa secuencias.fa  

El resultado será el siguiente:

 >QK4PW:00814:02853    >Marker_1    7    256   GTTGTGTCTTTAGCTCCCCTGAGCTGAACGACATCAGTACATCAGGTCCAACTTCTACAACAAGCTGGAGCTTTTCAGGTTTGACAGCAACCTGGGGAGTTTGTTGGATACACAGAATATGGAGTGAAACAAGCTGAATACAGGAACAACAACCCGTCATATATCGCATCACTGAGAGCTCGAGAGGGCAGACCGCCTGCCTGCACAACATTGGTATTGACTACCAGAACATTCTGACTAGATCAGGTGAGCCGAT  
 >QK4PW:03633:02817    >Marker_1    7    256   GTTGTGTCTTTAGCTCCCCTGAGCTGAAGGACATTCAGTACATCAACTCCTATTATTACAACAAGCTGGAATGGGCCAGGTTTGACAGCAACGTGGGTAGATATGTTGGATACACGAAGTTCGGAGTTTATAATGCAGAACGATGGAACAAAGACCCGTCAGAGATCGCTATGAGGATGGCTCAGAGGGAGACCTACTGCCTGCACAACATTGGTAACTGGTACCCAAACATGCTGACTAGATCAGGTGAGCCGAT  
 >QK4PW:01993:02583    >Marker_1    7    256   GTTGTGTCTTTAGCTCCCCTGAGCTCAAAGACATCCAGTACATCAGGTCCTATTATTACAACAAGCTGGAGTTCATCAGGTTTGACAGCAACGTGGGGGAGTTTGTTGGATACACGGAGCTGGGAGTGAAGAACGCTGAGCGGCTCAACAACGACCCGTCACAGATCGCTGGGTTGAAAGCTCAGAGGGAGGCCTACTGCATGAACCATGTTACTGCTTTCTACCCAAACGCTCTGGATATATCAGGTGAGCCGAT  
 >QK4PW:01336:01442    >Marker_2    7    314   ACTACAACCAGAGCGAGGGCGGCTCTCACACCATCCAGGAGATGTATGGCTGTGACGTGGGGTCGGACGGGCGCCTCCTCCGCGGGTACAGTCAGTTCGCCTACGATGGCCGCGATTACATCGCCCTGAACGAGGACCTGACGACATGGACGGCGGCGGACACGGCGGCGCAGATCTCCCAGCGCAAGTTGGAGCAGGCTGGTGCTGCTGAGATATACAAGGCCTACCTGGAGGACGAGTGCGTGCAGTGGCTGCGCAGACACCTGGAGAACGGGAAGGGACGAGCTGCTGCGCACAGTTCGCCCAAAGGCACA  
 >QK4PW:00383:00875    >Marker_2    7    314   ACTACAACCAGAGCGAGGGCGGCTCTCACACCATCCAGGAGATGTATGGCTGTGACGTGGGGTCGGACGGGCGCCTCCTCCACGGGTATTATCAGTTCGCCTACGACGGCCGCGATTACATCGCCCCTGAACGAGGACCTGAAGACGTGGACGGCAGCGGACACGGCGGCGCAGATCACCCAGCGCAAGTGGGAGCAGGCTGGTGATGCAGAGAGTCTGAAGGCCTTCCTGGAGGGCGAGTACACGCAGTGGCTGCGCAGATTCCTGGAGATCGGGAAGGACGCGCTGCTGCGCACAGTTCCCCCAAAGGCACA  

Cada columna indica lo siguiente:

1. Nombre de la secuencia que unen los primers.
2. Nombre de la pareja de primers.
3. Posición de inicio de la amplificación.
4. Longitud de la cadena amplificada.
5. La secuencia amplificada.

Este sencillo programita permite buscar 600000 secuencias amplificadas por los primers del ejemplo en un archivo de 5 millones de secuencias en un tiempo récord de 7 minutos en un ordenador normal con 4GB de RAM.

ACTUALIZACIÓN:
Probando un código similar escrito en Perl (gracias Bruno por el apunte) descubrimos que Perl es mucho más rápido que AWK y sólo tarda 1 minuto medio!!!

perl pcr.pl primers.fa secuencias.fa

Código en Perl en el fichero 'pcr.pl':

 open(PRIMERFILE, $ARGV[0]);     
 while(<PRIMERFILE>){  
      chomp;  
      if (/>/){  
           $primer_name=$_;  
      } elsif ($_) {  
           $primers{$primer_name}=$_;   
      }  
 }  
 close PRIMERFILE;  
 open(SEQSFILE, $ARGV[1]);     
 while(<SEQSFILE>){  
      chomp;  
      if (/>/){  
           $query_name=$_;  
      } elsif ($_) {  
           foreach my $primer_name (keys %primers) {  
                if (/$primers{$primer_name}/) {  
                     printf("%s\t%s\t%s\t%s\t%s\n",$query_name,$primer_name,$-[0]+1,length($&),$&);  
                }  
           }  
      }  
 }  
 close SEQSFILE;  

Para finalizar siento decir que esto no es nada nuevo, existen numerosos programas y utilidades web para realizar simulaciones de PCR a partir de un conjunto de primers y un fichero FASTA de secuencias, pero creo que ninguno tan sencillo y rápido como el código de AWK para procesar millones de secuencias:

• UCSC In-Silico PCR: Además de realizar la PCR permite controlar el tamaño de las secuencias amplificadas y el número mínimo de bases exactas en 3' para cada primer. El enlace nos lleva a su versión de servidor web, se puede descargar aquí. Debo decir que su versión en línea de comandos no me ha dado buenos resultados.
• e-PCR: permite simular PCRs y encontrar secuencias de DNA que serían reconocidas de forma específica y no específica por un conjunto de primers. Posee versión web y descargable.
• EHU In-Silico PCR: Permite con su interfaz web seleccionar un genoma, un par de primers y obtener los productos de los mismos.
• FastPCR: un programa comercial para el diseño de primers y su testeo in silico. Quizás sea el programa más completo si estamos interesados en el diseño de PCRs y nuestro laboratorio puede comprar el software.
• Primer-BLAST: sirve para diseñar primers y probarlos contra un fichero de secuencias por medio de BLAST para buscar todos los posibles productos amplificados (normalmente buscaremos subproductos generados por uniones no específicas de los primers).
• Simulate_PCR: un script en Perl que tomo como entrada un listado de primers y un fichero FASTA y devuelve los posibles productos de PCR usando alineamientos de BLAST.
• GASSST (Global Alignment Short Sequence Search Tool): genera alineamientos globales de un conjunto de primers contra un fichero de secuencias. Parseando sus resultados con un script podremos reconstruir productos amplificados por uniones específicas y no-específicas de los primers. La dificultad de tener que procesar los resultados es su principal problema, pero realiza alineamientos globales de secuencias muy cortas (como los primers) que no consigue casi ningún programa.

Artículo homenaje al inventor de la PCR, Kary Mullis.

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

kseq::klib to parse FASTQ/FASTA files

Buenas,
éste tiene pinta de ser el último artículo del año, no? Desde hace tiempo venimos usando en el laboratorio la librería kseq.h , integrante de [ https://github.com/attractivechaos/klib ], para las desagradecidas tareas de manipular grandes archivos en formato FASTQ. Esta librería es un único archivo .h en lenguaje C que de manera muy eficiente permite leer archivos FASTQ, y también FASTA, incluso comprimidos con GZIP, con zlib. La eficiencia es fundamental para estas aplicaciones, porque modificar archivos .fastq de varios GB puede consumir mucho tiempo.

Fuente: http://drive5.com/usearch/manual/fastq_files.html

Como resultado de esta experiencia hace unos meses colgamos de la web del labo un programita, que llamamos split_pairs. Este software facilita trabajos tediosos como pueden ser modificar las cabeceras de una archivo FASTQ con expresiones regulares para que un programa posterior no se queje, y de paso emparejar o interlazar los reads según convenga.  Está escrito en C++ y tiene asociado un script Perl.

De lo que quería hablar hoy es de que es posible utilizar kseq.h directamente desde nuestro lenguaje preferido. Para ellos es necesario emplear herramientas que permitan incluir, y compilar, código C mezclado con el de nuestro lenguaje interpretado. En el caso de Python estaríamos hablando por ejemplo de cython , aquí os muestro cómo hacerlo en Perl con el módulo Inline::CPP, dotando a kseq de una interfaz orientada a objetos para evitar pasar tipos de datos con tipos no estándar:

 package klib::kseq;  
 use Inline CPP => config => namespace => 'klib';    
 use Inline CPP => config => libs => '-lz';  
 use Inline CPP => <<'END_CPP';  
 /*  
 The Klib library [ https://github.com/attractivechaos/klib ], is a standalone and   
 lightweight C library distributed under MIT/X11 license.   
 Klib strives for efficiency and a small memory footprint.  
 This package contains C++ wrappers for components of Klib to make them easy to use from Perl scripts.  
 */  
 using namespace std;  
 #include <iostream>  
 #include <string>  
 #include <zlib.h>  
 #include "kseq.h"   
 /*   
 kseq is generic stream buffer and a FASTA/FASTQ format parser from Klib.  
 Wrapper written by Bruno Contreras Moreira.  
 */  
 KSEQ_INIT(gzFile, gzread)  
 class kseqreader
 {  
 public:  
    kseqreader()  
    {  
       fp = NULL;  
       fromstdin=false;  
       openfile=false;  
       filename = "";  
    }  
    ~kseqreader()  
    {  
       close_stream();  
    }  
    void rewind_stream()  
    {  
       if(openfile == true)  
       {  
          if(fromstdin == false) kseq_rewind(seq);  
          else  
             cerr<<"# reader: cannot rewind STDIN stream"<<endl;  
       }  
    }           
    void close_stream()  
    {  
       if(openfile == true)  
       {  
          kseq_destroy(seq);  
          if(fromstdin == false) gzclose(fp);  
          fp = NULL;  
          fromstdin=false;  
          openfile=false;  
          filename = "";  
       }     
    }  
    // accepts both FASTA and FASTQ files, which can be GZIP-compressed  
    // returns 1 if file was successfully opened, 0 otherwise  
    int open_stream(char *infilename)  
    {  
       filename = infilename;  
       if(filename.empty())  
       {   
          cerr<<"# reader: need a valid file name"<<endl;  
          filename = "";  
          return 0;  
       }     
       else  
       {  
          FILE *ifp = NULL;  
          if(filename == "stdin")  
          {  
             ifp = stdin;  
             fromstdin = true;  
          }  
          else ifp = fopen(filename.c_str(),"r");  
          if(ifp == NULL)     
          {  
             cerr<<"# reader: cannot open input file: "<<filename<<endl;  
             filename = "";  
             return 0;  
          }  
          else // open file and initialize stream  
          {  
             fp = gzdopen(fileno(ifp),"r");   
             seq = kseq_init(fp);   
             openfile=true;  
             return 1;  
          }  
       }     
    }  
     // returns length of current sequence, which can be 0,  
    // -1 in case of EOF and -2 in case of truncated quality string  
    int next_seq()  
    {   
       return kseq_read(seq);   
    }  
     const char *get_sequence()  
    {  
       return seq->seq.s;  
    }  
    // returns FASTA/FASTQ header up to first found space  
    const char *get_name()  
    {  
       return seq->name.s;  
    }  
    const char *get_quality()  
    {  
       if(seq->qual.l) return seq->qual.s;  
       else return "";   
    }  
    const char *get_comment()  
    {  
       if(seq->comment.l) return seq->comment.s;  
       else return "";   
    }  
    const char *get_full_header()  
    {  
       _header = seq->name.s;  
       if(seq->comment.l)  
       {   
          _header += " ";  
          _header + seq->comment.s;  
       }     
      return _header.c_str();  
    }  
 private:  
    gzFile fp;  
    kseq_t *seq;  
    bool fromstdin;  
    bool openfile;  
    string filename;  
    string _header;  
 };  
 END_CPP  
 1;  

Ahora estamos en disposición de probar a leer un archivo con este módulo:

 use strict;  
 use klib;  
 my ($length,$header,$sequence,$quality);  
 my $infile = "sample.fq.gz";  
 # my $infile = "stdin";  
 my $kseq = new klib::kseq::kseqreader();  
 if($kseq->open_stream($infile))  
 {  
   while( ($length = $kseq->next_seq()) != -1 )  
   {  
     ($header,$sequence,$quality) =  
       ($kseq->get_full_header(),$kseq->get_sequence(),$kseq->get_quality());   
   }  
 }  
 $kseq->close_stream()  

Qué ganancia en tiempo de procesamiento de archivos FASTQ se obtiene? Una prueba realizada con el módulo Benchmark , con diez repeticiones y un archivo .fastq.gz con 250K secuencias sugiere que kseq::klib es el doble de rápido:

            s/iter  perlio   kseq
perlio    1.32      --     -50%
kseq    0.660  100%     --

Hasta el próximo año, un saludo,
Bruno