计算化学公社

标题: 使用python一键绘制 (DCCM) Dynamical Cross-Correlation Matrix [打印本页]

作者
Author: casea 时间: 2022-9-2 11:25
标题: 使用python一键绘制 (DCCM) Dynamical Cross-Correlation Matrix
本帖最后由 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图。

附件即为脚本
(, 下载次数 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

效果图：
(, 下载次数 Times of downloads: 28)
脚本流程：

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

上传 Uploaded
点击下载
Click to download

作者
Author: 乐平 时间: 2022-9-4 10:35
赞！

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

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

作者
Author: casea 时间: 2022-9-4 10:53

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

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

作者
Author: 乐平 时间: 2022-9-11 20:27
本帖最后由乐平于 2022-9-12 11:09 编辑

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

作者
Author: youzi96 时间: 2022-12-5 11:13
老师好！我用了您的DCCM的代码，进行计算，但是出现了一些问题，麻烦您有空的时候帮我解答一下吧，谢谢！ (, 下载次数 Times of downloads: 18)

作者
Author: casea 时间: 2022-12-5 11:45

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

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

作者
Author: youzi96 时间: 2022-12-5 13:44

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

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

作者
Author: youzi96 时间: 2022-12-5 13:45
(, 下载次数 Times of downloads: 17)

作者
Author: youzi96 时间: 2022-12-5 15:41
非常感谢老师，100帧的出图成功了，但是X坐标轴叠起来了，请问要怎么修改呀？ (, 下载次数 Times of downloads: 19)

作者
Author: casea 时间: 2022-12-5 16:32

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即代表间隔，自己修改为合适地间隔就行

作者
Author: wsyfromwx 时间: 2025-3-20 11:01
请问有没有amber的py文件

作者
Author: 喵星大佬 时间: 2025-3-21 01:32
这玩意的意义和算协方差矩阵区别也不大吧

作者
Author: student0618 时间: 2025-3-21 12:39

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

Amber 的cpptraj可以直接算DCCM哦

欢迎光临计算化学公社 (http://bbs.keinsci.com/)