如何在Perl中计算蛋白质的几何中心？

ujv3wf0j 于 2023-01-09 发布在 Perl

关注(0)|答案(1)|浏览(137)

我有一个PDB文件，里面包含了特定蛋白质的信息，其中一个信息是不同原子的位置（XYZ坐标）。
文件是下面的https://files.rcsb.org/view/6U9D.pdb。用这个文件我想计算原子的几何中心。理论上我知道我需要做什么，但我写的脚本似乎不起作用。
程序的第一部分，在for循环之前，是任务的另一部分，它要求我读取序列，并将其从3个字母的命名法转换为1个字母的命名法。我感兴趣的部分是for循环，直到最后。我尝试模式匹配，以隔离XYZ坐标。然后我使用了一个计数器，我在开始时设置的变量$k，当我检查cygwin上的输出时，我得到的唯一输出是序列0 0 0，而不是每个维度的和除以$k。知道哪里出了问题吗？

$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";

perl

来源：https://stackoverflow.com/questions/75032662/how-can-i-calculate-the-geometric-center-of-a-protein-in-perl

1条答案

按热度按时间

pqwbnv8z1#

我不懂生物信息学，但看起来你应该这样做：

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";

这将给出以下输出：

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

赞(0）回复(0）举报 2023-01-09

我来回答

如何在Perl中计算蛋白质的几何中心？

1条答案

相关问题

热门标签

最新问答