Happy Scripting/PDBファイルを操作しよう

Parmalink: http://biokids.org/ece3b3.html

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

Contents

スクリプトでPDBファイルを操作しよう
- 対称操作した原子座標を取得する
- Scale matrixを格子定数から求める
  - Pythonによる実装

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

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

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

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

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

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

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

マトリックス変換やら対称操作やら言葉だけ聞くとややこしそうですが、実際の操作はさほど難しくありません。
まずは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.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_frac,Y_frac,Z_fracとなります。つまりX_frac=(1,1)+(1,2)+(1,3)ですね。(ここを勘違いしてました・・・)

対称操作

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

これを用いて、

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

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

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]は考慮に入れていません。なお、浮動小数点演算の丸め誤差により他のプログラムとは計算結果が微妙に異なることがあります(小数点の最後の桁)。

matrix2.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
# Symmetry operation implementation sample
# 2011-04-22 by Nob-rin
import re, math, numpy
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",
]

# Parsing SCALEn lines
sca = [re.split(" +", line)[1:4] for line in lines]
scale_matrix = numpy.array(sca, numpy.float64)
scale_inv_matrix = numpy.linalg.inv(scale_matrix)   # Compute the inverse

print "Scale matrix:"
print str(scale_matrix) + "\n"
print "Inverse matrix:"
print str(scale_inv_matrix) + "\n"

# Operating coordinate
orig = numpy.array([13.427,  8.085, 38.568], numpy.float64)
print "Original: (%7.3f, %7.3f, %7.3f)" % tuple(orig)
for idx in range(0, 3):
    frac = (scale_matrix * orig).sum(1)
    frac[idx] += 1     # Symmetry operation
    symm = (scale_inv_matrix * frac).sum(1)
    a = ["X", "Y", "Z"]
    a[idx] += "+1"
    print "%s : (%7.3f, %7.3f, %7.3f)" % ((",".join(a),) + tuple(symm))

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を格子定数から求める

上のスクリプトでは、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.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
# Scale matrix implementation sample
# Calculate scale matrix from cell constants
# 2011-04-22 by Nob-rin
# ref: http://en.wikipedia.org/wiki/Fractional_coordinates
"""
Origina PDB header
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
"""
import math, numpy
from decimal import Decimal

# Shortcuts
def COS(rad): return Decimal(str(math.cos(rad)))
def SIN(rad): return Decimal(str(math.sin(rad)))
def SQRT(v) : return Decimal(str(math.sqrt(v)))

# Definition and compute primitive values
a, b, c = (Decimal("38.996"), Decimal("62.743"), Decimal("65.724"))
alpha, beta, gamma = (Decimal("90.00"), Decimal("104.31"), Decimal("90.00"))
PI = Decimal(str(math.pi))
alpha = alpha / Decimal("180.0") * PI
beta  = beta  / Decimal("180.0") * PI
gamma = gamma / Decimal("180.0") * PI
V = SQRT(1 - COS(alpha) ** 2 - COS(beta) ** 2 - COS(gamma) ** 2 + 2 * COS(alpha) * COS(beta) * COS(gamma))

# Compute scale matrix
mat = numpy.zeros([3,3], numpy.float64)
mat[0][0] = 1 / a
mat[0][1] = -(COS(gamma) / (a * SIN(gamma)))
mat[0][2] = (COS(alpha) * COS(gamma) - COS(beta)) / (a * V * SIN(gamma))
mat[1][0] = Decimal(0)
mat[1][1] = 1 / (b * SIN(gamma))
mat[1][2] = (COS(beta) * COS(gamma) - COS(alpha)) / (b * V * SIN(gamma))
mat[2][0] = Decimal(0)
mat[2][1] = Decimal(0)
mat[2][2] = SIN(gamma) / (c * V)

# Rounding for 6 ditits
for y in range(3):
    for x in range(3):
        mat[x][y] = Decimal("%.6f" % mat[x][y])
print "Scale matrix:"
print mat

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]]

添付ファイル:

scale_matrix.py 2100件 [詳細]

matrix2.py 1992件 [詳細]

PDBファイルを操作しよう

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

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

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

直交軸→格子軸の操作

対称操作

格子軸→直交軸の操作

Pythonによる実装

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

Pythonによる実装

サイト内検索

クックブック

ナビゲーション

人気の20件

最新の10件

2022-07-22

2019-10-02

2017-12-11

2015-09-30

2014-11-14

2014-10-18

2014-10-08

2014-07-19

2014-07-05

2013-12-09

管理用