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