Curso de Python para biólogos - Lección 2. Nuevos tipos de datos: listas, diccionarios y matrices

Volvemos con una nueva lección...
por motivos de tiempo, esta vez voy a incluir las diapositivas en inglés, aunque traduciré las explicaciones ;)

También sería gratificante que si seguís el curso y os parece útil hiciérais algún comentario para ver si alguien lo usa o es un tutorial más de Python en internet...

En esta nueva lección introduciremos dos tipos más de datos: las listas y los diccionarios. Recordad que en la lección anterior vimos los números y las cadenas de texto, además de aprender a instalar Python en Windows y a ejecutar comandos en una consola

¿Qué es una lista?

Una lista es un listado de valores, y se introducen separados por comas y entre corchetes. Dichos valores suelen ser números o cadenas, pero también pueden ser otros tipos de datos como veremos posteriormente. En Perl las listas se llaman arrays o arreglos.

Para referirnos a un elemento de la lista, deberemos conocer su posición. Dicha posición se llama índice (index) y las posiciones de los elementos en las listas comienzan en 0 (como en las cadenas). Al igual que las cadenas, las listas pueden ser troceadas, pero además los elementos de una lista también pueden ser modificados.

¿Cómo modificar una lista?

Las listas poseen un conjunto de métodos o funciones que permiten hacer diversas cosas con ellas, por ejemplo insertar elementos, eliminarlos, extraerlos, ordenarlos, contarlos... Dichos métodos se aplican escribiendo el nombre de la lista sequido de un punto, el método y entre paréntesis sus parámetros si es que los requiere (el paréntesis es obligatorio pero puede estar vacío).

 Veamos unos ejemplos:

 Pero cuidado, a cada método hay que darle los parámetros adecuados entre paréntesis:

Veamos la diferencia entre el método 'pop' que extrae un elemento indicando su posición y lo borra de la otra forma de modificar el elemento sin borrarlo. También los métodos 'append' e 'insert':

Podemos borrar elementos de una lista con la función 'del', pero cuidado porque podemos borrar la lista completa... ('del' no es un método por ello no requiere punto ni paréntesis)

Copia por valor y por refencia

Creo que este concepto merece pararnos y explicarlo más en detalle. Normalmente estamos acostumbrados que cuando copiamos algo sea una copia por valor, es decir, que si modificamos la copia no se modifique el original. Pero éste no es el comportamiento por defecto de Python, por motivos de eficiencia Python siempre copia por referencia los datos si no se indica lo contrario. La copia por referencia es una copia de las posiciones en memoria de los datos, de forma que si modificamos posteriormente la copia también se mofica el original porque guardan los datos en el mismo lugar de la memoria.

Para hacer copia por valor de una lista deberemos utilizar el método 'copy'. Si usamos el símbolo '=' se realizará una copia por referencia.

¿Qué es un diccionario?

Un diccionario se parece a una lista, pero en vez de ordenar sus elementos o 'valores' e identificarlos por su posición o índice, un diccionario no ordena sus elementosy los identifica por la llamada 'clave'. La clave es normalmente un texto que nos ayuda a localizar su contenido asociado en el diccionario, este contenido es el valor. Los diccionarios en Perl se llaman hashes o arreglos asociativos.

Un diccionario se define entre llaves, y en su interior se incluye una lista de pares clave:valor separados por comas.

Los pares clave:valor de un diccionario NO están ordenados y  no podemos referirnos a ellos mediantes índices, la única forma de llamar a un valor del diccionario es conociendo su clave asociada. De la misma forma que con las listas, el comando 'del' permite borrar pares clave:valor del diccionario, o el diccionario completo.

¿Cómo modificar un diccionario?

Al igual que las listas, los diccionarios tienen un conjunto de métodos que nos permiten modificarlos.

Los métodos 'keys' y 'values' devuelven una lista de las claves o valores respectivamente, el método 'items' devuelve una lista de pares (clave, valor). Los métodos 'get', 'pop', 'clear' y 'copy' funcionan de forma similar a como lo hacen con las listas. Realicemos unos cuantos ejemplos:

Funciones para cadenas, listas y diccionarios

Existen un conjunto de funciones u operaciones comunes para estos tres tipos de datos, se listan a continuación:

 

Permiten conocer su longitud (o número de elementos o pares clave:valor), el tipo de datos guardado en una variable o convertir un tipo de datos en otro tipo diferente:

 

Listas anidadas o listas de listas

Una lista anidada, del inglés nested list, es una lista formada por varias listas, si éstas contienen sólo números podríamos pensar como si fuera una matriz:

Pero además de incluir en una lista otras listas, podemos incluir otros tipos de datos como un diccionario. En este caso no podremos acceder a los datos del diccionario usando índices.

  
Diccionarios anidados

Al igual que en el caso de las listas, podemos crear diccionarios de diccionarios (anidados). También podemos mezclar listas en diccionarios y otros tipos de datos.

Próxima lección

La cosa se va complicando, números, cadenas, listas, diccionarios... mix!!! y además cada tipo de datos tiene asociado unos métodos y operaciones para modificarlos. En el próximo capítulo la cosa se complica y aprenderemos a comparar los datos y poner condiciones al comportamiento de un programa.

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