告别手动！用Python脚本+Autodock Vina搞定多对多分子对接与热图绘制（附完整代码）

用Python+Autodock Vina实现全自动分子对接与热力分析

在药物发现和生物分子相互作用研究中，分子对接是预测小分子（配体）与生物大分子（受体）结合模式和亲和力的重要工具。传统手动操作不仅效率低下，还容易引入人为错误。本文将展示如何通过Python脚本实现从批量对接、结果解析到可视化呈现的全流程自动化。

1. 环境配置与文件组织

1.1 软件准备与依赖安装

确保系统已安装以下组件：

• Autodock Vina 1.2.3+
• Python 3.8+ 及以下关键库：

  pip install pandas numpy seaborn matplotlib

1.2 项目目录结构规范

建议采用以下目录结构保持工作流清晰：

project_root/ │── receptors/ # 存放受体.pdbqt文件 │ ├── receptor1.pdbqt │ ├── receptor1.txt # 对接参数文件 │── ligands/ # 存放配体.pdbqt文件 │── results/ # 脚本自动生成的结果目录 │── docking_script.py # 主程序脚本

每个受体需要配套的对接参数文件（.txt），示例内容：

center_x = 18.21 center_y = 12.45 center_z = -5.33 size_x = 25 size_y = 25 size_z = 25 num_modes = 10

2. 核心对接逻辑实现

2.1 多进程并行对接引擎

通过Python的subprocess模块调用Vina，并利用多进程加速：

import subprocess from multiprocessing import Pool def run_vina_docking(receptor, ligand, config): cmd = [ 'vina', '--receptor', receptor, '--ligand', ligand, '--config', config, '--out', f'results/{receptor.stem}_{ligand.stem}.pdbqt', '--log', f'results/{receptor.stem}_{ligand.stem}.log' ] try: subprocess.run(cmd, check=True, timeout=300) return (receptor.stem, ligand.stem, True) except Exception as e: print(f"Failed {receptor} - {ligand}: {str(e)}") return (receptor.stem, ligand.stem, False)

2.2 自动化任务调度器

批量生成对接任务队列并执行：

from pathlib import Path def batch_docking(): receptors = list(Path('receptors').glob('*.pdbqt')) ligands = list(Path('ligands').glob('*.pdbqt')) tasks = [] for rec in receptors: config = Path('receptors') / f"{rec.stem}.txt" for lig in ligands: tasks.append((rec, lig, config)) with Pool(processes=4) as pool: # 根据CPU核心数调整 results = pool.starmap(run_vina_docking, tasks) return results

3. 结果解析与数据清洗

3.1 亲和力数据提取优化

改进的日志解析函数，支持不同Vina版本：

def parse_affinity(log_file): with open(log_file, 'r') as f: lines = f.readlines() for line in reversed(lines[-20:]): # 从后向前搜索提高效率 if 'Affinity' in line: try: return float(line.split()[1]) except: continue return None

3.2 构建结构化结果矩阵

使用pandas构建交互数据框架：

import pandas as pd def build_results_df(): data = [] for log in Path('results').glob('*.log'): parts = log.stem.split('_') receptor, ligand = parts[0], '_'.join(parts[1:]) affinity = parse_affinity(log) if affinity is not None: data.append({ 'Receptor': receptor, 'Ligand': ligand, 'Affinity(kcal/mol)': affinity, 'PDBQT': f"{log.stem}.pdbqt" }) df = pd.DataFrame(data) df.to_csv('results/summary.csv', index=False) return df

4. 高级可视化与交互分析

4.1 动态热图生成技术

创建可交互的Heatmap并支持阈值筛选：

import seaborn as sns import matplotlib.pyplot as plt from ipywidgets import interact, FloatSlider def interactive_heatmap(df, threshold=-7.0): filtered = df[df['Affinity(kcal/mol)'] <= threshold] pivot_df = filtered.pivot(index='Ligand', columns='Receptor', values='Affinity(kcal/mol)') plt.figure(figsize=(12, 8)) ax = sns.heatmap( pivot_df, cmap='coolwarm', annot=True, fmt='.1f', linewidths=.5, cbar_kws={'label': 'Binding Affinity (kcal/mol)'} ) ax.set_title(f'Docking Results (Threshold ≤ {threshold} kcal/mol)') plt.tight_layout() plt.savefig(f'results/heatmap_threshold_{threshold}.png', dpi=300) plt.show() # 添加交互控件 def update_threshold(threshold=-7.0): df = pd.read_csv('results/summary.csv') interactive_heatmap(df, threshold) interact(update_threshold, threshold=FloatSlider(min=-15, max=0, step=0.5, value=-7.0))

4.2 结果筛选与导出功能

实现基于条件的自动筛选和报告生成：

def generate_report(affinity_cutoff=-7.0, top_n=5): df = pd.read_csv('results/summary.csv') # 生成各受体最佳配体 best_by_receptor = df.loc[df.groupby('Receptor')['Affinity(kcal/mol)'].idxmin()] # 全局最佳结果 global_best = df.nsmallest(top_n, 'Affinity(kcal/mol)') # 保存筛选结果 with pd.ExcelWriter('results/filtered_results.xlsx') as writer: df.to_excel(writer, sheet_name='All Results') best_by_receptor.to_excel(writer, sheet_name='Best per Receptor') global_best.to_excel(writer, sheet_name=f'Top {top_n}') print(f"Report generated with {len(df)} total results")

5. 错误处理与日志系统

5.1 健壮性增强设计

添加完善的异常处理机制：

import logging from datetime import datetime logging.basicConfig( filename=f'docking_log_{datetime.now().strftime("%Y%m%d")}.txt', level=logging.INFO, format='%(asctime)s - %(levelname)s - %(message)s' ) def safe_docking_wrapper(receptor, ligand, config): try: result = run_vina_docking(receptor, ligand, config) if result[2]: logging.info(f"Success: {result[0]} vs {result[1]}") else: logging.warning(f"Failed: {result[0]} vs {result[1]}") return result except Exception as e: logging.error(f"Critical error with {receptor} vs {ligand}: {str(e)}") return (receptor.stem, ligand.stem, False)

5.2 性能监控与优化

添加资源使用统计功能：

import time import psutil def monitor_resources(): start_time = time.time() cpu_percent = [] memory_usage = [] def callback(): cpu_percent.append(psutil.cpu_percent()) memory_usage.append(psutil.virtual_memory().percent) return (cpu_percent[-1], memory_usage[-1]) return callback, lambda: { 'total_time': time.time() - start_time, 'avg_cpu': sum(cpu_percent)/len(cpu_percent), 'max_memory': max(memory_usage) }

将核心对接函数修改为：

def run_vina_docking(receptor, ligand, config, monitor_cb=None): cmd = [...] # 同上 try: process = subprocess.Popen(cmd) while process.poll() is None: if monitor_cb: monitor_cb() time.sleep(1) return (receptor.stem, ligand.stem, True) except Exception as e: return (receptor.stem, ligand.stem, False)