20260910 144229 蛋白质结构预测与虚拟筛选全流程实现详解

靶点蛋白序列/PDB ID
        │
        ▼
┌──────────────────────────┐
│ Step 1: 结构获取           │  技能: bio-structural-biology-alphafold-predictions
│ · PDB 实验结构 → 直接使用  │
│ · 无实验结构 → AlphaFold DB│
│ · 需要复合物 → AF3/Boltz   │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 2: 结构质量评估       │  技能: bio-structural-biology-alphafold-predictions
│ · pLDDT 分段读取           │
│ · PAE 域间置信度           │
│ · 决策：AFDB vs 新预测     │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 3: 受体结构制备       │  技能: bio-structural-biology-structure-preparation
│ · PDBFixer 补缺失原子     │
│ · PROPKA 预测 pKa         │
│ · reduce 优化翻转/互变异构  │
│ · PDB2PQR 生成 PDBQT      │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 4: 结合口袋检测       │  技能: bio-structural-biology-binding-site-detection
│ · fpocket 几何枚举         │
│ · P2Rank ML 配体性评分     │
│ · 交叉验证，确定 docking 框 │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 5: 配体库制备         │  技能: bio-virtual-screening (ligand prep)
│ · RDKit SMILES → 3D       │
│ · ETKDGv3 构象嵌入         │
│ · MMFF94 能量最小化        │
│ · meeko → PDBQT           │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 6: 层次化虚拟筛选     │  技能: bio-virtual-screening
│ · Stage 1: Vina 筛全库     │
│ · Stage 2: GNINA CNN 重排  │
│ · Stage 3: top-1000 候选   │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 7: Pose 物理验证      │  技能: bio-pose-validation
│ · PoseBusters 全检         │
│ · 应变能量化               │
│ · 芳香环平面度             │
│ · 过滤 PB-invalid 构象     │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 8: ADMET 预测         │  技能: bio-admet-prediction
│ · ADMETlab 3.0 批量预测    │
│ · Lipinski/Veber 类药性    │
│ · hERG/CYP/AMES 毒性       │
│ · PAINS 结构警报            │
└──────────┬───────────────┘
           ▼
┌──────────────────────────┐
│ Step 9: 自由能精排序       │  技能: bio-free-energy-calculations
│ · RBFE 热力学循环          │
│ · MBAR 分析                │
│ · 闭环验证                 │
│ · 或 Boltz-2 快速替代      │
└──────────────────────────┘

import requests

def afdb_metadata(accession):
    """查询 AFDB 元数据接口，返回当前下载 URL"""
    url = f'https://alphafold.ebi.ac.uk/api/prediction/{accession}'
    r = requests.get(url)
    r.raise_for_status()
    return r.json()  # list；每个片段/异构体一个对象

def fetch_afdb(accession, out_dir='.'):
    """下载坐标文件和 PAE 矩阵"""
    entries = afdb_metadata(accession)
    if not entries:
        return None  # >2700aa 非人源蛋白或非 UniProt 序列常缺失
    entry = entries[0]
    cif_text = requests.get(entry['cifUrl']).text
    pae_json = requests.get(entry['paeDocUrl']).json()
    cif_path = f"{out_dir}/AF-{accession}.cif"
    with open(cif_path, 'w') as f:
        f.write(cif_text)
    return cif_path, pae_json

# BRD4 full-length UniProt: O60885
result = fetch_afdb('O60885')

from Bio.PDB import MMCIFParser

def extract_plddt(cif_file):
    """从 B-factor 列读取 pLDDT（注意：与热运动 B-factor 极性相反）"""
    parser = MMCIFParser(QUIET=True)
    structure = parser.get_structure('af', cif_file)
    plddt = {}
    for residue in structure[0].get_residues():
        if residue.id[0] == ' ' and 'CA' in residue:
            # pLDDT 存在 B-factor 列中，但是置信度（0-100，越高越好）
            # 与热运动 B-factor 极性相反——切勿当作热运动读取
            plddt[residue.id[1]] = residue['CA'].get_bfactor()
    return plddt

def plddt_band(score):
    if score > 90:
        return 'very_high'
    if score >= 70:
        return 'confident'
    if score >= 50:
        return 'low'
    return 'very_low'   # <50: 通常是无序区（IDR），不是建模错误

plddt = extract_plddt('AF-O60885.cif')
mean_plddt = sum(plddt.values()) / len(plddt)
disordered = [res for res, s in plddt.items() if s < 50]  # 候选 IDR

import numpy as np

def load_pae(pae_json):
    entry = pae_json[0] if isinstance(pae_json, list) else pae_json
    # 紧凑格式（2023+）: 2D num_res × num_res 数组
    # 遗留 1D 'distance'/'residue1'/'residue2' 字段已移除
    return np.array(entry['predicted_aligned_error'])

def interdomain_confidence(pae, domain_a, domain_b):
    """PAE 不对称，取两个 off-diagonal block 的平均"""
    block = np.concatenate([pae[np.ix_(domain_a, domain_b)].ravel(),
                            pae[np.ix_(domain_b, domain_a)].ravel()])
    return block.mean()  # 低（约 <5Å）= 域间相对位置可信

pae = load_pae(result[1])

from pdbfixer import PDBFixer
from openmm.app import PDBFile

# 注意：finder 调用顺序很重要——残基 → 非标准 → 原子
# hydrogens 必须最后加（重原子存在后才能加氢）
fixer = PDBFixer(filename='receptor.pdb')

fixer.findMissingResidues()           # 短内部缺口 + 末端
fixer.findNonstandardResidues()       # 如 MSE（硒代蛋氨酸）
fixer.replaceNonstandardResidues()    # 映射回标准母体
fixer.removeHeterogens(keepWater=False)  # 去缓冲离子/冷冻保护剂
fixer.findMissingAtoms()              # 截短侧链 + 上面找到的残基
fixer.addMissingAtoms()               # 补重原子；补建的 loop 是假设，需记录

# pH 7.0 是一个选择，不是安全默认——需声明并匹配实验/生物条件
# addMissingHydrogens 检测二硫键并将那些 Cys 留为 CYX
fixer.addMissingHydrogens(pH=7.0)

with open('receptor_prepared.pdb', 'w') as out:
    PDBFile.writeFile(fixer.topology, fixer.positions, out, keepIds=True)

# 记录 provenance
built = [(ci, pos, name) for (ci, pos), names in fixer.missingResidues.items() for name in names]
print(f'built {len(fixer.missingResidues)} missing-residue segment(s), {len(built)} residue(s); pH=7.0')

import subprocess

# -build 执行 -OH -ROTEXOH -HIS -FLIP：加氢并优化 OH/His 旋转 + Asn/Gln/His 翻转
# reduce 通过小探针全原子接触评分优化方向（Word 1999）
with open('receptor_reduced.pdb', 'w') as out:
    subprocess.run(['reduce', '-build', 'receptor_prepared.pdb'], stdout=out, check=True)

# Reduce2（CCTBX/Phenix）是维护中的继承者：mmtbx.reduce2 receptor_prepared.pdb

import subprocess

# propka3 写 receptor_prepared.pka；SUMMARY 列出预测 pKa vs 模型（标准）pKa
subprocess.run(['propka3', 'receptor_prepared.pdb'], check=True)

# PDB2PQR 通过 PROPKA 在 --with-ph 指定 pH 下指定质子化态，写出电荷/半径
# --ff 必须匹配下游力场；--with-ph 7.0 是声明的滴定条件
subprocess.run([
    'pdb2pqr30', '--ff=AMBER', '--with-ph', '7.0',
    '--titration-state-method', 'propka', '--keep-chain',
    'receptor_prepared.pdb', 'receptor.pqr',
], check=True)

# 再用 meeko mk_prepare_receptor.py 转 PDBQT
subprocess.run([
    'mk_prepare_receptor.py', '--read_pqr', 'receptor.pqr',
    '-o', 'receptor', '-p'
], check=True)

# 1) 先裁剪低 pLDDT 残基（<50-70 = 不可靠）
#    phenix.process_predicted_model 做这个 + PAE 域分割
# 2) 然后对裁剪后模型运行 PDBFixer 管线
fixer = PDBFixer(filename='af_model_trimmed.pdb')
fixer.findMissingAtoms()
fixer.addMissingAtoms()
fixer.addMissingHydrogens(pH=7.0)  # 预测模型从不携带实验氢

import subprocess
from pathlib import Path

def run_fpocket(pdb_path):
    subprocess.run(['fpocket', '-f', pdb_path], check=True)
    stem = Path(pdb_path).stem
    return Path(pdb_path).with_name(f'{stem}_out')

def parse_fpocket_info(out_dir):
    info = next(Path(out_dir).glob('*_info.txt'))
    pockets, current = [], None
    for raw in info.read_text().splitlines():
        line = raw.strip()
        if line.startswith('Pocket'):
            current = {'pocket': int(line.split()[1])}
            pockets.append(current)
        elif ':' in line and current is not None:
            key, val = (part.strip() for part in line.split(':', 1))
            try:
                current[key] = float(val)
            except ValueError:
                pass
    return pockets

out_dir = run_fpocket('receptor_prepared.pdb')
pockets = sorted(parse_fpocket_info(out_dir),
                 key=lambda p: p.get('Druggability Score', 0), reverse=True)
# 0.5 药性分是分诊线，不是定律
druggable = [p for p in pockets if p.get('Druggability Score', 0) >= 0.5]
for p in pockets[:5]:
    print(f"pocket {p['pocket']}: drug={p.get('Druggability Score', 0):.2f}")

import csv
import subprocess
from pathlib import Path

def run_p2rank(pdb_path, out_dir='p2rank_out'):
    subprocess.run(['prank', 'predict', '-f', pdb_path, '-o', out_dir], check=True)
    return Path(out_dir) / f'{Path(pdb_path).name}_predictions.csv'

def parse_p2rank(csv_path):
    with open(csv_path) as fh:
        reader = csv.DictReader(fh, skipinitialspace=True)  # P2Rank 列名有空格
        rows = [{k.strip(): v.strip() for k, v in row.items()} for row in reader]
    return [{'rank': int(r['rank']), 'score': float(r['score']),
             'probability': float(r['probability'])} for r in rows]

predictions = parse_p2rank(run_p2rank('receptor_prepared.pdb'))
top = predictions[0]
print(f"top pocket: probability={top['probability']:.2f} score={top['score']:.2f}")

from rdkit import Chem
from rdkit.Chem import AllChem
from meeko import MoleculePreparation, PDBQTWriterLegacy

def prepare_ligand(smiles, pdbqt_out):
    mol = Chem.MolFromSmiles(smiles)
    if mol is None:
        raise ValueError(f'invalid SMILES: {smiles}')
    mol = Chem.AddHs(mol)

    # ETKDGv3 3D 构象嵌入
    embed_status = AllChem.EmbedMolecule(mol, AllChem.ETKDGv3())
    if embed_status != 0:
        raise RuntimeError('ETKDGv3 failed to generate a ligand conformer')

    # MMFF94 能量最小化
    if not AllChem.MMFFHasAllMoleculeParams(mol):
        raise ValueError('MMFF94 parameters unavailable for this ligand')
    optimization_status = AllChem.MMFFOptimizeMolecule(mol)
    if optimization_status != 0:
        raise RuntimeError('MMFF94 optimization did not converge')

    # meeko 0.5+ API：prepare() 返回 MoleculeSetup 列表
    # 用 PDBQTWriterLegacy.write_string() 生成 PDBQT
    mk_prep = MoleculePreparation()
    setups = mk_prep.prepare(mol)
    pdbqt_text, is_ok, err = PDBQTWriterLegacy.write_string(setups[0])
    if not is_ok:
        raise RuntimeError(f'meeko PDBQT export failed: {err}')
    with open(pdbqt_out, 'w') as f:
        f.write(pdbqt_text)
    return pdbqt_out

# 批量制备
library = [
    'CC(=O)c1ccccc1',      # acetophenone
    'O=C(O)c1ccccc1O',     # salicylic acid
    # ... 更多 SMILES
]
for i, smi in enumerate(library):
    prepare_ligand(smi, f'ligand_{i:04d}.pdbqt')

from vina import Vina

def dock_single(receptor_pdbqt, ligand_pdbqt, center, box_size,
                exhaustiveness=8, n_poses=10):
    """
    center: (x, y, z) 口袋中心坐标
    box_size: (x, y, z) 搜索框大小
    """
    v = Vina(sf_name='vina')
    v.set_receptor(receptor_pdbqt)
    v.set_ligand_from_file(ligand_pdbqt)
    v.compute_vina_maps(center=center, box_size=box_size)
    v.dock(exhaustiveness=exhaustiveness, n_poses=n_poses)
    return v.energies(), v.poses()

gnina -r receptor.pdb -l ligand.sdf \
      --autobox_ligand reference_ligand.sdf \
      --cnn_scoring rescore \
      -o poses.sdf.gz \
      --num_modes 9 --exhaustiveness 8

import pandas as pd
from concurrent.futures import ProcessPoolExecutor
from functools import partial

def vs_pipeline(library_smi, receptor_pdbqt, center, box, output_dir, n_workers=16):
    # Stage 0: 读取化合物库
    df = pd.read_csv(library_smi)

    # Stage 1: 类药性过滤（Lipinski + Veber + PAINS）
    df_stage1 = drug_like_filter(df)  # 委托给 admet-prediction 技能

    # Stage 2: Vina 对接全过滤后库
    worker = partial(vina_dock, receptor_pdbqt=receptor_pdbqt,
                     center=center, box=box)
    with ProcessPoolExecutor(max_workers=n_workers) as ex:
        affinities = list(ex.map(worker, df_stage1['smiles']))
    df_stage1['vina_affinity'] = affinities

    # 取 Vina top 1%
    df_stage2 = df_stage1.nsmallest(int(len(df_stage1) * 0.01), 'vina_affinity')

    # Stage 3: GNINA CNN 重打分 top 1%
    df_stage2['gnina_affinity'] = df_stage2['smiles'].apply(
        lambda smi: gnina_rescore(smi, receptor_pdbqt, center, box))

    # 取 GNINA top 1000
    df_stage3 = df_stage2.nsmallest(1000, 'gnina_affinity')

    # Stage 4: Pose 物理验证
    return pose_validate(df_stage3)  # 委托给 pose-validation 技能

from posebusters import PoseBusters
import pandas as pd

bust = PoseBusters(config='dock')  # 'dock' = 全检除 RMSD 参考

results = bust.bust(
    mol_pred='/path/to/docked_poses.sdf',
    mol_cond='/path/to/receptor.pdb',
)

# 聚合所有 bool 列为 pb_valid 标志
check_cols = [
    col for col in results.select_dtypes(include='bool').columns
    if not col.lower().startswith('rmsd')
]
results['pb_valid'] = results[check_cols].all(axis=1)
valid = results[results['pb_valid']]
print(f'{len(valid)} / {len(results)} poses are PB-valid')

from rdkit import Chem
from rdkit.Chem import AllChem

def ligand_strain(docked_sdf, n_ref=20):
    """量化 docked pose 偏离最低能量自由构象的程度"""
    suppl = Chem.SDMolSupplier(docked_sdf, removeHs=False)
    strains = []
    for docked in suppl:
        if docked is None:
            continue
        smi = Chem.MolToSmiles(docked)
        ref = Chem.MolFromSmiles(smi)
        if ref is None:
            strains.append({'strain': None, 'note': 'reference_parse_failed'})
            continue
        ref = Chem.AddHs(ref)
        props_ref = AllChem.MMFFGetMoleculeProperties(ref)
        if props_ref is None:
            strains.append({'strain': None, 'note': 'no_mmff_params'})
            continue
        conf_ids = list(AllChem.EmbedMultipleConfs(
            ref, numConfs=n_ref, params=AllChem.ETKDGv3()))
        if not conf_ids:
            strains.append({'strain': None, 'note': 'embedding_failed'})
            continue
        AllChem.MMFFOptimizeMoleculeConfs(ref)

        ref_energies = []
        for c in conf_ids:
            ff = AllChem.MMFFGetMoleculeForceField(ref, props_ref, confId=c)
            if ff is not None:
                ref_energies.append(ff.CalcEnergy())
        min_ref = min(ref_energies) if ref_energies else None

        # 加 H 坐标，弛豫 H 同时固定重原子
        docked_h = Chem.AddHs(Chem.Mol(docked), addCoords=True)
        if docked_h.GetNumAtoms() != ref.GetNumAtoms():
            strains.append({'strain': None, 'note': 'atom_mismatch'})
            continue
        props_docked = AllChem.MMFFGetMoleculeProperties(docked_h)
        docked_ff = AllChem.MMFFGetMoleculeForceField(docked_h, props_docked) \
            if props_docked is not None else None
        if docked_ff is not None:
            for atom in docked_h.GetAtoms():
                if atom.GetAtomicNum() != 1:
                    docked_ff.AddFixedPoint(atom.GetIdx())
            docked_ff.Minimize(maxIts=200)
            docked_e = docked_ff.CalcEnergy()
            strains.append({
                'min_ref_energy': min_ref,
                'docked_energy': docked_e,
                'strain': docked_e - min_ref,
                'note': 'ok'
            })
        else:
            strains.append({'strain': None, 'note': 'docked_ff_failed'})
    return strains

from rdkit.Chem import Descriptors, Lipinski, QED
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams

def druglike_score(mol):
    """计算类药性指标"""
    mw = Descriptors.MolWt(mol)
    logp = Descriptors.MolLogP(mol)
    hbd = Lipinski.NumHDonors(mol)
    hba = Lipinski.NumHAcceptors(mol)
    tpsa = Descriptors.TPSA(mol)
    rotbonds = Lipinski.NumRotatableBonds(mol)
    qed = QED.qed(mol)

    lipinski_violations = sum([mw > 500, logp > 5, hbd > 5, hba > 10])
    veber_pass = rotbonds <= 10 and tpsa <= 140
    bbb_simple_screen = tpsa <= 90 and mw <= 500 and hbd <= 3  # BBB 启发式

    return {'MW': mw, 'LogP': logp, 'HBD': hbd, 'HBA': hba,
            'TPSA': tpsa, 'RotBonds': rotbonds, 'QED': qed,
            'Lipinski_violations': lipinski_violations,
            'Veber_pass': veber_pass,
            'BBB_simple_screen': bbb_simple_screen}

def alerts(mol, catalogs=('PAINS_A', 'BRENK', 'ZINC')):
    """PAINS/BRENK/REOS 结构警报"""
    params = FilterCatalogParams()
    for cat in catalogs:
        params.AddCatalog(getattr(FilterCatalogParams.FilterCatalogs, cat))
    catalog = FilterCatalog(params)
    hits = catalog.GetMatches(mol)
    return [h.GetDescription() for h in hits]

ΔΔG_binding = ΔG_lig1→lig2(bound) - ΔG_lig1→lig2(solvent)

from openfe import SmallMoleculeComponent, ProteinComponent, SolventComponent
from openfe.protocols.openmm_rfe import RelativeHybridTopologyProtocol

protein = ProteinComponent.from_pdb_file('receptor_prepared.pdb')
ligA = SmallMoleculeComponent.from_sdf_file('ligand_A.sdf')
ligB = SmallMoleculeComponent.from_sdf_file('ligand_B.sdf')
solvent = SolventComponent()

protocol = RelativeHybridTopologyProtocol(
    RelativeHybridTopologyProtocol.default_settings()
)
# 需创建/检查原子映射（Kartograf 是 OpenFE 1.7 CLI 默认映射器）
# 构建 bound 和 solvent Transformation 对象
# 通过 openfe quickrun 或 Python 执行接口运行

from alchemlyb import concat
from alchemlyb.parsing import gmx
from alchemlyb.estimators import MBAR
from alchemlyb.postprocessors.units import to_kcalmol

u_nks = []
for window in range(12):
    df = gmx.extract_u_nk(f'window_{window}.xvg', T=300)
    u_nks.append(df)

u_nk = concat(u_nks)
mbar = MBAR().fit(u_nk)
delta_g = to_kcalmol(mbar.delta_f_).iloc[0, -1]
d_delta_g = to_kcalmol(mbar.d_delta_f_).iloc[0, -1]
print(f'ΔG: {delta_g:.2f} ± {d_delta_g:.2f} kcal/mol')

def cycle_closure_residual(cycle):
    """cycle = [(lig_i, lig_j, delta_g, sd), ...]"""
    total = sum(d_g for _, _, d_g, _ in cycle)
    total_var = sum(sd**2 for _, _, _, sd in cycle)
    return total, total_var ** 0.5

# 残差应与传播不确定性一致（在 ~2σ 内）

靶点蛋白 BRD4-BD1 (O60885)
    │
    ├─ PDB 实验结构 (h3MXF) ──────────────────────────────┐
    │                                                      │
    └─ AlphaFold DB (O60885) ─┐                           │
         │                     │                           │
         ▼                     │                           │
    pLDDT/PAE 评估              │                           │
    裁剪低置信区域               │                           │
         │                     │                           │
         ▼                     ▼                           ▼
    ┌────────────────────────────────────────────────────────┐
    │           受体结构制备 (PDBFixer + PROPKA + reduce)     │
    │           补缺失原子 → 加氢 → 质子化 → PDBQT            │
    └────────────────────────┬───────────────────────────────┘
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  口袋检测 (fpocket + P2Rank)  ← 如果是 apo             │
    │  或从 holo 共晶配体中心定义搜索框                       │
    └────────────────────────┬───────────────────────────────┘
                             │
                  ┌──────────┴──────────┐
                  │                     │
            配体库制备               受体 PDBQT
          (RDKit + meeko)              │
                  │                     │
                  └──────────┬──────────┘
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  Stage 1: 类药性过滤 (Lipinski + Veber + PAINS)        │
    │  Stage 2: Vina 对接 (exhaustiveness=8, 并行)           │
    │  Stage 3: GNINA CNN rescore top 1%                    │
    │  Stage 4: PoseBusters 物理验证 top 1000               │
    │  Stage 5: ADMET 预测 (ADMETlab 3.0 + hERG + CYP)      │
    └────────────────────────┬───────────────────────────────┘
                             │
                    top ~50 候选
                             │
    ┌────────────────────────▼───────────────────────────────┐
    │  自由能精排序                                          │
    │  · 同系列: RBFE (OpenFE, MBAR, 闭环验证)              │
    │  · 跨骨架: ABFE 或 Boltz-2 快速分诊                    │
    └────────────────────────┬───────────────────────────────┘
                             │
                    top 5-10 候选
                             │
                    ▼
              实验验证 (SPR / ITC)