分子对接AutoDock实战¶

分子对接 = 找最佳"插入"姿势

蛋白质（受体）= 一把锁
小分子（配体）= 一把钥匙

对接过程：
1. 把"钥匙"放在"锁"旁边
2. 让"钥匙"旋转、翻转、弯曲
3. 试各种姿势（构象采样）
4. 用打分函数评估每种姿势的"契合程度"
5. 找到最佳姿势（最低结合能）

结合能（kcal/mol）：
  -3 ~ -5  → 弱结合（一般不行）
  -5 ~ -7  → 中等结合（值得关注）
  -7 ~ -9  → 强结合（好苗子）
  < -9      → 极强结合（可能的药物候选）

# 方法一：conda安装（推荐）
conda create -n docking python=3.10  # 创建对接专用环境
conda activate docking  # 激活环境
conda install -c conda-forge autodock-vina  # 安装AutoDock Vina

# 方法二：pip安装Python接口
pip install vina  # 安装Vina的Python绑定
pip install meeko  # 安装分子准备工具Meeko

# 方法三：直接下载二进制文件
wget https://github.com/ccsb-scripps/AutoDock-Vina/releases/download/v1.2.7/vina_1.2.7_linux_x86_64  # 下载Linux版（最新v1.2.7）
chmod +x vina_1.2.7_linux_x86_64  # 添加执行权限
mv vina_1.2.7_linux_x86_64 ~/bin/vina  # 移动到PATH目录

# 验证安装
vina --version  # 查看版本信息

# 安装MGLTools（分子准备工具）
conda install -c conda-forge mgltools  # 安装MGLTools

# 安装OpenBabel（分子格式转换）
conda install -c conda-forge openbabel  # 安装OpenBabel

# 安装RDKit（化学信息学工具包）
conda install -c conda-forge rdkit  # 安装RDKit

# 安装PLIP（蛋白质-配体互作分析）
pip install plip  # 安装互作分析工具

#!/usr/bin/env python3
"""准备蛋白质受体文件"""

from vina import Vina  # 导入Vina模块
import subprocess  # 用于调用外部命令

# ========== 第一步：获取蛋白质结构 ==========
# 从PDB数据库下载蛋白质结构
import urllib.request  # HTTP请求模块
pdb_id = "6LU7"  # 新冠病毒主蛋白酶
url = f"https://files.rcsb.org/download/{pdb_id}.pdb"  # PDB下载地址
urllib.request.urlretrieve(url, f"{pdb_id}.pdb")  # 下载PDB文件
print(f"已下载 {pdb_id}.pdb")

# ========== 第二步：清理蛋白质 ==========
def clean_pdb(input_pdb, output_pdb):
    """清理PDB文件：去水、去配体、只保留蛋白质"""
    with open(input_pdb, 'r') as f:  # 读取PDB文件
        lines = f.readlines()

    clean_lines = []
    for line in lines:  # 遍历每一行
        # 只保留ATOM记录（蛋白质原子）
        if line.startswith("ATOM"):  # ATOM=蛋白质原子
            clean_lines.append(line)
        elif line.startswith("END"):  # 保留END标记
            clean_lines.append(line)
        # 跳过HETATM（配体、水、离子等）

    with open(output_pdb, 'w') as f:  # 写入清理后的文件
        f.writelines(clean_lines)

    print(f"清理完成: {output_pdb}")

clean_pdb("6LU7.pdb", "6LU7_clean.pdb")  # 清理蛋白质

# ========== 第三步：添加氢原子并转为PDBQT ==========
# 使用OpenBabel添加氢原子
subprocess.run([
    "obabel", "6LU7_clean.pdb",  # 输入文件
    "-O", "6LU7_receptor.pdbqt",  # 输出PDBQT格式
    "-xr",  # 只输出受体（不含配体）
    "--addhydrogens"  # 添加氢原子
], check=True)
print("受体PDBQT文件已准备好")

#!/usr/bin/env python3
"""准备小分子配体文件"""

from rdkit import Chem  # RDKit化学模块
from rdkit.Chem import AllChem  # 分子力场模块
import subprocess  # 调用外部命令

# ========== 方法一：从SMILES生成配体 ==========
smiles = "CC(=O)Oc1ccccc1C(=O)O"  # 阿司匹林的SMILES式

mol = Chem.MolFromSmiles(smiles)  # 从SMILES创建分子对象
mol = Chem.AddHs(mol)  # 添加氢原子

# 生成3D构象
AllChem.EmbedMolecule(mol, randomSeed=42)  # 生成3D坐标
AllChem.MMFFOptimizeMolecule(mol)  # 用MMFF力场优化分子构型

# 保存为SDF文件
writer = Chem.SDWriter("aspirin.sdf")  # 创建SDF写入器
writer.write(mol)  # 写入分子
writer.close()
print("配体SDF文件已生成")

# ========== 转为PDBQT格式 ==========
# 使用Meeko准备配体（推荐，支持柔性）
subprocess.run([
    "mk_prepare_ligand.py",  # Meeko的配体准备脚本
    "-i", "aspirin.sdf",  # 输入SDF文件
    "-o", "aspirin.pdbqt"  # 输出PDBQT文件
], check=True)
print("配体PDBQT文件已准备好")

# ========== 方法二：从PubChem下载 ==========
import urllib.request
cid = 2244  # 阿司匹林的PubChem CID
url = f"https://pubchem.ncbi.nlm.nih.gov/rest/pug/compound/CID/{cid}/SDF"  # PubChem下载地址
urllib.request.urlretrieve(url, "aspirin_pubchem.sdf")  # 下载SDF文件
print("从PubChem下载配体完成")

#!/usr/bin/env python3
"""使用AutoDock Vina执行分子对接"""

from vina import Vina  # 导入Vina

# ========== 创建Vina对象 ==========
v = Vina(sf_name="vina")  # 使用Vina打分函数

# ========== 设置受体 ==========
v.set_receptor("6LU7_receptor.pdbqt")  # 加载受体PDBQT文件

# ========== 设置配体 ==========
v.set_ligand_from_file("aspirin.pdbqt")  # 加载配体PDBQT文件

# ========== 设置对接盒子 ==========
# 对接盒子参数需要根据活性位点位置设定
# 可以用PyMOL查看活性位点坐标来确定
v.compute_vina_maps(
    center=[-10.0, 12.0, 68.0],  # 盒子中心坐标（x, y, z）
    box_size=[25, 25, 25]  # 盒子大小（单位：埃）
)

# ========== 执行对接 ==========
v.dock(
    exhaustiveness=32,  # 搜索穷举度（越大越准，越慢），默认8
    n_poses=10  # 输出前10个最佳构象
)

# ========== 查看结果 ==========
energies = v.energies()  # 获取所有构象的能量
print("\n===== 对接结果 =====")
print(f"{'排名':<6}{'结合能(kcal/mol)':<20}{'RMSD_lb':<10}{'RMSD_ub':<10}")
print("-" * 46)
for i, e in enumerate(energies):  # 遍历每个构象
    print(f"{i+1:<6}{e[0]:<20.2f}{e[1]:<10.2f}{e[2]:<10.2f}")

# ========== 保存最佳构象 ==========
v.write_poses("docking_results.pdbqt", n_poses=10, overwrite=True)  # 保存所有构象
print("\n对接结果已保存到 docking_results.pdbqt")

# ========== 命令行方式运行AutoDock Vina ==========

# 准备配置文件 config.txt
cat > config.txt << 'EOF'
receptor = 6LU7_receptor.pdbqt
ligand = aspirin.pdbqt

center_x = -10.0
center_y = 12.0
center_z = 68.0

size_x = 25
size_y = 25
size_z = 25

exhaustiveness = 32
num_modes = 10
energy_range = 3
EOF

# 运行对接
vina --config config.txt --out docking_results.pdbqt --log docking_log.txt  # 执行对接

# 查看结果
cat docking_log.txt  # 查看对接日志
grep "^   " docking_log.txt  # 只看结果行（结合能排名）

# ========== PyMOL查看对接结果 ==========
# 在PyMOL命令行中运行

load 6LU7_clean.pdb, receptor  # 加载受体蛋白质
load docking_results.pdbqt, ligand  # 加载对接结果

# 显示设置
hide everything  # 隐藏所有
show cartoon, receptor  # 受体卡通显示
color marine, receptor  # 受体蓝色

show sticks, ligand  # 配体棍状显示
color yellow, ligand  # 配体黄色

# 显示结合口袋
select pocket, byres (ligand around 5)  # 配体5埃内的残基
show sticks, pocket  # 口袋棍状显示
color cyan, pocket  # 口袋青色

# 显示氢键
distance hbonds, ligand, pocket, 3.5, mode=2  # 分析氢键
set dash_color, red  # 氢键红色虚线

# 导出图片
orient  # 调整视角
ray 2400, 2000  # 渲染
png docking_result.png, dpi=300  # 保存图片

#!/usr/bin/env python3
"""批量虚拟筛选：对多个化合物进行对接"""

from vina import Vina  # 导入Vina
import os  # 文件操作
import pandas as pd  # 数据处理
from glob import glob  # 文件匹配

# ========== 准备 ==========
v = Vina(sf_name="vina")  # 创建Vina对象
v.set_receptor("6LU7_receptor.pdbqt")  # 设置受体（只需一次）

# 设置对接盒子（只需一次）
v.compute_vina_maps(
    center=[-10.0, 12.0, 68.0],  # 活性位点中心
    box_size=[25, 25, 25]  # 盒子大小
)

# ========== 批量对接 ==========
ligand_dir = "ligands/"  # 配体文件目录
results = []  # 存储结果

ligand_files = glob(f"{ligand_dir}/*.pdbqt")  # 获取所有配体文件
print(f"共找到 {len(ligand_files)} 个配体")

for lig_file in ligand_files:  # 遍历每个配体
    lig_name = os.path.basename(lig_file).replace(".pdbqt", "")  # 配体名称
    print(f"正在对接: {lig_name}...", end=" ")

    try:
        v.set_ligand_from_file(lig_file)  # 加载配体
        v.dock(exhaustiveness=16, n_poses=5)  # 对接（虚拟筛选用16就够）

        energy = v.energies()[0][0]  # 获取最佳构象的结合能
        results.append({
            "ligand": lig_name,  # 配体名
            "binding_energy": energy  # 结合能
        })

        # 保存最佳构象
        v.write_poses(f"results/{lig_name}_docked.pdbqt", n_poses=1, overwrite=True)
        print(f"结合能 = {energy:.2f} kcal/mol")

    except Exception as e:  # 捕获错误
        print(f"失败: {e}")
        results.append({"ligand": lig_name, "binding_energy": None})

# ========== 结果汇总 ==========
df_results = pd.DataFrame(results)  # 转为表格
df_results = df_results.sort_values("binding_energy")  # 按结合能排序
df_results.to_csv("virtual_screening_results.csv", index=False)  # 保存结果

print("\n===== 虚拟筛选结果TOP10 =====")
print(df_results.head(10).to_string(index=False))  # 显示结合能最低的前10个

========================================
分子对接AutoDock Vina 速查表
========================================

【安装】
conda安装            → conda install -c conda-forge autodock-vina
pip安装              → pip install vina meeko
验证                 → vina --version

【文件格式】
PDB → PDBQT（受体）  → obabel protein.pdb -O protein.pdbqt -xr
SDF → PDBQT（配体）  → mk_prepare_ligand.py -i ligand.sdf -o ligand.pdbqt

【对接参数】
exhaustiveness       → 搜索穷举度，默认8，建议16-32
n_poses / num_modes  → 输出构象数，默认9
energy_range         → 能量范围(kcal/mol)，默认3

【结合能判读】
> -5  kcal/mol       → 弱结合
-5 ~ -7              → 中等结合
-7 ~ -9              → 强结合
< -9                  → 极强结合

【关键步骤】
1. 下载蛋白质         → RCSB PDB
2. 清理受体           → 去水、去配体、加氢
3. 准备配体           → 3D构象、加氢、PDBQT
4. 设置盒子           → 覆盖活性位点
5. 执行对接           → vina --config config.txt
6. 分析结果           → PyMOL可视化

【面试考点】
Q: 分子对接的原理是什么？
A: 在受体活性位点附近搜索配体最佳构象，用打分函数评估结合亲和力

Q: AutoDock Vina的打分函数包含什么？
A: 疏水效应、氢键、范德华力、扭转惩罚等

Q: 虚拟筛选和分子对接的区别？
A: 虚拟筛选=批量对接，从化合物库中筛选可能的药物候选
========================================