analyze

REINVENT4 Molecule Analysis

Evaluate generated molecules using Python/RDKit/pandas directly. No MCP server needed — requires RDKit, pandas, and optionally umap-learn.

Workflow

1. Resolve paths

readlink -f <relative_path>

2. Identify inputs

• smiles_file (required): sampling CSV or .smi file. Use $0 if provided. SMILES column auto-detected from CSV headers (first column containing "smiles", case-insensitive).
• ref_smiles_file (optional): reference / known-active SMILES for novelty and similarity. Use $1 if provided.
• output_dir: default <stem>_analysis/ next to the input file.

3. Load molecules

import pandas as pd
from rdkit import Chem
from pathlib import Path

path = Path("/absolute/path/to/sampling.csv")

if path.suffix == ".csv":
    df = pd.read_csv(path)
    smi_col = next((c for c in df.columns if "smiles" in c.lower()), df.columns[0])
    raw_smiles = df[smi_col].dropna().astype(str).tolist()
else:
    raw_smiles = [l.split()[0] for l in path.read_text().splitlines() if l.strip()]

mols, valid_smiles = [], []
for smi in raw_smiles:
    mol = Chem.MolFromSmiles(smi)
    if mol:
        mols.append(mol)
        valid_smiles.append(Chem.MolToSmiles(mol))

print(f"Total: {len(raw_smiles)}, Valid: {len(mols)} ({100*len(mols)/max(len(raw_smiles),1):.1f}%)")

4. Physicochemical properties

from rdkit.Chem import Descriptors, rdMolDescriptors
import numpy as np

props = []
for mol in mols:
    props.append({
        "mw":      Descriptors.MolWt(mol),
        "logp":    Descriptors.MolLogP(mol),
        "tpsa":    Descriptors.TPSA(mol),
        "hbd":     rdMolDescriptors.CalcNumHBD(mol),
        "hba":     rdMolDescriptors.CalcNumHBA(mol),
        "rotbonds":rdMolDescriptors.CalcNumRotatableBonds(mol),
        "rings":   rdMolDescriptors.CalcNumRings(mol),
    })

df_props = pd.DataFrame(props)
print(df_props.describe().round(2))

5. Druglikeness

from rdkit.Chem import QED

qeds = [QED.qed(mol) for mol in mols]

# Lipinski RO5
def lipinski(mol):
    return (Descriptors.MolWt(mol) <= 500 and Descriptors.MolLogP(mol) <= 5 and
            rdMolDescriptors.CalcNumHBD(mol) <= 5 and rdMolDescriptors.CalcNumHBA(mol) <= 10)

ro5_pass = sum(1 for mol in mols if lipinski(mol))
print(f"QED — mean: {np.mean(qeds):.3f}, median: {np.median(qeds):.3f}")
print(f"Lipinski RO5 pass: {ro5_pass}/{len(mols)} ({100*ro5_pass/max(len(mols),1):.1f}%)")