PDBファイルを操作しよう

Last modified: 2010-04-08 (木) 13:31:19 (149d)

• [ トップ ]
• [ 差分 | 一覧 | 最終更新 | ヘルプ ]

スクリプトでPDBファイルを操作しよう  

対称操作した原子座標を取得する  

PDBファイル中のATOMレコードの座標は直交座標系(orthogonal coordinate system)なので、例えば(X-1, Y, Z)の対称操作を行おうするとちょっぴり煩雑な操作が必要になります。概要は以下の通り。

1. 直交座標系を格子軸座標系にマトリックス変換
2. Fractional座標に対して対称操作
3. 得られた座標を直交座標系に変換

このような手順で求める座標を得ることができます。手順を確認した後、Pythonで実装してみます 。

直交座標系と格子軸座標系  

マトリックス変換やら対称操作やら言葉だけ聞くとややこしそうですが、実際の操作はさほど難しくありません。
まずはPDB仕様書(v3.2)による直交座標系と格子軸座標系の関係を調べてみましょう。PDB仕様書のCrystallographic and Coordinate Transformation Section、SCALEnの解説には以下のような記述があります。

• http://www.wwpdb.org/documentation/format32/sect8.html -- Crystallographic and Coordinate Transformation Section

SCALEn

COLUMNS         DATA  TYPE    FIELD              DEFINITION
------------------------------------------------------------------
 1 -  6         Record name   "SCALEn" n=1,  2, or 3
11 - 20         Real(10.6)    s[n][1]            Sn1
21 - 30         Real(10.6)    s[n][2]            Sn2
31 - 40         Real(10.6)    s[n][3]            Sn3
46 - 55         Real(10.5)    u[n]               Un

使用されているオングストローム標準直交座標系はCRYST1レコードで示される格子軸座標系と以下の定義で関係づけられる。
ベクトルa、b、cが結晶学的な格子辺を表し、ベクトルA、B、Cがオングストローム標準直交座標系の単位格子ベクトルを表すとすると、A、B、Cおよびa、b、cは同じ原点を持ち、Aはaと平行でBはC×A(外積)と平行かつCはa×bと平行となる。オングストローム直交座標を(X,Y,Z)、Fractional格子座標を(X_frac,Y_frac,Z_frac)とすると、

• X_frac = S11X + S12Y + S13Z + U1
• Y_frac = S21X + S22Y + S23Z + U2
• Z_frac = S31X + S32Y + S33Z + U3
  となる。

直交軸→格子軸の操作  

さて、定義はわかりましたので、続いてこれらを用いて実際の変換を行ってみます。以下のようなPDBに対して操作を行うものとします。

         1         2         3         4         5         6         7         8
12345678901234567890123456789012345678901234567890123456789012345678901234567890
CRYST1   38.996   62.743   65.724  90.00 104.31  90.00 P 1 21 1      8
SCALE1      0.025644  0.000000  0.006541        0.00000
SCALE2      0.000000  0.015938  0.000000        0.00000
SCALE3      0.000000  0.000000  0.015702        0.00000
ATOM      4  O   GLY A   3      15.323  15.013  43.462  1.00 44.59           O  <-- これを変換

直交座標系から格子軸座標系への変換には上の解説の通りSCALE1,2,3を使用します。簡単に言うとATOMの座標に対し、行列による変換を行えばいいわけです。変換のための行列はSCALE1,2,3より、

scale_matrix = (0.025644, 0.000000, 0.000000)
               (0.000000, 0.015938, 0.000000)
               (0.006541, 0.000000, 0.015702)

となります。ここで気をつけなければならないことは、SCALEnレコードの行が行列の行に対応していないということです。SCALEnの計算式をよく見るとわかりますが、行と列が入れ替わった形になっています。おそらく計算しやすいようにという目的だと思いますが、格子軸→直交軸の変換のために逆行列が必要になるので、結局行列として扱う必要があります(もっと効果的な方法があればお教え下さい )。
ここで示した行列と移動対象ATOMより求めることができます。

(ATOM座標) × (scale_matrix)

で求められる行列の(1,1),(2,2),(3,3)がそれぞれ、X_frac,Y_frac,Z_fracとなります。

対称操作  

続いて、求められた座標に対してようやく対称操作を行うことができます。これは簡単で、対称操作の表記に従えばいいだけです。例えば(-X, Y - 1/2, -Z + 1)なら(-X_frac, Y_frac - 1/2, -Z_frac + 1)で求めることができます。

格子軸→直交軸の操作  

最後に必要になるのが格子軸→直交軸の操作です。これは直交軸→格子軸の逆の操作を行うだけなのでscale_matrixの逆行列を求めてやればいいわけです。逆行列とは、

(scale_matrix) × (scale_inv_matrix) = (unit matrix)

となる行列のことです。unit matrixは単位行列で(1,1),(2,2)...が1でそれ以外が0の行列のことですね。単位行列を求める方法はちゃんとあると思いますが、私はPythonのモジュールを使って求めました。それを求めると以下のようになります。小数第七位以下を丸めてあります。

scale_inv_matrix = ( 38.995477, 0.000000, 0.000000)
                   (  0.000000,62.743130, 0.000000)
                   (-16.244390, 0.000000,63.686155)

これを用いて、

(対称操作後ATOMfrac) × (scale_inv_matrix)

の計算を行うと直交座標系での原子座標を求めることができます。

Pythonによる実装  

以上の内容を踏まえて、ある原子に対して対称操作(X-1, Y, Z)を行う処理をPythonで実装してみました。Pythonの算術演算モジュールであるNumPyを使用しています。Vine 5.1 + Python 2.5.4で動作確認しました。
この例ではSCALEnのu[n]は考慮に入れていません。

• matrix.py

  1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

  #!/usr/bin/python # Symmetry operation implementation sample # 2010-04-08 by Nob-rin import sys, re, math, Numeric, LinearAlgebra from decimal import Decimal lines = [ "SCALE1 0.025644 0.000000 0.006541 0.00000", "SCALE2 0.000000 0.015938 0.000000 0.00000", "SCALE3 0.000000 0.000000 0.015702 0.00000", ] # Calculating matrices scale_matrix = Numeric.zeros([3,3],"decimal") sca = [] for line in lines: sca.append(re.split(" +", line)[1:]) for y in range(3): for x in range(3): scale_matrix[y][x] = Decimal(sca[x][y],6) print "Scale matrix:" print str(scale_matrix) + "\n" scale_inv_matrix = LinearAlgebra.inverse(scale_matrix) print "Inverse matrix:" print str(scale_inv_matrix) + "\n" # Converting coordinate orig = Numeric.array([15.323,15.013,43.462],"d") frac = orig * scale_matrix frac[0][0] -= 1 # Symmetry operation for (X-1,Y,Z) symm_mat = frac * scale_inv_matrix symm = Numeric.array([symm_mat[0][0], symm_mat[1][1], symm_mat[2][2]],"d") print "Original: (%7.3f, %7.3f, %7.3f)" % tuple(orig) print "Symmetry: (%7.3f, %7.3f, %7.3f)" % tuple(symm)

• 13-17行目：scale_matrixを求める
• 20行目：LinearAlgebra.inverse()でscale_matrixの逆行列を求める
• 25行目：移動させる原子座標(直交座標系)
• 26行目：格子軸座標系へ変換
• 27行目：X-1,Y,Zの対称操作
• 28-29行目：格子軸→直交軸の変換を行い対称操作後の原子座標を得る

添付ファイル: matrix.py 135件 [詳細]