研究メモ。
参考
CPPTRAJ Combined Cluster Analysis Example Tutorial
ambermd.org/doc12/Amber14.pdf
クラスタリングとは、ざっくり言うと構造のグループ分けのこと。
それではAmber付属のCPPTRAJでクラスタリングする例の日本語解説をしてみる。
ステップ１：トポロジーとトラジェクトリを読み込む
parm rGACC.nowat.parm7
trajin rGACC.MREMD1.nowat.nc.40
処理するデータを軽くしたければイオンとかを消す（この例では最初から水を消してある。）
strip :Na+ outprefix noions
ステップ２：クラスタリングのコマンド
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分の１から初めて後で追加する作戦。 - 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ではクラスタリングアルゴリズムは２種類（hieraggloとdbscan）しかないようだ。（Amber16ではk-meansが使える）
仕方がない。dbscanでやるか。hieraggloは階層的クラスタリング法なので。
メトリックはCα-RMSDで
rms :1-128@CA
最終的なcpptrajインプットはこんな感じ。
parm xxx.prmtop
trajin xxx.nccluster 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 netcdfrun
*sieve 10 randomはフレーム数をあらかじめ間引いていたので消した。フレーム数が少なすぎるとクラスタリングできない。
DBSCANについては
（追記）
Amber17のCPPTRAJを使えるようにしたのでk-meansのクラスタリングをやってみた。
インプットファイルはこんな感じ。
# for cpptraj17
parm xxx.prmtop
trajin xxx.nccluster 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 netcdfrun
これでn=5のk-meansクラスタリングができる。
