• Publicidad

Problema con análisis de PDB

Perl aplicado a la bioinformática

Problema con análisis de PDB

Notapor Diego3D » 2012-06-22 22:26 @976

Hola a todos, espero que estén bien :)

Mi problema es el siguiente:

Tengo que medir la calidad de un conjunto de PDB. En palabras simples, tengo que obtener el porcentaje de completitud de cada uno y corroborar si la numeración de los ATOM es consistente con la numeración del SEQRES.

Concretamente la tarea solicita construir dos subrutinas, una que traduzca los aminoácidos de tres letras a una letra, y otra subrutina que realice la reconstrucción de la proteína a partir de los ATOM, para obtener la completitud y la consistencia.

Con los requerimientos no he tenido problemas, pero al momento de ejecutar me aparece el típico warning "Use of uninitialized value in string eq at", pero resulta que yo he revisado por todos lados el arreglo y no tiene ningún valor nulo.

Sinceramente ya es algo que se me escapa de las manos. Por eso acudí al foro. Aquí adjunto el código que tengo hasta el momento, por si alguien me pudiese echar una mano :D

De antemano, ¡gracias!

Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1. #!/usr/bin/perl
  2. use strict;
  3. use warnings;
  4.  
  5. my %code=(
  6.              'ALA' => 'A',
  7.              'ASX' => 'B',
  8.              'CYS' => 'C',
  9.              'ASP' => 'D',
  10.              'GLU' => 'E',
  11.              'PHE' => 'F',
  12.              'GLY' => 'G',
  13.              'HIS' => 'H',
  14.              'ILE' => 'I',
  15.              'LYS' => 'K',
  16.              'LEU' => 'L',
  17.              'MET' => 'M',
  18.              'ASN' => 'N',
  19.              'PRO' => 'P',
  20.              'GLN' => 'Q',
  21.              'ARG' => 'R',
  22.              'SER' => 'S',
  23.              'THR' => 'T',
  24.              'VAL' => 'V',
  25.              'TRP' => 'W',
  26.              'XXX' => 'X',
  27.              'TYR' => 'Y',
  28.              'GLX' => 'Z',
  29.              'SEC' => 'U',
  30.              'PYL' => 'O',
  31.              'MSE' => 'M',
  32.              'AGM' => 'R',
  33.              'MNG' => 'Q',
  34.              'GL3' => 'G',
  35.              'SMC' => 'C',
  36.              'CMT' => 'C',
  37.              'CSO' => 'C',
  38.              'TPO' => 'T',
  39.              'TPR' => 'Y',
  40.              'ACE' => 'X',
  41.              'SEP' => 'S',
  42.              'MHO' => 'M',
  43.              '...' => '.',
  44.              'END' => '*'
  45.              );
  46.  
  47. opendir(my $dh, "/home/Diego/Universidad/bio2/tarea_pdbs/kinasas") || die "Error: can't opendir";
  48. open(RES, ">res.csv");
  49. my @osqr;
  50.  
  51. #while(readdir $dh) {
  52.  
  53.         my @pdb= ();
  54.         my @tsqr;
  55.         my $file = "1A06.pdb";#$_;
  56.         open (PRO, "<$file") || die "Error"; #, "/home/Diego/Universidad/bio2/tarea_pdbs/kinasas/$file") || die "Error: can't open file $file\n";
  57.         while(my $line = <PRO>){
  58.  
  59.                 push(@pdb, $line);
  60.         }
  61.         close (PRO);
  62.  
  63.         my @seqres=grep(/^SEQRES/, @pdb);
  64.         my @atoms = grep(/^ATOM/, @pdb);
  65.  
  66.         foreach my $line (@seqres){
  67.                 my @tmp= split(/\s+/, substr($line, 19, 51));
  68.                 push(@tsqr,@tmp);
  69.         }
  70.  
  71.         @osqr=translate(@tsqr);
  72.         rebuild(@atoms);
  73. #}
  74.  
  75. close(RES);
  76.  
  77. # subrutina que traduce aminoácidos de 3 letras a 1 letra
  78. sub translate{
  79.         my @tseq= @_;
  80.         my @oseq;
  81.         my $i=0;
  82.  
  83.         foreach my $res (@tseq){
  84.                 foreach my $key (keys %code)
  85.                 {
  86.                         if($res eq $key){
  87.                                 $oseq[$i]= $code{$key};
  88.                                 $i++;
  89.                         }
  90.                 }
  91.         }
  92.         return @oseq;
  93. }
  94.  
  95. # subrutina que reconstruye la proteína con los aminoácidos de ATOM, tomando como referencia seqres
  96. sub rebuild{
  97.  
  98.         my @atoms = @_;
  99.         my @tstruct_seq=();
  100.         my @ostruct_seq=();my @rebuild=();my @pos=();
  101.         my $i=0;my $j=0;my $check=0;my $cont=0;
  102.         my $tmp1;my $tmp2;
  103.  
  104.         #obtener el aminoácido del carbono alfa y su posición, a partir de ATOM.
  105.         foreach my $atom(@atoms){
  106.                 if($atom =~m/CA/){
  107.                         my @tmp= split(/\s+/,$atom);
  108.                         $tstruct_seq[$i]=$tmp[3];
  109.                         $pos[$i]=$tmp[5];
  110.                         $i++;
  111.                 }
  112.         }
  113.        
  114.         #traducir de proteína de código de 3 letras a 1 letra.
  115.         @ostruct_seq = translate(@tstruct_seq);
  116.         $i=0;
  117.         $j=0;
  118.  
  119.         #reconstruir la proteína con los aminoácidos de ATOM, tomando como referencia seqres
  120.         while($i<@ostruct_seq){
  121.  
  122.                 while($ostruct_seq[$i] ne $osqr[$j]){
  123.                         unless(defined($rebuild[$j])){
  124.                                 $rebuild[$j]= 'X';
  125.                                 $j++;
  126.                         }
  127.                 }
  128.              
  129.                 $tmp1=$i;
  130.                 $tmp2=$j;
  131.                 $cont=0;
  132.  
  133.                 while($ostruct_seq[$tmp1] eq $osqr[$tmp2]){
  134.                           $tmp1++;
  135.                           $tmp2++;
  136.                           $cont++;
  137.                 }
  138.  
  139.                 if($cont>3){
  140.                         while($ostruct_seq[$i] eq $osqr[$j]){
  141.                                 $rebuild[$j]=$ostruct_seq[$i];
  142.                                 $i++;
  143.                                 $j++;
  144.                         }
  145.                 }
  146.                
  147.                         else{
  148.                                 $tmp1=$i;
  149.                                 $tmp2=$j;
  150.                                 while($ostruct_seq[$tmp1] eq $osqr[$tmp2]){
  151.                                         $rebuild[$j]='X';
  152.                                         $tmp1++;
  153.                                         $tmp2++;
  154.                                 }
  155.                                 $j++;
  156.                         }
  157.         }
  158.  
  159.         #si la estructura no está completa hasta el final de la proteína, se rellena con X
  160.         if($j<@osqr){
  161.                 while($j<@osqr){
  162.                         $rebuild[$j] = 'X';
  163.                         $j++;
  164.                 }
  165.         }
  166.        
  167.         #comprobar que las posiciones de ATOM coinciden con las de seqres
  168.         $j=0;
  169.         for(my $i=0; $i < @rebuild; $i++){
  170.              
  171.                   next if($rebuild[$i] eq 'X');
  172.                   my $tmp=$i+1;
  173.                   if($tmp == $pos[$j]){
  174.                         $check++;
  175.                   }
  176.                   $j++;
  177.         }
  178.         print "\nPorcentaje de completitud:  ".((@ostruct_seq / @osqr)*100)." %\n";
  179.         if( $check == @ostruct_seq){
  180.                 print "La proteína tiene consistencia en la numeración.\n\n";
  181.         }
  182. }
Coloreado en 0.006 segundos, usando GeSHi 1.0.8.4
Última edición por Diego3D el 2012-06-22 22:57 @998, editado 1 vez en total
Diego3D
Perlero nuevo
Perlero nuevo
 
Mensajes: 9
Registrado: 2011-12-08 10:58 @498

Publicidad

Re: Problema con análisis de PDB

Notapor explorer » 2012-06-23 07:15 @343

¿Podrías crear un archivo de entrada, lo suficientemente pequeño como para que 1) se reproduzca el error, y 2) puedas publicarlo aquí?

Así podremos probarlo.
JF^D Perl programming & Raku programming. Grupo en Telegram: https://t.me/Perl_ES
Avatar de Usuario
explorer
Administrador
Administrador
 
Mensajes: 14476
Registrado: 2005-07-24 18:12 @800
Ubicación: Valladolid, España

Re: Problema con análisis de PDB

Notapor explorer » 2012-06-24 07:26 @351

Hay algunas cosas que se pueden abreviar bastante...
Líneas 56 a 61:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         #$file = "/home/Diego/Universidad/bio2/tarea_pdbs/kinasas";
  2.         open    PRO, "<$file"  or  die "Error: can't open file $file: $!\n";
  3.         @pdb = <PRO>;
  4.         close   PRO;
Coloreado en 0.001 segundos, usando GeSHi 1.0.8.4

Líneas 83 a 91:
Sintáxis: [ Descargar ] [ Ocultar ]
Using perl Syntax Highlighting
  1.         for my $res (@tseq){
  2.             next if not $code{$res};
  3.             push @oseq, $code{$res};
  4.         }
Coloreado en 0.001 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: 14476
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 1 invitado

cron