278 lines
8.3 KiB
Python
Executable File
278 lines
8.3 KiB
Python
Executable File
#!/usr/bin/env python
|
||
# -*- encoding: utf-8 -*-
|
||
'''
|
||
@file :simemacrocycle_repair.py
|
||
@Description: :SIME工具生成的16元环大环内酯化合物的自动化价态修复工具
|
||
@Date :2025/03/29 18:29:52
|
||
@Author :lyzeng
|
||
@Email :pylyzeng@gmail.com
|
||
@version :1.0
|
||
# 安装依赖
|
||
pip install rdkit pandas swifter joblib tqdm matplotlib
|
||
|
||
# 运行脚本
|
||
python simemacrocycle_repair.py
|
||
'''
|
||
|
||
##############################
|
||
# 模块说明
|
||
##############################
|
||
"""
|
||
主要功能：
|
||
1. 批量修复SIME生成的含双键16元环大环内酯的价态错误
|
||
2. 自动检测并处理以下问题：
|
||
- 碳原子显式价态超限（如5价碳）
|
||
- 不合理的显式氢配置
|
||
- 双键立体化学冲突
|
||
3. 提供修复结果统计和可视化分析
|
||
|
||
输入输出：
|
||
- 输入：包含SMILES的文本文件（每行一个分子）
|
||
- 输出：
|
||
- 修复后的CSV文件（含原始/修正SMILES和状态）
|
||
- 修复统计图表（PNG）
|
||
- 摘要报告（TXT）
|
||
|
||
依赖环境：
|
||
- Python >= 3.7
|
||
- RDKit >= 2022.03
|
||
- pandas >= 1.3
|
||
- swifter >= 1.3
|
||
"""
|
||
|
||
import pandas as pd
|
||
from pathlib import Path
|
||
import swifter
|
||
from joblib import Parallel, delayed
|
||
from tqdm.auto import tqdm
|
||
import matplotlib.pyplot as plt
|
||
from collections import Counter
|
||
import re
|
||
from rdkit import Chem
|
||
import warnings
|
||
warnings.filterwarnings('ignore', category=UserWarning)
|
||
|
||
##############################
|
||
# 核心修复函数
|
||
##############################
|
||
|
||
def safe_sanitize(mol):
|
||
"""
|
||
安全标准化分子结构
|
||
|
||
Parameters:
|
||
mol (rdkit.Chem.Mol): 待检测的分子对象
|
||
|
||
Returns:
|
||
int/None: 返回错误原子索引（如存在价态错误），无错误则返回None
|
||
"""
|
||
try:
|
||
Chem.SanitizeMol(mol)
|
||
return None
|
||
except Chem.AtomValenceException as e:
|
||
match = re.search(r'atom # (\d+)', str(e))
|
||
return int(match.group(1)) if match else None
|
||
except:
|
||
return None
|
||
|
||
|
||
def fix_valence_error(smi):
|
||
"""
|
||
修复单个SMILES的价态错误
|
||
|
||
Parameters:
|
||
smi (str): 输入SMILES字符串
|
||
|
||
Returns:
|
||
tuple: (修正后SMILES, 状态, 描述信息)
|
||
状态可能值:
|
||
- 'valid': 原始有效
|
||
- 'corrected': 修复成功
|
||
- 'kekule': 需Kekule形式
|
||
- 'failed': 修复失败
|
||
"""
|
||
try:
|
||
mol = Chem.MolFromSmiles(smi, sanitize=False)
|
||
if not mol:
|
||
return smi, "invalid", "无法解析SMILES"
|
||
|
||
error_atom_idx = safe_sanitize(Chem.Mol(mol))
|
||
if error_atom_idx is None:
|
||
return Chem.MolToSmiles(mol), "valid", "原始有效"
|
||
|
||
rw_mol = Chem.RWMol(mol)
|
||
atom = rw_mol.GetAtomWithIdx(error_atom_idx)
|
||
|
||
# 修复策略序列
|
||
repair_actions = [
|
||
lambda: atom.SetNumExplicitHs(0),
|
||
lambda: (atom.SetFormalCharge(0), atom.SetNumRadicalElectrons(0)),
|
||
lambda: rw_mol.RemoveBond(
|
||
list(atom.GetBonds())[-1].GetBeginAtomIdx(),
|
||
list(atom.GetBonds())[-1].GetEndAtomIdx()
|
||
) if atom.GetDegree() > 1 else None
|
||
]
|
||
|
||
for action in repair_actions:
|
||
action()
|
||
if safe_sanitize(rw_mol.GetMol()) is None:
|
||
return Chem.MolToSmiles(rw_mol.GetMol()), "corrected", f"修复原子 {error_atom_idx}"
|
||
|
||
Chem.Kekulize(rw_mol)
|
||
return Chem.MolToSmiles(rw_mol.GetMol(), kekuleSmiles=True), "kekule", "返回Kekule形式"
|
||
|
||
except Exception as e:
|
||
return smi, "failed", str(e)
|
||
|
||
##############################
|
||
# 并行处理模块
|
||
##############################
|
||
|
||
def batch_process(smi_chunk):
|
||
"""
|
||
分块处理SMILES列表（兼容并行化）
|
||
|
||
Parameters:
|
||
smi_chunk (list): SMILES字符串列表
|
||
|
||
Returns:
|
||
list: 包含(修正SMILES, 状态, 信息)的元组列表
|
||
"""
|
||
return [fix_valence_error(smi) for smi in smi_chunk]
|
||
|
||
|
||
def process_in_chunks(smi_list, chunk_size=50000, n_jobs=4):
|
||
"""
|
||
分块并行处理大规模SMILES数据
|
||
|
||
Parameters:
|
||
smi_list (list): 原始SMILES列表
|
||
chunk_size (int): 每块处理量
|
||
n_jobs (int): 并行进程数
|
||
|
||
Returns:
|
||
tuple: (修正SMILES列表, 状态列表, 信息列表)
|
||
"""
|
||
results = []
|
||
for i in tqdm(range(0, len(smi_list), chunk_size),
|
||
desc=f"Processing {len(smi_list):,} molecules"):
|
||
chunk = smi_list[i:i + chunk_size]
|
||
chunk_results = Parallel(n_jobs=n_jobs)(
|
||
delayed(batch_process)(chunk[i:i+1000])
|
||
for i in range(0, len(chunk), 1000)
|
||
)
|
||
results.extend([item for sublist in chunk_results for item in sublist])
|
||
|
||
return list(zip(*results)) if results else ([], [], [])
|
||
|
||
##############################
|
||
# 统计分析模块
|
||
##############################
|
||
|
||
def analyze_results(df):
|
||
"""
|
||
生成修复结果统计分析报告
|
||
|
||
Parameters:
|
||
df (pd.DataFrame): 包含修复结果的DataFrame
|
||
|
||
Returns:
|
||
dict: 包含关键统计指标的字典
|
||
"""
|
||
# 计算基本统计量
|
||
stats = {
|
||
'total_molecules': len(df),
|
||
'valid_count': len(df[df['status'] == 'valid']),
|
||
'corrected_count': len(df[df['status'] == 'corrected']),
|
||
'kekule_count': len(df[df['status'] == 'kekule']),
|
||
'failed_count': len(df[df['status'] == 'failed']),
|
||
'success_rate': (len(df[df['status'].isin(['valid', 'corrected'])]) / len(df))
|
||
}
|
||
|
||
# 错误分析（仅当存在失败时）
|
||
if stats['failed_count'] > 0:
|
||
stats['common_errors'] = dict(df[df['status'] == 'failed']['message'].value_counts().head(5))
|
||
else:
|
||
stats['common_errors'] = {}
|
||
|
||
# 可视化
|
||
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
|
||
|
||
# 状态分布饼图
|
||
status_counts = df['status'].value_counts()
|
||
status_counts.plot.pie(ax=ax1, autopct='%1.1f%%', startangle=90)
|
||
ax1.set_title('修复状态分布')
|
||
|
||
# 错误原因柱状图（仅当存在错误时）
|
||
if stats['common_errors']:
|
||
pd.Series(stats['common_errors']).plot.barh(ax=ax2)
|
||
ax2.set_title('Top 5错误原因')
|
||
else:
|
||
ax2.axis('off')
|
||
ax2.text(0.5, 0.5, '无修复失败记录',
|
||
ha='center', va='center', fontsize=12)
|
||
|
||
plt.tight_layout()
|
||
plt.savefig('repair_statistics.png', dpi=150)
|
||
plt.close()
|
||
|
||
return stats
|
||
|
||
##############################
|
||
# 主执行流程
|
||
##############################
|
||
|
||
def main(input_path, output_path="fixed_molecules.csv", n_jobs=4):
|
||
"""
|
||
主处理流程
|
||
|
||
Parameters:
|
||
input_path (str): 输入SMILES文件路径
|
||
output_path (str): 输出CSV文件路径
|
||
n_jobs (int): 并行进程数
|
||
"""
|
||
print(f"\n{' SIME大环内酯修复工具 ':=^50}\n")
|
||
|
||
# 数据加载
|
||
smi_list = [s.strip() for s in Path(input_path).read_text().splitlines() if s.strip()]
|
||
print(f"✅ 已加载 {len(smi_list):,} 个分子")
|
||
|
||
# 分子修复
|
||
fixed_smiles, statuses, messages = process_in_chunks(smi_list, n_jobs=n_jobs)
|
||
|
||
# 结果分析
|
||
df = pd.DataFrame({
|
||
'original_smiles': smi_list,
|
||
'fixed_smiles': fixed_smiles,
|
||
'status': statuses,
|
||
'message': messages
|
||
})
|
||
|
||
stats = analyze_results(df)
|
||
|
||
# 结果保存
|
||
df.to_csv(output_path, index=False)
|
||
print(f"\n{' 修复结果统计 ':=^50}")
|
||
print(f"总处理数: {stats['total_molecules']:,}")
|
||
print(f"成功率: {stats['success_rate']:.2%}")
|
||
print(f"\n输出文件已保存至: {output_path} 和 repair_statistics.png")
|
||
|
||
if __name__ == "__main__":
|
||
import argparse
|
||
parser = argparse.ArgumentParser(description="SIME大环内酯修复工具")
|
||
parser.add_argument('input', help="输入SMILES文件路径")
|
||
parser.add_argument('-o', '--output', default="fixed_molecules.csv", help="输出CSV路径")
|
||
parser.add_argument('-j', '--jobs', type=int, default=4, help="并行进程数")
|
||
args = parser.parse_args()
|
||
|
||
main(input_path=args.input, output_path=args.output, n_jobs=args.jobs)
|
||
|
||
'''
|
||
# 查看帮助
|
||
python simemacrocycle_repair.py -h
|
||
|
||
# 运行示例
|
||
python simemacrocycle_repair.py input.smi -o results.csv -j 8
|
||
|
||
python simemacrocycle_repair.py ../data/Macro16_SIME_Synthesis/2025-02-26-05-38-39_mcrl_1.smiles -o ../data/Macro16_SIME_Synthesis/fixed_macrolides_2025.csv -j 8
|
||
''' |