MDAnalysis通道孔径与HOLE操作

HOLE在蛋白通道中搜索能容纳的球形探针，给出沿通道方向的可用半径。本例采用上游公开的短杆菌肽1GRM结构，在鸢算科技集群实际运行；结果是静态几何孔径，不直接等于离子通透率、电导或自由能。

1 编译并运行#

下载HOLE孔径、结构与多构象示例，保存为channel-pore.zip：

python3 -m zipfile -e channel-pore.zip .
cd mdanalysis-channel-pore
bash setup.sh
bash run.sh
bash verify.sh
bash plot.sh

环境为MDAnalysis 2.9.0、mdahole2 0.5.0和HOLE 2.3.1 Apache版。首次setup.sh下载固定源码，在本例的tools/hole2中用GNU Fortran 13.3.0编译；使用make -j1，编译记录在tools/compile.log。所有依赖均留在本例目录，不修改共享安装。数值在results/FULL-FIELDS.npz，图片在plots/。

inputs/1grm_single.pdb有272个原子，本例选择全部原子，保留FOR、ETA等端基。不要盲目用protein选择后就假定保留了所有端基。该结构的通道沿y方向，本例明确使用CVECT 0 1 0；它不是所有蛋白通道的通用方向。

2 完整终端输入#

运行上述工作流后，results/single-coarse/simple2.rad已经生成。也可以直接运行原生HOLE：

cat > pore.inp <<'EOF'
COORD inputs/1grm_single.pdb
RADIUS results/single-coarse/simple2.rad
SPHPDB pore.sph
CPOINT -0.01784376 -0.01225000 4.21743749
CVECT 0 1 0
SAMPLE 0.2
ENDRAD 5.0
RASEED 335
IGNORE SOL WAT TIP HOH K   NA  CL 
SHORTO 0
EOF
tools/hole2/bin/hole < pore.inp > pore.out

CPOINT是本结构CA原子坐标的中心，作为通道内部起点；换结构时需要重新选择。SAMPLE=0.2是相邻搜索平面的间隔，单位Å，和MD时间步无关。ENDRAD=5.0表示搜索到较宽的通道出口后结束，不能把这个值当作孔的实际半径。RASEED固定随机搜索种子，便于比较结果。

HOLE会在已采样球心之间进一步计算中点的净空，因此0.2 Å搜索平面间隔的文本剖面中可以出现0.1 Å间隔；这些中点不是增加了一倍独立优化或新的MD帧。

3 提取剖面与最窄位置#

source env.sh
"$YS_PY" - <<'PYCODE'
from mdahole2.analysis.utils import collect_hole
import numpy as np
p = collect_hole('pore.out')[0]
i = p.radius.argmin()
print('minimum radius:', p.radius[i], 'angstrom')
print('channel y:', p.rxn_coord[i], 'angstrom')
np.savetxt('pore-profile.dat', np.c_[p.rxn_coord, p.radius, p.cen_line_D],
           header='channel_y_A radius_A centerline_distance_A')
PYCODE

本次0.2 Å步长、种子335的最小半径为1.19780 Å，对应直径约2.39560 Å。第三列是沿弯曲孔中心线的有符号距离，初始找到的中心点取零；它与y坐标不同，不应误当成半径。

4 检查步长、随机种子和坐标平移#

配套工作流实际比较0.2/0.1 Å步长、335/336种子，以及整个结构平移后的结果：

这些控制说明本例的数值对所选小改动稳定，并不保证任意复杂分叉通道都能找到全局最优路径。换体系时还要核对原子半径表、被忽略的水或离子、起点和方向。SPH文件中编号−888的球用于表面封口，不属于真实通道剖面，不能用它们扩大孔长度。

参考：HOLE Apache源码及安装说明、mdahole2源码与接口。