cluster探索
2025-10-31
目的是从10 us 的 CMD 轨迹中 聚类拿到 apo 的优势构象，看看 redock 结果是否会好于 apo 晶体结构
- 75.4 /home/databank/MD_traj_dataset/D3PM/analysis/cluster_com_top1_redock_suceess
轨迹是75.2 上处理好了后（轨迹处理）。同步到 75.4
根据乐乐老师的指导，先把 complex 体系 redock 的最优构象能复现 晶体结构的拿出来做
聚类的标准想用口袋rmsd，apo 的口袋index可用 pkt_rmsd.py 生成
- 75.2 /home/databank/MD_traj_dataset/D3PM/analysis/batch4/pkt_index
发现了问题：pkt_rmsd.py 是基于 complex 构象cpptraj 生成的口袋 index，这个index 很有可能不是apo 的口袋Index。
遂使用叠合（pymol super）后的apo结构，和配体组合，获得index，检查是否正确
其实不需要，先检查 complex 和 apo 的氨基酸序列是否一致，如果一致，直接从 mol-5A.dat 中提取氨基酸index 即可
一般MOL 都在最后或者最前，要注意检查是否是在最前，这样complex 生成的 mol-5A.dat 中氨基酸index 全部 减去1 才是 apo 口袋index
对于 batch4，因为是用 mdframe 进行的体系准备，文件的路径略有差异
以 4JVF_B_pair26 为例，pkt_rmsd.py 需要 sys_md_parameters 里面的 top 文件（gmx.top），但是mdframe 文件位置，名字都不一样
- cp /home/databank/MD_traj_dataset/D3PM/sys_md_parameters/batch4/complex/4JVF_B_pair26/topprep/com-sol-restrain.top /home/databank/MD_traj_dataset/D3PM/sys_md_parameters/batch4/4JVF_B_pair26/gmx.top
蛋白 rmsd 聚类
- 75.3 /home/data/MD_traj_dataset/D3PM/analysis/cluster_com_top1_redock_suceess/backbone_rmsd/4X1O_5BW4_pair15
调用环境
- # amber (cpptraj.OMP)
- # 75.3
- source /home/data/zyzhou/software/amber24/amber_bin/amber.sh
- # 75.2
- source /home/databank/zyzhou/software/amber24/amber_bin/amber.sh
- # 线程数根据实际情况调整
- export OMP_NUM_THREADS=30
提交任务：
对10条处理好了周期性的轨迹进行聚类
特征为蛋白的所有非氢原子
体系特征聚类参数
与师姐讨论后，师姐建议使用体系中关键氨基酸的二面角作为feather进行聚类，同时，将关键氨基酸的二面角统计分布计算后绘制频率分布直方图
cpptraj 计算特定二面角并聚类
下面两个路径下 聚类的参数 基本一致，是两个不同的二面角。
对于这个体系，个人认为 sys11_test2 里面的feather 更适合表征口袋构象特征：这个色氨酸W109的CA-CB-CG-CD2 二面角 决定了 是否与 配体存在碰撞
- 75.2 /home/databank/MD_traj_dataset/D3PM/analysis/cluster_com_top1_redock_suceess/pkt_feather/sys11_test2
- 75.2 /home/databank/MD_traj_dataset/D3PM/analysis/cluster_com_top1_redock_suceess/pkt_feather/4X1O_5BW4_pair15
聚类输入文件 cluster.in：
- parm ../pro.prmtop
- trajin ../../../../batch3/pbc_apo_10exp/4X1O_5BW4_pair15/dry-auto-apo-M.xtc
- dihedral dih_Trp109 :109@CA :109@CB :109@CG :109@CD2 out dih_Trp109_CA-CB-CG-CD2.dat
- cluster \
- hieragglo epsilon 30 clusters 20 averagelinkage \
- data dih_Trp109 \
- sieve 50 random \
- summary summary.dat \
- info info.dat \
- repout rep repfmt pdb \
- avgout Avg avgfmt restart
- run
因为10条轨迹，do.sh 可以帮助for循环都聚类
- for i in {1..10}
- do
- mkdir exp${i}
- sed "s/M/$i/g" cluster.in >exp${i}/cluster_exp${i}.in
- cd exp${i}
- cpptraj.OMP -i cluster_exp${i}.in
- cd ../
- done
gmx 二面角计算探索
- 75.2 /home/databank/MD_traj_dataset/D3PM/analysis/cluster_com_top1_redock_suceess/pkt_feather/gmx_angle_dihedral
mk_angndx 生成 计算 二面角所需的index文件
11-3 没来得及搞清楚其中的细节，dih.ndx里面的数据都不清楚是啥，明天看手册
- gmx_mpi mk_angndx -s /home/databank/MD_traj_dataset/D3PM/1um_cmd_12A/batch3/4X1O_A_pair15/4X1O_A_pair15.tpr -n dih.ndx
计算二面角
- gmx_mpi angle -f /home/databank/MD_traj_dataset/D3PM/analysis/batch3/pbc_apo_10exp/4X1O_5BW4_pair15/dry-auto-apo-1.xtc -n dih.ndx -type dihedral -ov -oc
算是能算，但是没搞懂算的是哪个二面角，具体的参数 -ov -oc 应该也要再摸索