OSCR

Drug screening for α-synuclein aggregation inhibitors via multimodal graph neural network.

Code ↔ Paper

6 matches between paragraphs of the paper and lines of its authors' code, computed by the harvester (lexical-v1). Click a colored paragraph or line to see its counterpart.

The 6 matches
  1. [1] § Results and discussion › Comparative analysis of diverse machine learning models ↔ pycaret.py, lines 243–302 · score 0.91 · gradient boosting, CatBoost, LightGBM, XGBoost, Extra Trees, linear regression
  2. [2] § Results and discussion › Interpretability analysis of drug molecules ↔ gnnexplainer.py, lines 115–232 · score 0.66 · pyridine rings, aromatic rings, Amide, methylene, amino, hydroxyl
  3. [3] § Results and discussion › Interpretability analysis of drug molecules ↔ gnnexplainer.py, lines 115–232 · score 0.63 · Aromatic rings, vinyl, Ether, amides, methylene, amino
  4. [4] § Results and discussion › Superiority of GraphSAGE over baseline Graph Neural Network (GNN) models ↔ statistical analysis.py, lines 135–232 · score 0.61 · odds ratio, log scaled, molecular descriptors odds, Forest
  5. [5] § Results and discussion › Distribution and differential analysis of molecular descriptors for high-/low-affinity ligands ↔ statistical analysis.py, lines 48–74 · score 0.59 · pChEMBL, aromatic rings, median, donors, affinity, bond
  6. [6] § Results and discussion › Interpretability analysis of drug molecules ↔ gnnexplainer.py, lines 1468–1512 · score 0.51 · stratified sampling, representative molecules, GNNExplainer, atom, predictive

Paper

Loaded from Europe PMC by your browser, not stored by OSCR: doi.org · Europe PMC

The paper is loaded when this pane is shown.

The authors' code

Python · 1,896 lines · 79 KB · no license · 3 matches

  1. import numpy as np
  2. import pandas as pd
  3. import torch
  4. import torch.nn as nn
  5. import torch.nn.functional as F
  6. from torch_geometric.data import Data, DataLoader
  7. from torch_geometric.explain import Explainer, GNNExplainer
  8. from torch_geometric.explain.config import ExplainerConfig, ModelConfig
  9. from torch_geometric.nn import GATConv, SAGEConv, global_max_pool
  10. import matplotlib.pyplot as plt
  11. import seaborn as sns
  12. import networkx as nx
  13. from rdkit import Chem
  14. from rdkit.Chem import Draw, AllChem, Descriptors, rdMolDescriptors
  15. from rdkit.Chem.Draw import rdDepictor
  16. import io
  17. import base64
  18. from PIL import Image
  19. from collections import defaultdict, Counter
  20. import warnings
  22. plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'Arial', 'sans-serif']
  23. plt.rcParams['axes.unicode_minus'] = False
  24. plt.rcParams['font.size'] = 14
  25. plt.rcParams['axes.titlesize'] = 18
  26. plt.rcParams['axes.labelsize'] = 16
  27. plt.rcParams['xtick.labelsize'] = 14
  28. plt.rcParams['ytick.labelsize'] = 14
  29. plt.rcParams['legend.fontsize'] = 14
  31. try:
  32. from rdkit.Chem.Draw import rdMolDraw2D
  33. except ImportError:
  34. try:
  35. from rdkit.Chem import rdMolDraw2D
  36. except ImportError:
  37. rdMolDraw2D = None
  38. print("Warning: rdMolDraw2D not available")
  40. warnings.filterwarnings('ignore')
  42. seed = 42
  43. torch.manual_seed(seed)
  44. np.random.seed(seed)
  47. def one_of_k_encoding_unk(x, valid_entries):
  48. if x not in valid_entries:
  49. x = 'Unknown'
  50. return [1 if entry == x else 0 for entry in valid_entries]
  53. class ModifiedGATLayer(nn.Module):
  54. def __init__(self, in_features, out_features):
  55. super(ModifiedGATLayer, self).__init__()
  56. self.query_transform = nn.Linear(in_features, out_features)
  57. self.key_transform = nn.Linear(in_features, out_features)
  58. self.value_transform = nn.Linear(in_features, out_features)
  59. self.conv3 = nn.Conv1d(in_channels=out_features, out_channels=out_features, kernel_size=3, padding=1)
  60. self.conv5 = nn.Conv1d(in_channels=out_features, out_channels=out_features, kernel_size=5, padding=2)
  61. self.linear_transform = nn.Linear(out_features * 2 + out_features, out_features)
  63. def forward(self, x):
  64. Q = self.query_transform(x)
  65. K = self.key_transform(x)
  66. V = self.value_transform(x)
  67. K = K.unsqueeze(2)
  68. K_conv3 = self.conv3(K)
  69. K_conv5 = self.conv5(K)
  70. K_concat = torch.cat((K_conv3, K_conv5, K), dim=1)
  71. K_new = self.linear_transform(K_concat.transpose(1, 2))
  72. attention_scores = torch.matmul(Q, K_new.transpose(1, 2)) / (K_new.size(-1) ** 0.5)
  73. attention_weights = F.softmax(attention_scores.squeeze(-1), dim=-1)
  74. output = torch.matmul(attention_weights, V) + V
  75. return output
  78. class GAT_GraphSAGE(nn.Module):
  79. def __init__(self, n_output=1, num_features_xd=35, output_dim=128, dropout=0.3):
  80. super(GAT_GraphSAGE, self).__init__()
  81. self.conv1 = ModifiedGATLayer(in_features=num_features_xd, out_features=num_features_xd)
  82. self.conv2 = SAGEConv(num_features_xd, num_features_xd)
  83. self.fc_g1 = nn.Linear(num_features_xd, 1500)
  84. self.fc_g2 = nn.Linear(1500, output_dim)
  85. self.relu = nn.ReLU()
  86. self.dropout = nn.Dropout(dropout)
  87. self.out = nn.Linear(output_dim, n_output)
  89. def forward(self, data):
  90. x, edge_index, batch = data.x, data.edge_index, data.batch
  91. x = self.conv1(x)
  92. x = self.relu(x)
  93. x = self.conv2(x, edge_index)
  94. x = self.relu(x)
  95. x = global_max_pool(x, batch)
  96. x = self.relu(self.fc_g1(x))
  97. x = self.dropout(x)
  98. x = self.fc_g2(x)
  99. out = self.out(x)
  100. return out
  103. class ExplainableGATGraphSAGE(nn.Module):
  104. def __init__(self, gat_graphsage_model):
  105. super(ExplainableGATGraphSAGE, self).__init__()
  106. self.gat_graphsage = gat_graphsage_model
  108. def forward(self, x, edge_index, batch=None, edge_attr=None):
  109. if batch is None:
  110. batch = torch.zeros(x.size(0), dtype=torch.long, device=x.device)
  111. data = Data(x=x, edge_index=edge_index, batch=batch)
  112. return self.gat_graphsage(data)
  115. class SubstructureIdentifier:
  116. def __init__(self):
  117. self.common_substructures = {
  118. 'hydroxyl': 'O',
  119. 'amino': 'N',
  120. 'carboxyl': 'C(=O)O',
  121. 'carbonyl': 'C=O',
  122. 'ester': 'C(=O)O[C,c]',
  123. 'amide': 'C(=O)N',
  124. 'ether': '[C,c]O[C,c]',
  125. 'nitro': 'N(=O)=O',
  126. 'sulfonyl': 'S(=O)(=O)',
  127. 'phosphate': 'P(=O)',
  129. 'benzene': 'c1ccccc1',
  130. 'pyridine': 'c1ccncc1',
  131. 'pyrimidine': 'c1cncnc1',
  132. 'imidazole': 'c1c[nH]cn1',
  133. 'thiophene': 'c1ccsc1',
  134. 'furan': 'c1ccoc1',
  135. 'indole': 'c1ccc2[nH]ccc2c1',
  136. 'quinoline': 'c1ccc2ncccc2c1',
  138. 'piperidine': 'C1CCNCC1',
  139. 'piperazine': 'C1CNCCN1',
  140. 'morpholine': 'C1COCCN1',
  141. 'pyrrolidine': 'C1CCNC1',
  142. 'tetrahydrofuran': 'C1CCOC1',
  144. 'methylene': 'CC',
  145. 'ethylene': 'CCC',
  146. 'propylene': 'CCCC',
  147. 'vinyl': 'C=C',
  148. 'acetylene': 'C#C',
  149. }
  151. def identify_substructures_in_molecule(self, mol):
  152. if mol is None:
  153. return {}
  155. found_substructures = {}
  157. for name, smarts in self.common_substructures.items():
  158. try:
  159. matches = mol.GetSubstructMatches(Chem.MolFromSmarts(smarts))
  160. if matches:
  161. found_substructures[name] = {
  162. 'pattern': smarts,
  163. 'matches': matches,
  164. 'count': len(matches)
  165. }
  166. except:
  167. continue
  169. return found_substructures
  171. def extract_important_substructures(self, mol, important_atom_indices, radius=2):
  172. if mol is None:
  173. return []
  175. important_substructures = []
  177. for atom_idx in important_atom_indices:
  178. if atom_idx >= mol.GetNumAtoms():
  179. continue
  181. env = Chem.FindAtomEnvironmentOfRadiusN(mol, radius, atom_idx)
  182. submol = Chem.PathToSubmol(mol, env)
  184. if submol:
  185. try:
  186. sub_smiles = Chem.MolToSmiles(submol)
  187. if sub_smiles and len(sub_smiles) > 1:
  188. important_substructures.append({
  189. 'center_atom': atom_idx,
  190. 'smiles': sub_smiles,
  191. 'num_atoms': submol.GetNumAtoms(),
  192. 'radius': radius
  193. })
  194. except:
  195. continue
  197. return important_substructures
  199. def get_functional_groups(self, mol):
  200. if mol is None:
  201. return []
  203. functional_groups = []
  205. try:
  206. from rdkit.Chem import Descriptors, Fragments
  208. groups_to_check = {
  209. 'Aromatic_Rings': Descriptors.NumAromaticRings(mol),
  210. 'Aliphatic_Rings': Descriptors.NumAliphaticRings(mol),
  211. 'Hydroxyl_Groups': Fragments.fr_Al_OH(mol),
  212. 'Carboxyl_Groups': Fragments.fr_COO(mol),
  213. 'Amino_Groups': Fragments.fr_NH2(mol),
  214. 'Ester_Groups': Fragments.fr_ester(mol),
  215. 'Ether_Groups': Fragments.fr_ether(mol),
  216. 'Amide_Groups': Fragments.fr_amide(mol),
  217. 'Nitro_Groups': Fragments.fr_nitro(mol),
  218. 'Benzene_Rings': Fragments.fr_benzene(mol),
  219. 'Pyridine_Rings': Fragments.fr_pyridine(mol),
  220. }
  222. for group_name, count in groups_to_check.items():
  223. if count > 0:
  224. functional_groups.append({
  225. 'name': group_name,
  226. 'count': count
  227. })
  229. except Exception as e:
  230. print(f"Error in functional group analysis: {e}")
  232. return functional_groups
  235. class SubstructureVisualizer:
  236. def __init__(self):
  237. plt.rcParams['font.sans-serif'] = ['DejaVu Sans', 'Arial', 'sans-serif']
  238. plt.rcParams['axes.unicode_minus'] = False
  240. def visualize_substructure_importance_summary(self, important_substructures, save_path=None):
  241. all_substructure_types = {}
  242. all_functional_groups = Counter()
  244. for struct in important_substructures:
  245. for sub_name, sub_info in struct['known_substructures'].items():
  246. if sub_name not in all_substructure_types:
  247. all_substructure_types[sub_name] = {
  248. 'count': 0,
  249. 'total_importance': 0,
  250. 'molecules': []
  251. }
  252. all_substructure_types[sub_name]['count'] += len(sub_info['matches'])
  253. all_substructure_types[sub_name]['total_importance'] += sum(sub_info['importance_scores'])
  254. all_substructure_types[sub_name]['molecules'].append(struct['smiles'][:30])
  256. for fg in struct['functional_groups']:
  257. all_functional_groups[fg['name']] += fg['count']
  259. if not all_substructure_types:
  260. print("No substructure data available for visualization")
  261. return
  263. sorted_substructures = sorted(all_substructure_types.items(),
  264. key=lambda x: x[1]['count'], reverse=True)[:15]
  266. sub_names = [item[0] for item in sorted_substructures]
  267. sub_counts = [item[1]['count'] for item in sorted_substructures]
  269. plt.figure(figsize=(14, 10))
  270. bars1 = plt.barh(sub_names, sub_counts, color='skyblue', alpha=0.8)
  271. plt.title('Substructure Frequency (Top 15)', fontsize=20, fontweight='bold', pad=25)
  272. plt.xlabel('Frequency Count', fontsize=18, fontweight='bold')
  273. plt.ylabel('Substructure Type', fontsize=18, fontweight='bold')
  274. plt.xticks(fontsize=16)
  275. plt.yticks(fontsize=16)
  276. plt.grid(axis='x', alpha=0.3)
  278. for bar, count in zip(bars1, sub_counts):
  279. plt.text(bar.get_width() + max(sub_counts) * 0.01, bar.get_y() + bar.get_height() / 2,
  280. str(count), ha='left', va='center', fontsize=14, fontweight='bold')
  282. plt.subplots_adjust(left=0.25, right=0.95, top=0.9, bottom=0.15)
  283. if save_path:
  284. plt.savefig(f"{save_path}_subplot1_frequency.png", dpi=300, bbox_inches='tight')
  285. print(f"Substructure subplot 1 saved: {save_path}_subplot1_frequency.png")
  286. plt.show()
  288. avg_importance = []
  289. for item in sorted_substructures:
  290. count = item[1]['count']
  291. total_imp = item[1]['total_importance']
  292. avg_imp = total_imp / count if count > 0 else 0
  293. avg_importance.append(avg_imp)
  295. plt.figure(figsize=(14, 10))
  296. bars2 = plt.barh(sub_names, avg_importance, color='lightcoral', alpha=0.8)
  297. plt.title('Average Substructure Importance (Top 15)', fontsize=20, fontweight='bold', pad=25)
  298. plt.xlabel('Average Importance Score', fontsize=18, fontweight='bold')
  299. plt.ylabel('Substructure Type', fontsize=18, fontweight='bold')
  300. plt.xticks(fontsize=16)
  301. plt.yticks(fontsize=16)
  302. plt.grid(axis='x', alpha=0.3)
  304. for bar, imp in zip(bars2, avg_importance):
  305. plt.text(bar.get_width() + max(avg_importance) * 0.01, bar.get_y() + bar.get_height() / 2,
  306. f'{imp:.3f}', ha='left', va='center', fontsize=14, fontweight='bold')
  308. plt.subplots_adjust(left=0.25, right=0.95, top=0.9, bottom=0.15)
  309. if save_path:
  310. plt.savefig(f"{save_path}_subplot2_importance.png", dpi=300, bbox_inches='tight')
  311. print(f"Substructure subplot 2 saved: {save_path}_subplot2_importance.png")
  312. plt.show()
  314. top_fg = all_functional_groups.most_common(12)
  315. fg_names = [item[0].replace('_', ' ') for item in top_fg]
  316. fg_counts = [item[1] for item in top_fg]
  318. plt.figure(figsize=(14, 12))
  319. colors = plt.cm.Set3(np.linspace(0, 1, len(fg_names)))
  321. def autopct_func_fg(pct, idx=0):
  322. if idx < len(fg_counts) - 4:
  323. return f'{pct:.1f}%'
  324. else:
  325. return ''
  327. counter = [0]
  329. def autopct_with_counter(pct):
  330. result = autopct_func_fg(pct, counter[0])
  331. counter[0] += 1
  332. return result
  334. wedges, texts, autotexts = plt.pie(fg_counts, autopct=autopct_with_counter,
  335. colors=colors, startangle=90,
  336. textprops={'fontsize': 14, 'fontweight': 'bold'})
  337. plt.title('Functional Group Distribution', fontsize=20, fontweight='bold', pad=25)
  339. for i, autotext in enumerate(autotexts):
  340. autotext.set_fontsize(14)
  341. autotext.set_color('black')
  342. autotext.set_fontweight('bold')
  344. plt.legend(wedges, fg_names, title="Functional Groups",
  345. title_fontsize=16, fontsize=14,
  346. loc="center left", bbox_to_anchor=(1.1, 0.5))
  348. plt.subplots_adjust(left=0.1, right=0.75, top=0.9, bottom=0.1)
  349. if save_path:
  350. plt.savefig(f"{save_path}_subplot3_functional_groups.png", dpi=300, bbox_inches='tight')
  351. print(f"Substructure subplot 3 saved: {save_path}_subplot3_functional_groups.png")
  352. plt.show()
  354. plt.figure(figsize=(14, 10))
  356. x_counts = [item[1]['count'] for item in sorted_substructures[:20]]
  357. y_importance = []
  358. colors_scatter = []
  359. labels = []
  361. for item in sorted_substructures[:20]:
  362. count = item[1]['count']
  363. total_imp = item[1]['total_importance']
  364. avg_imp = total_imp / count if count > 0 else 0
  365. y_importance.append(avg_imp)
  367. if avg_imp > 0.6:
  368. colors_scatter.append('red')
  369. elif avg_imp > 0.4:
  370. colors_scatter.append('orange')
  371. else:
  372. colors_scatter.append('blue')
  374. labels.append(item[0])
  376. scatter = plt.scatter(x_counts, y_importance, c=colors_scatter,
  377. s=150, alpha=0.7, edgecolors='black')
  379. for i, label in enumerate(labels):
  380. if y_importance[i] > 0.5 or x_counts[i] > 200:
  381. plt.annotate(label, (x_counts[i], y_importance[i]),
  382. xytext=(5, 5), textcoords='offset points',
  383. fontsize=12, alpha=0.8, fontweight='bold')
  385. plt.title('Substructure Importance vs Frequency', fontsize=20, fontweight='bold', pad=25)
  386. plt.xlabel('Frequency Count', fontsize=18, fontweight='bold')
  387. plt.ylabel('Average Importance', fontsize=18, fontweight='bold')
  388. plt.xticks(fontsize=16)
  389. plt.yticks(fontsize=16)
  390. plt.grid(True, alpha=0.3)
  392. from matplotlib.patches import Patch
  393. legend_elements = [Patch(facecolor='red', alpha=0.7, label='High Importance (>0.6)'),
  394. Patch(facecolor='orange', alpha=0.7, label='Medium Importance (0.4-0.6)'),
  395. Patch(facecolor='blue', alpha=0.7, label='Low Importance (<0.4)')]
  396. plt.legend(handles=legend_elements, loc='upper right', fontsize=14)
  398. plt.subplots_adjust(left=0.15, right=0.95, top=0.9, bottom=0.15)
  399. if save_path:
  400. plt.savefig(f"{save_path}_subplot4_scatter.png", dpi=300, bbox_inches='tight')
  401. print(f"Substructure subplot 4 saved: {save_path}_subplot4_scatter.png")
  402. plt.show()
  404. print(f"\nAll four substructure subplots have been generated and saved separately!")
  406. def visualize_substructure_in_molecules(self, full_dataset_substructures, num_examples=6, save_path=None):
  407. if not full_dataset_substructures:
  408. print("No important substructures available for visualization")
  409. return
  411. print(f"Searching for qualified molecules from {len(full_dataset_substructures)} molecules in full dataset...")
  413. selected_molecules = []
  415. for struct in full_dataset_substructures:
  416. if 'target_value' not in struct or struct['target_value'] <= 6:
  417. continue
  419. has_high_importance_substructure = False
  420. for sub_name, sub_info in struct['known_substructures'].items():
  421. if sub_info['importance_scores'] and max(sub_info['importance_scores']) > 0.5:
  422. has_high_importance_substructure = True
  423. break
  425. if has_high_importance_substructure:
  426. selected_molecules.append(struct)
  427. if len(selected_molecules) >= num_examples:
  428. break
  430. print(f"Found {len(selected_molecules)} molecules with y > 6 and high importance substructures (> 0.5)")
  432. if not selected_molecules:
  433. print("No molecules found with y > 6 and high importance substructures (> 0.5)")
  434. print("Falling back to molecules with any high importance substructures...")
  435. for struct in full_dataset_substructures:
  436. if len(selected_molecules) >= num_examples:
  437. break
  439. has_high_importance_substructure = False
  440. for sub_name, sub_info in struct['known_substructures'].items():
  441. if sub_info['importance_scores'] and max(sub_info['importance_scores']) > 0.5:
  442. has_high_importance_substructure = True
  443. break
  445. if has_high_importance_substructure:
  446. selected_molecules.append(struct)
  448. if not selected_molecules:
  449. print("No molecules with high importance substructures found")
  450. selected_molecules = full_dataset_substructures[:num_examples]
  452. fig, axes = plt.subplots(2, 3, figsize=(20, 14))
  453. axes = axes.flatten()
  455. for i, struct in enumerate(selected_molecules[:num_examples]):
  456. try:
  457. mol = Chem.MolFromSmiles(struct['smiles'])
  458. if mol is None:
  459. continue
  461. best_substructure = None
  462. best_importance = 0
  464. for sub_name, sub_info in struct['known_substructures'].items():
  465. if sub_info['importance_scores']:
  466. avg_importance = np.mean(sub_info['importance_scores'])
  467. if avg_importance > best_importance:
  468. best_importance = avg_importance
  469. best_substructure = (sub_name, sub_info)
  471. if best_substructure is None:
  472. img = Draw.MolToImage(mol, size=(400, 400))
  473. axes[i].imshow(img)
  474. y_value = struct.get('target_value', 'N/A')
  475. axes[i].set_title(f'Molecule {i + 1}\nTrue y: {y_value}' +
  476. (f' (y={y_value:.3f})' if isinstance(y_value, (int, float)) else ''),
  477. fontsize=16, fontweight='bold', pad=10)
  478. else:
  479. sub_name, sub_info = best_substructure
  480. highlight_atoms = []
  482. for match in sub_info['matches']:
  483. highlight_atoms.extend(match)
  485. highlight_atoms = list(set(highlight_atoms))
  487. img = Draw.MolToImage(mol, size=(400, 400),
  488. highlightAtoms=highlight_atoms,
  489. highlightColors={atom: (1, 0.8, 0.8) for atom in highlight_atoms})
  491. axes[i].imshow(img)
  492. y_value = struct.get('target_value', 'N/A')
  493. title_text = f'Molecule {i + 1}: {sub_name}\n'
  494. if isinstance(y_value, (int, float)):
  495. title_text += f'True y: {y_value:.3f}\n'
  496. else:
  497. title_text += f'True y: {y_value}\n'
  498. title_text += f'Importance: {best_importance:.3f}'
  499. axes[i].set_title(title_text, fontsize=16, fontweight='bold', pad=10)
  501. axes[i].axis('off')
  503. except Exception as e:
  504. print(f"Error drawing molecule {i + 1}: {e}")
  505. axes[i].text(0.5, 0.5, f'Molecule {i + 1}\nDrawing Failed',
  506. ha='center', va='center', transform=axes[i].transAxes,
  507. fontsize=16, fontweight='bold')
  508. axes[i].axis('off')
  510. for i in range(len(selected_molecules), len(axes)):
  511. axes[i].axis('off')
  513. plt.suptitle('Important Substructures Highlighted in Molecules (Full Dataset: y > 6, Importance > 0.5)',
  514. fontsize=20, fontweight='bold', y=0.95)
  515. plt.subplots_adjust(top=0.85, bottom=0.05, left=0.05, right=0.95, hspace=0.3, wspace=0.2)
  517. if save_path:
  518. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  519. print(f"Highlighted molecules plot saved: {save_path}")
  521. plt.show()
  523. return selected_molecules
  525. def create_substructure_heatmap(self, important_substructures, save_path=None):
  526. if not important_substructures:
  527. print("No data available for heatmap")
  528. return
  530. all_substructures = set()
  531. for struct in important_substructures:
  532. all_substructures.update(struct['known_substructures'].keys())
  534. all_substructures = sorted(list(all_substructures))
  536. if len(all_substructures) == 0:
  537. print("No substructures found")
  538. return
  540. matrix_data = []
  541. molecule_labels = []
  543. for i, struct in enumerate(important_substructures[:40]):
  544. row = []
  545. molecule_labels.append(f"Mol{i + 1}")
  547. for sub_name in all_substructures:
  548. if sub_name in struct['known_substructures']:
  549. importance_scores = struct['known_substructures'][sub_name]['importance_scores']
  550. avg_importance = np.mean(importance_scores) if importance_scores else 0
  551. row.append(avg_importance)
  552. else:
  553. row.append(0)
  555. matrix_data.append(row)
  557. df_heatmap = pd.DataFrame(matrix_data,
  558. columns=all_substructures,
  559. index=molecule_labels)
  561. df_heatmap = df_heatmap.loc[:, (df_heatmap != 0).any(axis=0)]
  563. if df_heatmap.empty:
  564. print("Insufficient data for heatmap")
  565. return
  567. fig_width = max(18, len(df_heatmap.columns) * 1.2)
  568. fig_height = max(14, len(df_heatmap.index) * 0.45)
  570. plt.figure(figsize=(fig_width, fig_height))
  572. try:
  573. ax = sns.heatmap(df_heatmap,
  574. cmap='Blues',
  575. annot=False,
  576. fmt='.2f',
  577. cbar_kws={'label': 'Importance Score', 'shrink': 0.8},
  578. xticklabels=True,
  579. yticklabels=True,
  580. square=False)
  582. plt.title('Molecule-Substructure Importance Heatmap (Top 40 Molecules)',
  583. fontsize=22, fontweight='bold', pad=30)
  584. plt.xlabel('Substructure Type', fontsize=18, fontweight='bold', labelpad=15)
  585. plt.ylabel('Molecules (Top 40)', fontsize=18, fontweight='bold', labelpad=15)
  587. plt.xticks(rotation=45, ha='right', fontsize=14, fontweight='normal')
  588. plt.yticks(rotation=0, fontsize=14, fontweight='normal')
  590. cbar = ax.collections[0].colorbar
  591. cbar.ax.tick_params(labelsize=16)
  592. cbar.set_label('Importance Score', fontsize=18, fontweight='bold', labelpad=20)
  594. plt.subplots_adjust(left=0.15, bottom=0.25, right=0.85, top=0.85)
  596. except Exception as e:
  597. print(f"Error creating heatmap: {e}")
  599. if save_path:
  600. plt.savefig(save_path, dpi=300, bbox_inches='tight',
  601. facecolor='white', edgecolor='none')
  602. print(f"Substructure heatmap (Top 40 molecules) saved: {save_path}")
  604. plt.show()
  607. class MolecularExplainer:
  608. def __init__(self, model, device='cpu'):
  609. self.device = device
  610. self.model = model.to(device)
  611. self.model.eval()
  613. self.atom_symbols = ['C', 'N', 'O', 'S', 'F', 'P', 'Cl', 'Br', 'I', 'Unknown']
  615. self.substructure_identifier = SubstructureIdentifier()
  617. self.substructure_visualizer = SubstructureVisualizer()
  619. try:
  620. self.explainer = Explainer(
  621. model=self.model,
  622. algorithm=GNNExplainer(epochs=100, lr=0.01),
  623. explanation_type='model',
  624. node_mask_type='attributes',
  625. edge_mask_type='object',
  626. model_config=ModelConfig(
  627. mode='regression',
  628. task_level='graph',
  629. return_type='raw'
  630. )
  631. )
  632. except Exception as e:
  633. print(f"Warning: Could not initialize GNNExplainer: {e}")
  634. self.explainer = None
  636. def get_atom_symbol_from_features(self, atom_features):
  637. atom_idx = torch.argmax(atom_features[:10]).item()
  638. return self.atom_symbols[atom_idx]
  640. def simple_gradient_explanation(self, data):
  641. self.model.eval()
  642. data = data.to(self.device)
  644. if not hasattr(data, 'batch') or data.batch is None:
  645. data.batch = torch.zeros(data.x.size(0), dtype=torch.long, device=self.device)
  647. data.x.requires_grad_(True)
  648. prediction = self.model(data.x, data.edge_index, data.batch)
  650. prediction.backward()
  652. node_importance = torch.norm(data.x.grad, dim=1)
  654. return {
  655. 'node_mask': node_importance.cpu().detach(),
  656. 'edge_mask': None,
  657. 'prediction': prediction.cpu().detach(),
  658. 'data': data.cpu()
  659. }
  661. def explain_molecule(self, data):
  662. try:
  663. if self.explainer is not None:
  664. data = data.to(self.device)
  666. if not hasattr(data, 'batch') or data.batch is None:
  667. data.batch = torch.zeros(data.x.size(0), dtype=torch.long, device=self.device)
  669. explanation = self.explainer(
  670. x=data.x,
  671. edge_index=data.edge_index,
  672. batch=data.batch
  673. )
  675. return {
  676. 'node_mask': explanation.node_mask.cpu() if explanation.node_mask is not None else None,
  677. 'edge_mask': explanation.edge_mask.cpu() if explanation.edge_mask is not None else None,
  678. 'prediction': explanation.prediction.cpu() if hasattr(explanation, 'prediction') else None,
  679. 'data': data.cpu()
  680. }
  681. else:
  682. return self.simple_gradient_explanation(data)
  684. except Exception as e:
  685. print(f"Error in explanation: {e}")
  686. try:
  687. return self.simple_gradient_explanation(data)
  688. except Exception as e2:
  689. print(f"Error in gradient explanation: {e2}")
  690. return None
  692. def process_node_importance(self, node_mask, num_nodes):
  693. if node_mask is None:
  694. return np.full(num_nodes, 0.5)
  696. node_mask_np = node_mask.numpy()
  698. if len(node_mask_np.shape) > 1:
  699. if node_mask_np.shape[0] == num_nodes:
  700. node_colors = np.linalg.norm(node_mask_np, axis=1)
  701. elif node_mask_np.shape[1] == num_nodes:
  702. node_colors = np.linalg.norm(node_mask_np, axis=0)
  703. else:
  704. node_colors = np.full(num_nodes, np.mean(node_mask_np))
  705. else:
  706. node_colors = node_mask_np
  708. if len(node_colors) > num_nodes:
  709. node_colors = node_colors[:num_nodes]
  710. elif len(node_colors) < num_nodes:
  711. avg_importance = node_colors.mean() if len(node_colors) > 0 else 0.5
  712. padded_colors = np.full(num_nodes, avg_importance)
  713. padded_colors[:len(node_colors)] = node_colors
  714. node_colors = padded_colors
  716. if len(node_colors) > 0 and node_colors.max() > node_colors.min():
  717. node_colors = (node_colors - node_colors.min()) / (node_colors.max() - node_colors.min())
  718. else:
  719. node_colors = np.full(num_nodes, 0.5)
  721. return node_colors
  723. def visualize_molecule_explanation(self, explanation, smiles, target_value=None, title="Molecule Explanation",
  724. save_path=None):
  725. if explanation is None:
  726. print("No explanation to visualize")
  727. return
  729. data = explanation['data']
  730. node_mask = explanation['node_mask']
  731. edge_mask = explanation['edge_mask']
  732. prediction = explanation['prediction']
  734. fig, axes = plt.subplots(1, 2, figsize=(18, 8))
  736. try:
  737. mol = Chem.MolFromSmiles(smiles)
  738. if mol is not None:
  739. img = Draw.MolToImage(mol, size=(400, 400))
  740. axes[0].imshow(img)
  741. axes[0].set_title(f'Original Molecule\nSMILES: {smiles[:50]}...',
  742. fontsize=16, fontweight='bold', pad=15)
  743. axes[0].axis('off')
  744. except Exception as e:
  745. print(f"Warning: Could not draw molecule with RDKit: {e}")
  746. axes[0].text(0.5, 0.5, 'Could not draw molecule', ha='center', va='center',
  747. fontsize=16, fontweight='bold')
  748. axes[0].set_title('Original Molecule', fontsize=16, fontweight='bold')
  749. axes[0].axis('off')
  751. num_nodes = data.x.size(0)
  753. node_colors = self.process_node_importance(node_mask, num_nodes)
  755. G = nx.Graph()
  757. for i in range(num_nodes):
  758. atom_symbol = self.get_atom_symbol_from_features(data.x[i])
  759. G.add_node(i, symbol=atom_symbol)
  761. edge_index = data.edge_index.cpu().numpy()
  762. for i in range(edge_index.shape[1]):
  763. G.add_edge(edge_index[0, i], edge_index[1, i])
  765. try:
  766. pos = nx.spring_layout(G, seed=42)
  767. except:
  768. pos = nx.circular_layout(G)
  770. node_x = [pos[i][0] for i in range(num_nodes)]
  771. node_y = [pos[i][1] for i in range(num_nodes)]
  773. try:
  774. scatter = axes[1].scatter(
  775. node_x, node_y,
  776. c=node_colors,
  777. cmap='RdYlBu_r',
  778. s=600,
  779. alpha=0.8,
  780. edgecolors='black',
  781. linewidth=1.5
  782. )
  784. cbar = plt.colorbar(scatter, ax=axes[1], label='Node Importance', shrink=0.8)
  785. cbar.ax.tick_params(labelsize=14)
  786. cbar.set_label('Node Importance', fontsize=16, fontweight='bold')
  788. except Exception as e:
  789. print(f"Error in scatter plot: {e}")
  790. axes[1].scatter(node_x, node_y, c='lightblue', s=600, alpha=0.8, edgecolors='black')
  792. try:
  793. for edge in G.edges():
  794. x_coords = [pos[edge[0]][0], pos[edge[1]][0]]
  795. y_coords = [pos[edge[0]][1], pos[edge[1]][1]]
  796. axes[1].plot(x_coords, y_coords, 'k-', alpha=0.5, linewidth=2)
  797. except Exception as e:
  798. print(f"Warning: Could not draw edges: {e}")
  800. try:
  801. for i in range(num_nodes):
  802. atom_symbol = self.get_atom_symbol_from_features(data.x[i])
  803. axes[1].text(pos[i][0], pos[i][1], atom_symbol,
  804. ha='center', va='center', fontsize=14, fontweight='bold')
  805. except Exception as e:
  806. print(f"Warning: Could not add labels: {e}")
  808. pred_val = prediction.item() if prediction is not None else "N/A"
  809. if target_value is not None:
  810. title_text = f'Node Importance\nTrue y: {target_value:.4f}\nPrediction: {pred_val:.4f}'
  811. else:
  812. title_text = f'Node Importance\nPrediction: {pred_val:.4f}'
  814. axes[1].set_title(title_text, fontsize=16, fontweight='bold', pad=15)
  815. axes[1].axis('off')
  816. axes[1].set_aspect('equal')
  818. plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1, wspace=0.3)
  819. if save_path:
  820. try:
  821. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  822. print(f"Image saved: {save_path}")
  823. except Exception as e:
  824. print(f"Failed to save image: {e}")
  826. plt.show()
  828. def visualize_selected_molecule(self, explanation, smiles, target_value=None, title="Selected Molecule",
  829. save_path=None):
  830. if explanation is None:
  831. print("No explanation to visualize")
  832. return
  834. data = explanation['data']
  835. node_mask = explanation['node_mask']
  836. edge_mask = explanation['edge_mask']
  837. prediction = explanation['prediction']
  839. fig, axes = plt.subplots(1, 2, figsize=(18, 8))
  841. try:
  842. mol = Chem.MolFromSmiles(smiles)
  843. if mol is not None:
  844. img = Draw.MolToImage(mol, size=(400, 400))
  845. axes[0].imshow(img)
  846. axes[0].set_title(f'Original Molecule\nSMILES: {smiles[:50]}...',
  847. fontsize=16, fontweight='bold', pad=15)
  848. axes[0].axis('off')
  849. except Exception as e:
  850. print(f"Warning: Could not draw molecule with RDKit: {e}")
  851. axes[0].text(0.5, 0.5, 'Could not draw molecule', ha='center', va='center',
  852. fontsize=16, fontweight='bold')
  853. axes[0].set_title('Original Molecule', fontsize=16, fontweight='bold')
  854. axes[0].axis('off')
  856. num_nodes = data.x.size(0)
  858. node_colors = self.process_node_importance(node_mask, num_nodes)
  860. G = nx.Graph()
  862. for i in range(num_nodes):
  863. atom_symbol = self.get_atom_symbol_from_features(data.x[i])
  864. G.add_node(i, symbol=atom_symbol)
  866. edge_index = data.edge_index.cpu().numpy()
  867. for i in range(edge_index.shape[1]):
  868. G.add_edge(edge_index[0, i], edge_index[1, i])
  870. try:
  871. pos = nx.spring_layout(G, seed=42)
  872. except:
  873. pos = nx.circular_layout(G)
  875. node_x = [pos[i][0] for i in range(num_nodes)]
  876. node_y = [pos[i][1] for i in range(num_nodes)]
  878. try:
  879. scatter = axes[1].scatter(
  880. node_x, node_y,
  881. c=node_colors,
  882. cmap='RdYlBu_r',
  883. s=600,
  884. alpha=0.8,
  885. edgecolors='black',
  886. linewidth=1.5
  887. )
  889. cbar = plt.colorbar(scatter, ax=axes[1], label='Importance Score', shrink=0.8)
  890. cbar.ax.tick_params(labelsize=14)
  891. cbar.set_label('Importance Score', fontsize=16, fontweight='bold')
  893. except Exception as e:
  894. print(f"Error in scatter plot: {e}")
  895. axes[1].scatter(node_x, node_y, c='lightblue', s=600, alpha=0.8, edgecolors='black')
  897. try:
  898. for edge in G.edges():
  899. x_coords = [pos[edge[0]][0], pos[edge[1]][0]]
  900. y_coords = [pos[edge[0]][1], pos[edge[1]][1]]
  901. axes[1].plot(x_coords, y_coords, 'k-', alpha=0.5, linewidth=2)
  902. except Exception as e:
  903. print(f"Warning: Could not draw edges: {e}")
  905. try:
  906. for i in range(num_nodes):
  907. atom_symbol = self.get_atom_symbol_from_features(data.x[i])
  908. axes[1].text(pos[i][0], pos[i][1], atom_symbol,
  909. ha='center', va='center', fontsize=14, fontweight='bold')
  910. except Exception as e:
  911. print(f"Warning: Could not add labels: {e}")
  913. if target_value is not None:
  914. title_text = f'True y: {target_value:.4f}'
  915. else:
  916. title_text = 'True y: N/A'
  918. axes[1].set_title(title_text, fontsize=16, fontweight='bold', pad=15)
  919. axes[1].axis('off')
  920. axes[1].set_aspect('equal')
  922. plt.subplots_adjust(left=0.05, right=0.95, top=0.9, bottom=0.1, wspace=0.3)
  923. if save_path:
  924. try:
  925. plt.savefig(save_path, dpi=300, bbox_inches='tight')
  926. print(f"Image saved: {save_path}")
  927. except Exception as e:
  928. print(f"Failed to save image: {e}")
  930. plt.show()
  932. def analyze_feature_importance(self, explanations_data):
  933. all_node_importances = []
  934. all_atom_types = []
  935. all_predictions = []
  937. for exp_data in explanations_data:
  938. explanation = exp_data['explanation']
  939. if explanation is None or explanation['node_mask'] is None:
  940. continue
  942. node_mask = explanation['node_mask']
  943. data = explanation['data']
  944. prediction = explanation['prediction']
  945. num_nodes = data.x.size(0)
  947. node_importances = self.process_node_importance(node_mask, num_nodes)
  949. for i in range(num_nodes):
  950. atom_features = data.x[i]
  951. atom_symbol = self.get_atom_symbol_from_features(atom_features)
  953. all_node_importances.append(node_importances[i])
  954. all_atom_types.append(atom_symbol)
  955. all_predictions.append(prediction.item() if prediction is not None else 0)
  957. importance_df = pd.DataFrame({
  958. 'atom_type': all_atom_types,
  959. 'importance': all_node_importances,
  960. 'prediction': all_predictions
  961. })
  963. return importance_df
  965. def find_important_substructures(self, explanations_data, importance_threshold=0.5):
  966. important_substructures = []
  968. for exp_data in explanations_data:
  969. explanation = exp_data['explanation']
  970. smiles = exp_data['smiles']
  971. target_value = exp_data.get('original_target', 0)
  973. if explanation is None or explanation['node_mask'] is None:
  974. continue
  976. data = explanation['data']
  977. node_mask = explanation['node_mask']
  978. prediction = explanation['prediction']
  979. num_nodes = data.x.size(0)
  981. node_importances = self.process_node_importance(node_mask, num_nodes)
  983. important_atoms = []
  984. important_atom_indices = []
  986. for i in range(num_nodes):
  987. if node_importances[i] > importance_threshold:
  988. atom_symbol = self.get_atom_symbol_from_features(data.x[i])
  989. important_atoms.append((i, atom_symbol, node_importances[i]))
  990. important_atom_indices.append(i)
  992. if not important_atoms:
  993. continue
  995. try:
  996. mol = Chem.MolFromSmiles(smiles)
  997. if mol is None:
  998. continue
  1000. all_substructures = self.substructure_identifier.identify_substructures_in_molecule(mol)
  1002. functional_groups = self.substructure_identifier.get_functional_groups(mol)
  1004. local_substructures = self.substructure_identifier.extract_important_substructures(
  1005. mol, important_atom_indices, radius=2
  1006. )
  1008. relevant_substructures = {}
  1009. for sub_name, sub_info in all_substructures.items():
  1010. for match in sub_info['matches']:
  1011. if any(atom_idx in match for atom_idx in important_atom_indices):
  1012. if sub_name not in relevant_substructures:
  1013. relevant_substructures[sub_name] = {
  1014. 'pattern': sub_info['pattern'],
  1015. 'matches': [],
  1016. 'importance_scores': []
  1017. }
  1019. match_importance = []
  1020. for atom_idx in match:
  1021. if atom_idx in important_atom_indices:
  1022. idx_in_list = important_atom_indices.index(atom_idx)
  1023. match_importance.append(important_atoms[idx_in_list][2])
  1025. if match_importance:
  1026. avg_importance = np.mean(match_importance)
  1027. relevant_substructures[sub_name]['matches'].append(match)
  1028. relevant_substructures[sub_name]['importance_scores'].append(avg_importance)
  1030. edge_index = data.edge_index.cpu().numpy()
  1031. important_edges = []
  1033. for i in range(edge_index.shape[1]):
  1034. if (edge_index[0, i] in important_atom_indices and
  1035. edge_index[1, i] in important_atom_indices):
  1036. important_edges.append((edge_index[0, i], edge_index[1, i]))
  1038. important_substructures.append({
  1039. 'smiles': smiles,
  1040. 'important_atoms': important_atoms,
  1041. 'important_edges': important_edges,
  1042. 'prediction': prediction.item() if prediction is not None else 0,
  1043. 'target_value': target_value,
  1044. 'num_important_atoms': len(important_atoms),
  1046. 'known_substructures': relevant_substructures,
  1047. 'functional_groups': functional_groups,
  1048. 'local_substructures': local_substructures,
  1050. 'substructure_summary': {
  1051. 'num_known_substructures': len(relevant_substructures),
  1052. 'num_functional_groups': len(functional_groups),
  1053. 'num_local_fragments': len(local_substructures)
  1054. }
  1055. })
  1057. except Exception as e:
  1058. print(f"Error analyzing substructures for {smiles}: {e}")
  1059. important_substructures.append({
  1060. 'smiles': smiles,
  1061. 'important_atoms': important_atoms,
  1062. 'important_edges': important_edges,
  1063. 'prediction': prediction.item() if prediction is not None else 0,
  1064. 'target_value': target_value,
  1065. 'num_important_atoms': len(important_atoms),
  1066. 'known_substructures': {},
  1067. 'functional_groups': [],
  1068. 'local_substructures': [],
  1069. 'substructure_summary': {
  1070. 'num_known_substructures': 0,
  1071. 'num_functional_groups': 0,
  1072. 'num_local_fragments': 0
  1073. }
  1074. })
  1076. return important_substructures
  1078. def analyze_full_dataset_substructures(self, test_csv_file, importance_threshold=0.3):
  1079. print("\n" + "=" * 50)
  1080. print("Full Dataset Substructure Analysis")
  1081. print("=" * 50)
  1083. test_df = pd.read_csv(test_csv_file)
  1084. full_dataset_substructures = []
  1085. successful_count = 0
  1087. print(f"Analyzing substructures for all {len(test_df)} molecules in the dataset...")
  1089. for i, row in test_df.iterrows():
  1090. try:
  1091. smiles = str(row['Smiles'])
  1092. target_value = row.get('pchembl', 0)
  1094. atom_features, edge_index = smiles_to_graph(smiles)
  1095. data = Data(x=atom_features, edge_index=edge_index)
  1097. explanation = self.explain_molecule(data)
  1099. if explanation is not None:
  1100. node_mask = explanation['node_mask']
  1101. num_nodes = data.x.size(0)
  1102. node_importances = self.process_node_importance(node_mask, num_nodes)
  1104. important_atoms = []
  1105. important_atom_indices = []
  1107. for j in range(num_nodes):
  1108. if node_importances[j] > importance_threshold:
  1109. atom_symbol = self.get_atom_symbol_from_features(data.x[j])
  1110. important_atoms.append((j, atom_symbol, node_importances[j]))
  1111. important_atom_indices.append(j)
  1113. if important_atoms:
  1114. mol = Chem.MolFromSmiles(smiles)
  1115. if mol is not None:
  1116. all_substructures = self.substructure_identifier.identify_substructures_in_molecule(mol)
  1117. functional_groups = self.substructure_identifier.get_functional_groups(mol)
  1118. local_substructures = self.substructure_identifier.extract_important_substructures(
  1119. mol, important_atom_indices, radius=2
  1120. )
  1122. relevant_substructures = {}
  1123. for sub_name, sub_info in all_substructures.items():
  1124. for match in sub_info['matches']:
  1125. if any(atom_idx in match for atom_idx in important_atom_indices):
  1126. if sub_name not in relevant_substructures:
  1127. relevant_substructures[sub_name] = {
  1128. 'pattern': sub_info['pattern'],
  1129. 'matches': [],
  1130. 'importance_scores': []
  1131. }
  1133. match_importance = []
  1134. for atom_idx in match:
  1135. if atom_idx in important_atom_indices:
  1136. idx_in_list = important_atom_indices.index(atom_idx)
  1137. match_importance.append(important_atoms[idx_in_list][2])
  1139. if match_importance:
  1140. avg_importance = np.mean(match_importance)
  1141. relevant_substructures[sub_name]['matches'].append(match)
  1142. relevant_substructures[sub_name]['importance_scores'].append(avg_importance)
  1144. edge_index = data.edge_index.cpu().numpy()
  1145. important_edges = []
  1147. for j in range(edge_index.shape[1]):
  1148. if (edge_index[0, j] in important_atom_indices and
  1149. edge_index[1, j] in important_atom_indices):
  1150. important_edges.append((edge_index[0, j], edge_index[1, j]))
  1152. full_dataset_substructures.append({
  1153. 'smiles': smiles,
  1154. 'important_atoms': important_atoms,
  1155. 'important_edges': important_edges,
  1156. 'prediction': explanation['prediction'].item() if explanation[
  1157. 'prediction'] is not None else 0,
  1158. 'target_value': target_value,
  1159. 'num_important_atoms': len(important_atoms),
  1160. 'known_substructures': relevant_substructures,
  1161. 'functional_groups': functional_groups,
  1162. 'local_substructures': local_substructures,
  1163. 'substructure_summary': {
  1164. 'num_known_substructures': len(relevant_substructures),
  1165. 'num_functional_groups': len(functional_groups),
  1166. 'num_local_fragments': len(local_substructures)
  1167. }
  1168. })
  1169. successful_count += 1
  1171. except Exception as e:
  1172. continue
  1174. if (i + 1) % 100 == 0:
  1175. print(f" Progress: {i + 1}/{len(test_df)}, successful: {successful_count}")
  1177. print(f"Full dataset analysis completed: {successful_count} molecules with important substructures")
  1178. return full_dataset_substructures
  1180. def plot_feature_importance_summary(self, importance_df, save_path=None):
  1181. if importance_df.empty:
  1182. print("No importance data to plot")
  1183. return
  1185. custom_colors = ['#98CFE6', '#ADE7A8', '#F39F4E', '#EEB7D3', '#DBDAD3', '#FFDF97']
  1187. try:
  1188. atom_importance = importance_df.groupby('atom_type')['importance'].agg(['mean', 'std', 'count'])
  1189. atom_importance = atom_importance.sort_values('mean', ascending=False)
  1190. atom_counts = importance_df['atom_type'].value_counts()
  1192. plt.figure(figsize=(14, 10))
  1193. for i, (atom_type, stats) in enumerate(atom_importance.iterrows()):
  1194. color_idx = i % len(custom_colors)
  1195. plt.bar(atom_type, stats['mean'], yerr=stats['std'], capsize=5,
  1196. alpha=0.8, color=custom_colors[color_idx], edgecolor='white', linewidth=1)
  1198. plt.title('Average Atom Importance', fontsize=20, fontweight='bold', pad=25)
  1199. plt.xlabel('Atom Type', fontsize=18, fontweight='bold')
  1200. plt.ylabel('Average Importance', fontsize=18, fontweight='bold')
  1201. plt.xticks(rotation=45, fontsize=16)
  1202. plt.yticks(fontsize=16)
  1203. plt.grid(axis='y', alpha=0.3)
  1205. for i, (atom_type, stats) in enumerate(atom_importance.iterrows()):
  1206. plt.text(i, stats['mean'] + stats['std'] + 0.01, f'{stats["mean"]:.3f}',
  1207. ha='center', va='bottom', fontsize=12, fontweight='bold')
  1209. plt.subplots_adjust(left=0.15, right=0.95, top=0.9, bottom=0.2)
  1210. if save_path:
  1211. plt.savefig(f"{save_path}_subplot1_atom_importance.png", dpi=300, bbox_inches='tight')
  1212. print(f"Subplot 1 saved: {save_path}_subplot1_atom_importance.png")
  1213. plt.show()
  1215. plt.figure(figsize=(14, 10))
  1217. total_contribution = importance_df.groupby('atom_type')['importance'].sum().sort_values(ascending=False)
  1218. cumulative_contribution = total_contribution.cumsum() / total_contribution.sum() * 100
  1220. plt.plot(range(1, len(cumulative_contribution) + 1), cumulative_contribution.values,
  1221. 'o-', linewidth=6, markersize=12, color=custom_colors[0], markerfacecolor='white',
  1222. markeredgewidth=3, markeredgecolor=custom_colors[0])
  1224. plt.axhline(y=80, color='red', linestyle='--', alpha=0.7, label='80% Contribution', linewidth=4)
  1226. plt.xlabel('Atom Type Rank', fontsize=18, fontweight='bold')
  1227. plt.ylabel('Cumulative Contribution (%)', fontsize=18, fontweight='bold')
  1228. plt.title('Cumulative Importance Contribution', fontsize=20, fontweight='bold', pad=25)
  1229. plt.xticks(fontsize=16)
  1230. plt.yticks(fontsize=16)
  1231. plt.grid(True, alpha=0.3)
  1232. plt.legend(fontsize=16, loc='lower right')
  1234. for i, (atom_type, contrib) in enumerate(cumulative_contribution.items()):
  1235. if i < 8:
  1236. plt.annotate(atom_type, (i + 1, contrib), xytext=(0, 15),
  1237. textcoords='offset points', fontsize=14, fontweight='bold',
  1238. ha='center', va='bottom')
  1240. plt.subplots_adjust(left=0.15, right=0.95, top=0.9, bottom=0.15)
  1241. if save_path:
  1242. plt.savefig(f"{save_path}_subplot2_cumulative_contribution.png", dpi=300, bbox_inches='tight')
  1243. print(f"Subplot 2 saved: {save_path}_subplot2_cumulative_contribution.png")
  1244. plt.show()
  1246. plt.figure(figsize=(14, 12))
  1247. pie_colors = [custom_colors[i % len(custom_colors)] for i in range(len(atom_counts))]
  1249. def autopct_func_atoms(pct, idx=0):
  1250. if idx < len(atom_counts) - 5:
  1251. return f'{pct:.1f}%'
  1252. else:
  1253. return ''
  1255. counter = [0]
  1257. def autopct_with_counter_atoms(pct):
  1258. result = autopct_func_atoms(pct, counter[0])
  1259. counter[0] += 1
  1260. return result
  1262. wedges, texts, autotexts = plt.pie(atom_counts.values, autopct=autopct_with_counter_atoms,
  1263. colors=pie_colors, startangle=90,
  1264. textprops={'fontsize': 16, 'fontweight': 'bold'})
  1265. plt.title('Atom Type Distribution', fontsize=20, fontweight='bold', pad=25)
  1267. for autotext in autotexts:
  1268. autotext.set_fontsize(16)
  1269. autotext.set_color('black')
  1270. autotext.set_fontweight('bold')
  1272. plt.legend(wedges, atom_counts.index, title="Atom Types",
  1273. title_fontsize=16, fontsize=14,
  1274. loc="center left", bbox_to_anchor=(1.1, 0.5))
  1276. plt.subplots_adjust(left=0.1, right=0.75, top=0.9, bottom=0.1)
  1277. if save_path:
  1278. plt.savefig(f"{save_path}_subplot3_atom_distribution.png", dpi=300, bbox_inches='tight')
  1279. print(f"Subplot 3 saved: {save_path}_subplot3_atom_distribution.png")
  1280. plt.show()
  1282. plt.figure(figsize=(14, 10))
  1284. atom_types_for_box = []
  1285. importance_values_for_box = []
  1287. for atom_type in atom_importance.index[:10]:
  1288. atom_data = importance_df[importance_df['atom_type'] == atom_type]['importance']
  1289. atom_types_for_box.extend([atom_type] * len(atom_data))
  1290. importance_values_for_box.extend(atom_data.values)
  1292. box_df = pd.DataFrame({
  1293. 'atom_type': atom_types_for_box,
  1294. 'importance': importance_values_for_box
  1295. })
  1297. import seaborn as sns
  1299. sns.set_context("paper", font_scale=1.5)
  1301. sns.boxplot(data=box_df, x='atom_type', y='importance',
  1302. palette=custom_colors[:len(atom_importance.index[:10])])
  1304. plt.title('Importance Distribution by Atom Type', fontsize=20, fontweight='bold', pad=25)
  1305. plt.xlabel('Atom Type', fontsize=18, fontweight='bold')
  1306. plt.ylabel('Importance Score', fontsize=18, fontweight='bold')
  1307. plt.xticks(rotation=45, fontsize=16)
  1308. plt.yticks(fontsize=16)
  1309. plt.grid(axis='y', alpha=0.3)
  1311. plt.subplots_adjust(left=0.15, right=0.95, top=0.9, bottom=0.2)
  1312. if save_path:
  1313. plt.savefig(f"{save_path}_subplot4_importance_distribution.png", dpi=300, bbox_inches='tight')
  1314. print(f"Subplot 4 saved: {save_path}_subplot4_importance_distribution.png")
  1315. plt.show()
  1317. sns.reset_defaults()
  1319. print(f"\nAll four subplots have been generated and saved separately!")
  1321. except Exception as e:
  1322. print(f"Error plotting feature importance: {e}")
  1324. def visualize_all_substructures(self, important_substructures, full_dataset_substructures,
  1325. base_save_path="substructure_analysis"):
  1326. print("\nGenerating substructure visualizations...")
  1328. print(" - Generating substructure importance summary (4 separate plots)...")
  1329. self.substructure_visualizer.visualize_substructure_importance_summary(
  1330. important_substructures,
  1331. save_path=f"{base_save_path}_importance_summary"
  1332. )
  1334. print(" - Generating highlighted molecules from full dataset...")
  1335. selected_molecules = self.substructure_visualizer.visualize_substructure_in_molecules(
  1336. full_dataset_substructures,
  1337. num_examples=6,
  1338. save_path=f"{base_save_path}_highlighted_molecules.png"
  1339. )
  1341. print(" - Generating substructure heatmap...")
  1342. self.substructure_visualizer.create_substructure_heatmap(
  1343. important_substructures,
  1344. save_path=f"{base_save_path}_heatmap.png"
  1345. )
  1347. print("Substructure visualization completed!")
  1349. return selected_molecules
  1352. def load_best_model(model_path='best_model.pth'):
  1353. try:
  1354. checkpoint = torch.load(model_path, map_location='cpu')
  1356. gat_graphsage_model = GAT_GraphSAGE(n_output=1, num_features_xd=35)
  1358. gat_graphsage_model.load_state_dict(checkpoint['gat_graphsage_model_state_dict'])
  1360. scaler = checkpoint.get('scaler', None)
  1362. print("Model loaded successfully!")
  1363. return gat_graphsage_model, scaler
  1364. except Exception as e:
  1365. print(f"Error loading model: {e}")
  1366. return None, None
  1369. def smiles_to_graph(smiles):
  1370. mol = Chem.MolFromSmiles(smiles)
  1371. if mol is None:
  1372. raise ValueError(f"Invalid SMILES string: {smiles}")
  1374. num_atoms = mol.GetNumAtoms()
  1375. atom_features = []
  1377. for atom in mol.GetAtoms():
  1378. results = one_of_k_encoding_unk(atom.GetSymbol(), ['C', 'N', 'O', 'S', 'F', 'P', 'Cl', 'Br', 'I', 'Unknown']) + \
  1379. one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4, 5, 6]) + \
  1380. one_of_k_encoding_unk(atom.GetImplicitValence(), [0, 1, 2, 3, 4, 5, 6]) + \
  1381. one_of_k_encoding_unk(atom.GetHybridization(), [
  1382. Chem.rdchem.HybridizationType.SP, Chem.rdchem.HybridizationType.SP2,
  1383. Chem.rdchem.HybridizationType.SP3, Chem.rdchem.HybridizationType.SP3D,
  1384. Chem.rdchem.HybridizationType.SP3D2
  1385. ]) + [atom.GetIsAromatic()] + \
  1386. one_of_k_encoding_unk(atom.GetTotalNumHs(), [0, 1, 2, 3, 4])
  1387. atom_feats = np.array(results).astype(np.float32)
  1388. atom_features.append(atom_feats)
  1390. atom_features = torch.tensor(atom_features, dtype=torch.float)
  1391. adj_matrix = torch.zeros((num_atoms, num_atoms), dtype=torch.float)
  1393. for bond in mol.GetBonds():
  1394. start, end = bond.GetBeginAtomIdx(), bond.GetEndAtomIdx()
  1395. adj_matrix[start, end] = 1.0
  1396. adj_matrix[end, start] = 1.0
  1398. edge_index = adj_matrix.nonzero(as_tuple=False).t().long()
  1399. return atom_features, edge_index
  1402. def quick_importance_analysis_all(test_csv_file, model_path):
  1403. print("Stage 1: Quick analysis of all molecules...")
  1404. print("-" * 50)
  1406. test_df = pd.read_csv(test_csv_file)
  1407. gat_model, scaler = load_best_model(model_path)
  1408. explainable_model = ExplainableGATGraphSAGE(gat_model)
  1409. explainer = MolecularExplainer(explainable_model)
  1411. molecule_info = []
  1412. successful_count = 0
  1414. for i, row in test_df.iterrows():
  1415. try:
  1416. smiles = str(row['Smiles'])
  1417. atom_features, edge_index = smiles_to_graph(smiles)
  1418. data = Data(x=atom_features, edge_index=edge_index)
  1420. explanation = explainer.simple_gradient_explanation(data)
  1422. if explanation:
  1423. molecule_info.append({
  1424. 'index': i,
  1425. 'smiles': smiles,
  1426. 'prediction': explanation['prediction'].item(),
  1427. 'avg_importance': explanation['node_mask'].mean().item(),
  1428. 'max_importance': explanation['node_mask'].max().item(),
  1429. 'std_importance': explanation['node_mask'].std().item(),
  1430. 'num_atoms': data.x.size(0),
  1431. 'target': row.get('pchembl', 0)
  1432. })
  1433. successful_count += 1
  1435. except Exception as e:
  1436. continue
  1438. if (i + 1) % 100 == 0:
  1439. print(f" Quick analysis progress: {i + 1}/961, successful: {successful_count}")
  1441. print(f"Stage 1 completed: Successfully analyzed {successful_count} molecules")
  1442. return molecule_info
  1445. def stratified_sample_by_column(df, column, target_count):
  1446. try:
  1447. df_copy = df.copy()
  1448. df_copy['quartile'] = pd.qcut(df_copy[column], q=5, labels=False, duplicates='drop')
  1450. sampled_indices = []
  1451. samples_per_quartile = target_count // 5
  1453. for q in df_copy['quartile'].unique():
  1454. if pd.isna(q):
  1455. continue
  1456. quartile_data = df_copy[df_copy['quartile'] == q]
  1457. if len(quartile_data) > 0:
  1458. sample_count = min(samples_per_quartile, len(quartile_data))
  1459. sampled = quartile_data.sample(n=sample_count, random_state=42)
  1460. sampled_indices.extend(sampled['index'].tolist())
  1462. return sampled_indices
  1463. except Exception as e:
  1464. print(f"Stratified sampling failed, using random sampling: {e}")
  1465. return df.sample(n=min(target_count, len(df)), random_state=42)['index'].tolist()
  1468. def select_representative_molecules(molecule_info, target_count=200):
  1469. print("\nStage 2: Selecting representative molecules for detailed analysis...")
  1470. print("-" * 50)
  1472. df = pd.DataFrame(molecule_info)
  1474. if len(df) < target_count:
  1475. print(f"Available molecules ({len(df)}) < target count ({target_count}), will analyze all available")
  1476. return df['index'].tolist()
  1478. selected_indices = []
  1480. print(" - Stratified sampling by prediction values...")
  1481. pred_samples = stratified_sample_by_column(df, 'prediction', int(target_count * 0.4))
  1482. selected_indices.extend(pred_samples)
  1484. print(" - Stratified sampling by average importance...")
  1485. remaining_df = df[~df['index'].isin(selected_indices)]
  1486. if len(remaining_df) > 0:
  1487. imp_samples = stratified_sample_by_column(remaining_df, 'avg_importance', int(target_count * 0.3))
  1488. selected_indices.extend(imp_samples)
  1490. print(" - Stratified sampling by molecule size...")
  1491. remaining_df = df[~df['index'].isin(selected_indices)]
  1492. if len(remaining_df) > 0:
  1493. size_samples = stratified_sample_by_column(remaining_df, 'num_atoms', int(target_count * 0.2))
  1494. selected_indices.extend(size_samples)
  1496. print(" - Random sampling for remaining molecules...")
  1497. remaining_df = df[~df['index'].isin(selected_indices)]
  1498. remaining_count = target_count - len(selected_indices)
  1499. if remaining_count > 0 and len(remaining_df) > 0:
  1500. random_samples = remaining_df.sample(n=min(remaining_count, len(remaining_df)),
  1501. random_state=42)['index'].tolist()
  1502. selected_indices.extend(random_samples)
  1504. print(f"Selected {len(selected_indices)} representative molecules for detailed analysis")
  1506. selected_df = df[df['index'].isin(selected_indices)]
  1507. print(f"\nSampling statistics:")
  1508. print(f" Prediction range: {selected_df['prediction'].min():.3f} - {selected_df['prediction'].max():.3f}")
  1509. print(f" Importance range: {selected_df['avg_importance'].min():.3f} - {selected_df['avg_importance'].max():.3f}")
  1510. print(f" Molecule size range: {selected_df['num_atoms'].min()} - {selected_df['num_atoms'].max()} atoms")
  1512. return selected_indices
  1515. def perform_detailed_analysis(test_csv_file, model_path, selected_indices, importance_threshold=0.3):
  1516. print("\nStage 3: Detailed substructure analysis of selected molecules...")
  1517. print("-" * 50)
  1519. test_df = pd.read_csv(test_csv_file)
  1520. gat_model, scaler = load_best_model(model_path)
  1521. explainable_model = ExplainableGATGraphSAGE(gat_model)
  1522. explainer = MolecularExplainer(explainable_model)
  1524. explanations_data = []
  1525. successful_count = 0
  1527. for i, idx in enumerate(selected_indices):
  1528. row = test_df.iloc[idx]
  1529. smiles = str(row['Smiles'])
  1531. try:
  1532. atom_features, edge_index = smiles_to_graph(smiles)
  1533. data = Data(x=atom_features, edge_index=edge_index)
  1534. explanation = explainer.explain_molecule(data)
  1536. if explanation is not None:
  1537. explanations_data.append({
  1538. 'smiles': smiles,
  1539. 'explanation': explanation,
  1540. 'original_target': row.get('pchembl', 0),
  1541. 'original_index': idx
  1542. })
  1543. successful_count += 1
  1545. if (i + 1) % 20 == 0:
  1546. print(f" Detailed analysis progress: {i + 1}/{len(selected_indices)}, successful: {successful_count}")
  1548. except Exception as e:
  1549. print(f" Molecule {idx} processing failed: {e}")
  1550. continue
  1552. print(f"Stage 3 completed: Successfully explained {successful_count} molecules")
  1554. print("\nPerforming feature importance analysis...")
  1555. importance_df = explainer.analyze_feature_importance(explanations_data)
  1557. print("Identifying important substructures from sampled data...")
  1558. important_substructures = explainer.find_important_substructures(
  1559. explanations_data, importance_threshold
  1560. )
  1562. print("Analyzing full dataset for qualified molecules...")
  1563. full_dataset_substructures = explainer.analyze_full_dataset_substructures(
  1564. test_csv_file, importance_threshold
  1565. )
  1567. return {
  1568. 'explanations_data': explanations_data,
  1569. 'importance_df': importance_df,
  1570. 'important_substructures': important_substructures,
  1571. 'full_dataset_substructures': full_dataset_substructures,
  1572. 'explainer': explainer
  1573. }
  1576. def combine_quick_and_detailed_results(quick_results, detailed_results):
  1577. print("\nStage 4: Combining and summarizing analysis results...")
  1578. print("-" * 50)
  1580. quick_df = pd.DataFrame(quick_results)
  1582. global_stats = {
  1583. 'total_molecules_analyzed': len(quick_df),
  1584. 'prediction_range': (quick_df['prediction'].min(), quick_df['prediction'].max()),
  1585. 'prediction_mean': quick_df['prediction'].mean(),
  1586. 'prediction_std': quick_df['prediction'].std(),
  1587. 'avg_importance_range': (quick_df['avg_importance'].min(), quick_df['avg_importance'].max()),
  1588. 'avg_importance_mean': quick_df['avg_importance'].mean(),
  1589. 'molecule_size_range': (quick_df['num_atoms'].min(), quick_df['num_atoms'].max()),
  1590. 'molecule_size_mean': quick_df['num_atoms'].mean()
  1591. }
  1593. combined_results = {
  1594. 'global_statistics': global_stats,
  1595. 'quick_analysis_results': quick_results,
  1596. 'detailed_analysis_results': detailed_results,
  1597. 'summary': {
  1598. 'total_molecules': len(quick_df),
  1599. 'detailed_molecules': len(detailed_results['explanations_data']),
  1600. 'identified_substructures': len(detailed_results['important_substructures']),
  1601. 'full_dataset_substructures': len(detailed_results['full_dataset_substructures']),
  1602. 'analysis_completeness': len(detailed_results['explanations_data']) / len(quick_df) * 100
  1603. }
  1604. }
  1606. return combined_results
  1609. def hybrid_analysis_strategy(test_csv_file, model_path, target_detailed_count=200, importance_threshold=0.3):
  1610. print("=" * 70)
  1611. print("Hybrid Analysis Strategy: Molecular Model Explainability Analysis")
  1612. print("=" * 70)
  1613. print(f"Data file: {test_csv_file}")
  1614. print(f"Model file: {model_path}")
  1615. print(f"Target detailed analysis count: {target_detailed_count}")
  1616. print(f"Importance threshold: {importance_threshold}")
  1617. print("=" * 70)
  1619. quick_results = quick_importance_analysis_all(test_csv_file, model_path)
  1621. if not quick_results:
  1622. print("Quick analysis failed, cannot continue")
  1623. return None
  1625. representative_indices = select_representative_molecules(
  1626. quick_results,
  1627. target_count=target_detailed_count
  1628. )
  1630. detailed_results = perform_detailed_analysis(
  1631. test_csv_file,
  1632. model_path,
  1633. representative_indices,
  1634. importance_threshold
  1635. )
  1637. final_results = combine_quick_and_detailed_results(quick_results, detailed_results)
  1639. generate_comprehensive_report(final_results)
  1641. return final_results
  1644. def generate_comprehensive_report(results):
  1645. print("\n" + "=" * 70)
  1646. print("Comprehensive Analysis Report")
  1647. print("=" * 70)
  1649. global_stats = results['global_statistics']
  1650. summary = results['summary']
  1652. print(f"\n[Global Statistics - Based on {global_stats['total_molecules_analyzed']} molecules]")
  1653. print("-" * 50)
  1654. print(f"Prediction distribution:")
  1655. print(f" Range: {global_stats['prediction_range'][0]:.3f} - {global_stats['prediction_range'][1]:.3f}")
  1656. print(f" Mean: {global_stats['prediction_mean']:.3f} ± {global_stats['prediction_std']:.3f}")
  1658. print(f"\nImportance distribution:")
  1659. print(f" Range: {global_stats['avg_importance_range'][0]:.3f} - {global_stats['avg_importance_range'][1]:.3f}")
  1660. print(f" Mean: {global_stats['avg_importance_mean']:.3f}")
  1662. print(f"\nMolecule size distribution:")
  1663. print(f" Range: {global_stats['molecule_size_range'][0]} - {global_stats['molecule_size_range'][1]} atoms")
  1664. print(f" Mean: {global_stats['molecule_size_mean']:.1f} atoms")
  1666. detailed = results['detailed_analysis_results']
  1668. print(f"\n[Detailed Analysis Results - Based on {summary['detailed_molecules']} representative molecules]")
  1669. print("-" * 50)
  1671. if not detailed['importance_df'].empty:
  1672. print("Atom importance statistics:")
  1673. atom_stats = detailed['importance_df'].groupby('atom_type')['importance'].agg(['mean', 'std', 'count'])
  1674. atom_stats = atom_stats.sort_values('mean', ascending=False)
  1676. for atom_type, stats in atom_stats.head(8).iterrows():
  1677. print(
  1678. f" {atom_type}: Average importance {stats['mean']:.3f} (±{stats['std']:.3f}), appeared {stats['count']} times")
  1680. important_substructures = detailed['important_substructures']
  1681. full_dataset_substructures = detailed['full_dataset_substructures']
  1683. if important_substructures:
  1684. print(f"\nFound {len(important_substructures)} molecules containing important substructures (from sample)")
  1685. print(
  1686. f"Found {len(full_dataset_substructures)} molecules containing important substructures (from full dataset)")
  1688. all_substructure_types = {}
  1689. all_functional_groups = Counter()
  1691. for struct in important_substructures:
  1692. for sub_name, sub_info in struct['known_substructures'].items():
  1693. if sub_name not in all_substructure_types:
  1694. all_substructure_types[sub_name] = {
  1695. 'count': 0,
  1696. 'total_importance': 0,
  1697. 'molecules': []
  1698. }
  1699. all_substructure_types[sub_name]['count'] += len(sub_info['matches'])
  1700. all_substructure_types[sub_name]['total_importance'] += sum(sub_info['importance_scores'])
  1701. all_substructure_types[sub_name]['molecules'].append(struct['smiles'][:30])
  1703. for fg in struct['functional_groups']:
  1704. all_functional_groups[fg['name']] += fg['count']
  1706. print(f"\nMost common important substructures (Top 10):")
  1707. sorted_substructures = sorted(all_substructure_types.items(),
  1708. key=lambda x: x[1]['count'], reverse=True)
  1710. for i, (sub_name, info) in enumerate(sorted_substructures[:10], 1):
  1711. avg_importance = info['total_importance'] / info['count'] if info['count'] > 0 else 0
  1712. print(f" {i:2d}. {sub_name:15s}: appeared {info['count']:3d} times, avg importance {avg_importance:.3f}")
  1714. print(f"\nMost common functional groups (Top 10):")
  1715. for i, (fg_name, count) in enumerate(all_functional_groups.most_common(10), 1):
  1716. print(f" {i:2d}. {fg_name:20s}: {count:3d} times")
  1718. print(f"\n[Analysis Completeness]")
  1719. print("-" * 30)
  1720. print(f"Total molecules: {summary['total_molecules']}")
  1721. print(f"Detailed analysis molecules: {summary['detailed_molecules']}")
  1722. print(f"Full dataset substructure analysis: {summary['full_dataset_substructures']}")
  1723. print(f"Analysis coverage: {summary['analysis_completeness']:.1f}%")
  1724. print(f"Identified important substructures: {summary['identified_substructures']}")
  1726. if not detailed['importance_df'].empty:
  1727. print(f"\nGenerating visualizations...")
  1728. try:
  1729. detailed['explainer'].plot_feature_importance_summary(
  1730. detailed['importance_df'],
  1731. 'enhanced_feature_importance_analysis'
  1732. )
  1734. selected_highlighted_molecules = detailed['explainer'].visualize_all_substructures(
  1735. detailed['important_substructures'],
  1736. detailed['full_dataset_substructures'],
  1737. base_save_path="substructure_analysis"
  1738. )
  1740. print(" - Generating selected molecule visualizations from full dataset (True y only)...")
  1742. qualified_molecules = []
  1743. for struct in full_dataset_substructures:
  1744. if (struct.get('target_value', 0) > 6 and
  1745. any(max(sub_info['importance_scores']) > 0.5
  1746. for sub_info in struct['known_substructures'].values()
  1747. if sub_info['importance_scores'])):
  1748. qualified_molecules.append(struct)
  1749. if len(qualified_molecules) >= 6:
  1750. break
  1752. print(f" Found {len(qualified_molecules)} qualified molecules from full dataset")
  1754. for i, struct_data in enumerate(qualified_molecules):
  1755. try:
  1756. smiles = struct_data['smiles']
  1757. target_value = struct_data.get('target_value')
  1759. atom_features, edge_index = smiles_to_graph(smiles)
  1760. data = Data(x=atom_features, edge_index=edge_index)
  1761. explanation = detailed['explainer'].explain_molecule(data)
  1763. if explanation is not None:
  1764. print(f" Generating visualization for selected molecule {i + 1}...")
  1765. detailed['explainer'].visualize_selected_molecule(
  1766. explanation, smiles,
  1767. target_value=target_value,
  1768. title=f'Selected Molecule {i + 1}',
  1769. save_path=f'selected_molecule_{i + 1}.png'
  1770. )
  1772. except Exception as e:
  1773. print(f" Error generating visualization for molecule {i + 1}: {e}")
  1775. except Exception as e:
  1776. print(f"Error generating visualizations: {e}")
  1778. print(f"\n" + "=" * 70)
  1779. print("Analysis completed!")
  1780. print("=" * 70)
  1781. print(f"Generated files:")
  1782. print(f"- enhanced_feature_importance_analysis_subplot1_atom_importance.png: Average Atom Importance")
  1783. print(f"- enhanced_feature_importance_analysis_subplot2_cumulative_contribution.png: Cumulative contribution")
  1784. print(f"- enhanced_feature_importance_analysis_subplot3_atom_distribution.png: Atom type distribution")
  1785. print(f"- enhanced_feature_importance_analysis_subplot4_importance_distribution.png: Importance distribution")
  1786. print(f"- substructure_analysis_importance_summary_subplot1_frequency.png: Substructure frequency")
  1787. print(f"- substructure_analysis_importance_summary_subplot2_importance.png: Substructure importance")
  1788. print(f"- substructure_analysis_importance_summary_subplot3_functional_groups.png: Functional group distribution")
  1789. print(f"- substructure_analysis_importance_summary_subplot4_scatter.png: Importance vs frequency scatter")
  1790. print(
  1791. f"- substructure_analysis_highlighted_molecules.png: Highlighted substructure molecules (Full Dataset: y > 6, Importance > 0.5)")
  1792. print(f"- substructure_analysis_heatmap.png: Substructure-molecule heatmap (Top 40 molecules)")
  1793. print(
  1794. f"- selected_molecule_1.png to selected_molecule_6.png: Selected molecules (Full Dataset: y > 6, Importance > 0.5, True y only, No predict values)")
  1797. def display_analysis_results(results):
  1798. if results is None:
  1799. print("No analysis results to display")
  1800. return
  1802. print("\n" + "=" * 60)
  1803. print("Detailed Analysis Results View")
  1804. print("=" * 60)
  1806. quick_df = pd.DataFrame(results['quick_analysis_results'])
  1807. print(f"\n=== Global prediction distribution (based on {len(quick_df)} molecules) ===")
  1808. print(quick_df['prediction'].describe())
  1810. print(f"\n=== Importance distribution ===")
  1811. print(quick_df['avg_importance'].describe())
  1813. substructures = results['detailed_analysis_results']['important_substructures']
  1814. full_dataset_substructures = results['detailed_analysis_results']['full_dataset_substructures']
  1816. print(f'\n=== Important substructure analysis ===')
  1817. print(f'Found {len(substructures)} molecules containing important substructures (from sample)')
  1818. print(f'Found {len(full_dataset_substructures)} molecules containing important substructures (from full dataset)')
  1820. all_substructure_types = {}
  1821. for struct in substructures:
  1822. for sub_name, sub_info in struct['known_substructures'].items():
  1823. if sub_name not in all_substructure_types:
  1824. all_substructure_types[sub_name] = 0
  1825. all_substructure_types[sub_name] += len(sub_info['matches'])
  1827. if all_substructure_types:
  1828. print(f"\nMost common chemical substructures:")
  1829. sorted_subs = sorted(all_substructure_types.items(), key=lambda x: x[1], reverse=True)
  1830. for sub_name, count in sorted_subs[:10]:
  1831. print(f" {sub_name}: {count} times")
  1833. importance_df = results['detailed_analysis_results']['importance_df']
  1834. if not importance_df.empty:
  1835. print(f"\n=== Atom importance analysis ===")
  1836. atom_stats = importance_df.groupby('atom_type')['importance'].agg(['mean', 'std', 'count'])
  1837. atom_stats = atom_stats.sort_values('mean', ascending=False)
  1838. print("Atom type importance ranking:")
  1839. for atom_type, stats in atom_stats.head(8).iterrows():
  1840. print(
  1841. f" {atom_type}: Average importance {stats['mean']:.3f} (±{stats['std']:.3f}), appeared {stats['count']} times")
  1843. print(f"\n=== Specific molecule examples (from full dataset) ===")
  1844. qualified_count = 0
  1845. for struct in full_dataset_substructures[:10]:
  1846. if (struct.get('target_value', 0) > 6 and
  1847. any(max(sub_info['importance_scores']) > 0.5
  1848. for sub_info in struct['known_substructures'].values()
  1849. if sub_info['importance_scores'])):
  1850. qualified_count += 1
  1851. print(f"\nQualified Molecule {qualified_count}:")
  1852. print(f" SMILES: {struct['smiles'][:60]}...")
  1853. print(f" Prediction: {struct['prediction']:.4f}")
  1854. print(f" True y value: {struct.get('target_value', 'N/A')}")
  1855. print(f" Important atoms: {struct['num_important_atoms']}")
  1857. if struct['known_substructures']:
  1858. print(f" Found substructures:")
  1859. for sub_name, sub_info in list(struct['known_substructures'].items())[:3]:
  1860. avg_imp = np.mean(sub_info['importance_scores']) if sub_info['importance_scores'] else 0
  1861. print(f" • {sub_name}: average importance {avg_imp:.3f}")
  1863. if qualified_count >= 3:
  1864. break
  1866. print(
  1867. f"\nTotal qualified molecules found in full dataset: {sum(1 for struct in full_dataset_substructures if struct.get('target_value', 0) > 6 and any(max(sub_info['importance_scores']) > 0.5 for sub_info in struct['known_substructures'].values() if sub_info['importance_scores']))}")
  1870. if __name__ == "__main__":
  1871. results = hybrid_analysis_strategy(
  1872. test_csv_file='D:\\pycharm\\gutingle\\pythonProject2\\突出核蛋白\\test_data.csv',
  1873. model_path='D:\\pycharm\\gutingle\\pythonProject2\\突出核蛋白\\best_model.pth',
  1874. target_detailed_count=200,
  1875. importance_threshold=0.3
  1876. )
  1878. if results is not None:
  1879. print("\nAnalysis results saved in 'results' variable")
  1880. print("You can access them through:")
  1881. print("- results['global_statistics']: Global statistics")
  1882. print("- results['quick_analysis_results']: Quick analysis results")
  1883. print("- results['detailed_analysis_results']: Detailed analysis results")
  1884. print("- results['summary']: Analysis summary")
  1886. display_analysis_results(results)
  1888. print("\nYou can use the following commands for further analysis:")
  1889. print("# View specific atom type importance")
  1890. print("importance_df = results['detailed_analysis_results']['importance_df']")
  1891. print("print(importance_df[importance_df['atom_type'] == 'N'].describe())")
  1893. print("\n# View molecules with highest predictions")
  1894. print("quick_df = pd.DataFrame(results['quick_analysis_results'])")
  1895. print("top_predictions = quick_df.nlargest(10, 'prediction')")
  1896. print("print(top_predictions[['smiles', 'prediction', 'avg_importance']])")

gnnexplainer.py at commit 76ec685, no license · at the source

Overview

Authors: Tingle Gu1, Zixu Ran2, Wenyin Li3, Xudong Guo2, Bo Li4, Fuyi Li2,5, Cangzhi Jia1
ORCID iDs: Fuyi Li, Cangzhi Jia
  1. School of Science, Dalian Maritime University, No. 1 Linghai Road, Dalian 116026, Liao Ning, China
  2. College of Information Engineering, Northwest A&F University, No. 3 Taicheng Road, Yangling 712100, Shanxi, China
  3. The First Clinical College, Liaoning University of Traditional Chinese Medicine, No. 79, Chongshen East Road, Huanggu District, Shenyang 110847, Liaoning, China
  4. Department of Dermatology, Dalian Dermatosis Hospital, ChangJiang Road 788, Dalian 116021, Liaoning, China
  5. South Australian immunoGENomics Cancer Institute (SAiGENCI), College of Health, Adelaide University, Adelaide 5005, Australia
Institutions: Dalian Maritime University (China); Northwest A&F University (China); Liaoning University of Traditional Chinese Medicine (China); Affiliated Hospital of Liaoning University of Traditional Chinese Medicine (China); Dalian Dermatosis Hospital (China); Adelaide University (Australia); The University of Adelaide (Australia)
Journal: Briefings in bioinformatics, volume 27, issue 2, article bbag118
Dates: received 5 November 2025; accepted 20 February 2026; published online 23 March 2026; in print March 2026
Type: Research article · Language: English
License: CC BY
Identifiers: DOI 10.1093/bib/bbag118 · PMID 41870129 · PMCID PMC13006971 · OpenAlex W7140079345
Open access: gold, a free copy (OpenAlex)
Status: code verified
Categories: computational modeling (no new data) (modality), human (organism), Parkinson's (population), cellular / molecular (subfield)
Methods: Machine learning
Keywords: α-synuclein, composite regularization, QSAR, graph contextual attention, dual-channel feature fusion
MeSH: alpha-Synuclein*, Drug Evaluation, Preclinical, Graph Neural Networks, Humans, Ligands, Molecular Docking Simulation (* major topic)
Topic: Computational Drug Discovery Methods (Computational Theory and Mathematics, Computer Science), according to OpenAlex
Funding: National Natural Science Foundation of China (62071079, 62202388); Hainan Normal University, Ministry of Education (JSKX202203)
Citations: cited by 1 paper (Europe PMC); 24 references in the paper

Abstract

The pathological aggregation of α-synuclein (α-syn) constitutes a pivotal hallmark in the progression of neurodegenerative disorders, including Parkinson’s disease, underscoring the imperative need for identifying site-specific ligands. This study presents, for the first time, an advanced deep learning framework specifically designed for the prediction of molecular properties associated with α-syn. The framework integrates graph-based contextual attention mechanisms, structural feature aggregation protocols, and dual-channel feature integration, complemented by a composite regularization strategy that synergizes mean squared error minimization, Kullback–Leibler divergence–induced latent space regularization, and L2 norm penalization, thereby delivering outstanding predictive accuracy on the independent test dataset with MSE of 0.1812. Mechanistic insights derived from GNNExplainer analysis and molecular docking studies (PDB: 6A6B) elucidated that aromatic ring systems (benzene ring significance: 0.737) and hydrogen bond donor groups (amino group significance: 0.438) play critical roles in mediating high-affinity ligand–receptor interactions through π–π stacking within the hydrophobic pocket formed by Val82 and Ala89 residues, as well as directed hydrogen bonding involving catalytic residues Ser42 and Lys45. These findings not only enhance the understanding of inhibitor mechanisms but also establish a novel framework for the preliminary screening of small-molecule therapeutics, thereby laying a rigorous groundwork for structure-guided drug optimization and rational molecular design.

Reproduced under the paper's license (CC BY), from the paper cited above.

Repository

Its files are read in the Code ↔ Paper reader above, with 6 matches between paragraphs and lines of code.

JiaCZ-Computational-Biology/M-GAT-GraphSAGE

License: none: the authors keep all their rights
State: the link answers, verified on 30 September 2026
Evidence: files inventoried
Commit: 76ec685148bad7b3a0362e816cc24b70f4946ddc, 4 November 2025
Languages: Python (43)
Size: 44 files, 43 scripts
Software Heritage: not archived
Found in: “Data availability”
Holds: README
Not found: license file, CITATION.cff, environment file, tests, continuous integration, documentation
Tools: NumPy (43 files), pandas (43 files), RDKit (43 files), PyTorch Geometric (41 files), PyTorch (41 files), scikit-learn (27 files), SciPy (17 files), Matplotlib (3 files), seaborn (2 files), LightGBM (1 file), NetworkX (1 file), Pillow (1 file), XGBoost (1 file)
Availability: 1 check, the latest on 30 September 2026: the link answers
  • 30 September 2026: the link answers
44 files
At the source: github.com/JiaCZ-Computational-Biology/M-GAT-GraphSAGE

The paper's code and data availability statement is in the Data section.

Tracing map

Proposed by the machine: these links were found in the paper and verified at the source, without human review. The map will receive a Zenodo DOI once one of the paper's authors has validated it with their ORCID.

What the map holds:

  • 1 repository of the authors' code, each at its verified commit, with its license and how the link was found in the paper;
  • 43 scripts, each with its path and the digest of its content;
  • 6 matches between paragraphs of the paper and lines of the code (method lexical-v1);
  • neither the text of the paper nor the code itself.

Its JSON (tracing-map.json) is deposited on Zenodo with its DOI once the map is validated.

Data

No dataset and no data link were found in the paper.

Data availability

The source code and the package are available at https://github.com/JiaCZ-Computational-Biology/M-GAT-GraphSAGE.

Reproduced under the paper's license (CC BY), from the paper cited above.

Versions

The history of this record: each version stored by the harvester or made by a correction of its authors or of the maintainers of its code, and what changed in its facts. The texts of the paper (its abstract, its availability statements) are not part of it; versions that changed only those are not listed.

Version 1, 30 September 2026: the first record

Recorded: type, language, journal, volume, issue, pages, dates, 7 authors, 5 keywords, 6 MeSH terms, 2 funders, 21 references.

Cite

This paper

Gu, T., Ran, Z., Li, W., Guo, X., Li, B., Li, F., & Jia, C. (2026). Drug screening for α-synuclein aggregation inhibitors via multimodal graph neural network. Briefings in bioinformatics, 27(2), bbag118. https://doi.org/10.1093/bib/bbag118

BibTeX

@article{gu2026drug,
author = {Gu, Tingle and Ran, Zixu and Li, Wenyin and Guo, Xudong and Li, Bo and Li, Fuyi and Jia, Cangzhi},
title = {{Drug screening for α-synuclein aggregation inhibitors via multimodal graph neural network}},
journal = {Briefings in bioinformatics},
year = {2026},
month = mar,
volume = {27},
number = {2},
pages = {bbag118},
publisher = {Oxford University Press},
issn = {1467-5463},
doi = {10.1093/bib/bbag118},
url = {https://doi.org/10.1093/bib/bbag118},
pmid = {41870129},
pmcid = {PMC13006971}
}

RIS

TY - JOUR
AU - Gu, Tingle
AU - Ran, Zixu
AU - Li, Wenyin
AU - Guo, Xudong
AU - Li, Bo
AU - Li, Fuyi
AU - Jia, Cangzhi
TI - Drug screening for α-synuclein aggregation inhibitors via multimodal graph neural network
T2 - Briefings in bioinformatics
J2 - Brief Bioinform
PY - 2026
DA - 2026/03/01
VL - 27
IS - 2
SP - bbag118
SN - 1467-5463
PB - Oxford University Press
DO - 10.1093/bib/bbag118
UR - https://doi.org/10.1093/bib/bbag118
LA - en
ER -

CSL-JSON

{
"id": "10.1093/bib/bbag118",
"type": "article-journal",
"title": "Drug screening for α-synuclein aggregation inhibitors via multimodal graph neural network",
"container-title": "Briefings in bioinformatics",
"author": [
{
"family": "Gu",
"given": "Tingle"
},
{
"family": "Ran",
"given": "Zixu"
},
{
"family": "Li",
"given": "Wenyin"
},
{
"family": "Guo",
"given": "Xudong"
},
{
"family": "Li",
"given": "Bo"
},
{
"family": "Li",
"given": "Fuyi"
},
{
"family": "Jia",
"given": "Cangzhi"
}
],
"container-title-short": "Brief Bioinform",
"volume": "27",
"issue": "2",
"page": "bbag118",
"DOI": "10.1093/bib/bbag118",
"PMID": "41870129",
"PMCID": "PMC13006971",
"ISSN": "1467-5463",
"publisher": "Oxford University Press",
"URL": "https://doi.org/10.1093/bib/bbag118",
"language": "en",
"issued": {
"date-parts": [
[
2026,
3,
1
]
]
}
}

The tracing map gets a citation of its own once an author has validated it and it has a DOI.

Similar papers

The papers with a page that share the most with this one: the tools found in their code, their categories, datasets, cited references and authors, the rarest counting most.

[1] doi:10.1038/s42003-026-10957-8 [code]
Brain defence by the extracellular matrix protein Cochlin.
Journal: Communications biology
In common: RDKit, XGBoost, NetworkX, 8 other tools, cellular / molecular
[2] doi:10.7554/elife.110588 [code]
Opening the black box toward a modular approach to spike sorting.
Journal: eLife
In common: LightGBM, XGBoost, NetworkX, 8 other tools
[3] doi:10.1016/j.isci.2026.116055 [code]
Mapping the transcriptional diversity of calcium signaling in the mouse and human brain.
Journal: iScience
In common: LightGBM, PyTorch Geometric, NetworkX, 7 other tools
[4] doi:10.1371/journal.pone.0345854 [code]
Shedding light on neural learning to rank models for anticancer drug prioritization.
Journal: PloS one
In common: RDKit, PyTorch Geometric, NetworkX, 7 other tools
[5] doi:10.1093/nar/gkag706 [code]
scDifformer: diffusion-based post-training for virtual cell modeling across large-scale single-cell data.
Journal: Nucleic acids research
In common: RDKit, PyTorch Geometric, NetworkX, 7 other tools
[6] doi:10.1186/s13244-026-02365-7 [code]
Super-resolution MRI and 2.5D deep learning for intratumoral-peritumoral radiomics in preoperative prediction of rectal cancer perineural invasion.
Journal: Insights into imaging
In common: LightGBM, XGBoost, Pillow, 7 other tools
[7] doi:10.1038/s41586-026-10670-w [code]
Zero-shot design of drug-binding proteins via neural iterative selection-expansion.
Journal: Nature
In common: RDKit, PyTorch Geometric, PyTorch, 6 other tools, computational modeling (no new data), cellular / molecular
[8] doi:10.1038/s41398-026-03965-z [code]
Disentangling individual heterogeneity reveals robust network and molecular signatures of major depressive disorder with suicidal ideation.
Journal: Translational psychiatry
In common: PyTorch Geometric, NetworkX, Pillow, 7 other tools, cellular / molecular
[9] doi:10.1002/hbm.70469 [code]
VarCoNet: A Variability-Aware Self-Supervised Framework for Functional Connectome Extraction From Resting-State fMRI.
Journal: Human brain mapping
In common: PyTorch Geometric, NetworkX, Pillow, 7 other tools
[10] doi:10.1371/journal.pcbi.1014555 [code]
Body surface potential driven personalisation of electrophysiological digital twins in hypertrophic cardiomyopathy.
Journal: PLoS computational biology
In common: XGBoost, NetworkX, Pillow, 7 other tools

Contribute

The authors of this paper can claim it, correct its record and validate its tracing map, and the maintainers of its code (its owner, or a public member of its organization) correct what it says of their repository; anyone signed in can ask for its removal. Every request goes to OSCR's own machine, which answers it; your account page follows them.

Sign in with ORCID to claim this paper as one of its authors, correct its record or validate its tracing map: when the paper's metadata lists your ORCID iD, you are recognized at once. Maintainers of its code: sign in with GitHub, then claim the repository on your account page.

Request its removal

To ask OSCR to remove this record, the copies of its authors' scripts or its tracing map, use the removal request page: signed in, you say who you are, what to remove and why, then review and confirm the request. Published rules decide every request (how).

Discussion, reproductions, activity

Discussion: questions and error reports about this paper and its code, from signed-in readers and its authors. It opens with sign-in.

Reproductions: reports from readers who ran the authors' code: what they reproduced, with which environment, commit and data. It opens with sign-in.

Activity: what happens around this paper: new versions of its record, its map's validation, discussions and reproductions. It opens with sign-in.