如何用Python寻找马尔可夫状态模型中4个亚稳态对应的结构

前提准备 #

你需要提前备好构建MSM时用到的轨迹原子坐标数据、拓扑文件，以及每个轨迹帧对应的MSM状态标签（聚类结果）。下面以常用的pyemma（MSM建模）和mdtraj（轨迹处理）库为例，给出具体实现步骤：

1. 提取亚稳态对应的轨迹帧 #

首先加载训练好的MSM模型，筛选出每个亚稳态（编号假设为0-3）对应的所有轨迹帧索引：

import pyemma
import numpy as np

# 加载已训练完成的MSM模型
msm = pyemma.msm.load('your_trained_msm.pkl')
# 获取所有轨迹的状态序列列表（每个元素对应一条轨迹的状态标签）
state_labels = msm.dtrajs

# 定义目标亚稳态编号
metastable_states = [0, 1, 2, 3]
# 存储每个亚稳态对应的（轨迹索引，帧索引）
state_frames = {}

for state in metastable_states:
    frames = []
    for traj_idx, traj_states in enumerate(state_labels):
        # 找到当前轨迹中属于该亚稳态的所有帧索引
        idx = np.where(traj_states == state)[0]
        if len(idx) > 0:
            frames.append((traj_idx, idx))
    state_frames[state] = frames

2. 提取亚稳态对应的结构 #

可以选择提取平均结构或中心代表性构象（该状态下出现频率最高、RMSD最小的构象）：

提取平均结构

import mdtraj as md

# 加载原始轨迹和拓扑文件
traj_files = ['traj1.xtc', 'traj2.xtc']  # 替换为你的轨迹文件路径
top_file = 'protein_topology.pdb'  # 替换为你的拓扑文件路径
trajs = [md.load(traj, top=top_file) for traj in traj_files]

for state in metastable_states:
    coords_list = []
    # 收集该亚稳态所有帧的原子坐标
    for traj_idx, frame_idx in state_frames[state]:
        coords_list.append(trajs[traj_idx].xyz[frame_idx])
    coords = np.concatenate(coords_list, axis=0)
    # 计算平均坐标
    avg_coords = np.mean(coords, axis=0)
    # 生成平均结构并保存
    avg_traj = md.Trajectory(avg_coords[np.newaxis, :, :], trajs[0].topology)
    avg_traj.save(f'metastable_state_{state}_average.pdb')

提取中心代表性构象

from sklearn.cluster import KMeans

for state in metastable_states:
    coords_list = []
    for traj_idx, frame_idx in state_frames[state]:
        coords_list.append(trajs[traj_idx].xyz[frame_idx])
    coords = np.concatenate(coords_list, axis=0).reshape(-1, coords_list[0].size)
    
    # 用KMeans聚类取中心构象
    kmeans = KMeans(n_clusters=1, random_state=42).fit(coords)
    center_pos = np.argmin(kmeans.transform(coords))
    
    # 找到中心构象对应的轨迹和帧索引
    traj_idx, frame_idx = None, None
    count = 0
    for t_idx, f_idx_list in state_frames[state]:
        if count + len(f_idx_list) > center_pos:
            traj_idx = t_idx
            frame_idx = f_idx_list[center_pos - count]
            break
        count += len(f_idx_list)
    
    # 保存中心构象
    center_traj = trajs[traj_idx][frame_idx]
    center_traj.save(f'metastable_state_{state}_center.pdb')

3. 验证结构一致性 #

计算亚稳态内部的RMSD分布，确认该状态下结构的聚集程度：

for state in metastable_states:
    coords_list = []
    for traj_idx, frame_idx in state_frames[state]:
        coords_list.append(trajs[traj_idx].xyz[frame_idx])
    coords = np.concatenate(coords_list, axis=0)
    avg_coords = np.mean(coords, axis=0)
    # 计算所有帧相对于平均结构的RMSD
    rmsd_values = md.rmsd(md.Trajectory(coords, trajs[0].topology), 
                         md.Trajectory(avg_coords[np.newaxis], trajs[0].topology))
    print(f"亚稳态{state} - RMSD均值: {np.mean(rmsd_values):.2f} nm，标准差: {np.std(rmsd_values):.2f} nm")

内容的提问来源于stack exchange，提问作者Giovanni Wences