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

ここでは原子座標の対称操作と格子定数からScale matrixを求めるスクリプトを紹介します。

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

2011-04-22以前に公開していたmatrix.pyスクリプトは不具合があり対称操作が行われていませんでした
行列の計算処理に勘違いがあり、プログラム的にちゃんと処理ができていませんでした。スミマセン・・・

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

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

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

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

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

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レコードで示される格子軸座標系と以下の定義で関係づけられる。
ベクトルabcが結晶学的な格子辺を表し、ベクトルABCがオングストローム標準直交座標系の単位格子ベクトルを表すとすると、ABCおよびabcは同じ原点を持ち、Aaと平行でBC×A(外積)と平行かつCa×bと平行となる。オングストローム直交座標を(X,Y,Z)、Fractional格子座標を(Xfrac,Yfrac,Zfrac)とすると、

  • Xfrac = S11X + S12Y + S13Z + U1
  • Yfrac = S21X + S22Y + S23Z + U2
  • Zfrac = 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.006541)
               (0.000000, 0.015938, 0.000000)
               (0.000000, 0.000000, 0.015702)

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

(scale_matrix) × (ATOM座標)

で求められる行列の(1,1),(2,2),(3,3)行を合計したものが、それぞれ、Xfrac,Yfrac,Zfracとなります。つまりXfrac=(1,1)+(1,2)+(1,3)ですね。(ここを勘違いしてました・・・)

対称操作  

続いて、求められた座標に対してようやく対称操作を行うことができます。これは簡単で、対称操作の表記に従えばいいだけです。例えば(-X, Y - 1/2, -Z + 1)なら(-Xfrac, Yfrac - 1/2, -Zfrac + 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,-16.244390)
                   (  0.000000, 62.743130,  0.000000)
                   (  0.000000,  0.000000, 63.686155)

これを用いて、

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

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

Pythonによる実装  

以上の内容を踏まえて、ある原子に対して対称操作(X+1, Y, Z), (X, Y+1, Z), (X, Y, Z+1)を行う処理をPythonで実装してみました。Pythonの算術演算モジュールであるNumPyを使用しています。Vine 5.2 + Python 2.5.4で動作確認しました。
この例ではSCALEnのu[n]は考慮に入れていません。なお、浮動小数点演算の丸め誤差により他のプログラムとは計算結果が微妙に異なることがあります(小数点の最後の桁)。

実行結果
Scale matrix:
[[ 0.025644  0.        0.006541]
 [ 0.        0.015938  0.      ]
 [ 0.        0.        0.015702]]

Inverse matrix:
[[ 38.99547652   0.         -16.24439001]
 [  0.          62.74312963   0.        ]
 [  0.           0.          63.68615463]]

Original: ( 13.427,   8.085,  38.568)
X+1,Y,Z : ( 52.422,   8.085,  38.568)
X,Y+1,Z : ( 13.427,  70.828,  38.568)
X,Y,Z+1 : ( -2.817,   8.085, 102.254)

Scale matrixを格子定数から求める  

上のスクリプトでは、SCALEnレコードを使って座標の対称移動を行う操作をまとめましたが、PDBファイルに常にSCALEnレコードがあるとは限りません。精密化中であるなどPDBに未登録のデータの場合はCRYST1レコードしかないかもしれません。
そこで格子定数よりScale matrixを求めるための方法をまとめてみます。つまり直交座標系から格子座標系への変換行列の求め方です。
座標の定義をきちんと追いかければイチから求めることができると思いますが、ここでは実践としてNumPyを使った求め方をスクリプトで解説します。
式はWikipediaに解説してありますので参考にしてください。

参考
WikipediaのFractional coorinates -- http://en.wikipedia.org/wiki/Fractional_coordinates

Pythonによる実装  

Wikipediaで示されている式を元にPythonで実装してみました。PDBファイルを使った場合だと、CRYST1レコードからSCALEnレコードを求めることができれば正解なので確認は楽ですね [smile]

実行結果
Scale matrix:
[[ 0.025644  0.        0.006541]
 [ 0.        0.015938  0.      ]
 [ 0.        0.        0.015702]]
添付ファイル: filescale_matrix.py 975件 [詳細] filematrix2.py 881件 [詳細]