PDBファイルの糖鎖位置合わせ
データベースからダウンロードしてきたとあるタンパク質の一部です。
タンパク質本体から、1) GlcNAc-GlcNAc-Manの3つの糖、2) アセチルグルコサミンが伸びています。
PDBファイルから該当箇所を一部抜き出します。
1)
ATOM 2751 C5 NAG F 2 122.975 143.965 148.767 1.00 75.59 C
ATOM 2752 C6 NAG F 2 122.358 143.644 150.123 1.00 75.59 C
ATOM 2753 C7 NAG F 2 122.707 145.556 143.723 1.00 75.59 C
ATOM 2754 C8 NAG F 2 123.484 146.265 142.655 1.00 75.59 C
ATOM 2755 N2 NAG F 2 123.396 145.238 144.813 1.00 75.59 N
ATOM 2756 O3 NAG F 2 123.503 147.164 146.913 1.00 75.59 O
ATOM 2757 O4 NAG F 2 123.848 146.064 149.510 1.00 75.59 O
ATOM 2758 O5 NAG F 2 122.224 143.325 147.739 1.00 75.59 O
ATOM 2759 O6 NAG F 2 121.040 144.201 150.188 1.00 75.59 O
ATOM 2760 O7 NAG F 2 121.525 145.287 143.597 1.00 75.59 O
ATOM 2761 C1 BMA F 3 123.097 147.020 150.274 1.00 97.22 C
ATOM 2762 C2 BMA F 3 123.656 147.051 151.689 1.00 97.22 C
ATOM 2763 C3 BMA F 3 122.850 148.034 152.525 1.00 97.22 C
ATOM 2772 C1 NAG A 404 88.555 146.492 149.278 1.00 57.75 C
ATOM 2773 C2 NAG A 404 87.820 145.854 150.451 1.00 57.75 C
ATOM 2774 C3 NAG A 404 88.459 146.305 151.753 1.00 57.75 C
ATOM 2775 C4 NAG A 404 89.933 145.940 151.735 1.00 57.75 C
ATOM 2776 C5 NAG A 404 90.594 146.574 150.520 1.00 57.75 C
ATOM 2777 C6 NAG A 404 92.060 146.177 150.432 1.00 57.75 C
ATOM 2778 C7 NAG A 404 85.470 145.354 150.143 1.00 57.75 C
ATOM 2779 C8 NAG A 404 84.107 145.932 149.903 1.00 57.75 C
ATOM 2780 N2 NAG A 404 86.421 146.229 150.450 1.00 57.75 N
ATOM 2781 O3 NAG A 404 87.819 145.642 152.846 1.00 57.75 O
ATOM 2782 O4 NAG A 404 90.553 146.423 152.931 1.00 57.75 O
ATOM 2783 O5 NAG A 404 89.936 146.142 149.332 1.00 57.75 O
ATOM 2784 O6 NAG A 404 92.145 144.807 150.026 1.00 57.75 O
ATOM 2785 O7 NAG A 404 85.691 144.158 150.064 1.00 57.75 O
これをdoGlycansで処理します。
B=P.74/HD22:-(ND2,C1,<Aa>)-4YB
A=P.143/HD22:-(ND2,C1,<Af>)-4YB-(O4,C1,<Ff>)-4YB-(O4,C1,<a>)-0MA
Aa=ND2,C1,[OD1 CG ND2 C1 250.0 CG ND2 C1 C2 180]
Af=ND2,C1,[OD1 CG ND2 C1 250.0 CG ND2 C1 C2 150]
Ff= O4,C1,[C3 C4 O4 C1 120.0 C4 O4 C1 C2 180.0]
a = O4,C1,[C3 C4 O4 C1 120.0 C4 O4 C1 C2 120.0]
二面角は適当です。すると該当糖鎖はこうつきました。
異常に近づいていますね。もうちょっとマシな配置にしたいですが、doglycansのインプットの2面角はこれだけでも8か所あります。2面角の決め方、動き方を完全に理解していても毎回元のPDBファイルから該当箇所の2面角を測定してインプットファイルに反映するのは面倒です。
ていうかdoglycansの実行中にDihedralとtorsionがこれくらいでRotationしました、ていう出力があるのですが、ここの入力の数値を変えても反映されてる気配がなくて。。。
もとのデータベースから落としてきたファイルには一応CNO原子の位置は正しく記録されているので、いっそこの座標を拾ってきて上書きするというのはどうかと思いつきます。そこで、doglycansで出力されたPDBファイルを見てみます。
ATOM 5540 O6 4YB A 337 124.110 143.027 148.551 1.00 0.00 O
ATOM 5541 H6O 4YB A 337 123.868 142.597 147.726 1.00 0.00 H
ATOM 5542 C4 4YB A 337 127.556 144.075 147.289 1.00 0.00 C
ATOM 5543 H4 4YB A 337 127.221 144.915 146.677 1.00 0.00 H
ATOM 5544 C3 4YB A 337 128.672 143.306 146.543 1.00 0.00 C
ATOM 5545 H3 4YB A 337 129.063 142.526 147.200 1.00 0.00 H
ATOM 5546 O3 4YB A 337 129.761 144.238 146.255 1.00 0.00 O
ATOM 5547 H3O 4YB A 337 130.532 143.769 145.918 1.00 0.00 H
ATOM 5548 C2 4YB A 337 128.117 142.636 145.254 1.00 0.00 C
ATOM 5549 H2 4YB A 337 127.779 143.428 144.581 1.00 0.00 H
ATOM 5550 N2 4YB A 337 129.165 141.828 144.584 1.00 0.00 N
ATOM 5551 H2N 4YB A 337 129.369 140.920 144.966 1.00 0.00 H
ATOM 5552 C2N 4YB A 337 129.695 142.154 143.401 1.00 0.00 C
ATOM 5553 O2N 4YB A 337 129.434 143.212 142.840 1.00 0.00 O
ATOM 5554 CME 4YB A 337 130.694 141.147 142.796 1.00 0.00 C
ATOM 5555 H3M 4YB A 337 130.375 140.122 142.994 1.00 0.00 H
ATOM 5556 H2M 4YB A 337 130.759 141.288 141.716 1.00 0.00 H
ATOM 5557 H1M 4YB A 337 131.680 141.308 143.234 1.00 0.00 H
ATOM 5558 O4 4YB A 337 128.099 144.585 148.526 1.00 0.00 O
ATOM 5559 C1 0MA A 338 127.996 146.001 148.539 1.00 0.00 C
ATOM 5560 H1 0MA A 338 128.972 146.480 148.616 1.00 0.00 H
ATOM 5561 C2 0MA A 338 127.119 146.408 149.757 1.00 0.00 C
ATOM 5562 H2 0MA A 338 127.517 145.943 150.661 1.00 0.00 H
ATOM 5563 C3 0MA A 338 125.658 145.949 149.552 1.00 0.00 C
ATOM 5564 H3 0MA A 338 125.637 144.857 149.531 1.00 0.00 H
水素原子が増えたうえ、CNO原子も順番と名前が変わってしまっているため、単純にPDBファイルの対応行をとってきて座標のカラムを置き換えるだけとはいきません。水素原子の座標を推定しつつ、きちんと化学結合をたどって化学的に等価な原子なら座標を置換するという手順を踏まないといけません。
そこでRDKitをフル活用してやってみます。
一応、DoGlycansで正しく原子が付き元のデータベースのPDBの原子を再構築できるような、普通の糖鎖であることを前提としています。
大まかな流れ
- RCSBからダウンロードしたpdbファイルを処理していき、DoGlycansで糖を付けたタンパク質のpdbファイルを得ます。
- もとのPDBファイルと生成したPDBファイル、あと比較するべき残基の対応を書いたファイルを作っておきます。
- もとのPDBファイルには水素がついていませんので、3次元構造を考慮して水素を付けます。
- 対応する残基＋αの原子を検索、それだけの分子を作ります。
- 分子構造の比較により、元のPDBファイルと生成したPDBファイルで対応する原子ペアを得ます。
- その原子ペア情報を使って、生成したPDBファイルの座標をもとのPDBファイルからとってきます。
- 結果を書き出し。
詳細
インターフェース
- CNOなど水素以外の原子の座標が正しく記録されたPDBファイルを用意する。これを
orig.pdbとする。 - DoGlycanで処理してできたPDBファイルを
ProtG.pdbとする。このファイル中の座標をorig.pdbのデータで上書きしたい。 -
orig.pdbとProtG.pdbでどの残基同士が対応するか、また一続きの糖鎖はどれかを指定したファイルを用意、residue_list.txtとします。
前述のように原子数も順番も名前も違うため、「３次元座標で差がある原子・残基を検出して自動的にフィットする残基を検出」ができないので、参照用のPDBファイル、修正したいPDBファイルの他に修正対象の残基名リストを渡してやる必要があります。
residule_list.txtは以下のようなフォーマットにします。
335 A 404 g1
336 F 1 g2
337 F 2 g2
338 F 3 g2
これは、
-
ProtG.pdbの 335番目の残基がorig.pdbの ChainAの404番目の残基に相当して１つのかたまりg1であること、 -
ProtG.pdbの 336～338番目の残基がorig.pdbのChainFの1～3番目の残基に相当してまた別の１つの塊g2であること
を書いています。
この準備をしたのち、
python adjust_residue.py ProtG.pdb orig.pdb residue_list.txt
として実行するようなpythonコードを書いていきます。
main部
引数を3つ受け取るスクリプトファイルです。 main 関数の書き方をしています。
あと必要なライブラリをインポートします。
import sys
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdFMCS # 分子構造から原子対応をとるためのモジュール
def main():
if len(sys.argv) != 4:
print(f"usage: {sys.argv[0]} protG.PDB orig.PDB residue_list")
return
"""
(処理を書くところ)
"""
if __name__ == "__main__":
main()
以下、「処理を書くところ」を埋めていきます。
ファイル読み込み
コマンドライン引数で渡したファイルを読み込みます。
まず、 sys.argv[1] のpdbファイルは座標を修正したい ProtG.pdb です。水素原子をそのままにしておきたいですし、残基同士が重なっているため結合情報も自動で検出してもらっては困ります。なのでいろいろなオプションを False にします。
一方、 sys.argv[2] のpdbファイルは座標の参考になる orig.pdb です。水素原子はありませんので、読み込み後 AddHs を適用します。そして参考用の座標データを取得します。
mod_mol = Chem.MolFromPDBFile(sys.argv[1], sanitize=False, removeHs=False, proximityBonding=False)
ref_mol = Chem.MolFromPDBFile(sys.argv[2])
ref_mol = Chem.AddHs(ref_mol, addCoords=True, addResidueInfo=True)
ref_coord = np.array(ref_mol.GetConformer().GetPositions())
続いて、残基の対応を記録したファイルを読み込みます。
res_list_dict = {}
with open(sys.argv[3], "r") as f:
for L in f.readlines():
cols = L.strip().split()
g = cols[3]
if not g in res_list_dict.keys():
res_list_dict[g] = []
res_list_dict[g].append((cols[0], cols[1], cols[2]))
例として、このようなデータが得られます。
{'g1': [('335', 'A', '404')],
'g2': [('336', 'F', '1'), ('337', 'F', '2'), ('338', 'F', '3')]}
処理全体の枠組み
糖鎖グループごとに処理をします。各糖鎖ごとに、orig.pdb と ProtG.pdb からRDkitのMolオブジェクトとして取り出し、分子構造のマッチング、 ProtG.pdb の座標を更新します。すべての糖鎖で処理ができれば、 modified.pdb に書き出します。まとめると、main関数は以下のような構造となります。
def main():
"""（先ほどの読み込み部）"""
for _, res_list in res_list_dict.items():
""" orig.pdb から糖鎖スニペット ref_glycomol の作成 """
""" ProtG.pdb から糖鎖スニペット mod_glycomol の作成 """
""" 分子構造マッチング """
""" mod_mol の座標の上書き """
Chem.MolToPDBFile(mod_mol, "modified.pdb")
これからループ内の各部分を深堀します。
元となる ref_mol から糖鎖スニペットの作成
まず原子インデックスを格納する配列と糖鎖部分だけの分子オブジェクトを用意します。
ref_atom_idx = []
ref_glycomol = Chem.RWMol()
ref_mol の該当する原子を分子構造マッチング用分子 ref_glycomol に追加しつつ、その ref_mol での AtmIdx をref_atom_idx に順番を保って追加します。これで、分子構造で対応をとりつつ、元の ref_mol でのインデックスを得て座標を得ることができます。
ref_mol での原子をひとつづつ確認していきます。
指定したChain ID、残基番号でなければスキップします。
for a in ref_mol.GetAtoms():
resinfo = a.GetPDBResidueInfo()
if not resinfo.GetChainId() in ref_chains:
continue
if not f"{resinfo.GetResidueNumber()}" in ref_resids:
continue
この if に該当しなければ対象の原子です。 ref_glycomol に追加する用の原子オブジェクトを生成し、追加します。その時同時に ref_atom_idx にも追加し、 ref_glycomol での AtmIdx をインデックスとして渡せば ref_mol での AtmIdx が得られるようにします。
new_atom = Chem.Atom(a.GetAtomicNum())
new_atom.SetMonomerInfo(a.GetMonomerInfo())
ref_glycomol.AddAtom(new_atom)
ref_atom_idx.append(a.GetIdx())
続いてref_glycomol に結合を追加します。いま処理している原子 a が関わる結合をループして舐めていきますが、その時に結合先原子もすでに ref_glycomol にある結合だけ処理するようにします。こうすれば同じ結合を2回追加することはないと思います。
for bond in a.GetBonds():
a1, a2 = bond.GetBeginAtom(), bond.GetEndAtom()
if (not a2.GetIdx() in ref_atom_idx) or (not a1.GetIdx() in ref_atom_idx):
continue
a1_newidx = ref_atom_idx.index(a1.GetIdx())
a2_newidx = ref_atom_idx.index(a2.GetIdx())
ref_glycomol.AddBond(a1_newidx, a2_newidx, bond.GetBondType())
これで ref_mol の原子チェックのループが終わりました。これを抜けたあと、 ref_glycomol の後処理をします。
分子スニペットの拡張
いろいろやって分かったのですが、doglycansの処理により、糖鎖の残基だけでなく、糖鎖が結合したアミノ酸残基の ND2 とさらにそれが付いたカルボニル基も回転しています。つまり、
(アミノ酸)-CO-NH-(糖鎖)
となったときのこの CO-NH も対象に加えないといけません。ということで、元の分子を見て結合3つ分拡張してやります。拡張する関数 expand_snippet を以下のように定義します。つまり、分子スニペットの原子を見て、結合先原子がまだスニペット（正確には index list）にいなければ追加するという処理を3回やります。
ただindex listにそのまま追加するとループ中に今しがた追加した原子について結合を調べ追加するという処理をしてしまい、1回のループで ref_mol の全体が入ってしまうので、ループはコピーしてから行います。あと、なんかアミノ酸中の C=O 結合がなんか orig.pdb だと単結合、 ProtG.pdb だと2重結合ということでもう単結合固定で追加していきます。
def expand_snippet(molobj, snippetobj, index_list, N):
for _ in range(N):
for a1_idx, ai in enumerate(index_list.copy()):
a = molobj.GetAtomWithIdx(ai)
for na in a.GetNeighbors():
if na.GetIdx() in index_list:
continue
nai = na.GetIdx()
new_atom = Chem.Atom(na.GetAtomicNum())
new_atom.SetMonomerInfo(na.GetMonomerInfo())
snippetobj.AddAtom(new_atom)
index_list.append(nai)
a2_idx = len(index_list) - 1
snippetobj.AddBond(a1_idx, a2_idx, Chem.BondType.SINGLE)
この関数を使って結合3つ分拡張します。続いて編集可能な RWMol のオブジェクトだった ref_glycomol をラスタライズします。
expand_snippet(ref_mol, ref_glycomol, ref_atom_idx, 3)
ref_glycomol = ref_glycomol.GetMol()
これで ref_mol から糖鎖部分を抜き出したオブジェクトを作れました。
この節は2重ループ（正確には3重ループ）中のコードをパートごとに分けて書いていますので、全体の流れが分かりにくいと思います。この記事の最後にプログラム全体をあげていますので、参考にしてください。
mod_mol から糖鎖スニペットの作成
ぶっちゃけ ref_mol と ref_glycomol と ref_atom_idx をそれぞれ mod_mol と mod_glycomol と mod_atom_idx に変えたらさっきとほとんど同じですので、省略します。
ChainIDによるチェックがないくらいでしょうか。
詳しくは最後のプログラム全体を参照してください。
あと、座標修正用に、 Conformer オブジェクトを取得しておきます。このオブジェクトの関数を使って mod_mol の原子座標を変更できます。
mod_coord = mod_mol.GetConformer()
分子構造マッチングと座標の上書き
マッチングに使うのは FindMCS 関数です。これは Mol オブジェクトを複数受け取り、共通部分の構造を得ることができます。
続いて、その構造をSMART表記で受け取り、SMART構造で検索できるように query オブジェクトとします。
これで元の mod_glycomol と ref_glycomol を検索すると当たり前にマッチしますが、そのとき原子の対応が SMARTを通して得られます。
mcs = rdFMCS.FindMCS([mod_glycomol, ref_glycomol])
query = Chem.MolFromSmarts(mcs.smartsString)
match1 = mod_glycomol.GetSubstructMatch(query)
match2 = ref_glycomol.GetSubstructMatch(query)
match* はどうもSMARTに対応する AtmIdx の配列のようです。ということは、この match* での数字は糖鎖スニペット中でのインデックスですので、 mod_atom_idx ref_atom_idx リストを使って元のmod_mol と ref_mol での AtmIdx を得ます。それを使って座標を上書きします。
for i in range(len(match1)):
atm1_idx = match1[i]
atm2_idx = match2[i]
mod_atom = mod_atom_idx[atm1_idx]
ref_atom = ref_atom_idx[atm2_idx]
ref_atom_coord = ref_coord[ref_atom, :]
mod_coord.SetAtomPosition(mod_atom, ref_atom_coord)
GetSubstructMatch を実行するの、元の mod_mol と ref_mol でもいいかな・・・でも同じ配列の糖鎖が複数あったりしたら破綻するし、検索に時間かかったりするか・・・
ともあれ、これで一通りやりたいことの処理ができました。
結果
以上の内容をまとめて実装した pythonコードを実行しました。
最初に出した、原子同士が重なってしまっていた糖鎖の部分も、こんな感じで参考PDBファイルに従った位置に持ってくることができました。
左側が orig.pdb の糖鎖の部分、右側が今回のプログラムで修正した pdb の該当部分です。
コード全体
import sys
import numpy as np
from rdkit import Chem
from rdkit.Chem import AllChem
from rdkit.Chem import rdFMCS
def expand_snippet(molobj, snippetobj, index_list, N):
for _ in range(N):
for a1_idx, ai in enumerate(index_list.copy()):
a = molobj.GetAtomWithIdx(ai)
for na in a.GetNeighbors():
if na.GetIdx() in index_list:
continue
nai = na.GetIdx()
new_atom = Chem.Atom(na.GetAtomicNum())
new_atom.SetMonomerInfo(na.GetMonomerInfo())
snippetobj.AddAtom(new_atom)
index_list.append(nai)
a2_idx = len(index_list) - 1
# snippetobj.AddBond(a1_idx, a2_idx, molobj.GetBondBetweenAtoms(ai, nai).GetBondType())
snippetobj.AddBond(a1_idx, a2_idx, Chem.BondType.SINGLE)
def main():
"""（ファイル読み込み部）"""
if len(sys.argv) != 4:
print(f"usage: {sys.argv[0]} protG.PDB orig.PDB residue_list")
return
mod_mol = Chem.MolFromPDBFile(sys.argv[1], sanitize=False, removeHs=False, proximityBonding=False)
ref_mol = Chem.MolFromPDBFile(sys.argv[2])
ref_mol = Chem.AddHs(ref_mol, addCoords=True, addResidueInfo=True)
# RDkitのバージョンが新しいとAddHsは以下のようにします：
# AddHsParam = Chem.AddHsParameters()
# AddHsParam.addCoords = True
# AddHsParam.addResidueInfo = True
# ref_mol = Chem.AddHs(ref_mol, AddHsParam)
ref_coord = np.array(ref_mol.GetConformer().GetPositions())
res_list_dict = {}
with open(sys.argv[3], "r") as f:
for L in f.readlines():
cols = L.strip().split()
g = cols[3]
if not g in res_list_dict.keys():
res_list_dict[g] = []
res_list_dict[g].append((cols[0], cols[1], cols[2]))
for _, res_list in res_list_dict.items():
mod_resids = [r[0] for r in res_list]
ref_chains = [r[1] for r in res_list]
ref_resids = [r[2] for r in res_list]
""" orig.pdb から糖鎖スニペット ref_glycomol の作成 """
ref_atom_idx = []
ref_glycomol = Chem.RWMol()
for a in ref_mol.GetAtoms():
resinfo = a.GetPDBResidueInfo()
if not resinfo.GetChainId() in ref_chains:
continue
if not f"{resinfo.GetResidueNumber()}" in ref_resids:
continue
new_atom = Chem.Atom(a.GetAtomicNum())
new_atom.SetMonomerInfo(a.GetMonomerInfo())
ref_glycomol.AddAtom(new_atom)
ref_atom_idx.append(a.GetIdx())
for bond in a.GetBonds():
a1, a2 = bond.GetBeginAtom(), bond.GetEndAtom()
if (not a2.GetIdx() in ref_atom_idx) or (not a1.GetIdx() in ref_atom_idx):
continue
a1_newidx = ref_atom_idx.index(a1.GetIdx())
a2_newidx = ref_atom_idx.index(a2.GetIdx())
ref_glycomol.AddBond(a1_newidx, a2_newidx, bond.GetBondType())
expand_snippet(ref_mol, ref_glycomol, ref_atom_idx, 3)
ref_glycomol = ref_glycomol.GetMol()
""" ProtG.pdb から糖鎖スニペット mod_glycomol の作成 """
mod_atom_idx = []
mod_first_atom = None
mod_glycomol = Chem.RWMol()
for a in mod_mol.GetAtoms():
resinfo = a.GetPDBResidueInfo()
if not f"{resinfo.GetResidueNumber()}" in mod_resids:
continue
new_atom = Chem.Atom(a.GetAtomicNum())
new_atom.SetMonomerInfo(a.GetMonomerInfo())
mod_glycomol.AddAtom(new_atom)
mod_atom_idx.append(a.GetIdx())
for bond in a.GetBonds():
a1, a2 = bond.GetBeginAtom(), bond.GetEndAtom()
if a1.GetIdx() > a2.GetIdx():
continue
if (not a2.GetIdx() in mod_atom_idx) or (not a1.GetIdx() in mod_atom_idx):
continue
a1_newidx = mod_atom_idx.index(a1.GetIdx())
a2_newidx = mod_atom_idx.index(a2.GetIdx())
mod_glycomol.AddBond(a1_newidx, a2_newidx, bond.GetBondType())
expand_snippet(mod_mol, mod_glycomol, mod_atom_idx, 3)
mod_glycomol = mod_glycomol.GetMol()
mod_coord = mod_mol.GetConformer()
""" 分子構造マッチング """
mcs = rdFMCS.FindMCS([mod_glycomol, ref_glycomol])
query = Chem.MolFromSmarts(mcs.smartsString)
match1 = mod_glycomol.GetSubstructMatch(query)
match2 = ref_glycomol.GetSubstructMatch(query)
""" mod_mol の座標の上書き """
for i in range(len(match1)):
atm1_idx = match1[i]
atm2_idx = match2[i]
mod_atom = mod_atom_idx[atm1_idx]
ref_atom = ref_atom_idx[atm2_idx]
ref_atom_coord = ref_coord[ref_atom, :]
mod_coord.SetAtomPosition(mod_atom, ref_atom_coord)
Chem.MolToPDBFile(mod_mol, "modified.pdb")
if __name__ == "__main__":
main()
さいごに
これでやりたいこと（糖鎖も含めて絡み合ったタンパク質複合体）ができるかも？
Discussion