Molgen — Computer-Aided Drug Design (CADD) Pipeline

End-to-end pipeline for automatic generation of new drug-like molecules, their chemical/pharmacological filtering, structural preparation, molecular docking, and generation of interactive result reports.

The project combines a deep generative model (GraphAF, via TorchDrug) with RDKit for chemical analysis and with the Schrödinger Suite (LigPrep/Glide) for ligand preparation and docking.

Typical paths used by the scripts: ~/CADD/<user>/molgen/{batch, filtered, ligprep, docking, reports, pubchem}

---

🔬 Process Overview

1. MOLECULE GENERATION (GraphAF / TorchDrug)
   generate_molecules*.py
        │  produces batches of valid SMILES (valid_smiles_50k_batch_N)
        ▼
2. BATCH COPY/ORGANIZATION
   my_script.py
        │  moves new batches into molgen/batch
        ▼
3. STRUCTURAL FILTERING / CONVERSION (Schrödinger structconvert)
   filter_mols.sh, filter_smiles.sh
        │  SMILES → MAE → SMILES (structure validation/normalization)
        ▼
4. (Optional) PUBCHEM MATCHING
   pubchem_filter.py, pubchempy.py
        │  matches PubChem CIDs to molecules via canonical SMILES
        ▼
5. LIGAND PREPARATION (Schrödinger LigPrep)
   prepare_mols.sh
        │  generates protonation states/tautomers (Epik) in .maegz format
        ▼
6. MOLECULAR DOCKING (Schrödinger Glide)
   docking_mols.sh
        │  SP docking against a target grid, computes r_i_glide_gscore
        ▼
7. DOCKING SCORE EXTRACTION AND ANALYSIS
   extract_docking_scores.ipynb
        │  aggregates and sorts docking scores across batches
        ▼
8. DRUG-LIKENESS + PAINS FILTERS + HTML REPORT
   filter_smarts.ipynb, batch_report.py
        │  descriptor calculation, filters (Lipinski/Ghose/Veber/Rule of 3/REOS/Drug-like),
        │  PAINS screening, similarity clustering, Murcko scaffold
        ▼
   Interactive HTML report (DataTables) with molecules, scores, and properties

---

📁 File Descriptions

File	Description
generate_molecules.py	Main script for molecule generation with GraphAF (autoregressive graph model, TorchDrug). Trains (or loads) the model on the train.csv dataset, generates SMILES in batches of 50,000, validates molecules with RDKit/TorchDrug, and avoids duplicates against both the training set and previously generated batches. Assigns sequential IDs (GM<n>) and saves the batches to CSV.
generate_molecules_drugbank_rf.py	Variant of the generation script (likely trained on/referring to a DrugBank-type dataset), without explicit filtering against the training set or SMILES canonicalization on output.
generate_molecules_parallel.py	Experimental variant that (incompletely) introduces multiprocessing (multiprocessing.Pool) to parallelize generation; does not filter out molecules already present in the training set.
my_script.py	Utility script that copies newly generated SMILES batches from Documents/molgen into the working folder molgen/batch, and sets up the ligprep folder for the next step.
filter_mols.sh / filter_smiles.sh	Bash scripts that use Schrödinger's structconvert to convert .smi files to .mae and back, as a structural validation/normalization step before docking. filter_mols.sh automatically manages the queue of unprocessed batches; filter_smiles.sh is a parametric version (input/output directories passed as command-line arguments).
pubchem_filter.py	Compares the generated canonical SMILES against a local PubChem dump (CID-SMILES) to retrieve, where available, the corresponding PubChem CID (exact match on canonical string, much faster than API queries).
pubchempy.py	Variant that queries PubChemPy (PubChem REST API) to retrieve the CID of each molecule via structural identity search; slower, but does not require a local database dump.
prepare_mols.sh	Prepares the filtered ligands with LigPrep (Schrödinger): generates protonation states at pH 7 (±2) using Epik, tautomers, and stereoisomers, producing .maegz files ready for docking.
docking_mols.sh	Runs molecular docking with Glide (SP precision, OPLS4 force field) of the prepared molecules against a target grid (gridfile.zip), processing molecules not yet docked.
extract_docking_scores.ipynb	Notebook that reads Glide results (r_i_glide_gscore) for all docked batches, sorts them, and aggregates them into a single DataFrame to identify the best hits.
filter_smarts.ipynb	Prototype notebook (later merged into batch_report.py) that: computes molecular descriptors, applies drug-likeness filters, performs PAINS screening via SMARTS patterns, computes fingerprint/similarity and hierarchical clustering, generates the Murcko scaffold image, and produces the final HTML report.
batch_report.py	Production/automated version of filter_smarts.ipynb: for each docked batch not yet processed, filters by docking score (< -6), computes descriptors and drug-likeness filters, performs PAINS screening, similarity-based clustering (RDKit fingerprints + hierarchical linkage), and generates a Murcko scaffold; produces an interactive .html report (a DataTables table with molecule images, filterable columns, and cluster-based ordering).

---

🧠 Models Used

GraphAF (Generative Flow-based Autoregressive Graph model)

A generative model for molecular graphs based on autoregressive normalizing flows, which generates molecules atom-by-atom/bond-by-bond sequentially, ensuring chemical validity at every step.

• Backbone: RGCN (Relational Graph Convolutional Network) — hidden_dims=[256, 256, 256]
• Prior: independent Gaussian distributions for nodes (atoms) and edges (bonds)
• Two separate flows: node_flow (atom types) and edge_flow (bond types), 12 layers each
• Training task: AutoregressiveGeneration with NLL (negative log-likelihood) criterion
• Training dataset: train.csv with a SMILES column (loaded as a custom MoleculeDataset class, AirPollutants)
• Reference paper: GraphAF: a Flow-based Autoregressive Model for Molecular Graph Generation (Shi et al., ICLR 2020) — https://arxiv.org/abs/2001.09382

"Classical" (non-neural) Scoring/Filtering Models

Applied downstream of generation, on filtered/docked molecules:

• Drug-likeness rules: Lipinski's Rule of Five, Ghose Filter, Veber Filter, Rule of 3, REOS, generic Drug-like filter
• PAINS (Pan-Assay Interference Compounds) — screening via SMARTS patterns from the WEHI_PAINS dataset (Baell & Holloway, 2010)
• Murcko scaffold — extraction of the core molecular skeleton (MurckoScaffold)
• Docking score — Glide SP r_i_glide_gscore (Schrödinger)

---

📚 Main Libraries and Tools

Library/Tool	Role in the Project	Link
TorchDrug	Deep learning framework for molecules/graphs; GraphAF model, molecular datasets, training engine	https://torchdrug.ai/ · https://github.com/DeepGraphLearning/torchdrug
PyTorch	Deep learning backend for TorchDrug	https://pytorch.org/
RDKit	Computational chemistry: SMILES parsing, molecular descriptors, fingerprints, SMARTS, molecule drawing, Murcko scaffold	https://www.rdkit.org/ · https://github.com/rdkit/rdkit
PubChemPy	Python interface to the PubChem REST API (CID lookup)	https://pubchempy.readthedocs.io/
pandas	Tabular data manipulation (SMILES, scores, descriptors)	https://pandas.pydata.org/
NumPy	Numerical computation (similarity matrices)	https://numpy.org/
SciPy (scipy.cluster.hierarchy)	Hierarchical clustering of molecules by similarity (linkage, dendrogram)	https://docs.scipy.org/doc/scipy/reference/cluster.hierarchy.html
Matplotlib	Supports visualizations (imported in batch_report.py)	https://matplotlib.org/
Schrödinger Suite (LigPrep, Glide, structconvert)	Ligand preparation (protonation/tautomers with Epik), molecular docking, structural file format conversion	https://www.schrodinger.com/
DataTables (jQuery)	Interactive, sortable, filterable HTML tables in the generated reports	https://datatables.net/
WEHI PAINS filters	Dataset of SMARTS patterns for identifying PAINS compounds (assay false positives)	http://cbligand.org/PAINS/

---

⚙️ Requirements

• Python ≥ 3.8
• RDKit
• TorchDrug + PyTorch (GPU support recommended for GraphAF training)
• pandas, numpy, scipy, matplotlib
• PubChemPy (only for pubchempy.py)
• Schrödinger Suite with a valid license (LigPrep, Glide, structconvert) for the .sh scripts
• Local wehi_pains.csv dataset (PAINS SMARTS patterns) for PAINS screening

⚠️ The scripts contain environment-specific absolute paths (e.g. ~/CADD/dsardina/molgen/..., /home/dssardina@Fondazionerimed.com/...) that need to be adapted to your own setup before use.

---

🚀 Typical Execution

# 1. Generate new molecules with GraphAF
python generate_molecules.py

# 2. Organize the generated batches
python my_script.py

# 3. Structural filtering/validation
./filter_mols.sh

# 4. (Optional) PubChem matching
python pubchem_filter.py

# 5. Ligand preparation (LigPrep)
./prepare_mols.sh

# 6. Molecular docking (Glide)
./docking_mols.sh path/to/gridfile.zip

# 7. Docking score analysis
jupyter notebook extract_docking_scores.ipynb

# 8. Final report generation (filters + PAINS + clustering)
python batch_report.py

---

📄 License

GNU General Public License v3.0