MDAnalysis中文教程-UserGuide

快速入门指南

本指南旨在作为 MDAnalysis 的基本介绍，以帮助您入门并运行。您可以在示例notebooks中看到更复杂的任务。本页概述了如何：

加载分子动力学结构或轨迹

使用 AtomGroups，这是 MDAnalysis 中的一个中心数据结构

使用轨迹

写出坐标

使用 MDAnalysis 中的分析算法

纠正和自动引用 MDAnalysis 和算法

import MDAnalysis as mda
from MDAnalysis.tests.datafiles import PSF, DCD, GRO, XTC
## 注意MDanalysisTests未与MDAnalysis捆绑下载，需另外安装

import warnings
## suppress some MDAnalysis warnings about PSF files
warnings.filterwarnings('ignore')

from matplotlib import pyplot as plt

print(mda.Universe(PSF, DCD))
print("Using MDAnalysis version", mda.__version__)

%matplotlib inline

正常无报错

概述

MDAnalysis 是一个 Python 包，它提供了访问和分析分子动力学轨迹中数据的工具。几个关键数据结构构成了 MDAnalysis 的支柱。

Atom：分子系统由粒子组成，除原子外，Atom对象也可表示粗粒化的珠子

AtomGroup：Atoms被分组到 AtomGroups 中。这可能是 MDAnalysis 中最重要的类，因为几乎所有内容都可以通过它访问。

Universe：包含分子系统AtomGroup中可通过.atoms属性访问的的所有粒子，并在.trajectory轨迹中将它们组合

MDAnalysis 中的一个基本概念是，在任何时候，只访问轨迹的一个时间范围。Universe的trajectory属性通常是文件读取器。将轨迹视为帧索引的函数，该帧索引仅使来自该特定帧的数据可用。此结构很重要，因为它允许 MDAnalysis 处理太大而无法放入计算机内存的轨迹文件。 MDAnalysis的内部单位是Å和ps。

载入结构或者轨迹

使用 MDAnalysis 通常首先将数据加载到 MDAnalysis 中的中心数据结构Universe中。
首要的参数就是拓扑文件和轨迹文件。

拓扑文件必选参数，接受pdf、pdb、crd、gro格式

可选参数，接受*任意*数量的轨迹文件，支持dcd、xtc、trr、xyz格式以及单帧的上述格式；未提供实参则仅载入拓扑文件的结构

psf = mda.Universe(PSF)
print(psf)
print(hasattr(psf, 'trajectory'))
#<Universe with 3341 atoms>
#False

psf文件不含坐标信息，所以没有trajectory属性

gro = mda.Universe(GRO)
print(gro)
print(len(gro.trajectory))
#<Universe with 47681 atoms>
#1

gro文件含坐标信息，所以有trajectory属性，且长度为1，即单帧
以下读取测试文件、创建一条轨迹用于后续分析。

u = mda.Universe(PSF, DCD)
print(u)
print(len(u.trajectory))
#<Universe with 3341 atoms>
#98

使用原子组

大多数分析都需要创建并操作一个“AtomGroup”，即一组Atom的集合。为了方便起见，您还可以处理具有化学意义的原子组，例如“Residue”或“Segment”。这些组都有与之对应的容器，如“ResidueGroup”和“SegmentGroup”。例如，一个“Universe”对象的“.residues”属性会返回一个“ResidueGroup”

1 2	print(u.residues) #<ResidueGroup [<Residue MET, 1>, <Residue ARG, 2>, <Residue ILE, 3>, ..., <Residue ILE, 212>, <Residue LEU, 213>, <Residue GLY, 214>]>

选择原子

获取Universe的原子的最简单方式是使用.atoms属性

1
2

print(u.atoms)
#<AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 2: HT1 of type 2 of resname MET, resid 1 and segid 4AKE>, <Atom 3: HT2 of type 2 of resname MET, resid 1 and segid 4AKE>, ..., <Atom 3339: C of type 32 of resname GLY, resid 214 and segid 4AKE>, <Atom 3340: OT1 of type 72 of resname GLY, resid 214 and segid 4AKE>, <Atom 3341: OT2 of type 72 of resname GLY, resid 214 and segid 4AKE>]>

返回一个AtomGroup，相当于Atom的list。大多数分析都涉及使用AtomGroup的atoms。AtomGroup也可通过另一个AtomGroup切片创建，如下：

1
2

print(u.atoms[-5:])
#<AtomGroup [<Atom 3337: HA1 of type 6 of resname GLY, resid 214 and segid 4AKE>, <Atom 3338: HA2 of type 6 of resname GLY, resid 214 and segid 4AKE>, <Atom 3339: C of type 32 of resname GLY, resid 214 and segid 4AKE>, <Atom 3340: OT1 of type 72 of resname GLY, resid 214 and segid 4AKE>, <Atom 3341: OT2 of type 72 of resname GLY, resid 214 and segid 4AKE>]>

也可通过索引进行创建

1
2

print(u.atoms[[0, 3, -1, 1, 3, 0]])
## <AtomGroup [<Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>, <Atom 4: HT3 of type 2 of resname MET, resid 1 and segid 4AKE>, <Atom 3341: OT2 of type 72 of resname GLY, resid 214 and segid 4AKE>, <Atom 2: HT1 of type 2 of resname MET, resid 1 and segid 4AKE>, <Atom 4: HT3 of type 2 of resname MET, resid 1 and segid 4AKE>, <Atom 1: N of type 56 of resname MET, resid 1 and segid 4AKE>]>

也可通过类似于现有VMD、PyMol等软件包的原子选择语句.select_atoms()方法进行选择

1
2

print(u.select_atoms('resname ASP or resname GLU'))
#<AtomGroup [<Atom 318: N of type 54 of resname GLU, resid 22 and segid 4AKE>, <Atom 319: HN of type 1 of resname GLU, resid 22 and segid 4AKE>, <Atom 320: CA of type 22 of resname GLU, resid 22 and segid 4AKE>, ..., <Atom 3271: OE2 of type 72 of resname GLU, resid 210 and segid 4AKE>, <Atom 3272: C of type 20 of resname GLU, resid 210 and segid 4AKE>, <Atom 3273: O of type 70 of resname GLU, resid 210 and segid 4AKE>]>

注意：原子选择语句中的数值范围是闭区间，包括’-‘和’:'，不同于list切片的左闭右开

print(u.select_atoms('resid 50-100').n_residues)
print(u.residues[50:100].n_residues)
#51
#50

select_atoms支持布尔运算、通配符

1 2	u.select_atoms("(resname GLU or resname HS*) and name CA and (resid 1:100)") #<AtomGroup with 6 atoms>

注意，选择创建的AtomGroup被排序和去重，而这可能不同于切片产生的AtomGroup
详见AtomGroup篇

从AtomGroup提取原子信息

AtomGroup对象包含许多属性，这些属性可以方便地提取原子信息。这些属性包括：names、masses、residues、segements、resnames等，详细属性见topology_system篇

print(list(zip(u.atoms[:20].names, u.atoms[:20].masses)))
## [('N', 14.007), ('HT1', 1.008), ('HT2', 1.008), ('HT3', 1.008), ('CA', 12.011), ('HA', 1.008), ('CB', 12.011), ('HB1', 1.008), ('HB2', 1.008), ('CG', 12.011), ('HG1', 1.008), ('HG2', 1.008), ('SD', 32.06), ('CE', 12.011), ('HE1', 1.008), ('HE2', 1.008), ('HE3', 1.008), ('C', 12.011), ('O', 15.999), ('N', 14.007)]

print(u.atoms[:20].residues)
print(u.atoms[-20:].segments)
## <ResidueGroup [<Residue MET, 1>, <Residue ARG, 2>]>
## <SegmentGroup [<Segment 4AKE>]>

print(u.atoms[:20].resnames)
## ['MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'MET', 'ARG']
## 注意.residues去重了，而resname没有

通过拓扑属性分组原子

near_met = u.select_atoms('not resname MET and (around 2 resname MET)')
sorted(near_met.groupby(['resnames', 'names']))
#[('ALA', 'C'),
## ('ALA', 'HN'),
## ('ARG', 'N'),
## ('ASN', 'O'),
## ('ASP', 'C'),
## ('ASP', 'N'),
## ('GLN', 'C'),
## ('GLU', 'N'),
## ('ILE', 'C'),
## ('LEU', 'N'),
## ('LYS', 'N'),
## ('THR', 'N')]

AtomGroup位置和方法

.position属性可能是AtomGroup中最有用的属性之一。它返回一个形状为(n_atoms, 3)的NumPy数组，其中n_atoms是AtomGroup中的原子数，3是每个原子的x、y、z坐标。该属性是只读的，但您可以使用它来计算原子组的位置，例如质心

ca = u.select_atoms('resid 1-5 and name CA')
print(ca.positions)
print(ca.positions.shape)
## [[11.664622    8.393473   -8.983231  ]
##  [11.414839    5.4344215  -6.5134845 ]
##  [ 8.959755    5.612923   -3.6132305 ]
##  [ 8.290068    3.075991   -0.79665166]
##  [ 5.011126    3.7638984   1.130355  ]]
## (5, 3)

AtomGroup对象还包含许多方法，这些方法可以方便地计算原子组的位置和运动。这些方法包括：
center_of_mass()：质心
center_of_geometry()：重心
total_mass()：总质量
total_charge()：净电荷
radius_of_gyration()：回转半径
bsphere()：重心和边界半径

MDAnalysis也可以创建键、键角、二面角、离面弯曲角(improper)，但要满足一定条件，如原子数量，这是很自然的。

1
2
3

nhh = u.atoms[:3]
print(nhh.names)
## ['N', 'HT1', 'HT2']

对于角度部分，输入原子的顺序非常关键。二面角是输入原子的∠123，2是顶点原子。

angle_nhh = nhh.angle
print(f"N-H-H angle: {angle_nhh.value():.2f}")
## N-H-H angle: 37.99
hnh = u.atoms[[1, 0, 2]]
angle_hnh = hnh.angle
print(f"H-N-H angle: {angle_hnh.value():.2f}")
## H-N-H angle: 106.20

使用轨迹

Universe的trajectory包含坐标变化信息。轨迹的长度就是帧数，即len(u.trajectory)

获取轨迹每帧的信息的标准方式就是迭代它。当时刻变化时universe仅包含当前时间步的相关信息。
*此处timestep被替代为“时刻”，timestep通常为MD的时间步长、时间间隔参数，如1fs，本文timestep应该表达时刻，如第10帧、第10ps

for ts in u.trajectory[:20]:
    time = u.trajectory.time
    rgyr = u.atoms.radius_of_gyration()
    print(f"Frame: {ts.frame:3d}, Time: {time:4.0f} ps, Rgyr: {rgyr:.4f} A")
## Frame:   0, Time:    1 ps, Rgyr: 16.6690 A
## Frame:   1, Time:    2 ps, Rgyr: 16.6732 A
## Frame:   2, Time:    3 ps, Rgyr: 16.7315 A
## Frame:   3, Time:    4 ps, Rgyr: 16.7223 A
## Frame:   4, Time:    5 ps, Rgyr: 16.7440 A
## Frame:   5, Time:    6 ps, Rgyr: 16.7185 A
## Frame:   6, Time:    7 ps, Rgyr: 16.7741 A
## Frame:   7, Time:    8 ps, Rgyr: 16.7764 A
## Frame:   8, Time:    9 ps, Rgyr: 16.7894 A
## Frame:   9, Time:   10 ps, Rgyr: 16.8289 A
## Frame:  10, Time:   11 ps, Rgyr: 16.8521 A
## Frame:  11, Time:   12 ps, Rgyr: 16.8549 A
## Frame:  12, Time:   13 ps, Rgyr: 16.8723 A
## Frame:  13, Time:   14 ps, Rgyr: 16.9108 A
## Frame:  14, Time:   15 ps, Rgyr: 16.9494 A
## Frame:  15, Time:   16 ps, Rgyr: 16.9810 A
## Frame:  16, Time:   17 ps, Rgyr: 17.0033 A
## Frame:  17, Time:   18 ps, Rgyr: 17.0196 A
## Frame:  18, Time:   19 ps, Rgyr: 17.0784 A
## Frame:  19, Time:   20 ps, Rgyr: 17.1265 A

迭代完后，轨迹的指针重新定位到第一帧，详见trajectories篇。

1 2	print(u.trajectory.frame) ## 0

除了迭代之后，也可通过帧索引获取结构，如u.trajectory[10]，或手动迭代，u.trajectory.next()
通常，轨迹分析首先在列表中收集逐帧数据。
然后可以将其转换为其他数据结构，例如 numpy 数组或 pandas DataFrame。它可以绘制（如下所示），或用于进一步分析。

rgyr = []
time = []
protein = u.select_atoms("protein")
for ts in u.trajectory:
    time.append(u.trajectory.time)
    rgyr.append(protein.radius_of_gyration())
    
import pandas as pd
rgyr_df = pd.DataFrame(rgyr, columns=['Radius of gyration (A)'], index=time)
rgyr_df.index.name = 'Time (ps)'

rgyr_df.head()

rgyr_df.plot(title='Radius of gyration')
plt.show()

动态选择

选择通常是在静态属性上定义的，这些属性在穿过轨迹时不会改变。但是一些选择包含依赖距离的查询，如around、point等，详见selection篇。这种情况下，可以通过设置“updating=True”使用动态选择对选择按时更新，这时select_atoms返回的是UpdatingAtomGroup而不是AtomGroup。

dynamic = u.select_atoms('around 2 resname ALA', updating=True)
print(type(dynamic))
dynamic
## <class 'MDAnalysis.core.groups.UpdatingAtomGroup'>
## <AtomGroup with 54 atoms, with selection 'around 2 resname ALA' on the entire Universe.>
static = u.select_atoms('around 2 resname ALA')
print(type(static))
static
## <class 'MDAnalysis.core.groups.AtomGroup'>
## <AtomGroup with 54 atoms>

#当调用universe的下一帧时，动态选择的原子会被更新，而静态选择的原子保持不变。
u.trajectory.next()
dynamic,static
## (<AtomGroup with 56 atoms, with selection 'around 2 resname ALA' on the entire Universe.>,<AtomGroup with 54 atoms>)

写出坐标

MDAnalysis 支持将数据写入一系列文件格式，包括单帧格式（例如 PDB、GRO）和轨迹写入器（例如 XTC、DCD 和多帧 PDB 文件）。（您还可以将选择写入其他程序，例如 VMD 宏，详见selection_exporters篇。)用户指南有一个完整的列表格式，每个格式都有自己的参考页面。

单帧

写入只能容纳单个帧的文件的最直接方法是使用任意AtomGroup的 write 的方法。MDAnalysis 使用文件扩展名来确定输出文件格式。

1 2	ca = u.select_atoms('name CA') ca.write('ca.gro')

轨迹

写出轨迹的标准方式是：

打开轨迹Writer，识别每帧将包含的原子数；
遍历轨迹，使用Wrier.write()将每帧写入Writer；
关闭Writer。

ca = u.select_atoms('name CA')
with mda.Writer('calphas.xtc', ca.n_atoms) as w:
    for ts in u.trajectory:
        w.write(ca)

分析

MDAnalysis 附带了一组不同的分析模块，以及用于实现您自己的构建块，详见analysis篇。
其中大多数遵循一个通用接口：

使用Universe和相关的参数初始化分析；
使用.run()运行分析。可选参数有起始帧索引、结束帧索引、步长、toggling verbose。默认分析整条轨迹。
分析结果存储在Analysis对象中
另外，函数可用于分析单帧
但是，不是所有分析都使用这个模型。检查每个分析的文档是很重要的。分析详见examples篇。

RMSD

不是所有的子模块都使用“import MDAnalysis”导入。大多数模块都必须显示导入，如“from MDAnalysis.analysis import rms”
MDAnalysis 提供了一个函数rmsd()，用于计算两个 numpy 坐标数组之间的 RMSD。

bb = u.select_atoms('backbone')

u.trajectory[0] ## first frame
first = bb.positions

u.trajectory[-1] #last frame
last = bb.positions

rms.rmsd(first, last)
## 6.852774844656239

类RMSD也提供了对轨迹的操作。

第一个参数是用于RMSD计算的AtomGroup/Universe
第二个参数是参考结构，可以是AtomGroup或Universe，默认为None，即使用第一个参数的当前帧作为参考结构
我们选择将轨迹对齐到骨架原子，然后计算select的骨架原子的RMSD。
然后，不用重新对齐轨迹，计算groupselections的RMSD。

u.trajectory[0] ## set to first frame
rmsd_analysis = rms.RMSD(u, select='backbone', groupselections=['name CA', 'protein'])
rmsd_analysis.run()
## <MDAnalysis.analysis.rms.RMSD at 0x18ee5c2d820>

结果被保存在属性中，即rmsd_analysis.results.rmsd，为shape为(n_frames, 3+len(groupselections))的numpy.ndarray，其中各列分别是

帧数
时间(ps)
RMSD(backbone)
RMSD(CA)
RMSD(protein)
我们可以将结果转换为DataFrame，并绘制出来。

import pandas as pd

rmsd_df = pd.DataFrame(rmsd_analysis.results.rmsd[:, 2:],
                       columns=['Backbone', 'C-alphas', 'Protein'],
                       index=rmsd_analysis.results.rmsd[:, 1])
rmsd_df.index.name = 'Time (ps)'
rmsd_df.head()

rmsd_df.plot(title='RMSD')
plt.show()

参考文献

在发表工作中使用MDAnalysis应当引用下列文献：

Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. J. Comput. Chem. 32 (2011), 2319–2327. doi:10.1002/jcc.21787

R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. MDAnalysis: A Python package for the rapid analysis of molecular dynamics simulations. In S. Benthall and S. Rostrup, editors, Proceedings of the 15th Python in Science Conference, pages 98-105, Austin, TX, 2016. SciPy. doi:10.25080/Majora-629e541a-00e

MDAnalysis包含许多应该单独引用的算法和模块,如MDAnalysis.analysis.rmsMDAnalysis.analysis.align需引用：

Douglas L. Theobald. Rapid calculation of RMSD using a quaternion-based characteristic polynomial. Acta Crystallographica A 61 (2005), 478-480.

Pu Liu, Dmitris K. Agrafiotis, and Douglas L. Theobald. Fast determination of the optimal rotational matrix for macromolecular superpositions. J. Comput. Chem. 31 (2010), 1561–1563.

每个模块的主要来源将会保存在他们的文档中。可以使用duecredit自动生成引用，详见references篇。

翻译https://userguide.mdanalysis.org/stable/examples/quickstart.html