Расстояние между одной точкой и всеми остальными в файле PDB

Question

У меня есть файл PDB. Теперь он состоит из двух частей, разделенных TER. Перед TER я называю это частью 1. Я хочу взять x, y, z из ATOM 1 первой части, т.е. до TER, и найти расстояние до всех x, y, z координат после TER, а затем второго 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.

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