#!/usr/bin/env python3 """ 碎片分析脚本 分析指定环大小的分子碎片： - 统计每个位置的碎片频率 - 分析碎片理化性质 - 找出裂解最多的位置（前5） - 找出变化最大的位置 """ import sys from pathlib import Path # 添加项目根目录到 Python 路径 script_dir = Path(__file__).resolve().parent project_root = script_dir.parent sys.path.insert(0, str(project_root)) import json import pandas as pd import numpy as np from collections import Counter, defaultdict from rdkit import Chem from rdkit.Chem import Descriptors, rdMolDescriptors from rdkit import DataStructs from rdkit.Chem import AllChem import matplotlib.pyplot as plt import seaborn as sns from scipy import stats as scipy_stats def load_fragments_from_ring_size(ring_size): """加载指定环大小的所有碎片""" output_dir = Path(f'output/ring{ring_size}_fragments') if not output_dir.exists(): print(f"❌ 未找到 {ring_size}元环的输出文件夹") return None all_fragments = [] for mol_dir in output_dir.iterdir(): if mol_dir.is_dir() and mol_dir.name.startswith(f'ring{ring_size}_mol_'): fragments_file = mol_dir / f"{mol_dir.name}_all_fragments.json" if fragments_file.exists(): with open(fragments_file, 'r', encoding='utf-8') as f: data = json.load(f) all_fragments.extend(data['fragments']) return all_fragments def calculate_fragment_properties(smiles): """计算碎片的理化性质""" mol = Chem.MolFromSmiles(smiles) if mol is None: return None properties = { 'smiles': smiles, 'molecular_weight': Descriptors.MolWt(mol), 'num_atoms': mol.GetNumAtoms(), 'num_heavy_atoms': mol.GetNumHeavyAtoms(), 'num_rotatable_bonds': Descriptors.NumRotatableBonds(mol), 'num_hbd': Descriptors.NumHDonors(mol), # 氢键供体 'num_hba': Descriptors.NumHAcceptors(mol), # 氢键受体 'logp': Descriptors.MolLogP(mol), 'tpsa': Descriptors.TPSA(mol), # 极性表面积 'num_rings': Descriptors.RingCount(mol), 'num_aromatic_rings': Descriptors.NumAromaticRings(mol), } return properties def calculate_shannon_entropy(fragment_list): """ 计算香农熵，用于衡量碎片多样性 熵越高，表示多样性越大 """ if not fragment_list: return 0 counts = Counter(fragment_list) total = len(fragment_list) entropy = 0 for count in counts.values(): p = count / total if p > 0: entropy -= p * np.log2(p) return entropy def calculate_structural_diversity(smiles_list): """ 计算结构多样性 使用Tanimoto相似度的平均值 相似度越低，多样性越高 """ if len(smiles_list) < 2: return 0 # 生成指纹 fps = [] for smiles in smiles_list: mol = Chem.MolFromSmiles(smiles) if mol: fp = AllChem.GetMorganFingerprintAsBitVect(mol, 2, nBits=1024) fps.append(fp) if len(fps) < 2: return 0 # 计算所有配对的Tanimoto相似度 similarities = [] for i in range(len(fps)): for j in range(i+1, len(fps)): sim = DataStructs.TanimotoSimilarity(fps[i], fps[j]) similarities.append(sim) # 返回平均相似度（1-相似度 = 多样性） avg_similarity = np.mean(similarities) diversity = 1 - avg_similarity return diversity def analyze_position_statistics(fragments, ring_size): """分析每个位置的统计信息""" # 按位置分组 position_data = defaultdict(list) for frag in fragments: pos = frag['cleavage_position'] position_data[pos].append(frag) # 统计信息 results = [] for pos in sorted(position_data.keys()): frags = position_data[pos] smiles_list = [f['fragment_smiles'] for f in frags] # 基础统计 total_count = len(smiles_list) unique_count = len(set(smiles_list)) most_common = Counter(smiles_list).most_common(3) # 多样性指标 shannon_entropy = calculate_shannon_entropy(smiles_list) structural_diversity = calculate_structural_diversity(list(set(smiles_list))) # 理化性质统计 props_list = [] for smiles in set(smiles_list): props = calculate_fragment_properties(smiles) if props: props_list.append(props) if props_list: mw_mean = np.mean([p['molecular_weight'] for p in props_list]) mw_std = np.std([p['molecular_weight'] for p in props_list]) atom_mean = np.mean([p['num_atoms'] for p in props_list]) atom_std = np.std([p['num_atoms'] for p in props_list]) else: mw_mean = mw_std = atom_mean = atom_std = 0 results.append({ 'position': pos, 'total_count': total_count, 'unique_count': unique_count, 'uniqueness_ratio': unique_count / total_count if total_count > 0 else 0, 'shannon_entropy': shannon_entropy, 'structural_diversity': structural_diversity, 'mw_mean': mw_mean, 'mw_std': mw_std, 'atom_mean': atom_mean, 'atom_std': atom_std, 'most_common': most_common }) return pd.DataFrame(results) def plot_position_analysis(df_stats, ring_size, output_dir): """绘制位置分析图表""" fig, axes = plt.subplots(2, 3, figsize=(18, 12)) fig.suptitle(f'{ring_size}元环碎片位置分析', fontsize=16, fontweight='bold') positions = df_stats['position'].values # 1. 碎片总数分布 ax = axes[0, 0] bars = ax.bar(positions, df_stats['total_count'], color='steelblue', edgecolor='black', alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('碎片总数', fontsize=12) ax.set_title('各位置碎片总数', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 标注前5 top5_indices = df_stats.nlargest(5, 'total_count').index for idx in top5_indices: pos = df_stats.loc[idx, 'position'] count = df_stats.loc[idx, 'total_count'] ax.text(pos, count, str(int(count)), ha='center', va='bottom', fontsize=9, fontweight='bold') # 2. 独特碎片数 ax = axes[0, 1] ax.bar(positions, df_stats['unique_count'], color='coral', edgecolor='black', alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('独特碎片数', fontsize=12) ax.set_title('各位置独特碎片数', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 3. 香农熵（多样性） ax = axes[0, 2] bars = ax.bar(positions, df_stats['shannon_entropy'], color='mediumseagreen', edgecolor='black', alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('香农熵', fontsize=12) ax.set_title('各位置碎片多样性（香农熵）', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 标注熵最高的位置 top_entropy_idx = df_stats['shannon_entropy'].idxmax() pos = df_stats.loc[top_entropy_idx, 'position'] entropy = df_stats.loc[top_entropy_idx, 'shannon_entropy'] ax.text(pos, entropy, f'{entropy:.2f}', ha='center', va='bottom', fontsize=9, fontweight='bold') # 4. 结构多样性 ax = axes[1, 0] ax.bar(positions, df_stats['structural_diversity'], color='mediumpurple', edgecolor='black', alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('结构多样性', fontsize=12) ax.set_title('各位置结构多样性（指纹相似度）', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 5. 分子量分布 ax = axes[1, 1] ax.errorbar(positions, df_stats['mw_mean'], yerr=df_stats['mw_std'], fmt='o-', color='darkred', ecolor='lightcoral', capsize=5, capthick=2, alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('分子量 (Da)', fontsize=12) ax.set_title('各位置碎片平均分子量', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 6. 原子数分布 ax = axes[1, 2] ax.errorbar(positions, df_stats['atom_mean'], yerr=df_stats['atom_std'], fmt='s-', color='darkgreen', ecolor='lightgreen', capsize=5, capthick=2, alpha=0.7) ax.set_xlabel('环位置', fontsize=12) ax.set_ylabel('原子数', fontsize=12) ax.set_title('各位置碎片平均原子数', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) plt.tight_layout() # 保存图表 plot_path = output_dir / f'ring{ring_size}_position_analysis.png' plt.savefig(plot_path, dpi=300, bbox_inches='tight') print(f"✓ 图表已保存: {plot_path}") plt.show() def plot_fragment_properties_distribution(fragments, ring_size, output_dir): """绘制碎片理化性质分布""" # 计算所有碎片的性质 all_props = [] unique_smiles = set([f['fragment_smiles'] for f in fragments]) for smiles in unique_smiles: props = calculate_fragment_properties(smiles) if props: all_props.append(props) df_props = pd.DataFrame(all_props) fig, axes = plt.subplots(2, 3, figsize=(18, 12)) fig.suptitle(f'{ring_size}元环碎片理化性质分布', fontsize=16, fontweight='bold') # 1. 分子量分布 ax = axes[0, 0] ax.hist(df_props['molecular_weight'], bins=30, color='steelblue', edgecolor='black', alpha=0.7) ax.set_xlabel('分子量 (Da)', fontsize=12) ax.set_ylabel('频数', fontsize=12) ax.set_title('分子量分布', fontsize=14, fontweight='bold') ax.axvline(df_props['molecular_weight'].median(), color='red', linestyle='--', label=f'中位数: {df_props["molecular_weight"].median():.1f}') ax.legend() ax.grid(axis='y', alpha=0.3) # 2. 原子数分布 ax = axes[0, 1] ax.hist(df_props['num_atoms'], bins=range(1, int(df_props['num_atoms'].max())+2), color='coral', edgecolor='black', alpha=0.7) ax.set_xlabel('原子数', fontsize=12) ax.set_ylabel('频数', fontsize=12) ax.set_title('原子数分布', fontsize=14, fontweight='bold') ax.axvline(df_props['num_atoms'].median(), color='red', linestyle='--', label=f'中位数: {df_props["num_atoms"].median():.0f}') ax.legend() ax.grid(axis='y', alpha=0.3) # 3. LogP分布 ax = axes[0, 2] ax.hist(df_props['logp'], bins=30, color='mediumseagreen', edgecolor='black', alpha=0.7) ax.set_xlabel('LogP', fontsize=12) ax.set_ylabel('频数', fontsize=12) ax.set_title('LogP分布（亲脂性）', fontsize=14, fontweight='bold') ax.axvline(df_props['logp'].median(), color='red', linestyle='--', label=f'中位数: {df_props["logp"].median():.2f}') ax.legend() ax.grid(axis='y', alpha=0.3) # 4. TPSA分布 ax = axes[1, 0] ax.hist(df_props['tpsa'], bins=30, color='mediumpurple', edgecolor='black', alpha=0.7) ax.set_xlabel('TPSA (Ų)', fontsize=12) ax.set_ylabel('频数', fontsize=12) ax.set_title('极性表面积分布', fontsize=14, fontweight='bold') ax.axvline(df_props['tpsa'].median(), color='red', linestyle='--', label=f'中位数: {df_props["tpsa"].median():.1f}') ax.legend() ax.grid(axis='y', alpha=0.3) # 5. 可旋转键数分布 ax = axes[1, 1] ax.hist(df_props['num_rotatable_bonds'], bins=range(int(df_props['num_rotatable_bonds'].max())+2), color='darkgoldenrod', edgecolor='black', alpha=0.7) ax.set_xlabel('可旋转键数', fontsize=12) ax.set_ylabel('频数', fontsize=12) ax.set_title('可旋转键数分布（柔性）', fontsize=14, fontweight='bold') ax.grid(axis='y', alpha=0.3) # 6. 氢键统计 ax = axes[1, 2] width = 0.35 x = np.arange(2) hbd_mean = df_props['num_hbd'].mean() hba_mean = df_props['num_hba'].mean() ax.bar(x[0]-width/2, hbd_mean, width, label='氢键供体 (HBD)', color='lightcoral', edgecolor='black') ax.bar(x[1]+width/2, hba_mean, width, label='氢键受体 (HBA)', color='lightblue', edgecolor='black') ax.set_ylabel('平均数量', fontsize=12) ax.set_title('氢键能力统计', fontsize=14, fontweight='bold') ax.set_xticks(x) ax.set_xticklabels(['HBD', 'HBA']) ax.legend() ax.grid(axis='y', alpha=0.3) plt.tight_layout() # 保存图表 plot_path = output_dir / f'ring{ring_size}_properties_distribution.png' plt.savefig(plot_path, dpi=300, bbox_inches='tight') print(f"✓ 图表已保存: {plot_path}") plt.show() return df_props def generate_analysis_report(ring_size): """生成完整的分析报告""" print("="*80) print(f"{ring_size}元环碎片分析") print("="*80) print() # 加载碎片 print("加载碎片数据...") fragments = load_fragments_from_ring_size(ring_size) if fragments is None or len(fragments) == 0: print(f"❌ 没有找到碎片数据") return print(f"✓ 加载了 {len(fragments)} 个碎片") print() # 创建输出文件夹 output_dir = Path(f'output/ring{ring_size}_analysis') output_dir.mkdir(parents=True, exist_ok=True) # 1. 位置统计分析 print("分析各位置统计信息...") df_stats = analyze_position_statistics(fragments, ring_size) # 保存统计表格 stats_csv = output_dir / f'ring{ring_size}_position_statistics.csv' df_stats.to_csv(stats_csv, index=False, encoding='utf-8') print(f"✓ 统计表格已保存: {stats_csv}") print() # 2. 打印关键结果 print("="*80) print("关键发现") print("="*80) print() # 裂解最多的前5个位置 top5_positions = df_stats.nlargest(5, 'total_count') print("📊 裂解最多的前5个位置:") for i, row in enumerate(top5_positions.itertuples(), 1): print(f" {i}. 位置 {row.position}: {row.total_count} 个碎片 " f"({row.unique_count} 种独特结构)") print() # 多样性最大的前5个位置（基于香农熵） top5_diversity = df_stats.nlargest(5, 'shannon_entropy') print("🌈 多样性最大的前5个位置（香农熵）:") for i, row in enumerate(top5_diversity.itertuples(), 1): print(f" {i}. 位置 {row.position}: 熵={row.shannon_entropy:.3f}, " f"独特比={row.uniqueness_ratio:.2%}, " f"{row.unique_count}/{row.total_count} 独特") print() # 结构多样性最大的前5个位置 top5_structural = df_stats.nlargest(5, 'structural_diversity') print("🔬 结构多样性最大的前5个位置（指纹相似度）:") for i, row in enumerate(top5_structural.itertuples(), 1): print(f" {i}. 位置 {row.position}: 多样性={row.structural_diversity:.3f}, " f"{row.unique_count} 种结构") print() # 碎片大小变化最大的位置 top5_mw_variance = df_stats.nlargest(5, 'mw_std') print("⚖️ 分子量变化最大的前5个位置:") for i, row in enumerate(top5_mw_variance.itertuples(), 1): print(f" {i}. 位置 {row.position}: 平均={row.mw_mean:.1f} Da, " f"标准差={row.mw_std:.1f} Da") print() # 3. 绘制图表 print("生成可视化图表...") plot_position_analysis(df_stats, ring_size, output_dir) print() # 4. 理化性质分析 print("分析碎片理化性质...") df_props = plot_fragment_properties_distribution(fragments, ring_size, output_dir) # 保存理化性质表格 props_csv = output_dir / f'ring{ring_size}_fragment_properties.csv' df_props.to_csv(props_csv, index=False, encoding='utf-8') print(f"✓ 理化性质表格已保存: {props_csv}") print() # 5. 总结统计 print("="*80) print("总体统计") print("="*80) print(f"总碎片数: {len(fragments)}") print(f"独特碎片数: {len(set([f['fragment_smiles'] for f in fragments]))}") print(f"平均分子量: {df_props['molecular_weight'].mean():.2f} ± {df_props['molecular_weight'].std():.2f} Da") print(f"平均原子数: {df_props['num_atoms'].mean():.2f} ± {df_props['num_atoms'].std():.2f}") print(f"分子量范围: {df_props['molecular_weight'].min():.1f} - {df_props['molecular_weight'].max():.1f} Da") print(f"原子数范围: {int(df_props['num_atoms'].min())} - {int(df_props['num_atoms'].max())}") print() # 6. 生成文本报告 report_path = output_dir / f'ring{ring_size}_analysis_report.txt' with open(report_path, 'w', encoding='utf-8') as f: f.write(f"{ring_size}元环碎片分析报告\n") f.write("="*80 + "\n\n") f.write(f"生成时间: {pd.Timestamp.now()}\n\n") f.write("位置统计摘要:\n") f.write(df_stats.to_string()) f.write("\n\n") f.write("理化性质摘要:\n") f.write(df_props.describe().to_string()) f.write("\n") print(f"✓ 完整报告已保存: {report_path}") print() print("="*80) print("分析完成！") print("="*80) if __name__ == '__main__': import argparse parser = argparse.ArgumentParser(description='分析大环内酯碎片') parser.add_argument('ring_size', type=int, help='环大小（如：16）') args = parser.parse_args() generate_analysis_report(args.ring_size)