MDAnalysis核酸骨架、糖环与翻转几何

核酸构象不只由碱基间距离描述。α、β、γ、δ、ε、ζ和χ反映骨架与碱基相对取向；糖环另有伪旋转描述，碱基翻转也需要特定几何定义。本文逐一调用原生函数，并给出真实结构中的数值。

1 下载并运行#

下载RNA与膜叶分类的完整输入、分析脚本、数值和结构图，保存为rna-membrane.zip：

module load pymatgen/2026.9.24
python -m zipfile -e rna-membrane.zip .
cd mdanalysis-rna-membrane
bash run.sh
source env.sh

输出在results/、plots/，请使用新的解压目录。环境是MDAnalysis 2.9.0和NetworkX 3.6.1；setup.sh在示例目录建立NetworkX的局部加载入口，使用集群已有依赖，不修改公共Python环境。env.sh供当前终端中的Python示例使用；没有调用qvasp或VASP。

数据来自MDAnalysisTests 2.9.0官方公开样例，来源、哈希及上游许可随包保留。RNA是PDB 1K5I派生数据中RNAA段的23个核苷酸、739个原子；膜是单个DPPC/胆固醇结构快照。这里在集群上实际执行几何和图分析，没有把这些上游快照说成新做的分子动力学。程序中的90°刚体旋转只检查几何量是否保持不变，不是第二个动力学时间点。

2 α至χ七种角#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis import nuclinfo
u = mda.Universe('inputs/rna.pdb')
print('alpha beta gamma delta epsilon zeta chi:')
print(nuclinfo.tors(u, 'RNAA', 5))
for name in ('alpha','beta','gamma','delta','eps','zeta','chi'):
    print(name, getattr(nuclinfo, 'tors_'+name)(u, 'RNAA', 5))
PYCODE

第5残基实际结果约为 302.2038、179.0431、35.2714、79.4997、201.0004、282.1432、210.7093°。组合tors与七个单项函数给出同样定义；脚本对第2–22残基的147个角全部保留。

需要相邻残基的角不能在缺失端点时硬补零。自己的拓扑若采用不同的糖原子命名，应先改选择，特别注意原子名里的撇号。角度跨0/360°时要用圆周差，不能把359°到1°算成变化358°。

3 糖环伪旋转与羟基角#

python - <<'PYCODE'
import MDAnalysis as mda
from MDAnalysis.analysis import nuclinfo
u = mda.Universe('inputs/rna.pdb')
print('CP phase:', nuclinfo.phase_cp(u, 'RNAA', 9))
print('AS phase:', nuclinfo.phase_as(u, 'RNAA', 11))
print('Hydroxyl:', nuclinfo.hydroxyl(u, 'RNAA', 5))
PYCODE

本例Cremer–Pople相位在第9残基约 3.16497°，Altona–Sundaralingam相位在第11残基约 171.71645°，第5残基羟基几何角约 156.62223°。两种伪旋转角的构造不同，图中分开呈现，不能直接混成一条时间序列。羟基函数在这个命名约定下使用C1′、C2′、O2′、H2′；缺少这些原子时不适用。

4 翻转拟二面角：先补齐真实质量#

该函数用四组原子的质量中心，所以质量来源也要正确。本包质量由对应原生PSF明确提取，不能只根据PDB名称猜测后就当作相同输入。

python - <<'PYCODE'
import json
import MDAnalysis as mda
from MDAnalysis.analysis import nuclinfo
u = mda.Universe('inputs/rna.pdb')
u.add_TopologyAttr('masses', json.load(open('inputs/rna-metadata.json'))['masses'])
for bp1,bp2,i in ((16,2,3),(8,9,10)):
    value = nuclinfo.pseudo_dihe_baseflip(u, bp1,bp2,i,
                seg1='RNAA',seg2='RNAA',seg3='RNAA')
    print(bp1,bp2,i,value)
PYCODE

两个选区的实际结果分别为 314.69804°和 34.50106°。四个点依次来自参考碱基对的质量中心、bp2糖环、目标i糖环和目标i碱基。bp1/bp2/i必须按研究结构定义，不能把任意三个残基都解释成发生了翻转。

这是同一个快照中的两个选区，不是前后两个翻转事件。若研究翻转动力学，应沿真实连续轨迹提取同一选区的角，结合周期边界、角度连续展开、采样长度和状态定义；单个角不能给出自由能势垒或速率。

5 保存完整数组#

python - <<'PYCODE'
import numpy as np
z = np.load('results/FULL-FIELDS.npz')
print(z['RNA_torsions_deg'][0].shape)
np.savetxt('my-rna-torsions.csv', z['RNA_torsions_deg'][0], delimiter=',',
    header='alpha,beta,gamma,delta,epsilon,zeta,chi', comments='')
print(z['RNA_pseudorotation_CP_AS_deg'][0])
PYCODE

[0]取原始快照；第二行控制构型是程序内的刚体旋转检查，不是动力学采样。更多接口与命名约定见MDAnalysis nuclinfo手册，原子对应距离见上一篇。