[Hiroki_Yamamura]

JSR株式会社の山村と申します。

CGモデルのdcd→pdb変換について、教えていただけないでしょうか。

genesis-1.7.1を用いて、AICG2+モデルのMD計算(1本のタンパク質のみ)を行っています。計算で得られたdcdファイルをpdb形式へ変換したいのですが、下記インプットファイルをcrd_convertで実行すると、TER文字がmodelの最後だけでなく1残基毎に記述されたpdbファイルが得られます(下記result1)。またcg_convertも使ってみたところ、こちらは全てのmodelにおいて1残基目の座標しか記録されませんでした(下記result2)。

コンバーターの使い方を何か間違っているのかと思いますが、対処法をご教示いただけないでしょうか。お手数おかけしますが、よろしくお願いいたします。

インプット———————————————————
[INPUT]
psffile = protein_cg.psf # protein structure file
reffile = protein_cg.pdb # PDB file

[OUTPUT]
trjfile = output.pdb # trajectory file

[TRAJECTORY]
trjfile1 = pro_md1.dcd # trajectory file
md_step1 = 50000000 # number of MD steps
mdout_period1 = 10000 # MD output period
ana_period1 = 10000 # analysis period
trj_format = DCD # (PDB/DCD)
trj_type = COOR+BOX # (COOR/COOR+BOX)
trj_natom = 0 # (0:uses reference PDB atom count)

[SELECTION]
group1 = all # select all atoms
group2 = an:CA

[FITTING]
fitting_method = NO # [NO,TR,TR+ROT,TR+ZROT,XYTR,XYTR+ZROT]

[OPTION]
check_only = NO
trjout_format = PDB # write PDB
trjout_type = COOR+BOX
trjout_atom = 1 # output all atoms
centering = YES # protein is centered at the origin
centering_atom = 2
center_coord = 0.0 0.0 0.0
pbc_correct = MOLECULE # molecules are wrapped into the unit cell
——————————————————————

result1———————————————————————-
CRYST1 180.000 180.000 180.000 90.00 90.00 90.00 P 1 1
MODEL 1
ATOM 1 CA ALA A 1 -6.262 72.645 -61.534 0.00 0.00 A
TER
ATOM 2 CA GLN A 2 -3.791 75.591 -61.523 0.00 0.00 A
TER
ATOM 3 CA GLY A 3 -3.271 75.270 -65.281 0.00 0.00 A
TER
・・・
—————————————————————————–

result2———————————————————————-
CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1
MODEL 1
ATOM 1 CA ALA A 1 -13.938 34.068 -23.143 0.00 0.00 PROT
TER
ENDMDL
CRYST1 0.000 0.000 0.000 90.00 90.00 90.00 P 1 1
MODEL 2
ATOM 1 CA ALA A 1 -13.921 34.070 -23.148 0.00 0.00 PROT
TER
ENDMDL
・・・
—————————————————————————–

[ckobayashi]

GENESIS開発グループ　小林です。

 

まず、最初にcg_convertは選択された原子グループ内の座標の平均を出力するツールですので、１粒子分の座標が出力されるのは仕様上正しいです。今回の目的には適さないと思います。

 

いただいた情報から “molecule”の設定が正しくできていないように思われます。このツールは分子ファイル(今回はpsf)からbondの情報を読み取り、bondがつながっている範囲を1分子とみなします。

crd_convertの標準出力ログを見ていただいて、[STEP2]の”Define Molecule>  Summary of Molecule”からの数行をご確認していただくと、”num_bonds  = “や”num_molecules = “という項目があるかと思います。こちらの数が、正しく表示されているかご確認いただければと思います。（AICG2+のタンパク質モデルならnum_bondsは残基数-1でしょうし、num_molecules は1ではないかと推測)

もしMoleculeの情報が正しい場合は、psf/pdbファイルを一度提示いただければと思います。

 

 

[Hiroki_Yamamura]

小林様

ご返信いただきありがとうございます。
cg_convertは用途が違うとのこと、承知いたしました。

また、ご指摘いただいたcrd_convertのログ[STEP2]を確認したところ、num_bonds=0、num_molecules=残基数となっておりました。
genesis_cg_tool(aa_2_cg.jl)を使ってインプットファイル群を作っていましたので、psfではなくgrotopを読み込む設定にしたところ、num_bonds=残基数-1、num_molecules=1となり、各modelの最後のみTERが付くpdbファイルを得られました。

今回crd_convertのサンプルページに合わせてpsfを読み込んでいたのですが、aa_2_cg.jlで出力できるpsfは他の解析に使用するためのものでコンバーターに使うものではなかったという理解でいいのでしょうか？
問題自体は解決しましたので勉強のためという形になりますが、ご教示いただけますと幸いです。
よろしくお願いいたします。

[ckobayashi]

ご確認ありがとうございます。

aa_2_cg.jlで作成されるpsfファイルはvmd/pymol等のGENESISのgrotop形式が使えない外部viewer、ツール等で利用するための物です。

crd_convertを始めとしたGENESISの解析ツール群ではMDと同様にgrotopファイルを利用していただくようにお願いします。

[Hiroki_Yamamura]

小林様

psfファイルは外部ソフト用とのこと承知いたしました。今後、注意いたします。
色々とご教示いただき、ありがとうございました。