PDB文件中一点到其他所有点之间的距离

Question

我有一个PDB文件。现在它有两个由TER分隔的部分。在TER之前，我把它称为第一部分。我想要把第一部分的ATOM 1（即TER之前）的x，y，z和TER之后的所有x，y，z坐标的距离，然后找到第一部分的第二个ATOM到所有ATOMS第二部分。对于第一部分的所有ATOMS =必须对第二部分的所有ATOMS重复。我必须将它自动化为20个文件。我的文件名称开始像1_0.pdb，2_0.pdb .... 20_0.pdb。 这是一个距离计算。我在PERL中尝试了一些东西，但非常粗糙。有人可以帮助一下。 文件看起来像：

----长文件（我截断它）----

ATOM 1279 C ALA 81  -1.925 -11.270 1.404 
ATOM 1280 O ALA 81  -0.279 9.355 15.557 
ATOM 1281 OXT ALA 81  -2.188 10.341 15.346 
TER 
ATOM 1282 N THR 82  29.632 5.205 5.525 
ATOM 1283 H1 THR 82  30.175 4.389 5.768 
ATOM 1284 H2 THR 82  28.816 4.910 5.008 

的代码是：到底它发现的最大距离及其共同坐标

my @points =(); 
open(IN, @ARGV[0]) or die "$!"; 
while (my $line = <IN>) { 

    chomp($line); 
    my @array = (split (/\s+/, $line))[5, 6, 7]; 
    print "@array\n"; 
    push @points, [ @array ]; 
} 
close(IN); 


$max=0; 
for my $i1 (0 .. $#points ) 

{ 
    my ($x1, $y1, $z1) = @{ $points[$i1] }; 
    my $dist = sqrt(($x1+1.925)**2 + ($y1+11.270)**2 + ($z1-1.404)**2); 
    print "distance from (-1.925 -11.270 1.404) to ($x1, $y1, $z1) is $dist\n"; 

    if ($dist > $max) 
    { $max = $dist; 
     $x=$x1; 
     $y=$y1; 
     $z=$z1; 
     }} 
    print "maximum value is : $max\n"; 
print "co ordinates are : $x $y $z\n";   

Answer 1

不知道我清楚地了解你想要什么，而是怎么样：

#!/usr/local/bin/perl 
use strict; 
use warnings; 

my (@refer, @points); 
my $part = 0; 
while (my $line = <DATA>) { 
    chomp($line); 
    if ($line =~ /^TER/) { 
     $part++; 
     next; 
    } 
    my @array = (split (/\s+/, $line))[5, 6, 7]; 
    if ($part == 0) { 
     push @refer, [ @array ]; 
    } else { 
     push @points, [ @array ]; 
    } 
} 
my %max = (val=>0, x=>0, y=>0, z=>0); 
foreach my $ref(@refer) { 
    my ($x1, $y1, $z1) = @{$ref}; 
    foreach my $atom(@points) { 
     my ($x, $y, $z) = @{$atom}; 
     my $dist = sqrt(($x-$x1)**2 + ($y-$y1)**2 + ($z-$z1)**2); 
     if ($dist > $max{val}) { 
      $max{val} = $dist; 
      $max{x} = $x; 
      $max{y} = $y; 
      $max{z} = $z; 
     } 
    } 
} 
print "max is $max{val}; coord: x=$max{x}, y=$max{y}, z=$max{z}\n"; 

__DATA__ 
ATOM 1279 C ALA 81  -1.925 -11.270 1.404 
ATOM 1280 O ALA 81  -0.279 9.355 15.557 
ATOM 1281 OXT ALA 81  -2.188 10.341 15.346 
TER 
ATOM 1282 N THR 82  29.632 5.205 5.525 
ATOM 1283 H1 THR 82  30.175 4.389 5.768 
ATOM 1284 H2 THR 82  28.816 4.910 5.008 

输出：

max is 35.9813670807545; coord: x=30.175, y=4.389, z=5.768 

Answer 2

这里的主要问题是读取数据。首先，请注意，由于字段是由位置定义的，而不是由分隔符定义的，因此不能使用拆分PDB文本文件。见Coordinate File Description (PDB Format)。

为了分离不同的聚合物的ATOM记录chains可以用一个简化的版本等

my $iblock = 0; 
my @atoms =(); 
while (my $line = <IN>) { 
    chomp($line); 

    # Switch blocks at TER lines 
    if ($line =~ /^TER/) { 
     $iblock++; 

    # Read ATOM lines 
    } elsif ($line =~ m/^ATOM/) { 
     my @xyz = (substr($line,7-1,9),substr($line,16-1,9),substr($line,25-1,9)); 
     printf "Block %d: atom at (%s)\n",$iblock,join (",",@xyz); 
     push @{$atoms[$iblock]},\@xyz; 

    # Parse additional line types (if needed) 
    } else { 
     ... 
    } 
} 

随之而来的是遍历所有对从不同的块，构成为坐标的开始如下：

# 1st block 
for my $iblock1 (0..$#atoms) { 

    # 2nd block 
    for my $iblock2 ($iblock1+1..$#atoms) { 

     # Compare all pairs of atoms 
     ... 
     my $xyz1 (@{$atoms[$iblock1]}) { 
     for my $xyz2 (@{$atoms[$iblock2]}) { 
      # Calculate distance and compare with $max_dist 
      ... 
     } 
     } 
     # Print the maximal distance between these two blocks 
     ... 
    } 
} 

当然，如果使用更复杂的数据结构或通过应用可用的PDB解析器之一（如Bioperl's），代码可能更通用。

Answer 3

通过适当的封装，这是非常简单，只需要轻微修改您的码。

ETA：增加了我手边的固定宽度解决方案。读取所有字段可能是最好的，而不是放弃前31个字符，然后将它们全部返回到散列引用中。这样，您可以使用相同的子程序处理所有行，并且只需在第一个字段变为TER时在零件之间切换即可。您应该很容易从给定的代码中推断出这一点。

你会注意到参考值是用循环读入的，因为我们需要在中断点处打开循环。其余的值是用map声明浑浊起来的。然后，我们只需将数据提供给我们从您的初始代码制作的子例程（有一些改进）。我为词法变量使用了相同的名称，以便于阅读代码。

use strict; 
use warnings; 

my @points; 
while (<DATA>) { 
    last if /^TER$/; 
    push @points, getpoints($_); 
} 
my @ref = map getpoints($_), <DATA>; 

for my $p (@points) { 
    getcoords($p, \@ref); 
} 

sub getpoints { 
    my $line = shift; 
    my @data = unpack "A31 A8 A8 A8", $line; 
    shift @data; 
    return \@data; 
} 
sub getcoords { 
    my ($p, $ref) = @_; 
    my ($p1,$p2,$p3) = @$p; 
    my $max=0; 
    my ($x,$y,$z); 
    for my $aref (@$ref) { 
     my ($x1, $y1, $z1) = @$aref; 
     my $dist = sqrt(
      ($x1-$p1)**2 + 
      ($y1-$p2)**2 + 
      ($z1-$p3)**2 
     ); 
     print "distance from ($p1 $p2 $p3) to ($x1, $y1, $z1) is $dist\n"; 

     if ($dist > $max) { 
      $max = $dist; 
      $x=$x1; 
      $y=$y1; 
      $z=$z1; 
     } 
    } 
    print "maximum value is : $max\n"; 
    print "co ordinates are : $x $y $z\n"; 
} 

__DATA__ 
ATOM 1279 C ALA 81  -1.925 -11.270 1.404 
ATOM 1280 O ALA 81  -0.279 9.355 15.557 
ATOM 1281 OXT ALA 81  -2.188 10.341 15.346 
TER 
ATOM 1282 N THR 82  29.632 5.205 5.525 
ATOM 1283 H1 THR 82  30.175 4.389 5.768 
ATOM 1284 H2 THR 82  28.816 4.910 5.008