How can I calculate the geometric center of a protein in Perl?-CodePudding

I have a PDB file which contains information about a specific protein. One of the information it holds is the positions of the different atoms (xyz coordinates).

The file is the following https://files.rcsb.org/view/6U9D.pdb . With this file I want to calculate the geometric center of the atoms. In theory I know what I need to do, but the script I wrote does not seem to work.

The first part of the program, before the for loop, is another part of the assignment which requires me to read the sequence and convert it from the 3 letter nomenclature to the 1 letter one. The part that interests me is the for loop until the end. I tried to pattern match in order to isolate the XYZ coordinates. Then I used a counter that I had set up in the beginning which is the $k variable. When I check the output on cygwin the only output I get is the sequence 0 0 0 instead of the sum of each dimension divided by $k. Any clue what has gone wrong?

$k=0;   
open (IN, '6U9D.pdb.txt');
%amino_acid_conversion = (
   ALA=>'A',TYR=>'Y',MET=>'M',LEU=>'L',CYS=>'C',
   GLY=>'G',ARG=>'R',ASN=>'N',ASP=>'D',GLN=>'Q',
   GLU=>'E',HIS=>'H',TRP=>'W',LYS=>'K',PHE=>'F',
   PRO=>'P',SER=>'S',THR=>'T',ILE=>'I',VAL=>'V'
);
while (<IN>) {
   if ($_=~m/HEADER\s (.*)/){
   print ">$1\n"; }
   if ($_=~m/^SEQRES\s \d \s \w \s \d \s (.*)/){
       $seq.=$1;
       $seq=~s/ //g;
   }
}

for ($i=0;$i<=length $seq; $i =3) {
   print "$amino_acid_conversion{substr($seq,$i,3)}";      
   if ($_=~m/^ATOM\s \d \s \w \s \w \s \w \s \d \s (\S )\s (\S )\s (\S )/) {
       $x =$1; $y =$2; $z =$3; $k  ;
   }
}
print "\n";
#print $k;
$xgk=($x/$k); $ygk=($y/$k); $zgk=($z/$k);
print "$xgk $ygk $zgk \n";

CodePudding user response：

I do not know bioinformatics but it seems like you should do something like this:

use feature qw(say);
use strict;
use warnings;

my $fn = '6U9D.pdb';
open ( my $IN, '<', $fn ) or die "Could not open file '$fn': $!";
my $seq = '';
my $x = 0;
my $y = 0;
my $z = 0;
my $k = 0;
while (<$IN>) {
   if ($_ =~ m/HEADER\s (.*)/) {
       say ">$1";
   }
   if ($_=~m/^SEQRES\s \d \s \w \s \d \s (.*)/){
       $seq .= $1;
   }
   if ($_ =~ m/^ATOM\s \d \s \w \s \w \s \w \s \d \s (\S )\s (\S )\s (\S )/) {
       $x =$1; $y =$2; $z =$3; $k  ;
   }
}
close $IN;
$seq =~ s/ //g;
my %amino_acid_conversion = (
   ALA=>'A',TYR=>'Y',MET=>'M',LEU=>'L',CYS=>'C',
   GLY=>'G',ARG=>'R',ASN=>'N',ASP=>'D',GLN=>'Q',
   GLU=>'E',HIS=>'H',TRP=>'W',LYS=>'K',PHE=>'F',
   PRO=>'P',SER=>'S',THR=>'T',ILE=>'I',VAL=>'V'
);
my %unknown_keys;
my $conversion = '';
say "Sequence length: ", length $seq;
for (my $i=0; $i < length $seq; $i  = 3) {
    my $key = substr $seq, $i, 3;
    if (exists $amino_acid_conversion{$key}) {
        $conversion.= $amino_acid_conversion{$key};
    }
    else {
        $unknown_keys{$key}  ;
    }
}
say "Conversion: $conversion";
say "Unknown keys: ", join ",", keys %unknown_keys;
say "Number of atoms: ", $k;
my $xgk=($x/$k);
my $ygk=($y/$k);
my $zgk=($z/$k);
say "Geometric center: $xgk $ygk $zgk";

This gives me the following output:

[...]
Number of atoms: 76015
Geometric center: 290.744642162734 69.196842162731 136.395842938893