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.

Escapar caracteres especiales en un texto antes de utilizar expresiones regulares

Hoy me he encontrado con un pequeño bug del script 'pfam_scan.pl' que sirve para buscar diferentes secuencias de proteínas en una base de datos de dominos de  Pfam (HMMs), todo ello usando la herramienta 'hmmscan' de HMMER3.

El error simplemente consistía en que una subrutina insertaba variables dentro de una expresión regular sin haber escapado los caracteres especiales previamente, con lo cual se insertaban directamente caracteres del tipo '(', '/', ... directamente dentro de la expresión regular y no funcionaba como era deseado.

Para evitarlo basta con escapar dichos caracteres especiales insertando el texto entre '\Q' y '\E' como bien nos ha sugerido Joaquín (más información en perlop):

 my $text_with_special_chars = 'abc$%()/def';  
 my $start_of_text_to_search = '$%';  
   
 $text_with_special_chars =~ /(\Q$start_of_text_to_search\E.+)/;  
   
 print $1; # $1 = '$%()/def'  

Buscar un elemento en un array de perl

En Perl no tenemos la útil función "in_array()" que posee PHP, por lo cual cuando queremos saber si un array contiene un determinado elemento tenemos que escribir un poquito más de código.

Aunque buscando en google "perl array search element" obtenemos muchas soluciones a este problema, aquí presentaré 2 muy sencillas.

Para saber únicamente si un elemento está o no está en el array:

 # Search an element inside an array and return 1 or 0 if the element is present or not, example: in_array(\@my_array, $element);  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my %items = map {$_ => 1} @$array; # create a hash out of the array values  
      my $exist = (exists($items{$search_for}))?1:0;  
      return $exist;
 }  

Si queremos además poder buscar expresiones regulares en todo el array y
conocer la o las posiciones de los resultados:

 # Search a regular expression inside an array and retrieve the positions of the hits, example: in_array(\@my_array, "/$expression/");  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my @pos;  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @pos = grep { $_ =~ /$search_for/ } @{$array};  
      } else {  
           @pos = grep { $_ eq $search_for } @{$array};  
      }  
      (scalar @pos)?return (1,\@pos):return (0);  
 }  

Otra opción sería usar el módulo Tie::IxHash que permite crear hashes manteniendo el orden de los elementos y buscar los índices de los mismos:

 # Search a regular expression inside an array and retrieve the positions of the hits, example: in_array(\@my_array, "/$expression/");  
 sub in_array {  
      my ($array,$search_for) = @_;  
      my @pos;  
      if ($search_for =~ /^\/.+\/$/){  
           $search_for =~ s/^\/|\/$//g;  
           @pos = grep { $_ =~ /$search_for/ } @{$array};  
      } else {  
           @pos = grep { $_ eq $search_for } @{$array};  
      }  
      (scalar @pos)?return (1,\@pos):return (0);  
 }