Foro - Perl en Español

por **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

#!usr/bin/perl
use strict;
 
 
 
#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.
#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.
 
#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.
 
#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".
 
#Pasos que dará el programa:
# 1) El programa nos pedirá el nombre del fichero que contiene la secuencia de ADN en formato FASTA.
# 2) El programa cargará el archivo de rebase "./basedatos.txt".
# 3) Cargar archivo que contiene la secuencia de ADN en formato FASTA.
# 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.
# 5) Se creará un archivo en la misma carpeta donde esté el programa llamado "Resultado.txt"
 
# 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.
print "Dame el archivo con las enzimas de restricción del MCS del plásmido: ";
my $archivo1 = <STDIN>;
chomp $_;
 
#Abre el archivo que le pedimos, y si no se encuentra el archivo que le pedimos da error.
my %ENZIMA;
my $contador;
open( fichero, $archivo1 ) || die "No se reconoce el archivo $archivo1\n";
 
#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
while (<fichero>) {
    chomp $_;
    $ENZIMA{$_} = $contador++;
 
    print "Control de errores: ¿Llama 1 a las enzimas de interés? $ENZIMA{$contador}\n ";
}
close fichero;
 
# 2) Cargar base de datos de las enzimas del archivo "./basedatos.txt"
my %enz;
my $nombre;
 
# 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".
# Codigos ambiguos:
# R = G or A
# Y = C or T
# M = A or C
# K = G or T
# S = G or C
# W = A or T
# B = not A (C or G or T)
# D = not C (A or G or T)
# H = not G (A or C or T)
# V = not T (A or C or G)
# N = A or C or G or T
 
my %amb = (
    "R" => "[AG]",
    "Y" => "[CT]",
    "M" => "[AC]",
    "K" => "[GT]",
    "S" => "[GC]",
    "W" => "[AT]",
    "B" => "[CGT]",
    "D" => "[AGT]",
    "H" => "[ACT]",
    "V" => "[ACG]",
    "N" => "[ACGT]"
);
 
open in, "./basedatos.txt";
while (<in>) {
    chomp $_;
 
    #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.
    if ( $_ =~ /^<1>/ ) {
        $nombre = $_;
        $nombre =~ s/<1>//;
    }
    elsif ( $_ =~ /^<3>/ ) {
        $_ =~ s/(<3>|\^)//g;
 
        #Se resuelven las ambigüedades
        foreach my $letra ( keys %amb ) {
            $_ =~ s/$letra/$amb{$letra}/g;
        }
 
        #Para que te reconozca las enzimas que antes le has dado valor ascendente
        if ($ENZIMA{$nombre} == (      #no se como poner que sea cada una de las enzimas que antes han sido numeradas#){
                $enz{$nombre} = $_;
 
                    #print "Control de errores: dame las secuencias de las enzimas que me interesan $enz{$nombre} \n";
            }
        }
    }
    close in;
 
    # 3) Cargar archivo que contiene la secuencia de ADN en formato FASTA
    print "Dame el archivo con la secuencia: ";
    my $archivo2 = <STDIN>;
    chomp $_;
    my $ADN;
    open( adn, $archivo2 ) || die "No se reconoce el archivo $archivo2\n";
    while (<adn>) {
        chomp $_;
        $ADN .= $_;
    }
    close adn;
 
    # 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".
    open OUTPUT, ">./Resultado.txt";
    my $longtotal = length($ADN);
    my $posini    = 0;
    foreach my $n ( keys %enz ) {
 
        #print OUTPUT "La enzima $n:\n";
        while ( substr( $ADN, $posini ) =~ /$enz{$n}/ ) {
            my $delante = $`;
            my $pos     = length($delante) + 1;
            my $long    = length( substr( $ADN, $posini ) );
            $posini = $longtotal - $long + $pos;
            print OUTPUT "Se encontró la enzima $n en la posición $posini\n";
            print "Se encontró la enzima $n en la posición $posini\n";
        }
    }
    close OUTPUT;
 
Coloreado en 0.005 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

AasIColoreado 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!

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!

por **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?

por **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.

por **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

#!usr/bin/perl
#
# Posiciones de enzimas. Joaquín Ferrero, 2013.
#
# Dada una base de datos de enzimas, una secuencia FASTA, y una serie de enzimas a buscar,
# muestra la posición de esas enzimas dentro de la secuencia.
#
# Versión 20130325
#
use v5.14;
use autodie;
use File::Slurp;
 
## Comprobación de los argumentos #############################################
@ARGV == 3  or  die "Uso: $0 <archivo de REBASE> <archivo con la secuencia> <archivo de enzimas a buscar>\n";
for (@ARGV) {
    -f $_   or  die "ERROR: No encuentro el archivo $_: $!\n";
}
my($archivo_reb, $archivo_seq, $archivo_enz) = @ARGV;
 
## Lectura del archivo REBASE #################################################
my %ENZIMAS;
{
    my @rebase = read_file( $archivo_reb, chomp => 1 );
 
    unless (grep {/REBASE version \d+/} @rebase) {
        die "ERROR: el archivo de REBASE indicado no está en ese formato.\n";
    }
 
    my %amb = (
        'R' => '(?:R|[AG])',
        'Y' => '(?:Y|[CT])',
        'M' => '(?:M|[AC])',
        'K' => '(?:K|[GT])',
        'S' => '(?:S|[GC])',
        'W' => '(?:W|[AT])',
        'B' => '(?:B|[^A])',
        'D' => '(?:D|[^C])',
        'H' => '(?:H|[^G])',
        'V' => '(?:V|[^T])',
        'N' => '.',
    );
 
    my($nombre, $seq);
 
    for (@rebase) {
 
        when (/^<1>(.+)/) {
            $nombre = $1;
        }
        when (/^<3>(.+)/) {
            $seq = $1;
            $seq =~ s/\^//g;
            $seq =~ s/\(.+?\)//g;
 
            while(my($l, $s) = each %amb) {
                $seq =~ s/$l/$s/g;
            }
            
            $ENZIMAS{ $nombre } = $seq;
        }
    }
}
 
## Lectura del archivo con la secuencia #######################################
my $ADN;
 
$ADN =  read_file( $archivo_seq ) ;
$ADN =~ s/\n//g;
$ADN =  uc $ADN;
 
$ADN =~ /^[ACTGRYMKSWBDHVN]+$/  or  die "ERROR: el archivo de secuencia de ADN no es correcto.\n";
 
## Lectura y procesamiento de las enzimas a extraer ###########################
open my $RESULTA, '>', 'Resultado.txt';
open my $ARCHIVO, '<', $archivo_enz;
 
while (my $nombre_enzima = <$ARCHIVO>) {
    chomp $nombre_enzima;
 
    if (exists $ENZIMAS{ $nombre_enzima }) {
        my $seq_enzima = $ENZIMAS{ $nombre_enzima };
        my @posiciones;
 
        while ($ADN =~ /$seq_enzima/pg) {
 
            pos($ADN) += 1 - length ${^MATCH};
 
            #say "$nombre_enzima => ", pos($ADN), " => ", substr($ADN, pos($ADN)-1, length ${^MATCH});
 
            push @posiciones, pos $ADN;
 
            # Aquí hay una "truco": pos() nos devuelve la posición de la enzima dentro de la cadena de $ADN, pero
            # basada en '0'. Por eso le sumamos '1' antes, para que se guarde en el informe el número de la posición
            # pero basado en '1'. Y, además, como debemos reiniciar la búsqueda de la enzima en *una* posición más
            # allá de donde fue encontrada, pos($ADN), al haber sido aumentada en '1', ya está realmente apuntando a
            # la nueva posición de búsqueda.
            # Sería así: pos($ADN) = pos($ADN) - 1 + 1; # que es lo mismo que no hacer nada.
        }
 
        if (@posiciones) {
            say $RESULTA "La enzima $nombre_enzima corta en la posición ", join ', ' => @posiciones;
        }
        else {
            say $RESULTA "La enzima $nombre_enzima no corta en ninguna posición.";
        }
    }
    else {
        say $RESULTA "La enzima $nombre_enzima no aparece en la secuencia.";
    }
}
 
close $ARCHIVO;
close $RESULTA;
 
__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, 133593Coloreado en 0.000 segundos,  usando GeSHi 1.0.8.4

por **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!

por **explorer** » 2013-03-26 08:00 @375

Si no sabes qué es lo que hace alguna parte, pregunta.

Foro - Perl en Español

Resolver programa con un contador

Resolver programa con un contador

Publicidad

Re: Resolver programa con un contador

Re: Resolver programa con un contador

Re: Resolver programa con un contador

Re: Resolver programa con un contador

Re: Resolver programa con un contador

Re: Resolver programa con un contador

Re: Resolver programa con un contador

¿Quién está conectado?