Extraer coordenadas de átomos en un fichero PDB

Los ficheros PDB (Protein Data Bank) contienen las coordenadas espaciales de los átomos de proteínas cuya estructura está resuelta por técnicas de rayos X o resonancia magnética nuclear (RMN). En una entrada de blog de Bruno podéis encontrar una muy buena explicación de estas técnicas realizada por el periódico El País.

Las estructuras de proteínas en formato PDB contienen mucha información que habitualmente no nos interesa para nuestros experimentos e interfiere con muchos programas que únicamente necesitan las coordenadas espaciales de los átomos. Además un fichero PDB normalmente contiene datos de múltiples moléculas y en diferentes posiciones en el cristal. Por ello es muy conveniente extraer de los ficheros PDB únicamente la información estructural que vamos a utilizar y descartar el resto.

En el post anterior vimos como extraer líneas de un texto o archivo. En la entrada de hoy veremos como reusar ese código para extraer los átomos deseados de un fichero PDB.

En el ejemplo se indica el fichero PDB original en la variable $pdb_file, en nuestro caso será la estructura de la famosa insulina que se puede descargar aquí (2INS.pdb). Y en el array @desired_coords le indicaremos el tipo de átomos que queremos extraer en un formato muy similar al usado por PyMOL.

   
 # File with 3D coords of insulin  
 my $pdb_file = "2INS.pdb";  
   
 # Extract coords of the alpha carbons of the chain A, and alpha and beta carbons of histidines of chain B  
 my @desired_coords = ('a/*/ca/', 'b/*/ca|cb/his');  
   
 open(PDBFILE,$pdb_file);  
 my $pdb_content = join('',<PDBFILE>);  
 close PDBFILE;  
   
 my $pdb_coords = join('',extract_pdb_coords($pdb_content,\@desired_coords));  
   
 open(PDBFILE,">$pdb_file.out");  
 print PDBFILE $pdb_coords;  
 close PDBFILE;  
   
 # Extract desired PDB coordinates from PDB entries  
 sub extract_pdb_coords {  
   
      my ($pdb_content, $types) = @_;  
   
      my $pdb_data;  
      my $patterns;  
      my @pattern_parts = ('chain','res_number','res_type','res_name');  
      foreach my $type (@{$types}) {  
           my $count = 0;  
           my %pattern;  
           while ($type =~ /([^\/]+)/g){  
                $type = $'; # Text to the right of the match  
                $pattern{$pattern_parts[$count]} = $1;  
                $count++;  
           }  
           push(@{$patterns},\%pattern);  
      }  
   
      foreach my $pattern (@{$patterns}){  
           my ($chain,$res_number,$res_type,$res_name);  
           if (!defined($pattern->{'chain'}) || !$pattern->{'chain'} || $pattern->{'chain'} eq '*'){  
                $chain = '\w{1}';  
           } else {  
                $chain = uc($pattern->{'chain'});  
           }  
           $pattern->{'res_number'} =~ s/\s+//;  
           if (!defined($pattern->{'res_number'}) || !$pattern->{'res_number'} || $pattern->{'res_number'} eq '*'){  
                $res_number = '\d+';  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)$/) {  
                $res_number = $1;  
           } elsif ($pattern->{'res_number'} =~ /^(\d+)-(\d+)$/) {  
                $res_number = '('.join('|', $1 .. $2).')';  
           } elsif ($pattern->{'res_number'} =~ /\d,/) {  
                $res_number = '('.join('|', split(",",$pattern->{'res_number'})).')';  
           }  
           if (!defined($pattern->{'res_type'}) || !$pattern->{'res_type'} || $pattern->{'res_type'} eq '*'){  
                $res_type = '[\w\d]+';  
           } else {  
                $res_type = uc($pattern->{'res_type'});  
           }  
           if (!defined($pattern->{'res_name'}) || !$pattern->{'res_name'} || $pattern->{'res_name'} eq '*'){  
                $res_name = '\w{3}';  
           } else {  
                $res_name = uc($pattern->{'res_name'});  
           }  
           $pattern = '/(ATOM|HETATM)\s+\d+\s+('.$res_type.')\s+('.$res_name.')\s('.$chain.')\s+('.$res_number.')\s+.+/';  
      }  
   
      my @pdb_data_lines = extract_lines_from_text($pdb_content, $patterns);  
      if (@pdb_data_lines){  
           $pdb_data = join("\n",@pdb_data_lines);  
      }  
   
      return $pdb_data;  
   
 }  
   
 # Extract lines from text with the desired patterns  
 sub extract_lines_from_text {  
   
      my ($text, $patterns) = @_;  
   
      my @data;  
      my @lines = split("\n",$text);  
   
      foreach my $line (@lines){  
           foreach my $pattern (@{$patterns}){  
                if ($pattern =~ /^\/(.+)\/$/){  
                     if ($line =~ /$1/){  
                          push(@data,$line);  
                          last;  
                     }  
                } else {  
                     if ($line =~ /\Q$pattern\E/){  
                          push(@data,$line);  
                          last;  
                     }  
                }  
           }  
      }  
   
      return @data;  
   
 }