Página 1 de 1

Problema con los frames y los ORF

NotaPublicado: 2012-04-14 01:14 @093
por Diego3D
Hola a todos :D

¡Espero que estén bien!

Resulta que estoy programando un script para anotar un genoma. Para eso busco los ORF en cada marco, los traduzco y los guardo en un archivo FASTA.

El problema es que el marco 1 y 2 me dan las mismas proteínas, y no entiendo por qué, ya que el marco 0 me funciona sin problemas. Lo mismo me pasa con los reversos complementarios.

Aquí está mi código:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4. use Bio::Seq;
  5. use Bio::SeqIO;
  6.  
  7.  
  8. my $obj_fasta = Bio::SeqIO->new(-file => "NC_008512.ffa",
  9.                                 -format => 'fasta');
  10. my $ob = $obj_fasta->next_seq;
  11.  
  12. my $seq_obj= Bio::Seq->new(-display_id  =>      $ob->id,
  13.                             -seq        =>      $ob->seq,
  14.                             -alphabet   =>      'dna'
  15.                       );
  16.  
  17. my $j=0;
  18. my $frame;
  19.  
  20. for(my $i=0; $i<6;$i++){ # Busca proteínas en los 6 marcos de lectura
  21.  
  22.         if($i==3){
  23.                 $j=0;
  24.         }
  25.  
  26.         if($i<=2){ # obtiene los tres primeros marcos directos
  27.                 $frame = substr($seq_obj->seq, $j,$seq_obj->length);
  28.                 open(MIFICH,">frame_$j.txt"); # se crea un archivo para cada frame
  29.         }
  30.  
  31.         else{ # obtiene los tres marcos reversos complementarios
  32.                 $frame = substr($seq_obj->revcom->seq, $j,$seq_obj->length);
  33.                 open(MIFICH,">frame_-$j.txt"); # se crea un archivo para cada frame
  34.         }
  35.  
  36.         $j++;
  37.         my @orfs = $frame =~ /(ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/g; # hace match con todos los ORF que encuentre
  38.        
  39.         foreach my $orf (@orfs) {
  40.                 if((length($orf))>=240){ # proteínas con largo >= a 80 aminoácidos
  41.  
  42.                         my $prot_ob= Bio::Seq->new( # se crea un objeto secuencia con el ORF
  43.                                 -seq            =>      $orf,
  44.                                 -alphabet       =>      'dna'
  45.                         );
  46.  
  47.                         if($i<=2){ # posición dentro del genoma del ORF directo
  48.                                 my $start = index($seq_obj->seq, $prot_ob->seq);
  49.                                 my $end = $start + length($prot_ob->seq);
  50.                                 print MIFICH "$start   $end \n";
  51.                         }
  52.  
  53.                         else{ # posición dentro del genoma del ORF reverso
  54.                                 my $start = index($seq_obj->seq, $prot_ob->revcom->seq);
  55.                                 my $end = $start + length($prot_ob->seq);
  56.                                 print MIFICH "$end   $start \n";
  57.                                
  58.                         }
  59.                         print MIFICH $prot_ob->seq->translate(-condotable=>11)->seq."\n\n"; # imprime la proteína en el archivo fasta generado para su respectivo frame
  60.                 }
  61.         }
  62. }
Coloreado en 0.004 segundos, usando GeSHi 1.0.8.4

Re: Problema con los frames y los ORF

NotaPublicado: 2012-04-14 10:15 @469
por explorer
Te falta un
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         close MIFICH;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Pero creo que el problema es otro...

Supongamos que tenemos esta secuencia (tomada de la página Marco abierto de lectura, de la Wikipedia): aactgcagtacgtaacgtca

Generamos las 6 combinaciones de secuencias posibles:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
  aactgcagtacgtaacgtca
 aactgcagtacgtaacgtca
aactgcagtacgtaacgtca
ttgacgtcatgcattgcagt
 ttgacgtcatgcattgcagt
  ttgacgtcatgcattgcagt
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

Estas secuencias, si los dividimos en codones, darían lugar a secuencias proteicas diferentes, ya que los codones serían distintos (cada secuencia es diferente, por estar desplazada una posición, luego se generan codones distintos).

Bueno... a cada uno de estos marcos le aplicas la expresión regular: /(ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/gi .

Y aquí está el problema: el patrón que tienes NO tiene en cuenta la posición desplazada de cada secuencia, sino que busca, de forma lineal, dentro de la secuencia, por lo que, es muy, muy posible, que estés encontrando los mismos marcos de lectura abierta.

Ejemplo: tenemos esta secuencia, con un ORF dentro de ella:
aactgcatgacgtaacgtca

Ese ORF es encontrado sin problemas por la expresión regular:
Sintáxis: [ Descargar ] [ Ocultar ]
Using bash Syntax Highlighting
  1. > perl -E '$seq = "aactgcatgacgtaacgtca"; say $1 while $seq =~ /(ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/gi;'
  2. atgacgtaa
Coloreado en 0.002 segundos, usando GeSHi 1.0.8.4


Ahora, desplazamos la secuencia una posición más a la derecha, y vemos qué ocurre:
_aactgcatgacgtaacgtca
Sintáxis: [ Descargar ] [ Ocultar ]
Using bash Syntax Highlighting
  1. > perl -E '$seq = " aactgcatgacgtaacgtca"; say $1 while $seq =~ /(ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/gi;'
  2. atgacgtaa
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4


Sale lo mismo... porque a la expresión regular le ha dado lo mismo encontrar la secuencia en la posición seis que en la siete.

Ese es el problema: la expresión regular sabe de codones y sus posiciones múltiples -el subpatrón (?:.{3})*?-, pero no tiene en cuenta la posición inicial del codón de arranque 'ATG'.

Entonces... si el problema es encontrar los ORF en secuencias desplazadas (y sus inversas), lo que hay que hacer es, después del desplazamiento (y reversión) de las secuencias, pasar la secuencia a codones, y a partir de ellos, localizar los ORF.

Aquí tienes un ejemplo sencillo de cómo se podría resolver (sin BioPerl, que ahora mismo no lo tengo instalado):
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use v5.14;
  3.  
  4. my $seq = "aactgcatgacgtaacgtccatgacgtagaactcgatcgatcgatcgactgactagcttatagcattacacg";
  5.  
  6. my $regex_orf_mala  = qr/            (ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/xi;
  7. my $regex_orf_buena = qr/\G(?:.{3})*?(ATG(?:.{3})*?(?:TAA|TAG|TGA|.{1,3}$))/i;
  8.  
  9. (my $seqrev = $seq) =~ tr/atcg/tagc/;
  10. $seqrev = reverse $seqrev;
  11.  
  12. say "Directa (mala):";
  13. find_orf($seq, $regex_orf_mala);
  14.  
  15. say "Inversa (mala):";
  16. find_orf($seqrev, $regex_orf_mala);
  17.  
  18. say;
  19.  
  20. say "Directa (buena):";
  21. find_orf($seq, $regex_orf_buena);
  22.  
  23. say "Inversa (buena):";
  24. find_orf($seqrev, $regex_orf_buena);
  25.  
  26. sub find_orf {
  27.     my($seq, $regex) = @_;
  28.    
  29.     #say "Analizando [$seq] con [$regex]";
  30.  
  31.     for my $d (0..2) {
  32.         my $s = ' ' x $d . $seq;
  33.         say "$d : [$1]" while $s =~ /$regex/g;
  34.     }
  35. }
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4
La salida es:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
Directa (mala):
0 : [atgacgtaa]
0 : [atgacgtag]
1 : [atgacgtaa]
1 : [atgacgtag]
2 : [atgacgtaa]
2 : [atgacgtag]
Inversa (mala):
0 : [atgctataa]
0 : [atggacgttacgtcatgcagtt]
1 : [atgctataa]
1 : [atggacgttacgtcatgcagtt]
2 : [atgctataa]
2 : [atggacgttacgtcatgcagtt]

Directa (buena):
0 : [atgacgtaa]
1 : [atgacgtag]
Inversa (buena):
0 : [atgctataa]
1 : [atggacgttacgtcatgcagtt]
2 : [atgcagtt]
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4
Observa cómo la primera regex localiza la misma secuencia, tres veces.

Un detalle importante, también: esta solución no tiene en cuenta el caso de que quieras buscar los ORF que están incluidos unos dentro de otros (el caso de encontrarse varios codones de arranque que terminan en el mismo codón de terminación). Si ese fuera también el caso... habría que hacer un ajuste más fino para la búsqueda de codones.

Re: Problema con los frames y los ORF

NotaPublicado: 2012-04-14 12:06 @546
por Diego3D
Muchas gracias, explorer, por la gran respuesta que me diste :D. Sospechaba que la expresión regular no discriminaba la distancia desde el inicio, cómo tu dices.

Lo que no entendí del código eso si era, ¿por qué hay una expresión regular buena y otra mala?

Re: Problema con los frames y los ORF

NotaPublicado: 2012-04-14 12:23 @557
por explorer
Es una forma de llamarlas, para distinguir entre el caso de un patrón que no tiene en cuenta la posición múltiplo de comienzo de los codones, de otra que sí.