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':
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':
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':
Fichero con el código AWK 'pcr.awk':
Y ejecutaremos el siguiente comando:
O en un sola línea podemos poner todo el código AWK:
El resultado será el siguiente:
Cada columna indica lo siguiente:
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!!!
Código en Perl en el fichero 'pcr.pl':
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:
- Nombre de la secuencia que unen los primers.
- Nombre de la pareja de primers.
- Posición de inicio de la amplificación.
- Longitud de la cadena amplificada.
- La secuencia amplificada.
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.faCó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:
Artículo homenaje al inventor de la PCR, Kary Mullis.
- 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.
