いつか博士になる人へ

あるポスドク研究員の思い出

AmberのCPPTRAJでクラスター解析する方法

研究メモ。

 

参考

CPPTRAJ Combined Cluster Analysis Example Tutorial

ambermd.org/doc12/Amber14.pdf

 

クラスタリングとは、ざっくり言うと構造のグループ分けのこと。

 

それではAmber付属のCPPTRAJでクラスタリングする例の日本語解説をしてみる。

 

 

ステップ1:トポロジーとトラジェクトリを読み込む

 

parm rGACC.nowat.parm7
trajin rGACC.MREMD1.nowat.nc.40

 

処理するデータを軽くしたければイオンとかを消す(この例では最初から水を消してある。)

 

strip :Na+ outprefix noions

 

 

ステップ2:クラスタリングのコマンド

 

cluster C0 \
dbscan minpoints 25 epsilon 0.9 sievetoframe \
rms :1@N2,O6,C1',P,:2@H2,N6,C1',P,:3@O2,H5,C1',P,:4@O2,H5,C1',P \
sieve 10 random \
out cnumvtime.dat \
sil Sil \
summary summary.dat \
info info.dat \
cpopvtime cpopvtime.agr normframe \
repout rep repfmt pdb \
singlerepout singlerep.nc singlerepfmt netcdf \
avgout Avg avgfmt restart

 

…ながっ

えーとなになに、

  • cluster C0
    clusterがメインのコマンドで、C0はアウトプットにつける名前(適当でよし)
  • dbscan minpoints 25 epsilon 0.9 sievetoframe
    dbscanって方法でクラスタリング。あとはそのオプション。
  • rms :1@N2,O6,C1',P,:2@H2,N6,C1',P,:3@O2,H5,C1',P,:4@O2,H5,C1',P
    RMSDを基準にクラスタリング。Pって何?Pはリン。rmsの後に書いてある原子を全部使って、任意のフレーム同士をBest fitした後にRMSDを計算。
  • sieve 10 random
    全部のMDスナップショットを使うと時間かかるから、最初10分の1から初めて後で追加する作戦。
  • out cnumvtime.dat
    メインのアウトプットファイル。
  • sil Sil
    まじなにこれ?Amber16以降のオプション。シルエットファイル(謎)の名前を与えている。
  • summary summary.dat
    結果のまとめはsummary.datに書いてちょ。
  • info info.dat
    結果の詳しい情報はinfo.datに書いてやんす
  • cpopvtime cpopvtime.agr normframe
    クラスターのポピュレーションの時系列データをcpopvtime.agrに書いて。frameの数で規格化しといてな
  • repout rep repfmt pdb
    クラスターの代表構造をPDBで書き出す(!)ちなみに代表構造とクラスターの中心構造(centroid)は同じでした。
  • singlerepout singlerep.nc singlerepfmt netcdf
    全代表構造を一つのファイルに書き出す
  • avgout Avg avgfmt restart
    各クラスタの平均構造をリスタートファイル形式で書き出す(Amber16以降のオプション)。

 

なるほど…実行は数分かかると。

チュートリアルにはこの後、二つのトラジェクトリを読み込んでクラスタリングした後、結果を各トラジェクトリごとに書き出す方法が書いてある。

これはクラスタリングの収束性を見るためやな。とりあえず今はいいか。

 

 それでここからは僕の今やりたいことの話。

k-meansでクラスタリングして、各cluster center構造をpdbで書き出したい。

 

が、manualを見たところAmber14のcpptrajではクラスタリングアルゴリズムは2種類(hieraggloとdbscan)しかないようだ。(Amber16ではk-meansが使える)

仕方がない。dbscanでやるか。hieraggloは階層的クラスタリング法なので。

 

メトリックはCα-RMSDで

rms :1-128@CA

 

最終的なcpptrajインプットはこんな感じ。

 

parm xxx.prmtop
trajin xxx.nc

cluster C0 \
dbscan minpoints 5 epsilon 0.9 sievetoframe \
rms :1-128@CA \
out cnumvtime.dat \
summary summary.dat \
info info.dat \
cpopvtime cpopvtime.agr normframe \
repout rep repfmt pdb \
singlerepout singlerep.nc singlerepfmt netcdf

run

 

*sieve 10 randomはフレーム数をあらかじめ間引いていたので消した。フレーム数が少なすぎるとクラスタリングできない。

 

DBSCANについては

DBSCAN - Wikipedia

f:id:kitos:20180615162349p:plain

 

(追記)

Amber17のCPPTRAJを使えるようにしたのでk-meansのクラスタリングをやってみた

 

インプットファイルはこんな感じ。

 

# for cpptraj17
parm xxx.prmtop
trajin xxx.nc

cluster C0 \
kmeans clusters 5 randompoint\
rms :1-128@CA \
out cnumvtime.dat \
summary summary.dat \
info info.dat \
cpopvtime cpopvtime.agr normframe \
repout rep repfmt pdb \
singlerepout singlerep.nc singlerepfmt netcdf

run

 

これでn=5のk-meansクラスタリングができる。