• Publicidad

Sitios de restricción de ADN

Perl aplicado a la bioinformática

Sitios de restricción de ADN

Notapor j_f_r_85 » 2010-12-06 13:16 @595

Compañeros:

Una vez más solicitando su ayuda. Deseo crear un script que me devuelva secuencias de ADN que se encuentra entre dos pequeñas secuencias del mismo patrón. Las secuencias de ADN se encuentran en un archivo .txt, y el usuario puede definir que patrón buscar y el archivo a introducir.

Es algo como esto: quiero extraer las secuencias que se encuentran entre el patrón ATC

CATCTAAAAGTGCGTAGATC
GTGTACCAACAGGAGGCTTG
TGTCAGATGAATACACTTGT
TGGGCGTGTTATTAATAAGA
ACTCGCATTCGCCTAGAGAA
C
ATCTAAAAGTGCGTAGATC
GTGTACCAACAGGAGGCTTG
TGTCAGATGAATACACTTGT


Deseo que se extraigan secuencias de toda la secuencia, y línea por línea. Además, deseo conocer la longitud de cada secuencia extraída. Algo así:

TAAAAGTGCGTAG 13
GTGTACCAACAGGAGGCTTGTGTCAGATGAATACACTTGTTGGGCGTGTTATTAATAAGAACTCGCATTCGCCTAGAGAAC 81
TAAAAGTGCGTAG 13



Lo más a lo que he llegado es a este script:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. # digest.pl
  3.  
  4. use strict; use warnings;
  5.  
  6. die "usage: digest.pl <pattern> <file>" unless @ARGV ==2;
  7.  
  8. my $pattern = $ARGV[0];
  9.  
  10. print "The pattern is: $pattern\n";
  11.  
  12. open (FILE, $ARGV[1]);
  13.  
  14. while (<FILE>) {
  15.         if (/$pattern/../$pattern/) {
  16.         print "$. $_\n";
  17.         }
  18. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4


Pero no me funciona, además soy malo con las regex. ¿Alguna recomendación?

¡Muchas gracias de antemano!

UPDATE: coloqué mal el patrón, el patrón es ATC.
Última edición por j_f_r_85 el 2010-12-06 15:45 @698, editado 1 vez en total
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Publicidad

Re: Sitios de restricción de ADN

Notapor explorer » 2010-12-06 16:17 @720

Aquí hay un problema: el operador rango viene bien cuando hay una expresión en una línea, y otra expresión unas líneas más abajo. Pero no te devuelve las partes de la secuencia entre los patrones. Y tampoco te vale si los patrones están en la misma línea.

La solución pasa por leer, antes, toda la secuencia, a un único escalar. Le quitas los caracteres de fin de línea y ya te quedas solo con las bases. Y ya puedes aplicar la expresión regular para buscar, que sería algo así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
while ($secuencia =~ /$patrón(.+?)(?=$patrón)/g) {
    print "$1\n";
}
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

La opción /g permite buscar la expresión regular de forma repetida a lo largo de la $secuencia. Y lo que buscamos es lo que hay entre dos patrones. El '(.+?)' quiere decir: "captura todo, pero solo antes del siguiente...". Y el (?=...) sirve para indicar que el segundo $patrón no formará parte de la captura, por lo que podemos detectar los casos en los que ese patrón forma tanto el final como el principio de la siguiente captura.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Sitios de restricción de ADN

Notapor j_f_r_85 » 2010-12-07 09:38 @443

Estimado explorer:

Buen día. Intenté modificar el script para que haga lo que deseo. Sé que aún hay trabajo por hacer. Este es el script que tengo al momento:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. # digest.pl
  3.  
  4. use strict; use warnings;
  5.  
  6. die "usage: digest.pl <pattern> <file>" unless @ARGV ==2;
  7.  
  8. my $pattern = $ARGV[0];
  9.  
  10. open (FILE, $ARGV[1]);
  11.  
  12. print "The pattern is: $pattern\n";
  13.  
  14. while (<FILE>) {
  15.         chomp;
  16.         while ($_ =~ /$pattern(.+?)(?=$pattern)/g) {
  17.     print "$1\n";
  18.         }
  19. }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Creo que me me da los sitios de restricción (fragmentos) para toda la secuencia al completo. ¿Qué opinas? Estoy usando un archivo de texto con 21 secuencias (una por línea) de 100 bases cada una.

Me imprime los fragmentos en pantalla. ¿Cómo podría hacer para que en un archivo me guardara solo los fragmentos y en otro solo las longitudes de cada fragmento?

Aún necesito avanzar para ver cómo podría hacer que haga esto pero línea por línea del archivo de texto. Ya que seguramente tendré archivos de texto con secuencias de 1000 bases por línea, y deseo conocer sitios de restricción y su longitud por cada secuencia. Espero poder seguir avanzando.

¡Gracias de antemano!
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Sitios de restricción de ADN

Notapor explorer » 2010-12-07 10:17 @470

Tal cual lo tienes puesto ahora, solo detectará las secuencias que estén completamente incluidas en una sola línea. No detectará las secuencias que ocupan más de una línea.

Fíjate que el bucle while() lo estás aplicando a la línea que has leído en el while() principal.

Esta es mi solución:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. #
  3. # Extracción de secuencias entre patrones
  4. #
  5. # JF 20101205.220008
  6. #
  7.  
  8. use common::sense;
  9. use File::Slurp;
  10.  
  11. ## Argumentos al programa
  12. @ARGV ==2  or  die "Uso: $0 <patrón> <fichero con la secuencia>\n";
  13.  
  14. my ($patrón, $fichero) = @ARGV;
  15.  
  16. $patrón =~ /^[ATCG]+$/  or  die "ERROR: el patrón no se compone de [ATCG]\n";
  17. -e $fichero             or  die "ERROR: no encuentro el fichero $fichero\n";
  18.  
  19.  
  20. ## leemos la secuencia
  21. my $secuencia = read_file($fichero);
  22.  
  23. $secuencia =~ s/\s+//g;         # quitamos blancos y finales de línea
  24.  
  25.  
  26. ## construimos la expresión regular
  27. # Al construirla fuera del bucle while(),
  28. # evitamos que sea recalculada en cada vuelta
  29. my $regex = qr/$patrón(.+?)(?=$patrón)/o;
  30.  
  31.  
  32. ## buscamos las secuencias entre patrones
  33. while ($secuencia =~ /$regex/g) {
  34.     say "[$1] ", length $1;
  35. }
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Sitios de restricción de ADN

Notapor j_f_r_85 » 2010-12-08 03:56 @206

Mi estimado explorer,

En serio que me has ayudado bastante, no solo corrigiendo los scripts. Si no también explicando los errores y soluciones. ¡Muchas gracias!

Seguiré practicando un poco más, por lo que creo que estaré preguntando más cosas.

Solo una cosa más por el momento; en mi caso, tenía deshabitada la característica "say". Entonces tuve que declarar su uso al principio del script, por lo que la primera parte de carga de bibliotecas me quedo así:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
use common::sense;
use File::Slurp;
use feature 'say';
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


¡Gracias nuevamente!
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Sitios de restricción de ADN

Notapor j_f_r_85 » 2010-12-08 06:06 @296

Estimado explorer,

Tengo una duda, ojalá me puedas guiar. Esta viendo las secuencias que me encuentra el script, así como los valores de su longitud. Y se me ocurría si es posible graficarlas en un histograma, algo así por ejemplo:

Tenemos esto:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
[TGTTGGGCGTGTTATTAATAAGA] 23
[CTACAATAATCTTTCT] 16
[TGTTGGGCGTGTTATTAATAAGA] 23
[CTACAATAATCTTTCT] 16
[CA] 2
[ATATAATGCATTCATATGTAATTATAAG] 28
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


¿Cómo se podría pedir al script que presente un histograma como el siguiente?

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
              23 |||||||||||||||||||||||||
              16 ||||||||||||||||
               2 ||
              28 ||||||||||||||||||||||||||||||
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


En este momento (tengo que pensar un poco más) se me ocurre que guardando los valores de longitud del fragmento en un hash, y con las siguiente instrucciones puede ser posible. ¿Qué opinas?

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
while (my ($keys, $values) = each %longitud) {
    print "$keys ", ('|' x $values), "\n";
}
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


También seria interesante hacer un pequeño y simple gráfico de la frecuencia de fragmentos encontrada en la salida después de los cortes de restricción. Se me ocurre que con una instrucción sort(), pero ¿cómo pedir por la frecuencia?

Bueno, por el momento es mucho preguntar, pensaré al respecto pero quisiera conocer tu opinión. ¡Gracias de antemano!
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Sitios de restricción de ADN

Notapor explorer » 2010-12-08 08:12 @383

Puedes generar el histograma "pintándolo" en una variable escalar, y luego imprimir todo. Algo así:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my $histograma;
  2.  
  3. while ($secuencia =~ /$regex/g) {
  4.     my $longitud = length $1;
  5.     $histograma .= sprintf "%3d %s\n", $longitud, ('|' x $longitud);
  6.  
  7.     say "[$1] ", $longitud;
  8. }
  9.  
  10. print $histograma;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Sitios de restricción de ADN

Notapor explorer » 2010-12-08 08:14 @384

j_f_r_85 escribiste:Solo una cosa más por el momento; en mi caso, tenía deshabitada la característica "say". Entonces tuve que declarar su uso al principio del script, por lo que la primera parte de carga de bibliotecas me quedo así:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
use common::sense;
use File::Slurp;
use feature 'say';
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


No hace falta. El módulo common::sense ya se encarga de eso.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Problema con sitios de restricción de ADN

Notapor j_f_r_85 » 2010-12-09 09:13 @426

Mi estimado explorer, he estado intentando ejecutar el script que me has ayudado a corregir, pero estoy obteniendo errores y una salida de ceros. Mira, según yo el script ha quedado así:

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2.  
  3. use common::sense;
  4. use File::Slurp;
  5.  
  6. @ARGV ==2  or  die "Uso: $0 <patrón> <fichero con la secuencia>\n";
  7.  
  8. my ($patron, $fichero) = @ARGV;
  9.  
  10. $patron =~ /^[ATCG]+$/  or  die "ERROR: el patrón no se compone de [ATCG]\n";
  11. -e $fichero             or  die "ERROR: no encuentro el fichero $fichero\n";
  12.  
  13. my $secuencia = read_file($fichero);
  14.  
  15. $secuencia =~ s/\s+//g;
  16.  
  17. my $regex = qr/$patron(.+?)(?=$patron)/o;
  18.  
  19. my $histograma;
  20.  
  21. while ($secuencia =~ /$regex/g) {
  22.     $histograma .= sprintf "%3d %s\n", $1, ('|' x $1);
  23.     say "[$1] ", length $1;
  24. }
  25.  
  26. print $histograma;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Pero la salida que obtengo es esta:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
Argument "ATCGTGATAACCGTGTCGAGGTGATCAGAGGCCGACGGAAATCAAACG" isn't numeric in sprintf at "la ruta donde se encuentra mi código" line 27.
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Este error una vez por línea, y después una columna de ceros así:

Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  0
  0
  0
  0
  0
  0
  0
  0
  0
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


¿Qué estoy haciendo mal?
¡Gracias por tu ayuda!
j_f_r_85
Perlero nuevo
Perlero nuevo
 
Mensajes: 10
Registrado: 2010-12-06 06:24 @309

Re: Sitios de restricción de ADN

Notapor wanako » 2010-12-09 11:58 @540

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. my $histograma;
  2. my $hist;
  3.  
  4. while ($secuencia =~ /$regex/g) {
  5.     $histograma = sprintf ("%s", '|' x length($1));
  6.     say "[$1]\t", length $1, "\t$histograma";
  7. }
  8.  
  9. say "\n==== O tal vez ====\n";
  10.  
  11. while ($secuencia =~ /$regex/g) {
  12.     $hist .= sprintf ("\t%d %s\n", length($1), '|' x length($1));
  13.     say "[$1]\t", length $1;
  14. }
  15. say $hist;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Puedo estar equivocado porque no conozco Bioinformática ni lo que se busca representar, para mejor formateo en un GNU/Linux debes usar 'tput' y 'cup' o crear otro archivo para el histograma, ¿gnuplot? a futuro es ideal.

Edito: sobraba un operador ".=", ahora está mejor, ahorramos 4 bytes :)
wanako
Perlero nuevo
Perlero nuevo
 
Mensajes: 27
Registrado: 2010-09-23 11:27 @519

Siguiente

Volver a Bioinformática

¿Quién está conectado?

Usuarios navegando por este Foro: No hay usuarios registrados visitando el Foro y 0 invitados

cron