Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Resolver programa con un contador

Perl aplicado a la bioinformática
Publicar una respuesta

Resolver programa con un contador

Notapor truku » 2013-03-24 04:46 @240

Hola a todos. Mirad: tengo un problema con un programa que me corre bastante prisa, porque lo tengo que tener listo hoy.

Este es el programa, y en él vienen explicados todos los pasos, y para lo que sirve (es muy cortito).
Sintáxis: (programa.pl) [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!usr/bin/perl
  2. use strict;
  3.  
  4.  
  5.  
  6. #Se quiere introducir la ruta de producción de baccatina III, un precursor del anticancerígeno taxol, en la microalga Chlamydiomonas. Para ello se quieren clonar los genes que nos interesan de la ruta de la baccatina III en un plásmido de este alga, con ayuda de las enzimas contenidas en el sitio de clonación múltiple. Este programa nos aportará información de la posición en la que se encuentran las dianas de  dichas enzimas en una secuencia de ADN en formato FASTA.
  7. #Es necesario saber la posición de inicio y fin, en la secuencia del organismo, del gen que queremos extraer para clonarlo en el plásmido de Chlamydomonas. Y en base a esa información se utilizarán las enzimas de restricción flanqueantes más próximas.
  8.  
  9. #Este programa se encarga de buscar, en una secuencia de ADN en formato FASTA que se encuentra en un fichero, las dianas de las enzimas contenidas en el sitio de clonación múltiple (MCS) del plásmido en el que queremos clonar nuestro gen de interés. EL programa comparara un fichero con el nombre de las enzimas en una columna, con las secuencias de las mismas localizadas en un archivo "./basedatos.txt". Y a su vez lo compara con el archivo que contiene la secuencia de ADN en FASTA e informa de la posición de inicio de la secuencia diana de dichas enzimas.
  10.  
  11. #El programa podría ser usado con otra base de datos de enzimas, siempre que tuviese el mismo formato, y que se modificase el nombre en el apartado "Cargar base de datos de las enzimas".
  12.  
  13. #Pasos que dará el programa:
  14. # 1) El programa nos pedirá el nombre del fichero que contiene la secuencia de ADN en formato FASTA.
  15. # 2) El programa cargará el archivo de rebase "./basedatos.txt".
  16. # 3) Cargar archivo que contiene la secuencia de ADN en formato FASTA.
  17. # 4) Compara las enzimas seleccionadas desde el rebase con la secuencia de ADN, y cuando encuentra la coincidencia, te dice la enzima que es, y la posición de inicio de la secuencia diana.
  18. # 5) Se creará un archivo en la misma carpeta donde esté el programa llamado "Resultado.txt"
  19.  
  20. # 1) Cargar nuestro archivo de enzimas: Inicialmente el programa nos pedirá el nombre del archivo en el que se encuentra el nombre de las enzimas contenidas en el MCS del plásmido.
  21. print "Dame el archivo con las enzimas de restricción del MCS del plásmido: ";
  22. my $archivo1 = <STDIN>;
  23. chomp $_;
  24.  
  25. #Abre el archivo que le pedimos, y si no se encuentra el archivo que le pedimos da error.
  26. my %ENZIMA;
  27. my $contador;
  28. open( fichero, $archivo1 ) || die "No se reconoce el archivo $archivo1\n";
  29.  
  30. #Se le asigna un valor ascendente a cada enzima del MCS para que al recorrer el rebase solo seleccione las enzimas que nos interesan cada una con su número correspondiente
  31. while (<fichero>) {
  32.     chomp $_;
  33.     $ENZIMA{$_} = $contador++;
  34.  
  35.     print "Control de errores: ¿Llama 1 a las enzimas de interés? $ENZIMA{$contador}\n ";
  36. }
  37. close fichero;
  38.  
  39. # 2) Cargar base de datos de las enzimas del archivo "./basedatos.txt"
  40. my %enz;
  41. my $nombre;
  42.  
  43. # 2.1) Array asociativo que soluciona el problema de las ambigüedades en las secuencias de nucleótidos de las enzimas del fichero "./basedatos.txt".
  44. # Codigos ambiguos:
  45. # R = G or A
  46. # Y = C or T
  47. # M = A or C
  48. # K = G or T
  49. # S = G or C
  50. # W = A or T
  51. # B = not A (C or G or T)
  52. # D = not C (A or G or T)
  53. # H = not G (A or C or T)
  54. # V = not T (A or C or G)
  55. # N = A or C or G or T
  56.  
  57. my %amb = (
  58.     "R" => "[AG]",
  59.     "Y" => "[CT]",
  60.     "M" => "[AC]",
  61.     "K" => "[GT]",
  62.     "S" => "[GC]",
  63.     "W" => "[AT]",
  64.     "B" => "[CGT]",
  65.     "D" => "[AGT]",
  66.     "H" => "[ACT]",
  67.     "V" => "[ACG]",
  68.     "N" => "[ACGT]"
  69. );
  70.  
  71. open in, "./basedatos.txt";
  72. while (<in>) {
  73.     chomp $_;
  74.  
  75.     #Eliminar "<1>" y "<3>", así como los "^" que indican el lugar de corte de la enzima para quedarnos unicamente con una línea que indica el nombre de la enzima --> $nombre, y otra linea que nos indica cual es la secuencia diana de esa enzima.
  76.     if ( $_ =~ /^<1>/ ) {
  77.         $nombre = $_;
  78.         $nombre =~ s/<1>//;
  79.     }
  80.     elsif ( $_ =~ /^<3>/ ) {
  81.         $_ =~ s/(<3>|\^)//g;
  82.  
  83.         #Se resuelven las ambigüedades
  84.         foreach my $letra ( keys %amb ) {
  85.             $_ =~ s/$letra/$amb{$letra}/g;
  86.         }
  87.  
  88.         #Para que te reconozca las enzimas que antes le has dado valor ascendente
  89.         if ($ENZIMA{$nombre} == (      #no se como poner que sea cada una de las enzimas que antes han sido numeradas#){
  90.                 $enz{$nombre} = $_;
  91.  
  92.                     #print "Control de errores: dame las secuencias de las enzimas que me interesan $enz{$nombre} \n";
  93.             }
  94.         }
  95.     }
  96.     close in;
  97.  
  98.     # 3) Cargar archivo que contiene la secuencia de ADN en formato FASTA
  99.     print "Dame el archivo con la secuencia: ";
  100.     my $archivo2 = <STDIN>;
  101.     chomp $_;
  102.     my $ADN;
  103.     open( adn, $archivo2 ) || die "No se reconoce el archivo $archivo2\n";
  104.     while (<adn>) {
  105.         chomp $_;
  106.         $ADN .= $_;
  107.     }
  108.     close adn;
  109.  
  110.     # 4) y 5) Compara las enzimas seleccionadas desde el rebase con la secuencia de ADN, y cuando encuentra la coincidencia, te dice la enzima que es, y la posición de inicio de la secuencia diana. El resultado se imprime en un archivo, que se localizará en la misma carpeta, llamado "Resultado.txt".
  111.     open OUTPUT, ">./Resultado.txt";
  112.     my $longtotal = length($ADN);
  113.     my $posini    = 0;
  114.     foreach my $n ( keys %enz ) {
  115.  
  116.         #print OUTPUT "La enzima $n:\n";
  117.         while ( substr( $ADN, $posini ) =~ /$enz{$n}/ ) {
  118.             my $delante = $`;
  119.             my $pos     = length($delante) + 1;
  120.             my $long    = length( substr( $ADN, $posini ) );
  121.             $posini = $longtotal - $long + $pos;
  122.             print OUTPUT "Se encontró la enzima $n en la posición $posini\n";
  123.             print "Se encontró la enzima $n en la posición $posini\n";
  124.         }
  125.     }
  126.     close OUTPUT;
  127.  
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4


Mi problema creo que está en el contador. Con él intento dar un valor numérico a cada una de las líneas de un fichero, para que después, en base a ese valor, se repita un bucle un número de veces.

Lo que quiero es que me diga cuántas veces corta y dónde un número de enzimas determinado.

Nota: Para que funcione hace falta un cuarto archivo que se llame, por ejemplo, "enzimas.txt" en el que introducís esto que pongo a continuación:
Sintáxis: (enzimas.txt) [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
HindIII
KpnI
BamHI
BstXI
EcoRI
BstXI
NotI
AaaaI
AagI
AanI
AaeI
AasI
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4


Así tal y como os lo pongo, deben ir en líneas separadas; son el nombre de las enzimas que el programa va a buscar en la base de datos, porque son las que nos interesan.

Si alguien pudiese echarme una mano estaría eternamente agradecido :)

Nota 2: El programa en una versión anterior, en lugar de llamar valores ascendentes a esas enzimas del archivo que debéis crear, los llamaba 1 a todas, pero eso no funcionaba, porque coge siempre la misma enzima, y lo hace todo para esa sola. Yo lo que quiero es que se repita el trabajo para cada una de esas enzimas del archivo "enzimas.txt".

Muchas gracias, ¡saludos!
Adjuntos
fastahomo.txt
Es una secuencia de nucleótidos en FASTA en la que queremos buscar los sitios de corte de diferentes enzimas (se te pedirá el nombre en el terminal)
(147.62 KiB) 150 veces
basedatos.txt
Base de datos que usa el programa para funcionar
(173.02 KiB) 146 veces
Última edición por truku el 2013-04-02 09:30 @437, editado 2 veces en total
truku
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2013-03-21 06:58 @332
Arriba

Publicidad

Re: Resolver programa con un contador

Notapor danimera » 2013-03-24 14:11 @632

Hummmm... Creo que no está muy claro qué es lo que quieres hacer... ¿A qué te refieres cuando dices:

Lo que quiero es que me diga cuántas veces corta y dónde un número de enzimas determinado.
?

"cuántas veces corta" ¿se refiere a cuántas veces aparece el archivo? Y "donde" ¿se refiere a en qué posición del archivo?

La verdad no tengo claro lo que necesitas.
100% Telch - Perl Web Programming
Cali PerlMongers: http://cali.pm.org
Avatar de Usuario
danimera
Perlero frecuente
Perlero frecuente
 
Mensajes: 871
Registrado: 2005-06-23 19:02 @834
Ubicación: Colombia
Arriba

Re: Resolver programa con un contador

Notapor truku » 2013-03-24 14:43 @655

Sí, te explico: el nombre de las enzimas que están en el archivo que te tienes que crear aparte, se compara con las enzimas de la base de datos del programa (donde está la secuencia de cada una), y lo que quiero que me diga el programa es:

La enzima "tal" corta en la posición 1, 23, 344, 2333, 4324, etc.
La enzima "cual" corta en la posición 34, 443, 3244, 5234 .

Y así con cada una de las enzimas.

Muchas gracias por intentar ayudar. ¡Saludos!
truku
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2013-03-21 06:58 @332
Arriba

Re: Resolver programa con un contador

Notapor danimera » 2013-03-24 15:30 @687

Pero, ¿qué es posición? ¿La línea del archivo donde está la encima o alguna otra cosa?
100% Telch - Perl Web Programming
Cali PerlMongers: http://cali.pm.org
Avatar de Usuario
danimera
Perlero frecuente
Perlero frecuente
 
Mensajes: 871
Registrado: 2005-06-23 19:02 @834
Ubicación: Colombia
Arriba

Re: Resolver programa con un contador

Notapor truku » 2013-03-24 15:54 @704

No, la posición es el número de carácter de la secuencia de caracteres. Ejemplo:

AAAAABAAAAAAA

B está en la posición 6.
truku
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2013-03-21 06:58 @332
Arriba

Re: Resolver programa con un contador

Notapor explorer » 2013-03-25 20:55 @913

Yo sigo sin entender la parte del contador.

«Se le asigna un valor ascendente a cada enzima del MCS para que al recorrer el rebase solo seleccione las enzimas que nos interesan cada una con su número correspondiente»

¿Para qué necesitamos dar un número a las enzimas? Entiendo que las enzimas están como claves, pero, ¿los números como valores de esas claves?

Las líneas 90 a 94 no las entiendo. Y por otro lado, se introduce por teclado algo, pero luego las encimas que nos interesan las pasamos por un cuarto archivo...

En el archivo FASTA, la secuencia de ADN no contiene solo las famosas 'A', 'C', 'G' y 'T', sino que tiene algunas 'Y' y 'R'.

¡Qué lío!

Esta es mi versión, siguiendo los pasos indicados, al menos, los que he entendido.
Los nombres de los archivos los paso como argumentos, en lugar de preguntarlos al usuario.
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!usr/bin/perl
  2. #
  3. # Posiciones de enzimas. Joaquín Ferrero, 2013.
  4. #
  5. # Dada una base de datos de enzimas, una secuencia FASTA, y una serie de enzimas a buscar,
  6. # muestra la posición de esas enzimas dentro de la secuencia.
  7. #
  8. # Versión 20130325
  9. #
  10. use v5.14;
  11. use autodie;
  12. use File::Slurp;
  13.  
  14. ## Comprobación de los argumentos #############################################
  15. @ARGV == 3  or  die "Uso: $0 <archivo de REBASE> <archivo con la secuencia> <archivo de enzimas a buscar>\n";
  16. for (@ARGV) {
  17.     -f $_   or  die "ERROR: No encuentro el archivo $_: $!\n";
  18. }
  19. my($archivo_reb, $archivo_seq, $archivo_enz) = @ARGV;
  20.  
  21. ## Lectura del archivo REBASE #################################################
  22. my %ENZIMAS;
  23. {
  24.     my @rebase = read_file( $archivo_reb, chomp => 1 );
  25.  
  26.     unless (grep {/REBASE version \d+/} @rebase) {
  27.         die "ERROR: el archivo de REBASE indicado no está en ese formato.\n";
  28.     }
  29.  
  30.     my %amb = (
  31.         'R' => '(?:R|[AG])',
  32.         'Y' => '(?:Y|[CT])',
  33.         'M' => '(?:M|[AC])',
  34.         'K' => '(?:K|[GT])',
  35.         'S' => '(?:S|[GC])',
  36.         'W' => '(?:W|[AT])',
  37.         'B' => '(?:B|[^A])',
  38.         'D' => '(?:D|[^C])',
  39.         'H' => '(?:H|[^G])',
  40.         'V' => '(?:V|[^T])',
  41.         'N' => '.',
  42.     );
  43.  
  44.     my($nombre, $seq);
  45.  
  46.     for (@rebase) {
  47.  
  48.         when (/^<1>(.+)/) {
  49.             $nombre = $1;
  50.         }
  51.         when (/^<3>(.+)/) {
  52.             $seq = $1;
  53.             $seq =~ s/\^//g;
  54.             $seq =~ s/\(.+?\)//g;
  55.  
  56.             while(my($l, $s) = each %amb) {
  57.                 $seq =~ s/$l/$s/g;
  58.             }
  59.            
  60.             $ENZIMAS{ $nombre } = $seq;
  61.         }
  62.     }
  63. }
  64.  
  65. ## Lectura del archivo con la secuencia #######################################
  66. my $ADN;
  67.  
  68. $ADN =  read_file( $archivo_seq ) ;
  69. $ADN =~ s/\n//g;
  70. $ADN =  uc $ADN;
  71.  
  72. $ADN =~ /^[ACTGRYMKSWBDHVN]+$/  or  die "ERROR: el archivo de secuencia de ADN no es correcto.\n";
  73.  
  74. ## Lectura y procesamiento de las enzimas a extraer ###########################
  75. open my $RESULTA, '>', 'Resultado.txt';
  76. open my $ARCHIVO, '<', $archivo_enz;
  77.  
  78. while (my $nombre_enzima = <$ARCHIVO>) {
  79.     chomp $nombre_enzima;
  80.  
  81.     if (exists $ENZIMAS{ $nombre_enzima }) {
  82.         my $seq_enzima = $ENZIMAS{ $nombre_enzima };
  83.         my @posiciones;
  84.  
  85.         while ($ADN =~ /$seq_enzima/pg) {
  86.  
  87.             pos($ADN) += 1 - length ${^MATCH};
  88.  
  89.             #say "$nombre_enzima => ", pos($ADN), " => ", substr($ADN, pos($ADN)-1, length ${^MATCH});
  90.  
  91.             push @posiciones, pos $ADN;
  92.  
  93.             # Aquí hay una "truco": pos() nos devuelve la posición de la enzima dentro de la cadena de $ADN, pero
  94.             # basada en '0'. Por eso le sumamos '1' antes, para que se guarde en el informe el número de la posición
  95.             # pero basado en '1'. Y, además, como debemos reiniciar la búsqueda de la enzima en *una* posición más
  96.             # allá de donde fue encontrada, pos($ADN), al haber sido aumentada en '1', ya está realmente apuntando a
  97.             # la nueva posición de búsqueda.
  98.             # Sería así: pos($ADN) = pos($ADN) - 1 + 1; # que es lo mismo que no hacer nada.
  99.         }
  100.  
  101.         if (@posiciones) {
  102.             say $RESULTA "La enzima $nombre_enzima corta en la posición ", join ', ' => @posiciones;
  103.         }
  104.         else {
  105.             say $RESULTA "La enzima $nombre_enzima no corta en ninguna posición.";
  106.         }
  107.     }
  108.     else {
  109.         say $RESULTA "La enzima $nombre_enzima no aparece en la secuencia.";
  110.     }
  111. }
  112.  
  113. close $ARCHIVO;
  114. close $RESULTA;
  115.  
  116. __END__
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4

La salida es (recortada):
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
La enzima HindIII corta en la posición 7376, 8788, 10477, 11832, 23062, 25992, 28975, 34535, 36539, 36991, 37315,
La enzima KpnI corta en la posición 11290, 53163, 55852, 84475, 110468, 111577, 113623, 115299, 131292, 134708
La enzima BamHI corta en la posición 7626, 14031, 22729, 27592, 32279, 35302, 35714, 37430, 49317, 49412, 90385, 90798,
La enzima BstXI corta en la posición 467, 840, 3178, 4470, 5239, 5809, 6622, 7199, 8830, 10238, 12587, 15538, 15992,
La enzima EcoRI corta en la posición 8614, 15138, 16708, 17712, 22066, 27213, 27549, 34185, 39203, 41046, 43655, 49166,
La enzima BstXI corta en la posición 467, 840, 3178, 4470, 5239, 5809, 6622, 7199, 8830, 10238, 12587, 15538, 15992,
La enzima NotI no corta en ninguna posición.
La enzima AaaaI no aparece en la secuencia.
La enzima AagI corta en la posición 9712, 50164, 128956
La enzima AanI corta en la posición 6424, 10293, 11220, 15069, 19806, 20073, 20445, 27133, 32086, 34468, 35676, 36515,
La enzima AaeI corta en la posición 7626, 14031, 22729, 27592, 32279, 35302, 35714, 37430, 49317, 49412, 90385, 90798,
La enzima AasI corta en la posición 1097, 15194, 26133, 45605, 47772, 61933, 83274, 93863, 105546, 114264, 121899, 133593
Coloreado en 0.000 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: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Resolver programa con un contador

Notapor truku » 2013-03-26 07:01 @334

Wouuuu muchas gracias, explorer :wink: ,
Tu versión es algo compleja para mi nivel de Perl actual, pero es verdad que hace su trabajo tal y como debe. Muchas gracias por el aporte, ¡Un saludo!
truku
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2013-03-21 06:58 @332
Arriba

Re: Resolver programa con un contador

Notapor explorer » 2013-03-26 08:00 @375

Si no sabes qué es lo que hace alguna parte, pregunta.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14480
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba
Publicar una respuesta

Volver a Bioinformática

¿Quién está conectado?

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

Powered by phpBB® Forum Software © phpBB Group
Traducción al Español Mexicano por nextgen en colaboración con phpbb-es.com
phpBB SEO
Advertisements by Advertisement Management
cron