基于Python与Pandas实现AutoDock Vina对接结果的自动化筛选与富集分析
1. 从零开始处理AutoDock Vina对接结果
刚接触AutoDock Vina的朋友可能会被一堆log.txt文件搞得头晕。我刚开始做分子对接时,每次跑完几十个化合物,看着满屏的文本文件就头疼——怎么才能快速找到那些结合能特别好的化合物呢?今天就跟大家分享我用Python和Pandas搭建的自动化处理流程。
这个流程特别适合需要处理大批量对接结果的情况。比如我最近做的一个项目,要筛选200多个天然产物化合物库,手动查看每个log文件根本不现实。通过Python脚本,不到5分钟就能完成所有数据的提取、清洗和筛选。
先说说我们需要准备什么:
- 一堆AutoDock Vina生成的log.txt文件(通常每个化合物一个文件夹)
- 一个记录化合物信息的Excel表格(包含CID、中英文名称等)
- 安装了Python和Pandas的环境
2. 批量读取log.txt文件
2.1 文件路径处理技巧
我建议把所有对接结果放在一个总文件夹下,每个化合物单独一个子文件夹。这样用os.walk遍历起来特别方便。这里有个小技巧:Windows路径中的反斜杠容易出问题,可以用raw string或者正斜杠:
import os
import pandas as pd
import numpy as np
# 推荐这种写法,避免转义字符问题
base_dir = r'C:\Users\yourname\dock_results'
# 或者用正斜杠
base_dir = 'C:/Users/yourname/dock_results'
2.2 高效提取结合能数据
AutoDock Vina的log文件有固定格式,结合能信息通常在特定行。我写了个更健壮的提取函数,能处理各种异常情况:
def extract_affinity(file_path):
try:
with open(file_path, 'r') as f:
lines = f.readlines()
# 通常第25行是第一个对接构象的结果
if len(lines) > 25:
line = lines[24].strip()
# 提取"-7.2"这样的数值
return float(line.split()[1])
except:
return None
这个改进版函数直接返回浮点数,方便后续处理。我还加了个try-except,防止文件损坏导致程序崩溃。
3. 构建并清洗数据框
3.1 创建结构化数据
用字典收集完所有数据后,转换成Pandas DataFrame时有几个注意事项:
# 收集数据的字典
results = {}
for root, dirs, files in os.walk(base_dir):
for dir_name in dirs:
log_path = os.path.join(root, dir_name, 'log.txt')
affinity = extract_affinity(log_path)
if affinity is not None: # 只保留有效数据
results[dir_name] = affinity
# 转换为DataFrame
df = pd.DataFrame.from_dict(results, orient='index', columns=['affinity'])
df = df.reset_index().rename(columns={'index': 'CID'})
3.2 数据清洗实战经验
清洗数据时我踩过不少坑,分享几个实用技巧:
- 处理异常值:有时对接失败会产生极大或极小的值
df = df[(df['affinity'] > -20) & (df['affinity'] < 0)] # 合理范围过滤
- 类型转换:确保affinity列是float类型
df['affinity'] = pd.to_numeric(df['affinity'], errors='coerce')
- 缺失值处理:
df = df.dropna(subset=['affinity'])
4. 高级筛选与排序
4.1 灵活设置筛选阈值
原始文章用-7.0作为阈值,但实际项目中可能需要动态调整:
threshold = float(input("请输入结合能阈值(如-7.0): ") or "-7.0")
filtered_df = df[df['affinity'] < threshold]
还可以筛选前N个最优结果:
top_n = filtered_df.nsmallest(10, 'affinity')
4.2 多条件复合筛选
比如同时要求结合能< -7.0且RMSD < 2.0:
filtered_df = df[(df['affinity'] < -7.0) & (df['rmsd'] < 2.0)]
5. 化合物信息关联匹配
5.1 处理CID格式问题
原始数据中的CID通常带有前缀,需要清洗:
# 更健壮的处理方式
filtered_df['CID'] = filtered_df['CID'].str.extract('(\d+)$')
5.2 数据库合并技巧
合并化合物信息时,要注意数据类型一致:
compound_db = pd.read_excel('compound_database.xlsx')
compound_db['CID'] = compound_db['CID'].astype(str) # 确保类型匹配
# 使用merge时指定后缀是个好习惯
final_df = pd.merge(
filtered_df,
compound_db,
on='CID',
how='left',
suffixes=('', '_db')
)
5.3 处理合并后的数据
合并后可能需要进一步清洗:
# 去除没有匹配到的化合物
final_df = final_df.dropna(subset=['compound_name'])
# 重命名列更直观
final_df = final_df.rename(columns={
'affinity': 'Binding_Energy',
'name_cn': 'Chinese_Name'
})
6. 结果输出与可视化
6.1 多样化输出格式
除了Excel,还可以输出CSV或JSON:
# Excel
final_df.to_excel('final_results.xlsx', index=False)
# CSV
final_df.to_csv('final_results.csv', index=False, encoding='utf-8-sig')
# JSON
final_df.to_json('final_results.json', orient='records', force_ascii=False)
6.2 简单可视化分析
用Matplotlib快速生成结合能分布图:
import matplotlib.pyplot as plt
plt.figure(figsize=(10,6))
plt.hist(final_df['Binding_Energy'], bins=20, color='skyblue')
plt.xlabel('Binding Energy (kcal/mol)')
plt.ylabel('Number of Compounds')
plt.title('Distribution of Binding Energies')
plt.grid(True)
plt.savefig('energy_distribution.png', dpi=300)
plt.close()
7. 脚本优化与批量处理
7.1 封装成可重用函数
把核心流程封装成函数,方便复用:
def process_vina_results(input_dir, db_path, threshold=-7.0):
# 实现上述所有步骤
...
return final_df
7.2 添加日志记录
添加日志功能方便调试:
import logging
logging.basicConfig(
filename='vina_processing.log',
level=logging.INFO,
format='%(asctime)s - %(levelname)s - %(message)s'
)
try:
results = process_vina_results(...)
logging.info(f"成功处理{len(results)}个化合物")
except Exception as e:
logging.error(f"处理失败: {str(e)}")
8. 实际应用中的经验分享
在真实项目中,我发现几个常见问题需要特别注意:
- 文件编码问题:有些log文件可能是gbk编码,需要指定:
with open(file_path, 'r', encoding='gbk') as f:
...
- 内存管理:处理上万化合物时,考虑分块处理:
chunksize = 1000
for chunk in pd.read_csv('huge_file.csv', chunksize=chunksize):
process(chunk)
- 并行处理加速:对于超大规模数据,可以用multiprocessing:
from multiprocessing import Pool
def process_file(file_path):
...
with Pool(processes=4) as pool:
results = pool.map(process_file, file_list)
最后保存结果时,我习惯添加处理日期和时间戳,避免覆盖之前的分析结果:
from datetime import datetime
timestamp = datetime.now().strftime("%Y%m%d_%H%M")
output_file = f"vina_results_{timestamp}.xlsx"
final_df.to_excel(output_file)
更多推荐
所有评论(0)