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

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

3 голосов
/ 06 марта 2012

У меня есть файл 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";

Ответы [ 3 ]

1 голос
/ 06 марта 2012

Не уверен, что я четко понимаю, что вы хотите, но как насчет: 

#!/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
1 голос
/ 06 марта 2012

Основная проблема здесь - чтение данных.Во-первых, обратите внимание, что нельзя использовать split с текстовыми файлами PDB, поскольку поля определяются положением, а не разделителями.См. Описание координатного файла (формат PDB) .

Чтобы разделить запись ATOM для различных полимеров цепей , вы можете начать с упрощенной версии, такой как 

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 .

0 голосов
/ 06 марта 2012

При правильной инкапсуляции это довольно просто и требует незначительных изменений вашего кода.

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
...