* スクリプトでPDBファイルを操作しよう [#t887d9b5]
#floatcontents
ここでは原子座標の対称操作と格子定数からScale matrixを求めるスクリプトを紹介します。
** 対称操作した原子座標を取得する [#p154e566]
:&color(red){2011-04-22以前に公開していたmatrix.pyスクリプトは不具合があり対称操作が行われていませんでした};|
行列の計算処理に勘違いがあり、プログラム的にちゃんと処理ができていませんでした。スミマセン・・・~
~
PDBファイル中のATOMレコードの座標は直交座標系(orthogonal coordinate system)なので、例えば(X-1, Y, Z)の対称操作を行おうするとちょっぴり煩雑な操作が必要になります。概要は以下の通り。
+ 直交座標系を格子軸座標系にマトリックス変換
+ Fractional座標に対して対称操作
+ 得られた座標を直交座標系に変換
このような手順で求める座標を得ることができます。手順を確認した後、Pythonで実装してみます⌣。
*** 直交座標系と格子軸座標系 [#q3f1a415]
マトリックス変換やら対称操作やら言葉だけ聞くとややこしそうですが、実際の操作はさほど難しくありません。~
まずはPDB仕様書(v3.2)による直交座標系と格子軸座標系の関係を調べてみましょう。PDB仕様書のCrystallographic and Coordinate Transformation Section、''SCALEn''の解説には以下のような記述があります。
- http://www.wwpdb.org/documentation/format32/sect8.html -- Crystallographic and Coordinate Transformation Section
>
&color(blue){''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&sub{frac};,Y&sub{frac};,Z&sub{frac};)とすると、
- X&sub{frac}; = S11X + S12Y + S13Z + U1
- Y&sub{frac}; = S21X + S22Y + S23Z + U2
- Z&sub{frac}; = S31X + S32Y + S33Z + U3
~
となる。
*** 直交軸→格子軸の操作 [#b53263c3]
さて、定義はわかりましたので、続いてこれらを用いて実際の変換を行ってみます。以下のような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の計算式をよく見るとわかりますが、''行と列が入れ替わった形''になっています。おそらく計算しやすいようにという目的だと思いますが、格子軸→直交軸の変換のために逆行列が必要になるので、結局行列として扱う必要があります(もっと効果的な方法があればお教え下さい⌣)。%%ここで示した行列と移動対象ATOMより求めることができます。
(scale_matrix) × (ATOM座標)
で求められる行列の%%(1,1),(2,2),(3,3)%%''行を合計したもの''が、それぞれ、X&sub{frac};,Y&sub{frac};,Z&sub{frac};となります。つまりX&sub{frac};=(1,1)+(1,2)+(1,3)ですね。(ここを勘違いしてました・・・)
*** 対称操作 [#mfa1c4b6]
続いて、求められた座標に対してようやく対称操作を行うことができます。これは簡単で、対称操作の表記に従えばいいだけです。例えば(-X, Y - 1/2, -Z + 1)なら(-X&sub{frac};, Y&sub{frac}; - 1/2, -Z&sub{frac}; + 1)で求めることができます。
*** 格子軸→直交軸の操作 [#l45d6b5e]
最後に必要になるのが格子軸→直交軸の操作です。これは直交軸→格子軸の逆の操作を行うだけなので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) × (対称操作後ATOM&sub{frac};)
の計算を行うと直交座標系での原子座標を求めることができます。
*** Pythonによる実装 [#yf6e894f]
以上の内容を踏まえて、ある原子に対して対称操作(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]は考慮に入れていません。なお、浮動小数点演算の丸め誤差により他のプログラムとは計算結果が微妙に異なることがあります(小数点の最後の桁)。
#script(matrix2.py)
- 13-15行目：SCALEnよりscale_matrixを作成する
- 16行目：numpy.linalg.inv()(旧：LinearAlgebra.inverse())でscale_matrixの逆行列を求める
- 24行目：移動させる原子座標(直交座標系)
- 27行目：格子軸座標系へ変換
- 28行目：対称操作(idx=0のときがX軸、1のときがY軸、2のときがZ軸です)
- 29行目：格子軸→直交軸の変換を行い対称操作後の原子座標を得る
- 使用したNumPyメソッド
-- numpy.array() -- 行列の構築
-- numpy.linalg.inv() -- 逆行列を求める
-- numpy.ndarray#sum() -- 合計を求める(引数をとらない場合は列を合計する)
:実行結果|
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を格子定数から求める [#r7905d38]
上のスクリプトでは、SCALEnレコードを使って座標の対称移動を行う操作をまとめましたが、PDBファイルに常にSCALEnレコードがあるとは限りません。精密化中であるなどPDBに未登録のデータの場合はCRYST1レコードしかないかもしれません。~
そこで格子定数よりScale matrixを求めるための方法をまとめてみます。つまり直交座標系から格子座標系への変換行列の求め方です。
~
座標の定義をきちんと追いかければイチから求めることができると思いますが、ここでは実践としてNumPyを使った求め方をスクリプトで解説します。
~
式はWikipediaに解説してありますので参考にしてください。
:参考|WikipediaのFractional coorinates -- http://en.wikipedia.org/wiki/Fractional_coordinates
*** Pythonによる実装 [#k5034d14]
Wikipediaで示されている式を元にPythonで実装してみました。PDBファイルを使った場合だと、CRYST1レコードからSCALEnレコードを求めることができれば正解なので確認は楽ですね⌣
#script(scale_matrix.py)
- 18-20行目：cos, sin, sqrtをDecimalで行うためのショートカット定義
- 22-28行目：格子定数の代入
- 29行目：'''v'''の計算(式はWikipediaを参照)
- 31-41行目：Scale matrixの各項を計算
:実行結果|
Scale matrix:
[[ 0.025644 0. 0.006541]
[ 0. 0.015938 0. ]
[ 0. 0. 0.015702]]
