• Publicidad

Detección de ORF

Perl aplicado a la bioinformática

Detección de ORF

Notapor 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
  1. $indica_orf = 0;                #Indicada cuando encuentra un orf
  2. $orf = "";
  3.  
  4.  
  5. my @var = ();                   #La usamos para guardar los orfn
  6. my $count = 0;                  #Contador de orf
  7. my $posi = 0;                   #Posicion inicial
  8. my $posf = 0;                   #Posicion final
  9.  
  10.  
  11. for (my $i = 0; $i < length $fasta; $i=$i+3) {
  12.   my $str ="";
  13.  
  14.   $a = substr ($fasta, $i, 3);
  15.  #Para que me vaya leyendo las letras de 3 en 3.
  16.  
  17.   if (($a eq 'ATG') and ($indica_orf == 0)) {           #Si $a es igual a ATG y el indicador esta a 0
  18.     $indica_orf = 1;                                    #Pon el indicador a uno ----->Ahora estamos dentro de un orf
  19.     $posi = $i+1;
  20.    
  21.    
  22.   }
  23.   if ($indica_orf == 1) {                               #Si estamos dentro de un orf
  24.     $str = modulo::codon($a);                           #modulo::codon -> llamamos a la funcion codon del paquete.
  25.     $orf = $orf.$str;                                   #Sustituye las 3 letras de a por su valor y va concatenando
  26.    
  27.  
  28.   }
  29.   if ( $a =~ /TA[AG]|TGA/i ) {                          #Cuando a sea igual a TAA TAG o TGA (codones de terminacion)
  30.                                                         #Y estemos dentro de un orf
  31.     if ($indica_orf == 1) {
  32.       $indica_orf = 0;                                  #Sal del orf
  33.       $orf=$orf."\n";                                   #Pon un salto de linea
  34.      
  35.      
  36.      
  37.       if ((length $orf)*3 >= $n) {              #$n = $ARGV[2] -> numero que se introduce como parametro.
  38.         push (@var, $orf);                      #Si la longitud del orf por 3 (porque un aminoacido vale tres codones) es mayor que n, guarda el orf
  39.         $count++;                               #Va contando los orf mayores q n
  40.         $posf = $i;
  41.        
  42.         print SALIDA "$cabecera_numerica"."_"."$count [$posi - $posf]$cabecera_infor";  #Cabecera (con su contador y sus cosas)
  43.         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
  44.       }
  45.     }  
  46.    
  47.     $orf = "";
  48.   }
  49.  
  50. }
  51.  
  52. #############REVERSE##############
  53.  
  54. #Exactamente lo mismo, pero al revés (mismo problema)
  55.  
  56. my $indica_reverse = 0;
  57. my $inicio = 0;
  58. my $fin =0;
  59. my $str_rvs = "";
  60. my $orf_rvs = "";
  61.  
  62.  
  63. $fasta_reverse = reverse $fasta;
  64.  
  65. for (my $i = 0; $i < length $fasta_reverse; $i = $i+3) {
  66.   my $a = substr ($fasta_reverse, $i, 3);
  67.  
  68.   if (($a eq 'ATG')and($indica_reverse == 0)) {
  69.     $indica_reverse = 1;
  70.     $inicio = $i +1;
  71.   }
  72.   if ($indica_reverse ==1) {
  73.     $str_rvs = modulo::codon($a);
  74.     $orf_rvs = $orf_rvs.$str_rvs;
  75.   }
  76.   if ($a =~ /TA[AG]|TGA/) {
  77.     if ($indica_reverse ==1) {
  78.       $indica_reverse = 0;                     
  79.       $orf_rvs = $orf_rvs."\n";
  80.         if ((length $orf_rvs)*3 > $n) {
  81.          push ($orf_rvs, @var);
  82.           $fin = $i;
  83.           $count++;
  84.      
  85.           print SALIDA "$cabecera_numerica"."_"."$count [$inicio - $fin]"."(REVERSE SENSE)"."$cabecera";
  86.           print SALIDA linea::modulo($orf_rvs);
  87.         }
  88.       }
  89.       $orf_rvs = "";
  90.     }
  91.   }
  92.  
  93. close FASTA;
  94. close SALIDA;
Coloreado en 0.008 segundos, usando GeSHi 1.0.8.4
Rodrigobiotec
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-01-28 16:11 @716

Publicidad

Re: Detección de ORF

Notapor pvaldes » 2011-01-28 17:35 @774

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.
Última edición por pvaldes el 2011-01-28 18:41 @820, editado 1 vez en total
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Detección de ORF

Notapor 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
Última edición por pvaldes el 2011-01-28 18:40 @819, editado 1 vez en total
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Detección de ORF

Notapor 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.
Rodrigobiotec
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-01-28 16:11 @716

Re: Detección de ORF

Notapor 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
Última edición por pvaldes el 2011-01-28 18:53 @828, editado 1 vez en total
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Detección de ORF

Notapor 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.
Rodrigobiotec
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-01-28 16:11 @716

Re: Detección de ORF

Notapor 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.
pvaldes
Perlero nuevo
Perlero nuevo
 
Mensajes: 129
Registrado: 2011-01-22 12:56 @580

Re: Detección de ORF

Notapor 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?
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14475
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Detección de ORF

Notapor 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
  1. #!/usr/bin/perl -w
  2.  
  3. use modulo;
  4.  
  5. #$ARGV[0];
  6. #$ARGV[1];
  7. #$ARGV[2];
  8. chomp $ARGV[2];
  9.  
  10. $NumARGV = $#ARGV;    #Ultima posicion de @ARGV (ultimo parametro))
  11.  
  12. if ( $NumARGV != 2 )
  13. {    #Si el numero de argumentos es distinto de 2 (son 3 pero el primero es 0)
  14.     print "ERROR !!!\nDebe intoducir tres parametros\n\n";
  15.     exit;
  16. }
  17.  
  18. if ( $ARGV[2] =~ /[^0-9]/ ) {    #Valor minimo del orf
  19.     print "El tercer parametro debe ser un numero!!\n";
  20.     exit;
  21. }
  22.  
  23. my ($fasta) = $ARGV[0];  #En el primer argumento pasamos el ficherp de entrada
  24.  
  25. open( FASTA, $fasta ) or die("No puedo abrir fichero $fasta: $! \n");
  26. open( SALIDA, ">$ARGV[1]" )
  27.     or die("No puedo abrir fichero de salida\n")
  28.     ;                    #El segundo argumento es un fichero de salida
  29.  
  30. @FASTA = <FASTA>;        #Convierto el fichero de entrada em un array
  31.  
  32. $cabecera = shift @FASTA;    #El primer elemnto del array va a ser la cabecera
  33.  
  34. $cabecera_numerica = modulo::obtener_numerico($cabecera)
  35.     ;    #Llamamos a la funcion definida en el otro paquete
  36. $cabecera_infor = modulo::obtener_infor($cabecera)
  37.     ;    #Llamamos a la funcion definida en el otro paquete
  38.  
  39. $fasta = join( "", @FASTA );    #Uno todos los elementos restantes del array
  40. $fasta =~ s/\s//g;              #Quito los espacios, saltos....
  41.  
  42. my ($n)
  43.     = $ARGV[2]
  44.     ; #$n es el parametro numero 3, Ejemplo: si $n = 300 todos los orf deben tener una longitud (en nucleotidos) mayor q 300
  45.  
  46. $indica_orf = 0;    #Indicada cuando encuentra un orf
  47. $orf        = "";
  48.  
  49. my @var   = ();     #La usamos para guardar los orfn
  50. my $count = 0;      #Contador de orf
  51. my $posi  = 0;      #Posicion inicial
  52. my $posf  = 0;      #Posicion final
  53. my $a     = "";
  54.  
  55. for ( my $i = 0; $i < length $fasta; $i = $i + 3 ) {
  56.     my $str = "";
  57.     $a = substr( $fasta, $i, 3 );
  58.  
  59.     #Para que me vaya leyendo las letras de 3 en 3.
  60.     if ( ( $a eq 'ATG' ) and ( $indica_orf == 0 ) )
  61.     {               #Si $a es igual a ATG y el indicador esta a 0
  62.         $indica_orf
  63.             = 1;  #Pon el indicador a uno ----->Ahora estamos dentro de un orf
  64.         $posi = $i + 1;
  65.     }
  66.     if ( $indica_orf == 1 ) {    #Si estamos dentro de un orf
  67.         $str = modulo::codon($a)
  68.             ;    #modulo::codon -> llamamos a la funcion codon del paquete.
  69.         $orf = $orf
  70.             . $str
  71.             ;    #Sustituye las 3 letras de a por su valor y va concatenando
  72.     }
  73.     if ( $a =~ /TA[AG]|TGA/i )
  74.     {            #Cuando a sea igual a TAA TAG o TGA (codones de terminacion)
  75.         if ( $indica_orf == 1 ) {    #Y estemos dentro de un orf
  76.             $indica_orf = 0;              #Sal del orf
  77.             $orf        = $orf . "\n";    #Pon un salto de linea
  78.  
  79.             if ( ( length $orf ) * 3 >= $n )
  80.             {    #$n = $ARGV[2] -> numero que se introduce como parametro.
  81.                 push( @var, $orf )
  82.                     ; #Si la longitud del orf por 3 (porque un aminoacido vale tres codones) es mayor que n, guarda el orf
  83.                 $count++;    #Va contando los orf mayores q n
  84.                 $posf = $i;
  85.  
  86.                 print SALIDA "$cabecera_numerica" . "_"
  87.                     . "$count [$posi - $posf]$cabecera_infor"
  88.                     ;        #Cabecera (con su contador y sus cosas)
  89.                 print SALIDA ($orf)
  90.                     ; #Llamada a una funcion para que me imprima un salto de linea cada 60 nucleotido (imprimi bonito)que no funciona
  91.             }
  92.         }
  93.         $orf = "";
  94.     }
  95. }
  96.  
  97. for ( my $i = 0; $i < length $fasta; $i = $i + 3 )
  98. {                     #Vuelo a hacer lo mismo para la T de ATG
  99.     my $str = "";
  100.     $a = substr( $fasta, $i, 3 );
  101.     if ( ( $a eq 'ATG' ) and ( $indica_orf == 0 ) ) {
  102.         $posi       = $i + 2;
  103.         $corte      = $i + 1;
  104.         $indica_orf = 1;
  105.  
  106.         #$cadena = modulo::codon_alt($a);
  107.         #$orf = $orf.$cadena;
  108.     }
  109.  
  110.     if ( $indica_orf == 1 ) {
  111.         $e = substr( $fasta, $corte, 3 );    #Aqui cojo la T o eso creo.
  112.         $str = modulo::codon($e)
  113.             ;    #Traduzco... (pero no me da Metionina al no ser ATG)
  114.         $orf = $orf . $str;
  115.     }
  116.  
  117.     if ( $e =~ /TA[AG]|TGA/i ) {
  118.         if ( $indica_orf == 1 ) {
  119.             $indica_orf = 0;
  120.             $orf        = $orf . "\n";
  121.             if ( ( length $orf ) * 3 >= $n ) {
  122.                 push( @var, $orf );
  123.                 $count++;
  124.                 $posf = $i;
  125.                 print SALIDA "$cabecera_numerica" . "_"
  126.                     . "$count [$posi - $posf]$cabecera_infor"
  127.                     ; #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.
  128.                 print SALIDA $orf;
  129.             }
  130.         }
  131.         $orf = "";
  132.     }
  133. }
  134.  
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2011-01-28 20:38 @901, editado 1 vez en total
Razón: Formateo del código con la ayuda de perltidy
Rodrigobiotec
Perlero nuevo
Perlero nuevo
 
Mensajes: 4
Registrado: 2011-01-28 16:11 @716

Re: Detección de ORF

Notapor 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.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14475
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España


Volver a Bioinformática

¿Quién está conectado?

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

cron