0

私は2つのチェーンを持つ300のpdbファイを持っています。最初の鎖の最初の原子から2番目の鎖のすべての原子までの距離を計算します。次に、最初の鎖の2番目の原子から2番目の鎖のすべての原子までの距離を計算します。これは、300ファイルに対して繰り返す必要があります。距離が5以上の場合にのみ原子ペアを印刷し、入力ファイル名を使用して出力を別のフォルダーに保存する必要があります。距離を求める式はsqrt（（x1-x2）^ 2 +（y1-y2）^ 2 +（z1-z2）^ 2）です。$ 5はチェーンIDで、$12はアトム名です。$ 7、$ 8、$ 9はx、y、z座標です。あなたの貴重な提案をいただければ幸いです！

ATOM      1  N   MET A   1     -16.220  53.312  36.564  1.00 32.19           N  
ATOM      2  CA  MET A   1     -15.722  52.290  37.522  1.00 28.47           C
ATOM      3  C   MET A   1     -14.451  51.635  37.011  1.00 26.82           C 
ATOM   2542  CG  ASN B  17      -1.077   9.776  13.155  1.00 18.23           C  
ATOM   2543  OD1 ASN B  17      -0.563   9.098  12.250  1.00 18.58           O 
ATOM   2544  ND2 ASN B  17      -0.632   9.746  14.418  1.00 14.82           N

必要な出力（距離の値が正しくありません）

N-C   8.90
N-O   10.3
N-N   7.62
C-C   12.45
C-O   9.0
C-N   9.89 
C-C   11.45
C-O   19.0
C-N   10.89
4

3 に答える 3

1

を使用する 1 つの方法を次に示しGNU awkます。次のように実行します。

awk -f script.awk file{,}

の内容script.awk

NR==1 {
    n = $5
}

FNR==NR && $5 != n {
    a[c++]=$0
}

FNR!=NR && $5 == n {
    for (i=0;i<=c-1;i++) {
        split (a[i],b)
        dist = sqrt (($7-b[7])^2 + ($8-b[8])^2 + ($9-b[9])^2)
        if (dist >= 5) {
            printf "%s-%s\t%.2f\n", $NF, b[NF], dist
        }
    }
}

タブ区切りの結果:

N-C 51.70
N-O 52.83
N-N 51.30
C-C 51.14
C-O 52.29
C-N 50.71
C-C 50.00
C-O 51.14
C-N 49.56

または、ここにワンライナーがあります：

awk 'NR==1 { n = $5 } FNR==NR && $5 != n { a[c++]=$0 } FNR!=NR && $5 == n { for (i=0;i<=c-1;i++) { split (a[i],b); dist = sqrt (($7-b[7])^2 + ($8-b[8])^2 + ($9-b[9])^2); if (dist >= 5) printf "%s-%s\t%.2f\n", $NF, b[NF], dist } }' file{,}

したがって、現在の作業ディレクトリ内の複数のファイルに対してこれを実行し、このディレクトリに目的のファイルしかないと仮定すると、ステートメントをforループで囲むことができます。明らかに、正しく機能させるためには、選択したパスawkに変更する必要があります。/path/to/folder/

for i in *; do awk 'NR==1 { n = $5 } FNR==NR && $5 != n { a[c++]=$0 } FNR!=NR && $5 == n { for (i=0;i<=c-1;i++) { split (a[i],b); dist = sqrt (($7-b[7])^2 + ($8-b[8])^2 + ($9-b[9])^2); if (dist >= 5) printf "%s-%s\t%.2f\n", $NF, b[NF], dist > "/path/to/folder/" FILENAME } }' "$i"{,}; done
于 2012-12-10T05:21:01.497 に答える
0

Python の MDAnalysis モジュールを見てみましょう -- http://code.google.com/p/mdanalysis/

http://code.google.com/p/mdanalysis/wiki/Examples#Extracting_a_chain_from_a_PDB

于 2012-12-10T03:45:23.760 に答える
0

サンプルデータを「テスト」として保存すると、次のようになります。

#! /usr/bin/python3.2

import re

class Atom:
    def __init__ (self, line):
        tokens = re.split (' +', line)
        self.chain = tokens [4]
        self.x = float (tokens [6] )
        self.y = float (tokens [7] )
        self.z = float (tokens [8] )
        self.name = tokens [11]

    def __repr__ (self):
        return self.name

    def distance (self, other):
        return ( (self.x - other.x) ** 2 + (self.y - other.y) ** 2 + (self.z - other.z) ** 2) ** .5

chains = [ [], [] ]

with open ('test', 'r') as pdb:
    for line in pdb.readlines ():
        atom = Atom (line.strip () )
        chains [0 if atom.chain == 'A' else 1].append (atom)

print ('\n'.join ( ['{} - {}\t{}'.format (a, b, a.distance (b) ) for a in chains [0] for b in chains [1] ] ) )

そのようなものを生成します：

N - C   51.697920905970676
N - O   52.8317143484858
N - N   51.29744063791097
C - C   51.14359109409506
C - O   52.28783920760161
C - N   50.709908814747436
C - C   50.001484907950484
C - O   51.140786403808846
C - N   49.559022700210704

.

編集：

あなたの >= 5.0 条件を忘れました: 最後の行を次のように置き換えてください:

print ('\n'.join ( ['{} - {}\t{}'.format (a, b, a.distance (b) ) for a in chains [0] for b in chains [1] if a.distance (b) >= 5.0] ) )
于 2012-12-10T04:17:27.737 に答える