Foro - Perl en Español

por **Rodrigobiotec** » 2011-01-28 16:20 @722

Hola a todos. Me presento: estoy estudiando 2º de grado de biotecnología y tengo que terminar este trabajo para aprobar la asignatura.

El trabajo consiste en, dado un fichero FASTA identificar todos los ORF tanto en sentido normal como inverso ("reverse"). El problema está en que no me encuentra todos los ORF que debería. Sin embargo, los que encuentra coinciden con la solución. Por ejemplo, mi ORF nº 1 equivale al ORF nº 7 del resultado.

Gracias de antemano; necesito ayuda.

Pego aquí el bucle que me lee el programa:

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

$indica_orf = 0;                #Indicada cuando encuentra un orf
$orf = "";
 
 
my @var = ();                   #La usamos para guardar los orfn
my $count = 0;                  #Contador de orf
my $posi = 0;                   #Posicion inicial
my $posf = 0;                   #Posicion final
 
 
for (my $i = 0; $i < length $fasta; $i=$i+3) {
  my $str ="";
 
  $a = substr ($fasta, $i, 3);
 #Para que me vaya leyendo las letras de 3 en 3.
 
  if (($a eq 'ATG') and ($indica_orf == 0)) {           #Si $a es igual a ATG y el indicador esta a 0
    $indica_orf = 1;                                    #Pon el indicador a uno ----->Ahora estamos dentro de un orf
    $posi = $i+1;
    
    
  }
  if ($indica_orf == 1) {                               #Si estamos dentro de un orf
    $str = modulo::codon($a);                           #modulo::codon -> llamamos a la funcion codon del paquete.
    $orf = $orf.$str;                                   #Sustituye las 3 letras de a por su valor y va concatenando
    
 
  }
  if ( $a =~ /TA[AG]|TGA/i ) {                          #Cuando a sea igual a TAA TAG o TGA (codones de terminacion)
                                                        #Y estemos dentro de un orf
    if ($indica_orf == 1) {
      $indica_orf = 0;                                  #Sal del orf
      $orf=$orf."\n";                                   #Pon un salto de linea
      
     
      
      if ((length $orf)*3 >= $n) {              #$n = $ARGV[2] -> numero que se introduce como parametro.
        push (@var, $orf);                      #Si la longitud del orf por 3 (porque un aminoacido vale tres codones) es mayor que n, guarda el orf
        $count++;                               #Va contando los orf mayores q n
        $posf = $i;
        
        print SALIDA "$cabecera_numerica"."_"."$count [$posi - $posf]$cabecera_infor";  #Cabecera (con su contador y sus cosas)
        print SALIDA  modulo::linea ($orf);             #Llamada a una funcion para que me imprima un salto de linea cada 60 nucleotido (imprimi bonito)que no funciona
      }
    }  
   
    $orf = "";
  }
 
}
 
#############REVERSE##############
 
#Exactamente lo mismo, pero al revés (mismo problema)
 
my $indica_reverse = 0;
my $inicio = 0;
my $fin =0;
my $str_rvs = "";
my $orf_rvs = "";
 
 
$fasta_reverse = reverse $fasta;
 
for (my $i = 0; $i < length $fasta_reverse; $i = $i+3) {
  my $a = substr ($fasta_reverse, $i, 3);
 
  if (($a eq 'ATG')and($indica_reverse == 0)) {
    $indica_reverse = 1;
    $inicio = $i +1;
  }
  if ($indica_reverse ==1) {
    $str_rvs = modulo::codon($a);
    $orf_rvs = $orf_rvs.$str_rvs;
  }
  if ($a =~ /TA[AG]|TGA/) {
    if ($indica_reverse ==1) {
      $indica_reverse = 0;                      
      $orf_rvs = $orf_rvs."\n";
        if ((length $orf_rvs)*3 > $n) {
         push ($orf_rvs, @var);
          $fin = $i;
          $count++;
      
          print SALIDA "$cabecera_numerica"."_"."$count [$inicio - $fin]"."(REVERSE SENSE)"."$cabecera";
          print SALIDA linea::modulo($orf_rvs);
        }
      }
      $orf_rvs = "";
    }
  }
 
close FASTA;
close SALIDA;
Coloreado en 0.006 segundos,  usando GeSHi 1.0.8.4

Lo primero de todo recordar que ORF significa Open Reading Frames, una secuencia que no contiene una señal de stop (TGA, TAA o TAG) digamos que cuando se traducen funcionan como un todo.

Buscaríamos fragmentos que empiecen con ATG (la señal de inicio habitual) y acaben con estas señales de stop, leídas al derecho o del revés.

por **pvaldes** » 2011-01-28 17:55 @788

Un posible enfoque

hacemos una copia del archivo original
obtenemos la posición de todas las ATG

Si el contador inicial vale 0
mientras la cadena tenga algún ATG (while)
buscamos el primer ATG de la lista,
descartamos todo lo anterior acortando la cadena
leemos desde ahí de 3 e 3 tomando primero la A hasta encontrar el siguiente codón de stop. Guardamos la cadena entre ambos
quitamos la A y repetimos el proceso
quitamos la T y repetimos el proceso

pasamos al siguiente ATG de la lista, acortamos lo anterior y repetimos el proceso

etc

cuando la cadena ya no tenga más ATG o no queden más codones (su tamaño sea <7 letras), salimos del bucle y anotamos una repetición a un contador con valor inicial 0

Si el contador vale 1

volvemos a leer la cadena original, invertimos la cadena, y sobreescribimos la variable, a la que seguimos llamando del mismo modo, como es una copia no importa que se machaque, eso evita repetir código

y repetimos el bucle while, lo reutilizamos entero tal cual

Si el contador vale 2 (hemos acabado la segunda vuelta)
imprimimos los resultados deseados
salimos del programa

por **Rodrigobiotec** » 2011-01-28 18:39 @818

Gracias por la ayuda y por la rapidez, no entiendo muy bien lo que me propones; pensaba que debía buscar desde ATG hasta los codones de terminación y guardarlo así con todos los que me salgan así que lo de contar desde la T de ATG y desde la G de ATG no lo entiendo, lo siento pero tampoco comprendo ni sé cómo realizar lo que me propones en el segundo mensaje.

Gracias por todo espero respuesta.

por **pvaldes** » 2011-01-28 18:44 @822

Tienes que hacer tres pasadas hacia adelante y otras 3 hacia atrás porque puede arrancar la secuencia de tripletes desde la A, la T o la G indistintamente, así que por cada sentido de lectura hay tres modos de agrupar las letras de 3 en 3

http://es.wikipedia.org/wiki/Marco_abierto_de_lectura

por **Rodrigobiotec** » 2011-01-28 18:52 @827

Entiendo, muchas gracias por el enlace y la explicación, entonces debería hacer otros 2 for() más para la "ida" y otros dos para la "vuelta" ¿correcto ? si no, ¿cómo sería?

Gracias.

por **pvaldes** » 2011-01-28 19:01 @834

Exacto, y por cierto, ahora que lo pienso, la cadena complementaria no es igual que darle la vuelta a la cadena inicial, habría que traducirla además antes de iniciar la segunda vuelta.

por **explorer** » 2011-01-28 19:33 @856

Bienvenido a los foros de Perl en Español, Rodrigobiotec

¿Podrías publicar un breve trozo de secuencia de ADN donde ocurre el error con tu código?

por **Rodrigobiotec** » 2011-01-28 20:03 @877

Estoy intentando hacer lo que pvaldes me ha explicado y no consigo hacerlo. Esto es lo que he escrito pero no avanzo (está abajo reflejado). Respecto a la parte de la secuencia del ADN que me da error no es así exactamente, si no que no me coge todos los ORF que debería. Además meto por pantalla la longitud mínima de cada ORF que ha de ser 300, la cadena entera consta de 300.000 nucleótidos. No sé qué parte pegar.

Gracias por todo.

Esto es lo que llevo hecho.

Sintáxis: [ Descargar ] [ Ocultar ]

Using perl Syntax Highlighting

#!/usr/bin/perl -w
 
use modulo;
 
#$ARGV[0];
#$ARGV[1];
#$ARGV[2];
chomp $ARGV[2];
 
$NumARGV = $#ARGV;    #Ultima posicion de @ARGV (ultimo parametro))
 
if ( $NumARGV != 2 )
{    #Si el numero de argumentos es distinto de 2 (son 3 pero el primero es 0)
    print "ERROR !!!\nDebe intoducir tres parametros\n\n";
    exit;
}
 
if ( $ARGV[2] =~ /[^0-9]/ ) {    #Valor minimo del orf
    print "El tercer parametro debe ser un numero!!\n";
    exit;
}
 
my ($fasta) = $ARGV[0];  #En el primer argumento pasamos el ficherp de entrada
 
open( FASTA, $fasta ) or die("No puedo abrir fichero $fasta: $! \n");
open( SALIDA, ">$ARGV[1]" )
    or die("No puedo abrir fichero de salida\n")
    ;                    #El segundo argumento es un fichero de salida
 
@FASTA = <FASTA>;        #Convierto el fichero de entrada em un array
 
$cabecera = shift @FASTA;    #El primer elemnto del array va a ser la cabecera
 
$cabecera_numerica = modulo::obtener_numerico($cabecera)
    ;    #Llamamos a la funcion definida en el otro paquete
$cabecera_infor = modulo::obtener_infor($cabecera)
    ;    #Llamamos a la funcion definida en el otro paquete
 
$fasta = join( "", @FASTA );    #Uno todos los elementos restantes del array
$fasta =~ s/\s//g;              #Quito los espacios, saltos....
 
my ($n)
    = $ARGV[2]
    ; #$n es el parametro numero 3, Ejemplo: si $n = 300 todos los orf deben tener una longitud (en nucleotidos) mayor q 300
 
$indica_orf = 0;    #Indicada cuando encuentra un orf
$orf        = "";
 
my @var   = ();     #La usamos para guardar los orfn
my $count = 0;      #Contador de orf
my $posi  = 0;      #Posicion inicial
my $posf  = 0;      #Posicion final
my $a     = "";
 
for ( my $i = 0; $i < length $fasta; $i = $i + 3 ) {
    my $str = "";
    $a = substr( $fasta, $i, 3 );
 
    #Para que me vaya leyendo las letras de 3 en 3.
    if ( ( $a eq 'ATG' ) and ( $indica_orf == 0 ) )
    {               #Si $a es igual a ATG y el indicador esta a 0
        $indica_orf
            = 1;  #Pon el indicador a uno ----->Ahora estamos dentro de un orf
        $posi = $i + 1;
    }
    if ( $indica_orf == 1 ) {    #Si estamos dentro de un orf
        $str = modulo::codon($a)
            ;    #modulo::codon -> llamamos a la funcion codon del paquete.
        $orf = $orf
            . $str
            ;    #Sustituye las 3 letras de a por su valor y va concatenando
    }
    if ( $a =~ /TA[AG]|TGA/i )
    {            #Cuando a sea igual a TAA TAG o TGA (codones de terminacion)
        if ( $indica_orf == 1 ) {    #Y estemos dentro de un orf
            $indica_orf = 0;              #Sal del orf
            $orf        = $orf . "\n";    #Pon un salto de linea
 
            if ( ( length $orf ) * 3 >= $n )
            {    #$n = $ARGV[2] -> numero que se introduce como parametro.
                push( @var, $orf )
                    ; #Si la longitud del orf por 3 (porque un aminoacido vale tres codones) es mayor que n, guarda el orf
                $count++;    #Va contando los orf mayores q n
                $posf = $i;
 
                print SALIDA "$cabecera_numerica" . "_"
                    . "$count [$posi - $posf]$cabecera_infor"
                    ;        #Cabecera (con su contador y sus cosas)
                print SALIDA ($orf)
                    ; #Llamada a una funcion para que me imprima un salto de linea cada 60 nucleotido (imprimi bonito)que no funciona
            }
        }
        $orf = "";
    }
}
 
for ( my $i = 0; $i < length $fasta; $i = $i + 3 )
{                     #Vuelo a hacer lo mismo para la T de ATG
    my $str = "";
    $a = substr( $fasta, $i, 3 );
    if ( ( $a eq 'ATG' ) and ( $indica_orf == 0 ) ) {
        $posi       = $i + 2;
        $corte      = $i + 1;
        $indica_orf = 1;
 
        #$cadena = modulo::codon_alt($a);
        #$orf = $orf.$cadena;
    }
 
    if ( $indica_orf == 1 ) {
        $e = substr( $fasta, $corte, 3 );    #Aqui cojo la T o eso creo.
        $str = modulo::codon($e)
            ;    #Traduzco... (pero no me da Metionina al no ser ATG)
        $orf = $orf . $str;
    }
 
    if ( $e =~ /TA[AG]|TGA/i ) {
        if ( $indica_orf == 1 ) {
            $indica_orf = 0;
            $orf        = $orf . "\n";
            if ( ( length $orf ) * 3 >= $n ) {
                push( @var, $orf );
                $count++;
                $posf = $i;
                print SALIDA "$cabecera_numerica" . "_"
                    . "$count [$posi - $posf]$cabecera_infor"
                    ; #Tengo dudas con los print..si lo pongo fuera de los for no saldría nada... pero dentro de ellos el contador se pone a cero.
                print SALIDA $orf;
            }
        }
        $orf = "";
    }
}
 
Coloreado en 0.004 segundos,  usando GeSHi 1.0.8.4

por **explorer** » 2011-01-28 20:45 @906

No necesitas hacer varios bucles para detectar ORF en posiciones de la secuencia desplazadas 0, 1 o 2 bases desde el comienzo... fíjate que lo único que tienes que hacer es cambiar el '0' del inicio del bucle for() en la búsqueda de los ORF.

Tienes un ejemplo en el último código del hilo Buscar ORF.

Hay más hilos al respecto, en este foro de Bioinformática (no eres el único a quien le han encargado el mismo trabajo, a lo largo de estos años).

En cuanto a tu código, a primera vista, parece que está bien.

Cuando te decía que publicaras el ADN es si podías localizar o inventarte un trozo pequeño de ADN que, pasado por tu programa, fallara a la hora de localizar los ORF. Aquí hay una secuencia de ejemplo: post23998.html#p23998

Un detalle importante... el último ORF de la secuencia puede estar abierto (no detectar un codón de terminación al finalizar la secuencia). En algunos casos, ese ORF también forma parte de la salida correcta que se pide.

Foro - Perl en Español

Detección de ORF

Detección de ORF

Publicidad

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

Re: Detección de ORF

¿Quién está conectado?