3

以下に示すようなテキストファイルがあります

ATOM    920  CA  GLN A 203      39.292 -13.354  17.416  1.00 55.76           C 
ATOM    929  CA  HIS A 204      38.546 -15.963  14.792  1.00 29.53           C
ATOM    939  CA  ASN A 205      39.443 -17.018  11.206  1.00 54.49           C  
ATOM    947  CA  GLU A 206      41.454 -13.901  10.155  1.00 26.32           C
ATOM    956  CA  VAL A 207      43.664 -14.041  13.279  1.00 40.65           C 
.
.
.

ATOM    963  CA  GLU A 208      45.403 -17.443  13.188  1.00 40.25           C

2 つのアルファ炭素原子間の距離を計算したいと思います。つまり、最初の原子と 2 番目の原子の間の距離、次に 2 番目と 3 番目の原子の間の距離を計算します。2 つの原子間の距離は次のように表すことができます。distance = sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2) .

列 7、8、および 9 は、それぞれ x、y、および z 座標を表します。以下に示すように、距離と対応する残基のペア (列 4) を出力する必要があります (距離の値は実数ではありません)。

GLN-HIS   4.5
HIS-ASN   3.2
ASN-GLU   2.5

この計算を perl や python で行うにはどうすればよいですか?

4

6 に答える 6

6

空白で分割しないでください

ここで与えられた他の答えは、座標がスペースで区切られるという誤った仮定をしています。のPDB 仕様にATOMよると、これは必ずしもそうではありません。PDB レコードの値は列インデックスによって指定され、互いに流れ込む可能性があります。たとえば、最初のATOMレコードは次のようになります。

ATOM    920  CA  GLN A 203      39.292 -13.354  17.416  1.00 55.76           C

しかし、これも完全に有効です。

ATOM    920  CA  GLN A 203      39.292-13.3540  17.416  1.00 55.76           C

より良いアプローチ

列指定のインデックス、および PDB ファイルで発生する可能性のあるその他の問題の数のため、独自のパーサーを作成しないでください。PDB 形式は扱いにくく、多くの特殊なケースや不適切な形式のファイルを処理する必要があります。代わりに、すでに作成されているパーサーを使用してください。

私はBiopythonが好きですPDB.PDBParser。便利な機能を備えた Python オブジェクトに構造を解析します。Perl の方が好きなら、BioPerlをチェックしてください。

PDB.Residueオブジェクトは、名前によるアトムへのキーによるアクセスを許可し、PDB.Atomオブジェクトはオペレーターをオーバーロードして、-2 つのアトム間の距離を返します。これを使用して、きれいで簡潔なコードを書くことができます。

コード

from Bio import PDB
parser = PDB.PDBParser()

# Parse the structure into a PDB.Structure object
pdb_code = "1exm"
pdb_path = "pdb1exm.ent"
struct = parser.get_structure(pdb_code, pdb_path)

# Grab the first two residues from the structure
residues = struct.get_residues()
res_one = residues.next()
res_two = residues.next()

try:
    alpha_dist = res_one['CA'] - res_two['CA']
except KeyError:
    print "Alpha carbon missing, computing distance impossible!"
于 2012-12-01T17:07:57.153 に答える
5

データが空白で区切られている場合、単純なsplit作業を行うことができます。行をバッファリングして、それらを順次比較します。

use strict;
use warnings;

my @line;
while (<>) {
    push @line, $_;            # add line to buffer
    next if @line < 2;         # skip unless buffer is full
    print proc(@line), "\n";   # process and print 
    shift @line;               # remove used line 
}

sub proc {
    my @a = split ' ', shift;   # line 1
    my @b = split ' ', shift;   # line 2
    my $x = ($a[6]-$b[6]);      # calculate the diffs
    my $y = ($a[7]-$b[7]);
    my $z = ($a[8]-$b[8]);
    my $dist = sprintf "%.1f",                # format the number
                   sqrt($x**2+$y**2+$z**2);   # do the calculation
    return "$a[3]-$b[3]\t$dist"; # return the string for printing
}

出力 (サンプル データを使用):

GLN-HIS 3.8
HIS-ASN 3.8
ASN-GLU 3.9
GLU-VAL 3.8

データがタブで区切られている場合は、/\t/の代わりに分割できます' '。

于 2012-11-30T13:50:58.007 に答える
4

データが「atoms.txt」にあると仮定すると、これはデータを1行ずつ読み取り、エントリをリストに分割します。

import csv

with open("atoms.txt") as f:
  reader = csv.reader(f)
  for line, in reader:
      entries = line.split()
      print entries

次に、リストごとに必要な列を抽出し、距離などを計算します（Pythonのリストはゼロベースであることに注意してください）。

于 2012-11-30T12:44:17.230 に答える
3

理想的には、原子とセグメントを「スライス」し、それらの間の距離測定値を計算する Pythonic な方法にMDAnalysis パッケージを使用する必要があります。実際、MDAnalysis は複数の MD シミュレーションおよび化学構造フォーマットをサポートしています。

もう少し詳細な例については、Biostars.org の次のエントリも参照してください。

于 2013-05-02T07:26:17.893 に答える
1

1 つのペアだけに関心がある場合は、bash で問題なく動作します。これは私が使用するスクリプトです。最後に再起動するように設定しています (必要に応じてオフにします)。どのアトムを指定するかを尋ねられます。PDB ファイルには異なる列が設定されている可能性があるため、awk 行では、列が一致していることを確認してください。新しい pdb ファイルで使用する前に、手動でテスト ケースを実行します。これは些細なことですが、スクリプトで my pdb file を yours に変更します。

#!/usr/bin/env bash

echo " "
echo "Use one letter abbreviations. Case doesn't matter." 
echo "Example: A 17 CA or n 162 cg"

echo " - - - - - - - - - - - - - - - - - -"
#------------First Selection------------

read -e -p "What first atom? " sel1

# echo $sel1
sel1caps=${sel1^^}
# echo "sel1caps="$sel1caps

arr1=($sel1caps)
# echo "arr1[0]= "${arr1[0]}
# echo "arr1[1]= "${arr1[1]}
# echo "arr1[2]= "${arr1[2]}

#To convert one to three letters

if [ ${arr1[0]} = A ] ; then
    SF1=ALA
elif [ ${arr1[0]} = H ] ; then
    SF1=HIS
elif [ ${arr1[0]} = R ] ; then
    SF1=ARG
elif [ ${arr1[0]} = K ] ; then
    SF1=LYS
elif [ ${arr1[0]} = I ] ; then
    SF1=ILE
elif [ ${arr1[0]} = F ] ; then
    SF1=PHE 
elif [ ${arr1[0]} = L ] ; then
    SF1=LEU
elif [ ${arr1[0]} = W ] ; then
    SF1=TRP
elif [ ${arr1[0]} = M ] ; then
    SF1=MET
elif [ ${arr1[0]} = P ] ; then
    SF1=PRO 
elif [ ${arr1[0]} = C ] ; then
    SF1=CYS 
elif [ ${arr1[0]} = N ] ; then
    SF1=ASN
elif [ ${arr1[0]} = V ] ; then
    SF1=VAL
elif [ ${arr1[0]} = G ] ; then
    SF1=GLY 
elif [ ${arr1[0]} = S ] ; then
    SF1=SER 
elif [ ${arr1[0]} = Q ] ; then
    SF1=GLN 
elif [ ${arr1[0]} = Y ] ; then
    SF1=TYR 
elif [ ${arr1[0]} = D ] ; then
    SF1=ASP
elif [ ${arr1[0]} = E ] ; then
    SF1=GLU 
elif [ ${arr1[0]} = T ] ; then
    SF1=THR 
else
    echo "use one letter codes"
    echo "exiting"
    exit
fi

# echo "SF1 ="$SF1

#If there is nothing printing for line 1, check the expression for your pdb file. The spaces may need adjustment at the end.
line1=$(grep -E "${arr1[2]} *?${SF1}(.*?) ${arr1[1]}     " 1A23.pdb)
# echo $line1

ar_l1=($line1)
# echo "ar_l1[1]="${ar_l1[1]}

echo " - - - - - - - - - - - - - - - - - -"
#------------Second Selection------------

read -e -p "What second atom? " sel2

# echo $sel2

sel2caps=${sel2^^}
# echo "sel2caps="$sel2caps

arr2=($sel2caps)
# echo "arr2[0]= "${arr2[0]}
# echo "arr2[1]= "${arr2[1]}
# echo "arr2[2]= "${arr2[2]}

#To convert one to three letters

if [ ${arr2[0]} = A ] ; then
    SF2=ALA
elif [ ${arr2[0]} = H ] ; then
    SF2=HIS
elif [ ${arr2[0]} = R ] ; then
    SF2=ARG
elif [ ${arr2[0]} = K ] ; then
    SF2=LYS
elif [ ${arr2[0]} = I ] ; then
    SF2=ILE
elif [ ${arr2[0]} = F ] ; then
    SF2=PHE 
elif [ ${arr2[0]} = L ] ; then
    SF2=LEU
elif [ ${arr2[0]} = W ] ; then
    SF2=TRP
elif [ ${arr2[0]} = M ] ; then
    SF2=MET
elif [ ${arr2[0]} = P ] ; then
    SF2=PRO 
elif [ ${arr2[0]} = C ] ; then
    SF2=CYS 
elif [ ${arr2[0]} = N ] ; then
    SF2=ASN
elif [ ${arr2[0]} = V ] ; then
    SF2=VAL
elif [ ${arr2[0]} = G ] ; then
    SF2=GLY 
elif [ ${arr2[0]} = S ] ; then
    SF2=SER 
elif [ ${arr2[0]} = Q ] ; then
    SF2=GLN 
elif [ ${arr2[0]} = Y ] ; then
    SF2=TYR 
elif [ ${arr2[0]} = D ] ; then
    SF2=ASP
elif [ ${arr2[0]} = E ] ; then
    SF2=GLU 
elif [ ${arr2[0]} = T ] ; then
    SF2=THR 
else
    echo "use one letter codes"
    echo "exiting"
    exit
fi

# echo "SF2 ="$SF2


line2=$(grep -E "${arr2[2]} *?${SF2}(.*?) ${arr2[1]}     " 1A23.pdb)
# echo $line2

ar_l2=($line2)
# echo "ar_l2[1]="${ar_l2[1]}
# echo "ar_l2[1]="${ar_l2[1]}

atom1=${ar_l1[1]}
atom2=${ar_l2[1]}
echo "==========================="
echo ${arr1[0]}${arr1[1]}${arr1[2]}" to "${arr2[0]}${arr2[1]}${arr2[2]}":"
# 6, 7, 8 are column numbers in the pdb file. 
# If there are multiple molecules it should be 7, 8, 9.
awk '$2=='$atom1'{x1=$7;y1=$8;z1=$9}                                 # get the ATOM 1
     $2=='$atom2'{x2=$7;y2=$8;z2=$9}                               # get the ATOM 2
     END{print sqrt((x1-x2)*(x1-x2) + (y1-y2)*(y1-y2) + (z1-z2)*(z1-z2))}' 1A23.pdb    # calculate the distance.

echo "Angstroms"
echo "==========================="
echo " "
echo "-_-_-_-_Running Script Again_-_-_-_-"
./distance_soln.sh
于 2014-09-03T00:52:35.847 に答える