[GROMACS]模拟数据分析前轨迹文件生成-轨迹预处理

前言在gromacs动力学中原始轨迹文件.xtc/.trr因周期性边界条件PBC和整体平动/转动的存在无法直接用于定量分析或可视化。本文介绍三类核心预处理方法分别针对连续轨迹、结构叠合和可视化需求并明确每类输出文件的适用分析场景与注意事项。为了得到一个用于可视化和适配于各种数据分析使用的轨迹我们需要对轨迹进行三种预处理用于不同的数据分析和可视化。预处理类型输出文件示例核心目标典型应用PBC跳跃消除md_nojump.xtc还原真实连续运动保留绝对坐标MSD、RDF、能量rerun整体运动叠合md_fit.xtc固定蛋白质消除平动/转动RMSD、RMSF、氢键、DCCM可视化优化md_center.xtc画面整洁分子完整居中轨迹动画展示注意尽量保留原始轨迹每次处理生成新文件不要覆盖保证预处理逻辑更清晰一、PBC 跳跃消除获得连续轨迹md_nojump.xtc周期性边界条件 (PBC)在分子动力学模拟中为了模拟一个无限大的体系我们使用一个有限的盒子并让这个盒子在三个维度上无限重复。当一个分子比如你的小分子配体从盒子的一侧穿出去时它会立刻从盒子的另一侧穿进来。这在计算上是连续的但在轨迹文件中记录的坐标会发生一个巨大的跳跃从盒子的一端到另一端。1.1命令与参数详解gmx trjconv-smd.tpr-fmd.xtc-nindex.ndx-omd_nojump.xtc-pbcnojump参数含义详细说明-s md.tpr结构拓扑/输入文件包含盒子信息、原子质量、键连接等。-s通常指.tpr文件提供参考坐标和盒子。-f md.xtc输入轨迹原始轨迹可以是.xtc压缩或.trr全精度。-n index.ndx索引文件可选推荐。定义原子组如Protein、Ligand、Backbone用于后续选择参考组。若没有gmx trjconv会交互式提示你从预定义组中选择。-o md_nojump.xtc输出轨迹处理后的连续轨迹。-pbc nojump核心参数“No Jump”检测并消除分子因PBC产生的跳跃。算法对于每个分子追踪其质心轨迹当质心跨过盒子边界时将该分子的所有原子坐标平移回连续位置。-pbc nojump的底层原理假设盒子长度为 10 nm一个分子在第 10 帧坐标x 9.9 nm第 11 帧从右侧穿出实际物理位置应为x 10.1 nm但由于PBC轨迹记录为x 0.1 nm。-pbc nojump会检测到跳跃 半个盒长默认阈值 0.45×盒长自动给该分子所有原子加上±L使其连续。最终第 11 帧被修正为10.1 nm。如图上图为消除跳跃之前的会发现超出PBC盒子以外的部分发生跳跃从盒子对侧出现而当我们消除之后如下图你会发现即便超出盒子区域仍然不会发生跳跃重要-pbc nojump依赖分子完整性。若分子本身在轨迹中被“切开”例如配体的一部分在盒子一边另一部分在对面则需要先用-pbc whole将分子拼接完整。1.2 前置步骤确保分子完整-pbc whole什么时候必须用-pbc whole当模拟中分子跨过了盒子边界且轨迹记录时没有保持分子完整性默认GROMACS输出的轨迹中原子坐标按照其在盒子中的真实位置记录一个分子可能被“劈”成两半分别位于盒子两端。这种情况常见于长链分子如聚合物、某些配体膜蛋白中的脂质分子使用-pbc nojump前没有做任何处理如图为被劈成两半的氨基酸分子。解决方法先执行一次-pbc wholegmx trjconv-smd.tpr-fmd.xtc-nindex.ndx-omd_whole.xtc-pbcwhole-pbc whole将每个分子的所有原子拼回同一盒子内通过搜索每个原子的周期性映像使得分子内原子间距离最小化即分子不断裂。之后再对md_whole.xtc执行-pbc nojump。一条命令同时处理 whole nojump管道避免中间文件gmx trjconv-smd.tpr-fmd.xtc-nindex.ndx-omd_nojump.xtc-pbcwhole-pbcnojump注意多个-pbc参数按顺序执行先 whole 后 nojump。1.4 PBC跳跃消除轨迹-适用分析分析是否可用备注与必要条件均方位移MSD必须需要绝对坐标的连续性。使用gmx msd时必须输入-pbc nojump处理后的轨迹。径向分布函数RDF必须执行gmx rdf时必须加上-pbc no因为轨迹已经消除了跳跃无需再搜索周期性映像。能量 rerun如相互作用能可用前提轨迹中分子完整已用-pbc whole。若跳过 whole能量计算可能错误原子位置错误。回旋半径Rg完美无需叠合连续轨迹即可。氢键寿命可用依赖连续时间序列但不依赖绝对位置。自相关函数ACF可用同上。扩散系数Einstein relation必须MSD 的子项。PCA主成分分析不可PCA 需要先去除整体平动/转动否则前几个主成分为整体运动。RMSD / RMSF不可除非后续用gmx rms的-fit选项叠合但最好单独做叠合轨迹。DCCM / 动力学相关性不可需要去除整体运动。可视化VMD/PyMOL勉强分子可能不在盒子中央且可能仍有残留跳跃如果-pbc nojump未能完全处理某些情况。二、整体平动与转动消除结构叠合轨迹md_fit.xtc此处是为了消除动力学模拟中发生的2.1 命令与参数详解gmx trjconv-smd.tpr-fmd_nojump.xtc-nindex.ndx-omd_fit.xtc-fitrottrans-center参数拆解参数含义详细说明-fit rottrans叠合方式rot旋转 trans平动。将每一帧的指定参考组通常是蛋白质骨架通过最小二乘拟合叠合到第一帧默认参考帧。-center居中将指定组通常是蛋白质平移到盒子中心。通常与-fit配合使用。注意-fit和-center都需要交互选择组Least-squares fit group选择叠合参考组如Backbone或C-alpha。这个组的原子坐标将被用来计算旋转和平动矩阵。Output group输出哪个组通常选择System所有原子这样整个体系包括配体、水、离子都随蛋白质一起旋转平移。Centering group如果用了-center选择居中的组一般也是Protein。先做-pbc nojump再做叠合原因如果直接用原始轨迹做叠合由于 PBC 跳跃第 10 帧蛋白质在盒子中央第 11 帧蛋白质可能“跳”到盒子角落叠合算法会错误地将整体平移很大距离导致配体位置错乱。必须先消除跳跃使轨迹在绝对坐标下连续再做叠合。2.2 叠合的底层原理最小二乘拟合每一帧的原子坐标集合 {r_i}参考帧 {r_i^ref}。寻找旋转矩阵 R 和平移向量 t使得下式最小∑ w_i |R·r_i t - r_i^ref|^2其中 w_i 通常是原子质量默认或等权重。GROMACS 使用 Kabsch 算法。只对 fit group 计算 R,t然后应用到整个 output group。2.4 适用分析分析是否可用备注RMSD均方根偏差完美直接gmx rms -s md.tpr -f md_fit.xtcRMSF涨落完美先做 RMSD 拟合轨迹或用gmx rmsf -res氢键数量gmx hbond完美依赖相对距离和角度不受整体运动影响盐桥gmx saltbr完美同上溶剂可及表面积gmx sasa完美依赖构象与整体运动无关动态交叉关联矩阵DCCM必须DCCM 要求去除整体运动否则虚假相关自由能形貌图FEL完美基于 PCA 或 RMSD 等内部坐标二级结构DSSP完美gmx do_dssp基于局部几何平均结构gmx covargmx anaeig完美必须先叠合轨迹主成分分析PCA必须必须去除整体运动MM/PBSA 结合自由能谨慎若使用gmx mmgbsa需要连续轨迹且无整体运动但需注意MM/PBSA 通常使用无溶剂轨迹仅蛋白配体且需要确保坐标与拓扑一致。建议对md_fit.xtc提取蛋白配体子集。MSD不可绝对坐标已被改变RDF不可相对位置虽保留但盒子信息可能被破坏-center平移后盒子不变GROMACS 中盒子矢量不变仅原子坐标平移RDF 理论上仍正确实际上-center不改变相对距离但很多用户担心所以建议 RDF 用md_nojump.xtc。能量 rerun危险坐标可能超出原始盒子边界虽然盒子矢量没变但平移后部分原子可能跑到盒子外PME 计算可能出错。2.6 错误问题错误1叠合时选错了输出组。如果输出组只选了蛋白质配体就没有被旋转平移配体会相对于蛋白质“飞走”。必须选System。错误2参考组包含了配体或柔性 loop。配体运动会导致叠合扭曲RMSD 虚高。错误3没有先做-pbc nojump直接叠合。蛋白质的跳跃会导致叠合失败轨迹中出现异常跳跃。错误4使用-fit rottrans但忘记-center。虽然没有原则错误但蛋白质可能不在盒子中央可视化不美观。不过对 RMSD 等无影响。三、可视化优化最优轨迹md_center.xtc3.1 命令与参数详解gmx trjconv-smd.tpr-fmd.xtc-nindex.ndx-omd_center.xtc-pbcmol-center-urcompact参数详解参数含义详细说明-pbc mol分子完整化将每个分子的所有原子通过周期性平移放到同一盒子内类似-pbc whole但不改变质心位置只确保分子内原子不断裂。-center居中将指定组通常是蛋白质平移到盒子中心。-ur compact使用紧凑盒子表示重新定义盒子矢量使得盒子刚好包裹所有原子最小化空白区域。这改变了盒子的形状和大小因此会改变原子间的相对距离注意-ur compact会修改盒子尺寸从而改变分子间的距离例如如果盒子原来是立方体 10×10×10 nm但所有原子聚集在中央一个 3×3×3 区域-ur compact会将盒子缩小到刚好包裹所有原子可能变成 3.2×3.2×3.2 nm。这将导致周期性的距离计算错误因为盒子变小周期性映像更近对于多分子体系分子间的相对距离跨盒子被破坏因此此轨迹最好不要用于任何定量分析3.2 可视化需要的操作原始轨迹中蛋白质可能不在中心且水分子散落各处看起来杂乱。-pbc mol消除了断裂的分子如某水分子一半在左一半在右。-center让蛋白质居中便于观察。-ur compact移除大片空白区域视频文件更小渲染更快。3.3 适用场景仅限视觉用途用途说明制作模拟动画VMD/PyMOL画面整洁聚焦蛋白质-配体。提取首尾帧做构象对比图例如论文中的 Figure 1模拟前后叠合图。粗略观察配体结合模式变化仅定性判断不可计算距离。3.4 验证与注意事项验证在 VMD 中加载打开“Periodic”显示如果 VMD 能显示盒子你会看到盒子缩小了。检查运行gmx check -f md_center.xtc查看盒子尺寸应该远小于原始盒子如果体系浓缩。重要如果你想在可视化中保留正确的相对距离例如配体与蛋白质的距离看起来不变形不要使用-ur compact只使用-pbc mol -centergmx trjconv-smd.tpr-fmd.xtc-nindex.ndx-omd_visual_noCompact.xtc-pbcmol-center这样盒子保持不变仅居中距离信息保留。但注意即使如此这个轨迹仍不适合定量分析因为没有去除整体运动蛋白质可能在转动。参考文献GROMACS 官方手册gmx trjconv章节https://manual.gromacs.org/current/onlinehelp/gmx-trjconv.html“Handling periodic boundary conditions” – GROMACS Tutorials by Justin Lemkul“Best practices for analysis of MD trajectories” – BioExcel CoE个人博客地址https://blog.huimy.top/44/578.html本文遵循 CC BY-NC 4.0 许可欢迎转载但请注明出处。