用 ProLIF 提取相互作用指纹：把 Pose 变成特征

对接分数是一个数字，相互作用指纹则告诉你「这个姿势到底形成了哪些相互作用」。用它筛选对接结果，比只看分数可靠得多——能复现共晶关键相互作用的姿势，才值得采纳。

安装与基本用法

pip install prolif

import prolif as plf
from rdkit import Chem
import MDAnalysis as mda

# 受体
u = mda.Universe("receptor_fixed.pdb")
protein = plf.Molecule.from_mda(u.select_atoms("protein"))

# 配体姿势（对接输出的 SDF）
poses = list(plf.sdf_supplier("docked_poses.sdf"))

fp = plf.Fingerprint([
    "HBAcceptor", "HBDonor", "Hydrophobic",
    "PiStacking", "Anionic", "Cationic", "CationPi",
    "XBDonor", "MetalDonor",
])
fp.run_from_iterable(poses, protein)

df = fp.to_dataframe()
print(df.shape)          # (姿势数, 相互作用数)
print(df.columns[:10])   # MultiIndex: (ligand, residue, interaction)

用途一：按参考指纹筛选对接姿势

这是最有价值的用法。先算出共晶结构的相互作用指纹作为参考，再看哪些对接姿势能复现它：

import numpy as np

# 1) 参考：共晶配体的指纹
ref_pose = list(plf.sdf_supplier("ref_ligand.sdf"))
fp_ref = plf.Fingerprint(fp.interactions.keys())
fp_ref.run_from_iterable(ref_pose, protein)
ref_vec = fp_ref.to_dataframe().iloc[0]

# 2) 对齐列（对接姿势可能缺少某些相互作用）
all_cols = ref_vec.index.union(df.columns)
df_al  = df.reindex(columns=all_cols, fill_value=False)
ref_al = ref_vec.reindex(all_cols, fill_value=False)

# 3) 计算与参考的 Tanimoto 相似度
def tanimoto(row, ref):
    inter = (row & ref).sum()
    union = (row | ref).sum()
    return inter / union if union else 0.0

df_al["similarity_to_ref"] = df_al.apply(lambda r: tanimoto(r, ref_al), axis=1)
print(df_al["similarity_to_ref"].describe())

# 4) 只保留能复现主要相互作用的姿势
good = df_al[df_al["similarity_to_ref"] > 0.5]
print(f"{len(good)}/{len(df_al)} 个姿势复现了参考的主要相互作用")

用途二：检查关键相互作用是否保留

# 激酶体系的例子：铰链区氢键是活性的必要条件
KEY_INTERACTIONS = [
    ("MET793", "HBAcceptor"),    # EGFR 铰链
    ("MET793", "HBDonor"),
    ("LYS745", "Anionic"),       # 保守赖氨酸盐桥
]

def has_key_interactions(row, keys, require=1):
    hits = 0
    for res, itype in keys:
        cols = [c for c in row.index
                if len(c) == 3 and res in str(c[1]) and c[2] == itype]
        if any(row[c] for c in cols):
            hits += 1
    return hits >= require

mask = df.apply(lambda r: has_key_interactions(r, KEY_INTERACTIONS), axis=1)
print(f"{mask.sum()} 个姿势形成了铰链区相互作用")

# 结合对接分数一起筛
# 分数好 + 复现关键相互作用 = 真正值得关注的姿势

这一步能挡掉大量「分数很高但结合模式完全不对」的假阳性。打分函数会被大分子、高疏水性分子欺骗，而相互作用核对是独立的检验维度。

用途三：指纹作为机器学习特征

# 把相互作用指纹拼进特征矩阵
import pandas as pd

X_ifp = df.astype(int).values                    # 相互作用指纹
X_ecfp = compute_ecfp(smiles_list)               # 分子指纹（见 330）
X_combined = np.hstack([X_ecfp, X_ifp])

# 相互作用指纹提供了「结构互补性」信息，
# 而 ECFP 只描述配体本身 —— 两者互补

from sklearn.ensemble import RandomForestClassifier
clf = RandomForestClassifier(n_estimators=500, n_jobs=-1)
clf.fit(X_combined, y_active)

# 特征重要性可解释：哪些相互作用与活性相关
ifp_importance = clf.feature_importances_[X_ecfp.shape[1]:]
top_idx = np.argsort(ifp_importance)[-10:]
for i in top_idx:
    print(df.columns[i], round(ifp_importance[i], 4))

用途四：分析 MD 轨迹中相互作用的稳定性

u = mda.Universe("topology.pdb", "traj.dcd")
lig = u.select_atoms("resname MOL")
prot = u.select_atoms("protein")

fp_md = plf.Fingerprint(["HBAcceptor", "HBDonor", "Hydrophobic", "PiStacking"])
fp_md.run(u.trajectory[::5], lig, prot)

df_md = fp_md.to_dataframe()
occupancy = df_md.mean() * 100                  # 各相互作用的占有率（%）

stable = occupancy[occupancy > 70].sort_values(ascending=False)
transient = occupancy[(occupancy > 30) & (occupancy <= 70)]
print("稳定相互作用（>70%）：")
print(stable)
print(f"\n间歇性接触（30~70%）：{len(transient)} 个")
print(f"偶发接触（<30%）：{(occupancy <= 30).sum()} 个 —— 不应写进结论")

# 可视化
fp_md.plot_barcode()          # 相互作用随时间的条码图

占有率是把 MD 轨迹转成可读结论最直接的方式：>70% 是稳定的结构特征，<30% 属于偶发接触，不该写进结合模式描述。

注意事项

• 依赖输入结构的质量：缺氢会严重影响氢键判定，分析前务必用 PDBFixer 补氢（见 340《用 PDBFixer 修复蛋白结构》）。
• 几何判据可调：默认阈值适用于多数情况，金属体系或非常规相互作用可能需要调整。
• PDBQT 丢失键级：对接输出要用 mk_export.py 转回 SDF 再分析（见 341《用 Meeko 准备 Docking 文件》），否则芳香性、键级判断会出错。
• 只描述不评价：指纹告诉你存在什么相互作用，不告诉你结合有多强。

延伸资源

• 工具详解：207《ProLIF》；对照：208《PLIP》、339《用 PLIP 生成蛋白-配体相互作用报告》；
• 上一步：336《用 AutoDock Vina 跑 Docking》、337《用 GNINA 做 Rescoring》。