Curso de Python para biólogos - Lección 3. Data comparison and conditional statements

Comenzamos la semana con una nueva lección de Python, todavía no quedan un par de lecciones de conceptos básicos antes de empezar a hacer algún pequeño programa para resolver problemas biológicos. En la presente lección explicaremos cómo se pueden comparar datos y condicionar la ejecución del código de acuerdo a dichas comparaciones.

Expresiones booleanas

Una expresión booleana es una comparación entre 2 valores o variables que devuelve un resultado de cierto (true) si ambos valores son iguales o ambas variables contienen el mismo valor y falso (false) si son diferentes. True y false son los dos únicos valores de otro tipo de datos llamado booleano (pero no queremos explicar más tipos de datos, con los ya vistos es suficiente: números, cadenas, listas y diccionarios).

Operadores de comparación y lógicos

Se denominan operadores de comparación a las expresiones mostradas en la imagen (==, !=, <, >, <= y >=), permiten realizar comparaciones de dos valores o variables con 2 únicos posibles resultados: true o false.

Los comparadores lógicos son 'and', 'or' and 'not' y permiten combinar diferences comparaciones. En la tabla se muestran los posibles resultados. 'and' da como resultado true si todas las comparaciones son ciertas y false si una de ellas es falsa. 'or' devuelve true si al menos una de las comparaciones es cierta y false si todas son falsas. 'not' cambia el resultado de una comparación de true a false o al contrario.

 Veamos unos ejemplos:

Sentencias condicionales

Una sentencia condicional decide si se ejecuta o no el código indentado que le sigue a continuación, si la comparación es cierta el código se ejecuta, si es falsa no sucede nada. Las sentencias condicionales se escriben con el operador 'if' seguido de una comparación y termina con dos puntos ':'. En las siguientes líneas se escribe el código indentado con espacios o tabulador (recomendado 4 espacios). El código indentado sólo se ejecutará si la comparación da como resultado true.

Las sentencias condicionales se pueden hacer más completas añadiendo condiciones con el operador 'elif'. Si la primera condición (con 'if') no se cumple, se pasará a comprobar si se cumple la segunda (con 'elif') y así sucesivamente hasta que se agoten las condiciones. Si todas las comparaciones con 'if' y 'elif' son falsas, entonces se ejecutará el código especificado por el operador 'else'.

Para usar las sentencias condicionales con listas, o con valores de diccionarios (recordar dict.values() devuelve una lista de los valores asociados a las claves y dict.keys() devuelve una lista de todas las claves) existe la expresión if in :. Dicha expresión permite buscar un valor en una lista y si lo encuentra ejecuta el código indentado. Es una expresión muy útil que usaremos en futuras lecciones.

Condicionales anidados

Para terminar por hoy veremos une ejemplo donde hay una combinación de comparaciones y operadores lógicos muy larga y que hace confuso el código (lo que también propicia errores al escribirlo, por ejemplo si nos dejamos un paréntesis). El mismo código se puede escribir introduciendo unas sentencias condicionales dentro de otras, el código es más largo, pero también más claro y conciso. Cada sentencia condicional que se escribe debe ser indentada, si escribimos una sentencia condicional dentro de otra deberemos indentar el doble de espacios.

Próxima lección

En la próxima lección veremos como pedir al usuario que introduzca datos y cómo realizar bucles 'while'.

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.