Foro - Perl en Español

La comunidad de programadores en Perl de habla hispana.

Obviar


  • Publicidad

Subrutina para parsear anotaciones GenBank

Perl aplicado a la bioinformática
Publicar una respuesta

Subrutina para parsear anotaciones GenBank

Notapor leire_12 » 2009-07-31 08:49 @409

¡¡Hola!!

Necesito una subrutina para parsear las anotaciones de un fichero GenBank. Tengo que pasar como argumento una variable donde se encuentran las anotaciones (que ya tengo), y lo que me tiene que devolver es un hash en el que las llaves sea los títulos de las diferentes secciones (locus, features...) y los valores los datos correspondientes a esa sección. ¿Cómo lo hago?
leire_12
Perlero nuevo
Perlero nuevo
 
Mensajes: 18
Registrado: 2009-07-28 07:12 @342
Arriba

Publicidad

Re: Subrutina para parsear anotaciones GenBank

Notapor explorer » 2009-07-31 09:18 @429

Te recuerdo que este es un foro de programación Perl, no de programación en Bioinformática. Deberías detallar todos los conceptos que expones. Piensa que nosotros no somos bioinformáticos y no sabemos de qué hablas. :)

De todas maneras, buscando por el foro, te recomiendo que visites los hilos en los que han participado los usuarios raquel y wampaier, pues ya se trató de ese tema en ellos.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Subrutina para parsear anotaciones GenBank

Notapor leire_12 » 2009-07-31 09:29 @437

Raquel os puso lo un ejemplo de fichero GenBank. Yo lo que tengo que hacer es abrir ese fichero y crear un hash que tenga como llaves: LOCUS, DEFINITION, ACCESION etc, y como valores:

LOCUS la información

DEFINITION la información

.
.
.

Tengo el código hecho (lo he cogido del libro "begining perl for bioinformatics") pero no me funciona. ¿¿Me lo podrías echar un vistazo?? Con la subrutina get_file_data() se me separa por una parte el ADN y por la otra las anotaciones. Pero creo que lo que falla es la subrutina parse_annotation() porque no me devuelve nada.

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!usr/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5.  
  6. my $fh;
  7. my $record;
  8. my $dna='';
  9. my $annotation='';
  10. my $offset;
  11. my %fields = ();
  12. my @features;
  13. my @CDS;
  14. my @codones_key;
  15. my %codones;
  16. my $codon;
  17. my %Aminoacidos;
  18. my $CDS;
  19. my $anotaciones;
  20. my $secuencia_valida;
  21. my $fichero='Escherichiacoli.gb';
  22. my $fichero2= 'Arabidopsisthaliana.gb';
  23.  
  24.  
  25. $fh = open_file($fichero2);
  26.  
  27. $record = get_next_record($fh);
  28.  
  29. ($annotation, $dna) = get_annotation_and_dna($record);
  30.  
  31.  
  32.  
  33.  
  34. %fields = parse_annotation ($annotation);
  35.  
  36.  
  37. foreach my $key (keys %fields) {
  38.         print "$key \n";
  39.         print $fields{$key};
  40. }
  41.  
  42. ##################################################################
  43.  
  44. sub open_file{
  45.        
  46.         my ($fichero)=@_;
  47.         my $fh;
  48.         unless(open($fh,$fichero))
  49.         {
  50.                 print "Error a la hora de abrir el fichero";
  51.                 exit;
  52.         }
  53.        
  54.         return $fh;
  55.        
  56. }
  57.  
  58. sub get_next_record {
  59.         #given the filehandle, get the record
  60.      #(we can get the offset by first calling "tell")
  61.     my($fh) = @_;
  62.  
  63.     my($offset);
  64.     my($record) = '';
  65.     my($save_input_separator) = $/;
  66.     $/ = "//\n";
  67.     $record = <$fh>;
  68.     $/ = $save_input_separator;
  69.    
  70.     return $record;
  71. }    
  72. sub get_annotation_and_dna
  73. {
  74.         my($record)=@_;
  75.        
  76.         my($annotation)= '';
  77.         my($dna)= '';
  78.         my($dnaAux) = '';
  79.         my($annotation_aux)='';
  80.         ($annotation, $dna) = ($record=~ /^(LOCUS.*ORIGIN\s*\n)(.*)\/\/\n/s);
  81.         $annotation_aux = $annotation;
  82.         $dnaAux=$dna;
  83.         #$dna=~ s/[\s\/]//g;
  84.        
  85.         return ($annotation_aux, $dnaAux);
  86.        
  87. }
  88.  
  89. sub parse_annotation {
  90.         my($annotation)= @_;
  91.         my (%results) = ();
  92.        
  93.         while ($annotation =~ /^[A-Z].*\n(^\s.*\n)*/gm) {
  94.                 my $value = $&;
  95.                 (my $key = $value) =~ s/^([A-Z]+).*/$1/s;
  96.                 $results {$key} = $value;}
  97.        
  98.        
  99.        
  100.         return %results;
  101. }
Coloreado en 0.003 segundos, usando GeSHi 1.0.8.4
Última edición por explorer el 2009-07-31 09:49 @450, editado 1 vez en total
Razón: Ortografía, bloques de código
leire_12
Perlero nuevo
Perlero nuevo
 
Mensajes: 18
Registrado: 2009-07-28 07:12 @342
Arriba

Re: Subrutina para parsear anotaciones GenBank

Notapor explorer » 2009-07-31 10:27 @477

Con este fichero de entrada
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
LOCUS       X72459                   375 bp    mRNA    linear   PRI 14-NOV-2006
DEFINITION  H.sapiens mRNA for rearranged Ig kappa light chain variable region
            (I.169).
ACCESSION   X72459
VERSION     X72459.1  GI:441386
KEYWORDS    immunoglobulin; J-segment; kappa light chain; V-region.
SOURCE      Homo sapiens (human)
  ORGANISM  Homo sapiens
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Hominidae; Homo.
REFERENCE   1  (bases 1 to 375)
  AUTHORS   Klein,R., Jaenichen,R. and Zachau,H.G.
  TITLE     Expressed human immunoglobulin kappa genes and their hypermutation
  JOURNAL   Eur. J. Immunol. 23 (12), 3248-3262 (1993)
   PUBMED   8258341
REFERENCE   2  (bases 1 to 375)
  AUTHORS   Zachau,H.G.
  TITLE     Direct Submission
  JOURNAL   Submitted (26-APR-1993) H.G. Zachau, Institut fuer Physiologische
            Chemie, der Universitaet Muenchen, Schillerstr 44, 8000 Muenchen 2,
            FRG
FEATURES             Location/Qualifiers
     source          1..375
                     /organism="Homo sapiens"
                     /mol_type="mRNA"
                     /isolate="M.L."
                     /db_xref="taxon:9606"
                     /chromosome="2"
                     /clone="I.169"
                     /tissue_type="spleen"
                     /clone_lib="lambda zap II phage library"
     CDS             1..375
                     /note="Protein sequence is in conflict with the conceptual
                     translation"
                     /codon_start=1
                     /product="Ig kappa light chain (VJ)"
                     /protein_id="CAA51127.1"
                     /db_xref="GI:441387"
                     /translation="PAQLLGLLLLWLPGARCAIQLTQSPSSLSASVGDRVTITCRASQ
                     GISSALAWYQQKPGKAPKLLIYDASSLESGVPSRFSGSGSGTDFTLTISSLQPEDFAT
                     YYCQQFNTYPLTFGGGTKVEIKR"
     V_region        1..375
                     /product="Ig kappa light chain (VJ)"
     V_segment       53..336
     J_segment       337..375
                     /note="J-Kappa 4"
ORIGIN    
        1 cccgctcagc tcctggggct tctgctgctc tggctcccag gtgccagatg tgccatccag
       61 ttgacccagt ctccatcctc cctgtctgca tctgtaggag acagagtcac catcacttgc
      121 cgggcaagtc agggcataag cagtgcttta gcctggtatc agcagaaacc agggaaagct
      181 cctaagctcc tgatctatga tgcctccagt ttggaaagtg gggtcccatc aaggttcagc
      241 ggcagtggat ctgggacaga tttcactctc accatcagca gcctgcagcc tgaagatttt
      301 gcaacttatt actgtcaaca gtttaatact tacccgctca ctttcggcgg agggaccaag
      361 gtggagatca aacga
 
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

Y este programa:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4. use diagnostics;
  5.  
  6. ## Leemos todo el fichero de golpe
  7. open  ADN, '<adn.txt' or die;
  8. my $adn = join '', <ADN>;
  9. close ADN;
  10.  
  11. my %genbank;
  12.  
  13. ## Sacamos las secciones
  14. while ($adn =~ /^([A-Z]+)\s*(.*\n(?:^\s.*\n)*)/gm ) {
  15.     my ($clave,$valor) = ($1,$2);
  16.     chomp $valor;
  17.     $genbank{$clave} = $valor;
  18. }
  19.  
  20. ## Aspecto del diccionario
  21. use Data::Dumper;
  22. print Dumper \%genbank;
  23.  
  24. __END__
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Sale:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
$VAR1 = {
          'ACCESSION' => 'X72459',
          'DEFINITION' => 'H.sapiens mRNA for rearranged Ig kappa light chain variable region
            (I.169).',
          'ORIGIN' => '1 cccgctcagc tcctggggct tctgctgctc tggctcccag gtgccagatg tgccatccag
       61 ttgacccagt ctccatcctc cctgtctgca tctgtaggag acagagtcac catcacttgc
      121 cgggcaagtc agggcataag cagtgcttta gcctggtatc agcagaaacc agggaaagct
      181 cctaagctcc tgatctatga tgcctccagt ttggaaagtg gggtcccatc aaggttcagc
      241 ggcagtggat ctgggacaga tttcactctc accatcagca gcctgcagcc tgaagatttt
      301 gcaacttatt actgtcaaca gtttaatact tacccgctca ctttcggcgg agggaccaag
      361 gtggagatca aacga',
          'LOCUS' => 'X72459                   375 bp    mRNA    linear   PRI 14-NOV-2006',
          'REFERENCE' => '2  (bases 1 to 375)
  AUTHORS   Zachau,H.G.
  TITLE     Direct Submission
  JOURNAL   Submitted (26-APR-1993) H.G. Zachau, Institut fuer Physiologische
            Chemie, der Universitaet Muenchen, Schillerstr 44, 8000 Muenchen 2,
            FRG',
          'FEATURES' => 'Location/Qualifiers
     source          1..375
                     /organism="Homo sapiens"
                     /mol_type="mRNA"
                     /isolate="M.L."
                     /db_xref="taxon:9606"
                     /chromosome="2"
                     /clone="I.169"
                     /tissue_type="spleen"
                     /clone_lib="lambda zap II phage library"
     CDS             1..375
                     /note="Protein sequence is in conflict with the conceptual
                     translation"
                     /codon_start=1
                     /product="Ig kappa light chain (VJ)"
                     /protein_id="CAA51127.1"
                     /db_xref="GI:441387"
                     /translation="PAQLLGLLLLWLPGARCAIQLTQSPSSLSASVGDRVTITCRASQ
                     GISSALAWYQQKPGKAPKLLIYDASSLESGVPSRFSGSGSGTDFTLTISSLQPEDFAT
                     YYCQQFNTYPLTFGGGTKVEIKR"
     V_region        1..375
                     /product="Ig kappa light chain (VJ)"
     V_segment       53..336
     J_segment       337..375
                     /note="J-Kappa 4"',
          'KEYWORDS' => 'immunoglobulin; J-segment; kappa light chain; V-region.',
          'VERSION' => 'X72459.1  GI:441386',
          'SOURCE' => 'Homo sapiens (human)
  ORGANISM  Homo sapiens
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Hominidae; Homo.'
        };
Coloreado en 0.000 segundos, usando GeSHi 1.0.8.4

La expresión regular busca secciones compuestas de: palabras que comienzan en la primera columna del fichero, seguida por cero o más líneas que comienzan por caracteres espacio.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Subrutina para parsear anotaciones GenBank

Notapor leire_12 » 2009-07-31 13:52 @619

El problema es que tengo que utilizar el código que he escrito yo arriba porque es el que hemos dado en clase. Yo no le veo ningún error y no sé por qué no funciona.
leire_12
Perlero nuevo
Perlero nuevo
 
Mensajes: 18
Registrado: 2009-07-28 07:12 @342
Arriba

Re: Subrutina para parsear anotaciones GenBank

Notapor explorer » 2009-07-31 14:43 @654

He intentado correr tu programa con el fichero de ejemplo que he puesto, y he descubierto que el problema está en la expresión regular de la línea 80: exige que el fichero GenBank termine en una línea que comience por '//'.

Tienes dos opciones: o te aseguras que los ficheros leídos tienen esa marca al final, y se la pones si no la tiene, o modificas la expresión regular de la línea 80.

Esa expresión está así porque, como comenta el usuario wampaier en el hilo 'Contar ocurrencias', esos dos caracteres se usan de separador, dentro de un mismo fichero, para diferenciar distintos registros genbank.

Si vas a tratar ficheros genbank compuestos de un solo registro, entonces creo que se podría cambiar a esto: /^(LOCUS.*ORIGIN\s*\n)(.*)\n/s Es decir, quitar el separador '//'.

Otro detalle, en la línea 1:
#!usr/bin/perl
Si estás en Windows, no es problema, pero si estás en Unix/Linux, sí.
Lo correcto es:
#!/usr/bin/perl

Ahora ya me sale:
Sintáxis: [ Descargar ] [ Ocultar ]
Using text Syntax Highlighting
ORIGIN
ORIGIN
ACCESSION
ACCESSION   X72459
LOCUS
LOCUS       X72459                   375 bp    mRNA    linear   PRI 14-NOV-2006
FEATURES
FEATURES             Location/Qualifiers
     source          1..375
                     /organism="Homo sapiens"
                     /mol_type="mRNA"
                     /isolate="M.L."
                     /db_xref="taxon:9606"
                     /chromosome="2"
                     /clone="I.169"
                     /tissue_type="spleen"
                     /clone_lib="lambda zap II phage library"
     CDS             1..375
                     /note="Protein sequence is in conflict with the conceptual
                     translation"
                     /codon_start=1
                     /product="Ig kappa light chain (VJ)"
                     /protein_id="CAA51127.1"
                     /db_xref="GI:441387"
                     /translation="PAQLLGLLLLWLPGARCAIQLTQSPSSLSASVGDRVTITCRASQ
                     GISSALAWYQQKPGKAPKLLIYDASSLESGVPSRFSGSGSGTDFTLTISSLQPEDFAT
                     YYCQQFNTYPLTFGGGTKVEIKR"
     V_region        1..375
                     /product="Ig kappa light chain (VJ)"
     V_segment       53..336
     J_segment       337..375
                     /note="J-Kappa 4"
REFERENCE
REFERENCE   2  (bases 1 to 375)
  AUTHORS   Zachau,H.G.
  TITLE     Direct Submission
  JOURNAL   Submitted (26-APR-1993) H.G. Zachau, Institut fuer Physiologische
            Chemie, der Universitaet Muenchen, Schillerstr 44, 8000 Muenchen 2,
            FRG
KEYWORDS
KEYWORDS    immunoglobulin; J-segment; kappa light chain; V-region.
VERSION
VERSION     X72459.1  GI:441386
SOURCE
SOURCE      Homo sapiens (human)
  ORGANISM  Homo sapiens
            Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;
            Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;
            Catarrhini; Hominidae; Homo.
DEFINITION
DEFINITION  H.sapiens mRNA for rearranged Ig kappa light chain variable region
            (I.169).
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: 14486
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España
Arriba

Re: Subrutina para parsear anotaciones GenBank

Notapor leire_12 » 2009-08-13 05:01 @251

¡Muchas gracias! Ya está resuelto.
leire_12
Perlero nuevo
Perlero nuevo
 
Mensajes: 18
Registrado: 2009-07-28 07:12 @342
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