使用CyFES对配体运动轨迹进行数据透视

技术背景

如果我们有一个蛋白质X和一个配体Y，那么可以对这个X+Y的体系跑一段长时间的分子动力学模拟，以观测这个体系在不同结合位点下的稳定性。类似于前面一篇博客中计算等高面的方法，我们可以计算轨迹的KDE函数，然后保存成Cube格式（高斯中用于保存电子轨道的一种格式）的文件。然后就可以通过可视化软件如VMD等，来加载蛋白和cube文件进行分析。

软件介绍与安装

CyFES是一款基于Python开发用户层，基于Cython开发链路层以及CUDA C++开发物理层Kernel函数的一个高性能开源FES计算软件。

目前功能还比较简单，仅支持三维分子运动轨迹的数据透视。安装支持源码和pip安装，pip安装方法如下：

$ python3 -m pip install cyfes --user --upgrade -i https://pypi.org/simple

之所以建议采取--user的策略，是因为很多环境中没有权限在site-packages路径下创建文件，会导致动态链接库缺失的问题，--user可以一定程度上规避这个问题。下载的时候国内的镜像源同步会有一定的延迟，如果需要及时的下载最新版本的CyFES，可以使用-i https://pypi.org/simple参数配置。安装成功后，可以使用如下指令确认CyFES是否安装成功：

$ python3 -m cyfes --help
python3 -m cyfes --help
usage: __main__.py [-h] [-i I] [-ic IC] [-ib IB] [-s S] [-e E] [-g G] [-o O]
                   [-no_bias NO_BIAS] [-f32 F32] [-sigma SIGMA]

optional arguments:
  -h, --help        show this help message and exit
  -i I              Set the input record file path.
  -ic IC            Set the cv index of input record file. Default: 0,1,2
  -ib IB            Set the bias index of input record file. Default: 3
  -s S              CV length. Default: None
  -e E              Edge length. Default: 1.0
  -g G              Grid numbers. Default: 10,10,10
  -o O              Set the output FES file path.
  -no_bias NO_BIAS  Do not use the bias from input file. Default: false
  -f32 F32          Use float32. Default: false
  -sigma SIGMA      Sigma value when calculating FES. Default: 0.3

当然，也可以从CyFES-Gitee主页下载源码进行源码安装。安装完成后，可以执行一个python脚本确认一下动态链接库是否缺失的问题：

# check_dynamics.py
import os
import site
from pathlib import Path

site_path = Path(site.getsitepackages()[0])
site_file_path = site_path.parent.parent.parent / 'cyfes' / 'libcufes.so'
site_dynamics_path = str(site_file_path)

user_site_path = Path(site.USER_SITE)
user_file_path = user_site_path.parent.parent.parent / 'cyfes' / 'libcufes.so'
user_dynamics_path = str(user_file_path)

if not os.path.exists(site_dynamics_path) and not os.path.exists(user_dynamics_path):
    print ('Check dynamics complete, no libcufes.so file founded!')
else:
    print ('Installation of CyFES success!')

确认安装成功之后，我们可以开始使用CyFES进行轨迹数据透视。

轨迹输入

CyFES支持读取这样的一个轨迹文件xyz_bias.txt：

23.5578 33.8817 37.8341 0.000000
23.4752 33.7842 37.8489 0.882319
23.4557 33.7728 37.8236 1.544485
23.4979 33.6253 37.8524 1.952011
23.5502 33.6256 37.8981 2.140049
23.6389 33.6791 37.9141 1.437173
...

其中前三列表示轨迹的x、y、z坐标，最后一列表示对应坐标位置的偏置势bias，用于计算权重。在CyFES-2.6之后的版本中，可以支持不给定bias，那么默认就全部都是0。

CyFES终端指令

最简单的场景，我们可以使用这样一行简单的代码执行CyFES的计算：

$ python3 -m cyfes -i xyz_bias.txt -o z.cub

这样就可以根据轨迹文件生成一个可以用于可视化的cube格式文件。如果需要更多的一些配置，常用的有：

$ python3 -m cyfes -i xyz_bias.txt -o z.cub -e 5.0

表示边缘增加5A（这里注意单位是埃，而最终保存的cube格式文件会转为波尔bohr长度）的空隙。

$ python3 -m cyfes -i xyz_bias.txt -o z.cub -g 20,30,40

表示x、y、z方向的格点数量分别为20，30和40个格点。

$ python3 -m cyfes -i xyz_bias.txt -o z.cub -sigma 0.2

表示禁带宽度band width设置为0.2。屏幕打印输出的样式大致是这样的：

$ python3 -m cyfes -i xyz_bias.txt -e 5.0 -g 10,20,30 -sigma 0.1 -o z_x.cub
2024-08-19 09:34:26,304 [CyFES] Start to initialize parameters
2024-08-19 09:34:26,318 [CyFES] CV (1000, 3)
2024-08-19 09:34:26,318 [CyFES] Bias (1000,)
2024-08-19 09:34:26,319 [CyFES] Origin crd [16.0653 28.2584 32.4803]
2024-08-19 09:34:26,319 [CyFES] Final crd [29.4832 41.7948 47.3887]
2024-08-19 09:34:26,320 [CyFES] Grids (6000, 3)
2024-08-19 09:34:26,320 [CyFES] BandWidth [0.1 0.1 0.1]
2024-08-19 09:34:26,320 [CyFES] Start to calculate FES
2024-08-19 09:34:27,285 [CyFES] Writting FES into file /home/cy-fes/tests/z_x.cub
2024-08-19 09:34:27,306 [CyFES] Task complete :)

这里主要是记录一些参数和运行日志，最终的计算结果会按照Cube文件的格式输出到z_x.cub文件中，文件格式可以查阅下一个章节的介绍。

Cube输出

简单的说明一下我们这里生成的cube文件格式内容，大体如下：

Generated by CyFES
Total	1000000	grids
1	30.3589	53.4004	61.3785
100	0.256121	0	0
100	0	0.258383	0
100	0	0	0.284572
1	1.000000	43.1649	66.3195	75.6072
34.6304	34.1533	33.6897	33.2396	32.8028	32.3795	
31.9696	31.5731	31.19	30.8203	30.464	30.121	
29.7914	29.4751	29.1721	28.8823	28.6057	28.3422	
28.0918	27.8543	27.6297	27.4177	27.2182	27.031	
...

第一行是标题，第二行是格点数声明，第三行是原子数和原点坐标，第四行到第六行是x、y、z方向的格点数和偏移矢量。第七行是原子信息，这里我们直接在轨迹盒子的中心位置放了一个氢原子，第二列是核电荷数，后面三列是坐标。第八列开始是每一个格点的势能数值，每一行最多6个格点数据。