使用python一键绘制 (DCCM) Dynamical Cross-Correlation Matrix

casea · 发表于 Post on 2022-9-2 11:25:08

本帖最后由 casea 于 2022-12-5 14:25 编辑

****2022/09/05 更新****

经过优化和调整后，257个原子200帧算了1min，同时出图和Bio3D的一致
目前计算速度卡在了轨迹读取(耗时20 s)以及位移矢量的点积(耗时28 s)
轨迹读取的代码是我之前结合libxtc和pandas写的，速度太慢，之后可能会重写
numpy数组的点积目前没有更好的办法进行优化，欢迎提出建议

**********************
前言：Dynamical Cross-Correlation Matrix，即DCCM，通过DCCM可以看出在模拟期间不同残基之间的共同进化关系。
目前，绘制DCCM可以通过gromacs计算covar然后使用脚本生成DCCM 脚本地址，与此同时，也有很多其他的工具可以绘制出DCCM，如：MD-TASK，Bio3d
相对而言，Bio3D绘制DCCM图是最为广泛以及出图最精致的选择。Bio3D是R语言程序包，gromacs跑出来的轨迹需要先使用VMD等工具转化为DCD格式，在传入到Bio3D中进行计算。因为个人比较懒所以就想用python写一个脚本一键生成DCCM图。

附件即为脚本

DCCM.py (8.03 KB, 下载次数 Times of downloads: 146)

使用方法：
python DCCM.py md.gro md.xtc

注意事项：
xtc轨迹最好不要超过200帧
我测试的体系：257个原子200帧算了3min
DCCM.py默认选择CA原子进行计算，如果想计算其他的原子可以打开文件在ressele处修改

所需包
pandas
matplotlib
numpy
sys
mdanalysis

效果图：

脚本流程：

读入gro和xtc文件
通过选择语法选择所有的CA原子，即CA原子代表了整个残基
通过点云算法将所有帧的结构叠合到第一帧
计算交叉相关系数
附：交叉相关系数计算公式引用链接：

202209021113029603..png (15.34 KB, 下载次数 Times of downloads: 21)

下载附件 Download

2022-9-2 11:13 上传 Uploaded

乐平 · 发表于 Post on 2022-9-4 10:35:43

赞！

xtc轨迹最好不要超过200帧
我测试的体系：257个原子200帧算了3min

粗略看了一眼脚本，跑得慢应该是里面的 for 循环嵌套 太多了，还有很大的提升空间。
既然用了 numpy，应该是能够用高级切片功能，以及 axis= 选项来省掉过多的 for 循环的。

casea · 发表于 Post on 2022-9-4 10:53:29

乐平发表于 2022-9-4 10:35
赞！

谢谢老师的指导，当时做的时候只是想复现一下原图和了解一下原理。再加上自己也是python的初学者，对于numpy的一些用法还不是很熟。我会根据老师的建议尝试优化一下。再次感谢

乐平 · 发表于 Post on 2022-9-11 20:27:38

本帖最后由乐平于 2022-9-12 11:09 编辑

如果你能把测试的文件上传，并讲解一下程序的思路，我也许能再优化一下。
另外，计算交叉相关系数也许可以用 scipy.signal.fftconvolve() 函数加快运算速度

youzi96 · 发表于 Post on 2022-12-5 11:13:11

老师好！我用了您的DCCM的代码，进行计算，但是出现了一些问题，麻烦您有空的时候帮我解答一下吧，谢谢！

casea · 发表于 Post on 2022-12-5 11:45:58

youzi96 发表于 2022-12-5 11:13
老师好！我用了您的DCCM的代码，进行计算，但是出现了一些问题，麻烦您有空的时候帮我解答一下吧，谢谢！

这应该是warning，我当时写的时候只考虑复现，并没有优化性能，因此对于大轨迹可能就会出现warning。你可以先用100帧计算试试有没有出现结果。如果出现了，就可以忽略。

youzi96 · 发表于 Post on 2022-12-5 13:44:47

casea 发表于 2022-12-5 11:45
这应该是warning，我当时写的时候只考虑复现，并没有优化性能，因此对于大轨迹可能就会出现warning。你可 ...

老师好！谢谢解答。我轨迹文件缩小到了101帧，还是出现了如上的问题，麻烦您再帮我看一下，谢谢！

youzi96 · 发表于 Post on 2022-12-5 13:45:12

youzi96 · 发表于 Post on 2022-12-5 15:41:55

非常感谢老师，100帧的出图成功了，但是X坐标轴叠起来了，请问要怎么修改呀？

casea · 发表于 Post on 2022-12-5 16:32:44

youzi96 发表于 2022-12-5 15:41
非常感谢老师，100帧的出图成功了，但是X坐标轴叠起来了，请问要怎么修改呀？

我用的matplotlib绘图。你在py文件中找到
ax.xaxis.set_major_locator(matplotlib.ticker.MultipleLocator(25))
ax.yaxis.set_major_locator(matplotlib.ticker.MultipleLocator(25))
其中25即代表间隔，自己修改为合适地间隔就行

wsyfromwx · 发表于 Post on 2025-3-20 11:01:26

请问有没有amber的py文件

喵星大佬 · 发表于 Post on 2025-3-21 01:32:50

这玩意的意义和算协方差矩阵区别也不大吧

student0618 · 发表于 Post on 2025-3-21 12:39:38

wsyfromwx 发表于 2025-3-20 11:01
请问有没有amber的py文件

Amber 的cpptraj可以直接算DCCM哦

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[GROMACS] 使用python一键绘制 (DCCM) Dynamical Cross-Correlation Matrix

评分 Rate

浏览过的版块