{
  "skill": "alterlab-medchem",
  "evals": [
    {
      "id": "library-triage-rules",
      "prompt": "Triage my 5000-compound library: keep only molecules that pass Rule of Five and Veber, and drop anything with a known structural alert. Give me back the clean set.",
      "expected_output": "Invokes alterlab-medchem: builds mc.rules.RuleFilters(rule_list=['rule_of_five','rule_of_veber']) and applies it with n_jobs=-1/progress=True, then runs mc.structural.CommonAlertsFilters to flag alerts, and combines pass-rules AND not-has-alerts into a filtered DataFrame for the clean set.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Applies medchem rule filters (Ro5, Veber) and structural alerts with parallel n_jobs, then intersects the pass conditions." }
      ]
    },
    {
      "id": "pains-flagging",
      "prompt": "Before I trust these screening hits, flag any PAINS or reactive Michael-acceptor groups so I don't chase assay artifacts.",
      "expected_output": "Invokes alterlab-medchem: applies PAINS filtering (CommonAlertsFilters with the PAINS set, or HASALERT(\"pains\") in the query DSL) and flags reactive/Michael-acceptor motifs via the structural-alert sets (common alerts / NIBR) or the electrophilic-warhead ChemicalGroup catalogs (there is no 'michael_acceptors' group), returning per-molecule flags so the user can deprioritize likely assay-interference artifacts.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "output_contains", "value": "PAINS" }
      ]
    },
    {
      "id": "lilly-nibr-lead-opt",
      "prompt": "For lead optimization, run the NIBR filter set and Lilly demerits on my candidate list and tell me which ones exceed 100 demerits.",
      "expected_output": "Invokes alterlab-medchem: applies mc.structural.NIBRFilters and the Lilly demerits filter (mc.functional.lilly_demerit_filter, which needs the Lilly MedChem Rules tools installed via `medchem install-lilly`) with max_demerits=100 to honor the user's threshold (noting medchem's default ceiling is 160), reporting which molecules are rejected, for stricter lead-optimization triage.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Uses NIBRFilters plus the Lilly demerit filter (via medchem.functional.lilly_demerit_filter, flagging the external-binary requirement) and applies the 100-demerit threshold the user asked for." }
      ]
    },
    {
      "id": "complexity-and-constraints",
      "prompt": "Filter my set to molecular weight 200-500, logP between -2 and 5, TPSA under 140, and also drop anything with high Bertz molecular complexity.",
      "expected_output": "Invokes alterlab-medchem: enforces the MW/logP/TPSA property windows via the query DSL HASPROP (or mc.rules.in_range over computed descriptors, since medchem has no all-in-one Constraints property object) and applies the Bertz complexity filter (mc.complexity.ComplexityFilter / mc.functional.complexity_filter with complexity_metric='bertz' and a percentile limit), returning molecules that satisfy all conditions.",
      "assertions": [
        { "type": "should_trigger", "value": true },
        { "type": "behavior", "value": "Enforces the MW/logP/TPSA windows via HASPROP or in_range (NOT a fabricated Constraints(mw_range=...) object) and applies the Bertz complexity filter via ComplexityFilter/complexity_filter with complexity_metric='bertz'." }
      ]
    },
    {
      "id": "near-miss-datamol",
      "prompt": "I just need to compute MW, logP, HBD, HBA and TPSA descriptors in bulk for my SDF and standardize the structures first. No filtering decisions, just the numbers.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-datamol. The user only wants descriptor computation and standardization, not medicinal-chemistry filtering / alert flagging, so datamol's batch descriptor and standardization utilities are the right tool.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-datamol" }
      ]
    },
    {
      "id": "near-miss-pytdc",
      "prompt": "I want a labeled hERG cardiotoxicity dataset with a proper scaffold train/test split to train a toxicity classifier.",
      "expected_output": "Does NOT invoke this skill; defers to alterlab-pytdc. The user wants a curated, AI-ready benchmark dataset (Tox hERG) with scaffold splits for model training, not rule-based filtering of their own compounds, which is Therapeutics Data Commons territory.",
      "assertions": [
        { "type": "should_not_trigger", "value": true },
        { "type": "output_contains", "value": "alterlab-pytdc" }
      ]
    }
  ]
}

# Medchem API Reference

Verified against `medchem==2.1.0`. For the authoritative, always-current API, run
`mc.rules.RuleFilters.list_available_rules()`,
`mc.structural.CommonAlertsFilters.list_default_available_alerts()`,
`mc.groups.list_default_chemical_groups()`, and
`mc.catalogs.list_named_catalogs()` — and read https://medchem-docs.datamol.io/.

## Cross-cutting conventions

- Inputs accept SMILES strings **or** RDKit mols (`Sequence[Union[str, Chem.Mol]]`).
- Filter **classes** (`medchem.rules.RuleFilters`, `medchem.structural.*`) return a **pandas DataFrame**, one row per input mol, with a `mol` column.
- Filter **functions** (`medchem.functional.*`) return a **NumPy boolean array** (`True` = keep); pass `return_idx=True` for indices instead.
- Common kwargs: `n_jobs` (default `-1` = all cores), `progress` (bool), `scheduler` (`"auto"`/`"threads"`/`"processes"`).

---

## medchem.rules

### RuleFilters

```python
RuleFilters(rule_list: List[Union[str, Callable]], rule_list_names: Optional[List[Optional[str]]] = None)
__call__(mols, n_jobs=-1, progress=False, keep_props=False, fail_if_invalid=True, ...) -> pandas.DataFrame
```

Returns a DataFrame with columns: `mol`, `pass_all`, `pass_any`, and one boolean column per rule.

- `RuleFilters.list_available_rules()` -> DataFrame of rule names + metadata.

### medchem.rules.basic_rules

Single-molecule rule functions returning `bool`:

```python
rule_of_five(mol, mw=None, clogp=None, n_lipinski_hbd=None, n_lipinski_hba=None, **kwargs) -> bool
```

Available (subset): `rule_of_five`, `rule_of_five_beyond`, `rule_of_four`, `rule_of_three`,
`rule_of_three_extended`, `rule_of_two`, `rule_of_ghose`, `rule_of_veber`, `rule_of_reos`,
`rule_of_egan`, `rule_of_pfizer_3_75`, `rule_of_gsk_4_400`, `rule_of_oprea`, `rule_of_xu`,
`rule_of_cns`, `rule_of_respiratory`, `rule_of_zinc`, `rule_of_leadlike_soft`,
`rule_of_druglike_soft`, `rule_of_generative_design`, `rule_of_chemaxon_druglikeness` (needs ChemAxon).

> No `rule_of_drug`, `rule_of_leadlike_strict`, `golden_triangle`, or `pains_filter`. Use the
> alert system for PAINS and `rule_of_leadlike_soft`/`rule_of_oprea` for lead-likeness.

### Helpers

- `in_range(x, min_val=-inf, max_val=inf) -> bool` — gate a numeric descriptor.
- `list_descriptors()` -> names usable for property windows: `mw`, `clogp`, `tpsa`, `fsp3`, `qed`,
  `sas`, `n_lipinski_hbd`, `n_lipinski_hba`, `n_rotatable_bonds`, `n_rings`, `n_aromatic_rings`,
  `n_heavy_atoms`, `n_hetero_atoms`, `formal_charge`, ...
- `n_fused_aromatic_rings`, `n_heavy_metals`, `fraction_atom_in_scaff`, `has_spider_chains`.

---

## medchem.structural

### CommonAlertsFilters

```python
CommonAlertsFilters(alerts_set: Union[str, List[str], None] = None, alerts_db_path=None)
__call__(mols, n_jobs=-1, progress=False, batch_size=None, keep_details=False) -> pandas.DataFrame
```

Returns DataFrame columns: `mol`, `pass_filter` (bool; `True` = clean), `status` (`"ok"`/`"exclude"`),
`reasons` (`;`-joined alert names, `NaN` when clean).

- `CommonAlertsFilters.list_default_available_alerts()` -> DataFrame of alert sets
  (`Glaxo`, `Dundee`, `BMS`, `PAINS`, `SureChEMBL`, ...). Pass a subset via `alerts_set=["PAINS"]`.

### NIBRFilters

```python
NIBRFilters()
__call__(mols, n_jobs=-1, progress=False, keep_details=False) -> pandas.DataFrame
```

Returns DataFrame including `mol`, `pass_filter`, `severity`, `status`, `reasons`,
`n_covalent_motif`, `special_mol`.

### Lilly demerits

The class is `medchem.structural.lilly_demerits.LillyDemeritsFilters` (not exported at
`medchem.structural` top level) and **requires the compiled Lilly MedChem Rules tools**
(`medchem install-lilly` in medchem ≥ 2.1 builds the pinned upstream release; needs `make`,
a C++ compiler and zlib; not native Windows). Running it without them raises `ImportError`.
Prefer the functional entry point (below).

---

## medchem.functional

All return a NumPy bool array (`True` = keep); `return_idx=True` returns passing indices.

```python
rules_filter(mols, rules: Union[List, RuleFilters], n_jobs=None, progress=False, ...) -> np.ndarray
alert_filter(mols, alerts: List[str], alerts_db=None, n_jobs=1, progress=False, ...) -> np.ndarray
nibr_filter(mols, n_jobs=None, max_severity=10, progress=False, ...) -> np.ndarray
lilly_demerit_filter(mols, max_demerits=160, n_jobs=None, progress=False, ...) -> np.ndarray  # needs binaries
complexity_filter(mols, complexity_metric="bertz", limit="99", threshold_stats_file="zinc_15_available", ...) -> np.ndarray
chemical_group_filter(mols, chemical_group: ChemicalGroup, exact_match=False, ...) -> np.ndarray
```

> `alert_filter` alert names are **lowercase** (`"pains"`, `"tox"`, `"nih"`, ...). There is no
> `common_alerts_filter` / `lilly_demerits_filter`; the names are `alert_filter` / `lilly_demerit_filter`.

Other helpers: `atom_list_filter`, `bredt_filter`, `catalog_filter`, `halogenicity_filter`,
`macrocycle_filter`, `molecular_graph_filter`, `num_atom_filter`, `num_stereo_center_filter`,
`protecting_groups_filter`, `ring_infraction_filter`, `symmetry_filter`.

---

## medchem.groups

### ChemicalGroup

```python
ChemicalGroup(groups: Union[str, List[str], None] = None, n_jobs=None, groups_db=None)
has_match(mol, exact_match=False, terminal_only=False) -> bool
get_matches(mol, use_smiles=True, exact_match=False, terminal_only=False)
filter(names: List[str], fuzzy=False)
# also: dataframe, get_catalog, list_groups, mols, smarts, smiles
```

- `list_default_chemical_groups()` -> catalog group names (e.g. `hinge_binders`,
  `electrophilic_warheads_for_kinases`, `common_warhead_covalent_inhibitors`,
  `privileged_kinase_inhibitor_scaffolds`, `aggregator`, `privileged_scaffolds`).
- `list_functional_group_names()` -> ~579 fine-grained functional-group names used by the
  query DSL `HASGROUP(...)`.

> `phosphate_binders`, `michael_acceptors`, and `reactive_groups` are **not** default group names.

---

## medchem.complexity

```python
ComplexityFilter(limit="99", complexity_metric="bertz", threshold_stats_file="zinc_15_available")
__call__(mol) -> bool        # per molecule; True = within the limit
```

Metric classes for raw scores: `BaroneCT`, `WhitlockCT`, `SMCM`, `TWC` (metric strings:
`"bertz"`, `"barone"`, `"whitlock"`, `"smcm"`, `"twc"`). There is no `calculate_complexity` and no
`max_complexity` argument; the filter thresholds against a percentile of a reference set.

---

## medchem.catalogs

```python
NamedCatalogs.pains()   # RDKit FilterCatalog; also .pains_a/.pains_b/.pains_c, .brenk, .nih, .bms, ...
list_named_catalogs() -> List[str]   # ['tox','pains','pains_a','pains_b','pains_c','nih', ...]
catalog_from_smarts(...)             # build a catalog from custom SMARTS
merge_catalogs(...)
```

---

## medchem.constraints

```python
Constraints(core: Chem.Mol, constraint_fns: Dict[str, Callable], prop_name="query")
has_match(...) / validate(...) / get_matches(...)
```

Enforces **substructure / R-group** constraints around a query `core` — **not** physchem
property windows. For MW/logP/TPSA/HBD windows use `medchem.rules.in_range` over computed
descriptors, or the query DSL `HASPROP` (below).

---

## medchem.query

```python
QueryFilter(query: str, grammar=None, parser="lalr")
__call__(mols, n_jobs=-1, progress=False) -> array of bool
```

Structured boolean grammar (not free-form text). Primitives:

- `MATCHRULE("rule_of_five")`
- `HASPROP("mw" < 500)` — comparison operators `< > <= >= == !=`; descriptor names from `list_descriptors()`
- `HASALERT("pains")` — alert-set names are lowercase
- `HASGROUP("<functional_group_name>")`
- `HASSUBSTRUCTURE(...)`, `HASSUPERSTRUCTURE(...)`, `LIKE(...)`
- Combine with `AND`, `OR`, `NOT`, parentheses.

```python
qf = mc.query.QueryFilter('MATCHRULE("rule_of_five") AND HASPROP("mw" < 500) AND NOT HASALERT("pains")')
keep = qf(mol_list, n_jobs=-1)
```

> There is no `mc.query.parse()`. Construct `QueryFilter(query_string)` and call it directly.

---

## Working with DataFrames

```python
import pandas as pd
import datamol as dm
import medchem as mc

df = pd.read_csv("molecules.csv")
mols = [dm.to_mol(smi) for smi in df["smiles"]]

res = mc.rules.RuleFilters(rule_list=["rule_of_five", "rule_of_cns"])(mols=mols, n_jobs=-1)
df["passes_ro5"] = res["rule_of_five"].to_numpy()
df["passes_cns"] = res["rule_of_cns"].to_numpy()
df["pass_all"]   = res["pass_all"].to_numpy()

filtered = df[df["pass_all"]]
```

# Medchem Rules and Filters Catalog

Background and selection guidance for medicinal-chemistry rules, structural alerts, and filters.
The rule criteria and literature references below are conceptual background; for the exact
implementation and the canonical name list use `mc.rules.RuleFilters.list_available_rules()`.

> **Verified against `medchem==2.1.0`.** Not every rule below is shipped as a named medchem
> function (e.g. Rule of Drug, strict lead-likeness, and Golden Triangle have no dedicated
> function in this version — see notes). Apply those criteria via `RuleFilters` of the available
> rules plus property windows (`mc.rules.in_range`) or the query DSL.

## Table of Contents

1. [Drug-Likeness Rules](#drug-likeness-rules)
2. [Lead-Likeness Rules](#lead-likeness-rules)
3. [Fragment Rules](#fragment-rules)
4. [CNS Rules](#cns-rules)
5. [Structural Alert Filters](#structural-alert-filters)
6. [Chemical Group Patterns](#chemical-group-patterns)

---

## Drug-Likeness Rules

### Rule of Five (Lipinski)

**Reference:** Lipinski et al., Adv Drug Deliv Rev (1997) 23:3-25

**Purpose:** Predict oral bioavailability

**Criteria:**
- Molecular Weight ≤ 500 Da
- LogP ≤ 5
- Hydrogen Bond Donors ≤ 5
- Hydrogen Bond Acceptors ≤ 10

**Usage:**
```python
mc.rules.basic_rules.rule_of_five(mol)
```

**Notes:**
- One of the most widely used filters in drug discovery
- About 90% of orally active drugs comply with these rules
- Exceptions exist, especially for natural products and antibiotics

---

### Rule of Veber

**Reference:** Veber et al., J Med Chem (2002) 45:2615-2623

**Purpose:** Additional criteria for oral bioavailability

**Criteria:**
- Rotatable Bonds ≤ 10
- Topological Polar Surface Area (TPSA) ≤ 140 Ų

**Usage:**
```python
mc.rules.basic_rules.rule_of_veber(mol)
```

**Notes:**
- Complements Rule of Five
- TPSA correlates with cell permeability
- Rotatable bonds affect molecular flexibility

---

### Rule of Drug (composite — no single function)

**Purpose:** Combined drug-likeness assessment

**Criteria:**
- Passes Rule of Five
- Passes Veber rules
- Does not contain PAINS substructures

**There is no `rule_of_drug` function in medchem 2.1.0.** Compose it explicitly, e.g. with the query DSL:
```python
qf = mc.query.QueryFilter(
    'MATCHRULE("rule_of_five") AND MATCHRULE("rule_of_veber") AND NOT HASALERT("pains")'
)
keep = qf(mol_list, n_jobs=-1)
```

---

### REOS (Rapid Elimination Of Swill)

**Reference:** Walters & Murcko, Adv Drug Deliv Rev (2002) 54:255-271

**Purpose:** Filter out compounds unlikely to be drugs

**Criteria:**
- Molecular Weight: 200-500 Da
- LogP: -5 to 5
- Hydrogen Bond Donors: 0-5
- Hydrogen Bond Acceptors: 0-10

**Usage:**
```python
mc.rules.basic_rules.rule_of_reos(mol)
```

---

### Golden Triangle (no single function)

**Reference:** Johnson et al., J Med Chem (2009) 52:5487-5500

**Purpose:** Balance lipophilicity and molecular weight

**Criteria:**
- 200 ≤ MW ≤ 50 × LogP + 400
- LogP: -2 to 5

**No `golden_triangle` function in medchem 2.1.0.** Implement with computed `mw`/`clogp`
descriptors and `mc.rules.in_range`, or a custom callable passed to `RuleFilters(rule_list=[...])`.

**Notes:**
- Defines an optimal physicochemical space (a triangle on the MW vs LogP plot).

---

## Lead-Likeness Rules

### Rule of Oprea

**Reference:** Oprea et al., J Chem Inf Comput Sci (2001) 41:1308-1315

**Purpose:** Identify lead-like compounds for optimization

**Criteria:**
- Molecular Weight: 200-350 Da
- LogP: -2 to 4
- Rotatable Bonds ≤ 7
- Number of Rings ≤ 4

**Usage:**
```python
mc.rules.basic_rules.rule_of_oprea(mol)
```

**Rationale:** Lead compounds should have "room to grow" during optimization

---

### Rule of Leadlike (Soft)

**Purpose:** Permissive lead-like criteria

**Criteria:**
- Molecular Weight: 250-450 Da
- LogP: -3 to 4
- Rotatable Bonds ≤ 10

**Usage:**
```python
mc.rules.basic_rules.rule_of_leadlike_soft(mol)
```

---

### Rule of Leadlike (Strict) — not shipped

**Purpose:** Restrictive lead-like criteria (more aggressive than `rule_of_leadlike_soft`)

**Criteria (conceptual):**
- Molecular Weight: 200-350 Da
- LogP: -2 to 3.5
- Rotatable Bonds ≤ 7
- Number of Rings: 1-3

**medchem 2.1.0 ships only `rule_of_leadlike_soft`** (there is no `rule_of_leadlike_strict`).
For stricter lead-likeness, combine `rule_of_oprea` with tighter property windows
(`mc.rules.in_range`) or `HASPROP` queries.

---

## Fragment Rules

### Rule of Three

**Reference:** Congreve et al., Drug Discov Today (2003) 8:876-877

**Purpose:** Screen fragment libraries for fragment-based drug discovery

**Criteria:**
- Molecular Weight ≤ 300 Da
- LogP ≤ 3
- Hydrogen Bond Donors ≤ 3
- Hydrogen Bond Acceptors ≤ 3
- Rotatable Bonds ≤ 3
- Polar Surface Area ≤ 60 Ų

**Usage:**
```python
mc.rules.basic_rules.rule_of_three(mol)
```

**Notes:**
- Fragments are grown into leads during optimization
- Lower complexity allows more starting points

---

## CNS Rules

### Rule of CNS

**Purpose:** Central nervous system drug-likeness

**Criteria:**
- Molecular Weight ≤ 450 Da
- LogP: -1 to 5
- Hydrogen Bond Donors ≤ 2
- TPSA ≤ 90 Ų

**Usage:**
```python
mc.rules.basic_rules.rule_of_cns(mol)
```

**Rationale:**
- Blood-brain barrier penetration requires specific properties
- Lower TPSA and HBD count improve BBB permeability
- Tight constraints reflect CNS challenges

---

## Structural Alert Filters

### PAINS (Pan Assay INterference compoundS)

**Reference:** Baell & Holloway, J Med Chem (2010) 53:2719-2740

**Purpose:** Identify compounds that interfere with assays

**Categories:**
- Catechols
- Quinones
- Rhodanines
- Hydroxyphenylhydrazones
- Alkyl/aryl aldehydes
- Michael acceptors (specific patterns)

**Usage** (PAINS lives in the alert system, not a `pains_filter` rule function):
```python
# Functional API (lowercase set name), True = no PAINS / keep:
mc.functional.alert_filter(mol_list, alerts=["pains"], n_jobs=-1)
# Or via the CommonAlerts class restricted to the PAINS set:
mc.structural.CommonAlertsFilters(alerts_set=["PAINS"])(mols=mol_list, n_jobs=-1)
# Or the RDKit catalog: mc.catalogs.NamedCatalogs.pains()
```

**Notes:**
- PAINS compounds show activity in multiple assays through non-specific mechanisms
- Common false positives in screening campaigns
- Should be deprioritized in lead selection

---

### Common Alerts Filters

**Source:** Derived from ChEMBL curation and medicinal chemistry literature

**Purpose:** Flag common problematic structural patterns

**Alert Categories:**
1. **Reactive Groups**
   - Epoxides
   - Aziridines
   - Acid halides
   - Isocyanates

2. **Metabolic Liabilities**
   - Hydrazines
   - Thioureas
   - Anilines (certain patterns)

3. **Aggregators**
   - Polyaromatic systems
   - Long aliphatic chains

4. **Toxicophores**
   - Nitro aromatics
   - Aromatic N-oxides
   - Certain heterocycles

**Usage:**
```python
alert_filter = mc.structural.CommonAlertsFilters()   # or alerts_set=["PAINS", "BMS"]
df = alert_filter(mols=mol_list, n_jobs=-1, progress=True)
```

**Return format** — a pandas DataFrame, one row per molecule:

| column        | meaning                                              |
|---------------|------------------------------------------------------|
| `mol`         | the RDKit molecule                                   |
| `pass_filter` | `True` if clean (no alert triggered)                 |
| `status`      | `"ok"` or `"exclude"`                                |
| `reasons`     | `;`-joined matched alert names (`NaN` when clean)    |

(There is no `check_mol` method in 2.1.0.) Discover the available alert sets with
`mc.structural.CommonAlertsFilters.list_default_available_alerts()`.

---

### NIBR Filters

**Source:** Novartis Institutes for BioMedical Research

**Purpose:** Industrial medicinal chemistry filtering rules

**Features:**
- Proprietary filter set developed from Novartis experience
- Balances drug-likeness with practical medicinal chemistry
- Includes both structural alerts and property filters

**Usage:**
```python
nibr = mc.structural.NIBRFilters()
df = nibr(mols=mol_list, n_jobs=-1)        # DataFrame: mol, pass_filter, severity, status, reasons, ...
# or the functional one-liner returning a NumPy bool array:
keep = mc.functional.nibr_filter(mol_list, max_severity=10, n_jobs=-1)
```

---

### Lilly Demerits Filter

**Reference:** Bruns & Watson, J Med Chem (2012) 55:9763-9772, doi 10.1021/jm301008n
(275 rules developed over ~18 years).

**Purpose:** Identify assay interference and problematic functionalities via a demerit score.

**Mechanism:**
- Each matched pattern adds demerits; molecules above a demerit ceiling are rejected.
- The original paper rejects at >100 demerits; medchem exposes this as the `max_demerits`
  argument (**default 160** in 2.1.0 — set `max_demerits=100` to match the paper).
- High-severity patterns can hard-reject; lower-severity patterns accumulate.

**Requires the Lilly MedChem Rules tools** (`medchem install-lilly`, medchem ≥ 2.1; the
conda-forge `lilly-medchem-rules` 1.0.1 build is obsolete); the call raises `ImportError` if they
are missing. The class lives at
`medchem.structural.lilly_demerits.LillyDemeritsFilters`; use the functional entry point:

```python
keep = mc.functional.lilly_demerit_filter(mol_list, max_demerits=160, n_jobs=-1)
# NumPy bool array: True = within the demerit ceiling (kept)
```

---

## Chemical Group Patterns

`ChemicalGroup(groups=[...])` matches curated catalog groups. **Validate names against
`mc.groups.list_default_chemical_groups()`** — verified default groups include:

| group name                              | use                                  |
|-----------------------------------------|--------------------------------------|
| `hinge_binders`                         | kinase hinge-binding motifs          |
| `electrophilic_warheads_for_kinases`    | covalent kinase warheads             |
| `common_warhead_covalent_inhibitors`    | covalent-inhibitor warheads          |
| `privileged_kinase_inhibitor_scaffolds` | privileged kinase scaffolds          |
| `privileged_scaffolds`                  | privileged drug scaffolds            |
| `aggregator`                            | aggregation-prone motifs             |

```python
group = mc.groups.ChemicalGroup(groups=["hinge_binders"])
group.has_match(mol)                                       # bool for ONE molecule
keep = mc.functional.chemical_group_filter(mol_list, chemical_group=group)  # batch -> bool array
```

> `phosphate_binders`, `michael_acceptors`, and `reactive_groups` are **not** default catalog
> names in 2.1.0. For Michael acceptors / reactive electrophiles, use the alert filters
> (`alert_filter`, `CommonAlertsFilters`) or a custom SMARTS catalog (below). For covalent
> warheads, use `electrophilic_warheads_for_kinases` / `common_warhead_covalent_inhibitors`.

---

## Custom SMARTS Patterns

`ChemicalGroup` has no `custom_smarts` argument. Build a catalog from your own SMARTS with
`mc.catalogs.catalog_from_smarts`, then match via `mc.functional.catalog_filter`:

```python
cat = mc.catalogs.catalog_from_smarts(
    smarts=["[CX3]=[CX3]C(=O)[NX3]", "[C;H0](=O)C(F)(F)F"],  # acrylamide, CF3-ketone warhead
    labels=["acrylamide", "tfm_ketone"],
)
keep = mc.functional.catalog_filter(mol_list, catalogs=[cat])  # True = no match / keep
```

---

## Filter Selection Guidelines

### Initial Screening (High-Throughput)

Recommended filters:
- Rule of Five
- PAINS filter
- Common Alerts (permissive settings)

```python
rfilter = mc.rules.RuleFilters(rule_list=["rule_of_five"])
alert_filter = mc.structural.CommonAlertsFilters(alerts_set=["PAINS"])  # PAINS via alerts, not a rule
```

---

### Hit-to-Lead

Recommended filters:
- Rule of Oprea or Leadlike (soft)
- NIBR filters
- Lilly Demerits (needs external binaries)

```python
rfilter = mc.rules.RuleFilters(rule_list=["rule_of_oprea"])
nibr = mc.structural.NIBRFilters()
# Lilly via the functional API (requires the tools from `medchem install-lilly`):
# keep_lilly = mc.functional.lilly_demerit_filter(mols, n_jobs=-1)
```

---

### Lead Optimization

Recommended filters:
- Rule of Five + Veber (compose "rule of drug")
- Lead-likeness (`rule_of_oprea`) + tighter property windows
- Full structural alert analysis
- Complexity filter

```python
rfilter = mc.rules.RuleFilters(rule_list=["rule_of_five", "rule_of_veber", "rule_of_oprea"])
alert_filter = mc.structural.CommonAlertsFilters()
complexity_filter = mc.complexity.ComplexityFilter(limit="95", complexity_metric="bertz")
```

---

### CNS Targets

Recommended filters:
- Rule of CNS
- PAINS / alert screening
- BBB-oriented property windows (low TPSA, low HBD)

```python
rfilter = mc.rules.RuleFilters(rule_list=["rule_of_cns"])
# Property windows via the query DSL (no all-in-one Constraints object):
qf = mc.query.QueryFilter('HASPROP("tpsa" <= 90) AND HASPROP("n_lipinski_hbd" <= 2) AND HASPROP("mw" <= 450)')
```

---

### Fragment-Based Drug Discovery

Recommended filters:
- Rule of Three
- Low complexity
- Reactive-group / alert check

```python
rfilter = mc.rules.RuleFilters(rule_list=["rule_of_three"])
complexity_filter = mc.complexity.ComplexityFilter(limit="90", complexity_metric="bertz")
```

---

## Important Considerations

### False Positives and False Negatives

**Filters are guidelines, not absolutes:**

1. **False Positives** (good drugs flagged):
   - ~10% of marketed drugs fail Rule of Five
   - Natural products often violate standard rules
   - Prodrugs intentionally break rules
   - Antibiotics and antivirals frequently non-compliant

2. **False Negatives** (bad compounds passing):
   - Passing filters doesn't guarantee success
   - Target-specific issues not captured
   - In vivo properties not fully predicted

### Context-Specific Application

**Different contexts require different criteria:**

- **Target Class:** Kinases vs GPCRs vs ion channels have different optimal spaces
- **Modality:** Small molecules vs PROTACs vs molecular glues
- **Administration Route:** Oral vs IV vs topical
- **Disease Area:** CNS vs oncology vs infectious disease
- **Stage:** Screening vs hit-to-lead vs lead optimization

### Complementing with Machine Learning

Modern approaches combine rules with ML:

```python
# Rule-based pre-filtering (RuleFilters returns a DataFrame with a pass_all column)
res = mc.rules.RuleFilters(rule_list=["rule_of_five"])(mols=mols, n_jobs=-1)
filtered_mols = [mol for mol, ok in zip(mols, res["pass_all"]) if ok]

# ML model scoring on the filtered set
ml_scores = ml_model.predict(filtered_mols)

# Combined decision
final_candidates = [
    mol for mol, score in zip(filtered_mols, ml_scores)
    if score > threshold
]
```

---

## References

1. Lipinski CA et al. Adv Drug Deliv Rev (1997) 23:3-25
2. Veber DF et al. J Med Chem (2002) 45:2615-2623
3. Oprea TI et al. J Chem Inf Comput Sci (2001) 41:1308-1315
4. Congreve M et al. Drug Discov Today (2003) 8:876-877
5. Baell JB & Holloway GA. J Med Chem (2010) 53:2719-2740
6. Johnson TW et al. J Med Chem (2009) 52:5487-5500
7. Walters WP & Murcko MA. Adv Drug Deliv Rev (2002) 54:255-271
8. Hann MM & Oprea TI. Curr Opin Chem Biol (2004) 8:255-263
9. Rishton GM. Drug Discov Today (1997) 2:382-384
10. Bruns RF & Watson IA. J Med Chem (2012) 55:9763-9772 (Lilly MedChem Rules)

#!/usr/bin/env python3
"""
Batch molecular filtering using medchem library.

This script provides a production-ready workflow for filtering compound libraries
using medchem rules, structural alerts, and custom constraints.

Verified against medchem==2.1.0.

Usage:
    python filter_molecules.py input.csv --rules rule_of_five,rule_of_cns --nibr --output filtered.csv
    python filter_molecules.py input.sdf --rules rule_of_oprea --lilly --complexity-metric bertz --output results.csv
    python filter_molecules.py smiles.txt --nibr --pains --n-jobs -1 --output clean.csv
"""

import argparse
import sys
from pathlib import Path
from typing import List, Tuple

try:
    import numpy as np
    import pandas as pd
    import datamol as dm
    import medchem as mc
    from rdkit import Chem
    from tqdm import tqdm
except ImportError as e:
    print(f"Error: Missing required package: {e}")
    print("Install dependencies: uv pip install medchem datamol pandas tqdm")
    sys.exit(1)


def load_molecules(input_file: Path, smiles_column: str = "smiles") -> Tuple[pd.DataFrame, List[Chem.Mol]]:
    """
    Load molecules from various file formats.

    Supports:
    - CSV/TSV with SMILES column
    - SDF files
    - Plain text files with one SMILES per line

    Returns:
        Tuple of (DataFrame with metadata, list of RDKit molecules)
    """
    suffix = input_file.suffix.lower()

    if suffix == ".sdf":
        print(f"Loading SDF file: {input_file}")
        supplier = Chem.SDMolSupplier(str(input_file))
        mols = [mol for mol in supplier if mol is not None]

        # Create DataFrame from SDF properties
        data = []
        for mol in mols:
            props = mol.GetPropsAsDict()
            props["smiles"] = Chem.MolToSmiles(mol)
            data.append(props)
        df = pd.DataFrame(data)

    elif suffix in [".csv", ".tsv"]:
        print(f"Loading CSV/TSV file: {input_file}")
        sep = "\t" if suffix == ".tsv" else ","
        df = pd.read_csv(input_file, sep=sep)

        if smiles_column not in df.columns:
            print(f"Error: Column '{smiles_column}' not found in file")
            print(f"Available columns: {', '.join(df.columns)}")
            sys.exit(1)

        print(f"Converting SMILES to molecules...")
        mols = [dm.to_mol(smi) for smi in tqdm(df[smiles_column], desc="Parsing")]

    elif suffix == ".txt":
        print(f"Loading text file: {input_file}")
        with open(input_file) as f:
            smiles_list = [line.strip() for line in f if line.strip()]

        df = pd.DataFrame({"smiles": smiles_list})
        print(f"Converting SMILES to molecules...")
        mols = [dm.to_mol(smi) for smi in tqdm(smiles_list, desc="Parsing")]

    else:
        print(f"Error: Unsupported file format: {suffix}")
        print("Supported formats: .csv, .tsv, .sdf, .txt")
        sys.exit(1)

    # Filter out invalid molecules
    valid_indices = [i for i, mol in enumerate(mols) if mol is not None]
    if len(valid_indices) < len(mols):
        n_invalid = len(mols) - len(valid_indices)
        print(f"Warning: {n_invalid} invalid molecules removed")
        df = df.iloc[valid_indices].reset_index(drop=True)
        mols = [mols[i] for i in valid_indices]

    print(f"Loaded {len(mols)} valid molecules")
    return df, mols


def apply_rule_filters(mols: List[Chem.Mol], rules: List[str], n_jobs: int) -> pd.DataFrame:
    """Apply medicinal chemistry rule filters.

    RuleFilters(...)(...) returns a DataFrame with a `mol` column, `pass_all`,
    `pass_any`, and one boolean column per rule. We keep the per-rule booleans
    plus `pass_all` (renamed for clarity).
    """
    print(f"\nApplying rule filters: {', '.join(rules)}")

    rfilter = mc.rules.RuleFilters(rule_list=rules)
    res = rfilter(mols=mols, n_jobs=n_jobs, progress=True)

    df_results = res[rules].copy().reset_index(drop=True)
    df_results["passes_all_rules"] = res["pass_all"].to_numpy()
    return df_results


def apply_structural_alerts(mols: List[Chem.Mol], alert_type: str, n_jobs: int) -> pd.DataFrame:
    """Apply structural alert filters."""
    print(f"\nApplying {alert_type} structural alerts...")

    if alert_type == "common":
        res = mc.structural.CommonAlertsFilters()(mols=mols, n_jobs=n_jobs, progress=True)
        # `pass_filter` is True when the molecule is clean (no alert triggered).
        df_results = pd.DataFrame({
            "passes_common_alerts": res["pass_filter"].to_numpy(),
            "common_alert_details": res["reasons"].fillna("").to_numpy(),
        })

    elif alert_type == "nibr":
        res = mc.structural.NIBRFilters()(mols=mols, n_jobs=n_jobs, progress=True)
        df_results = pd.DataFrame({
            "passes_nibr": res["pass_filter"].to_numpy(),
            "nibr_severity": res["severity"].to_numpy(),
        })

    elif alert_type == "lilly":
        # Requires the Lilly MedChem Rules tools (medchem >= 2.1: `medchem install-lilly`).
        try:
            keep = mc.functional.lilly_demerit_filter(mols, n_jobs=n_jobs, progress=True)
        except ImportError:
            print("Error: Lilly demerits need the Lilly MedChem Rules tools. "
                  "Install them with `medchem install-lilly` (needs make, a C++ compiler "
                  "and zlib), then re-run with --lilly.")
            sys.exit(1)
        df_results = pd.DataFrame({"passes_lilly": np.asarray(keep, dtype=bool)})

    elif alert_type == "pains":
        keep = mc.functional.alert_filter(mols, alerts=["pains"], n_jobs=n_jobs, progress=True)
        df_results = pd.DataFrame({"passes_pains": np.asarray(keep, dtype=bool)})

    else:
        raise ValueError(f"Unknown alert type: {alert_type}")

    return df_results


def apply_complexity_filter(mols: List[Chem.Mol], method: str, limit: str, n_jobs: int) -> pd.DataFrame:
    """Flag molecules exceeding the complexity percentile threshold.

    `complexity_filter` returns a boolean array (True = within the limit / keep).
    """
    print(f"\nApplying complexity filter (metric={method}, limit={limit} percentile)...")
    keep = mc.functional.complexity_filter(
        mols, complexity_metric=method, limit=limit, n_jobs=n_jobs, progress=True
    )
    return pd.DataFrame({"passes_complexity": np.asarray(keep, dtype=bool)})


def apply_chemical_groups(mols: List[Chem.Mol], groups: List[str]) -> pd.DataFrame:
    """Detect chemical groups (one boolean column per group)."""
    print(f"\nDetecting chemical groups: {', '.join(groups)}")
    df_results = pd.DataFrame()
    for group in groups:
        cg = mc.groups.ChemicalGroup(groups=[group])
        matched = mc.functional.chemical_group_filter(mols, chemical_group=cg)
        df_results[f"has_{group}"] = np.asarray(matched, dtype=bool)
    return df_results


def apply_property_windows(mols: List[Chem.Mol], args) -> pd.DataFrame:
    """Enforce physchem property windows via RDKit descriptors + medchem.rules.in_range.

    medchem has no all-in-one property-window object, so we compute the standard
    descriptors with RDKit and gate each with `mc.rules.in_range`.
    """
    from rdkit.Chem import Descriptors, Crippen
    from rdkit.Chem import rdMolDescriptors as rd

    print("\nApplying property windows...")
    in_range = mc.rules.in_range
    passes = []
    for mol in tqdm(mols, desc="Properties"):
        ok = True
        if args.mw_range:
            lo, hi = map(float, args.mw_range.split(","))
            ok &= in_range(Descriptors.MolWt(mol), min_val=lo, max_val=hi)
        if args.logp_range:
            lo, hi = map(float, args.logp_range.split(","))
            ok &= in_range(Crippen.MolLogP(mol), min_val=lo, max_val=hi)
        if args.tpsa_max is not None:
            ok &= in_range(rd.CalcTPSA(mol), max_val=args.tpsa_max)
        if args.hbd_max is not None:
            ok &= in_range(rd.CalcNumLipinskiHBD(mol), max_val=args.hbd_max)
        if args.hba_max is not None:
            ok &= in_range(rd.CalcNumLipinskiHBA(mol), max_val=args.hba_max)
        if args.rotatable_bonds_max is not None:
            ok &= in_range(rd.CalcNumRotatableBonds(mol), max_val=args.rotatable_bonds_max)
        passes.append(bool(ok))
    return pd.DataFrame({"passes_properties": passes})


def generate_summary(df: pd.DataFrame, output_file: Path):
    """Generate filtering summary report."""
    summary_file = output_file.parent / f"{output_file.stem}_summary.txt"

    with open(summary_file, "w") as f:
        f.write("=" * 80 + "\n")
        f.write("MEDCHEM FILTERING SUMMARY\n")
        f.write("=" * 80 + "\n\n")

        f.write(f"Total molecules processed: {len(df)}\n\n")

        # Rule results
        rule_cols = [col for col in df.columns if col.startswith("rule_") or col == "passes_all_rules"]
        if rule_cols:
            f.write("RULE FILTERS:\n")
            f.write("-" * 40 + "\n")
            for col in rule_cols:
                if col in df.columns and df[col].dtype == bool:
                    n_pass = df[col].sum()
                    pct = 100 * n_pass / len(df)
                    f.write(f"  {col}: {n_pass} passed ({pct:.1f}%)\n")
            f.write("\n")

        # Structural alerts
        alert_cols = [col for col in df.columns if "alert" in col.lower() or "nibr" in col.lower() or "lilly" in col.lower() or "pains" in col.lower()]
        if alert_cols:
            f.write("STRUCTURAL ALERTS:\n")
            f.write("-" * 40 + "\n")
            if "passes_common_alerts" in df.columns:
                n_clean = df["passes_common_alerts"].sum()
                pct = 100 * n_clean / len(df)
                f.write(f"  No common alerts: {n_clean} ({pct:.1f}%)\n")
            if "passes_nibr" in df.columns:
                n_pass = df["passes_nibr"].sum()
                pct = 100 * n_pass / len(df)
                f.write(f"  Passes NIBR: {n_pass} ({pct:.1f}%)\n")
            if "passes_lilly" in df.columns:
                n_pass = df["passes_lilly"].sum()
                pct = 100 * n_pass / len(df)
                f.write(f"  Passes Lilly: {n_pass} ({pct:.1f}%)\n")
            if "passes_pains" in df.columns:
                n_pass = df["passes_pains"].sum()
                pct = 100 * n_pass / len(df)
                f.write(f"  Passes PAINS: {n_pass} ({pct:.1f}%)\n")
            f.write("\n")

        # Complexity
        if "passes_complexity" in df.columns:
            f.write("COMPLEXITY:\n")
            f.write("-" * 40 + "\n")
            n_pass = df["passes_complexity"].sum()
            pct = 100 * n_pass / len(df)
            f.write(f"  Within complexity limit: {n_pass} ({pct:.1f}%)\n")
            f.write("\n")

        # Overall pass rate
        pass_cols = [col for col in df.columns if col.startswith("passes_")]
        if pass_cols:
            df["passes_all_filters"] = df[pass_cols].all(axis=1)
            n_pass = df["passes_all_filters"].sum()
            pct = 100 * n_pass / len(df)
            f.write("OVERALL:\n")
            f.write("-" * 40 + "\n")
            f.write(f"  Molecules passing all filters: {n_pass} ({pct:.1f}%)\n")

        f.write("\n" + "=" * 80 + "\n")

    print(f"\nSummary report saved to: {summary_file}")


def main():
    parser = argparse.ArgumentParser(
        description="Batch molecular filtering using medchem",
        formatter_class=argparse.RawDescriptionHelpFormatter,
        epilog=__doc__
    )

    # Input/Output
    parser.add_argument("input", type=Path, help="Input file (CSV, TSV, SDF, or TXT)")
    parser.add_argument("--output", "-o", type=Path, required=True, help="Output CSV file")
    parser.add_argument("--smiles-column", default="smiles", help="Name of SMILES column (default: smiles)")

    # Rule filters
    parser.add_argument("--rules", help="Comma-separated list of rules (e.g., rule_of_five,rule_of_cns)")

    # Structural alerts
    parser.add_argument("--common-alerts", action="store_true", help="Apply common structural alerts")
    parser.add_argument("--nibr", action="store_true", help="Apply NIBR filters")
    parser.add_argument("--lilly", action="store_true", help="Apply Lilly demerits filter")
    parser.add_argument("--pains", action="store_true", help="Apply PAINS filter")

    # Complexity (percentile threshold against a reference set)
    parser.add_argument("--complexity", action="store_true",
                        help="Apply the complexity filter (percentile threshold)")
    parser.add_argument("--complexity-method", default="bertz",
                        choices=["bertz", "whitlock", "barone", "smcm", "twc"],
                        help="Complexity metric (default: bertz)")
    parser.add_argument("--complexity-limit", default="99",
                        help="Percentile limit against the reference set (default: 99)")

    # Property windows (RDKit descriptors + medchem.rules.in_range)
    parser.add_argument("--mw-range", help="Molecular weight range (e.g., 200,500)")
    parser.add_argument("--logp-range", help="cLogP range (e.g., -2,5)")
    parser.add_argument("--tpsa-max", type=float, help="Maximum TPSA")
    parser.add_argument("--hbd-max", type=int, help="Maximum Lipinski H-bond donors")
    parser.add_argument("--hba-max", type=int, help="Maximum Lipinski H-bond acceptors")
    parser.add_argument("--rotatable-bonds-max", type=int, help="Maximum rotatable bonds")

    # Chemical groups
    parser.add_argument("--groups", help="Comma-separated chemical groups to detect")

    # Processing options
    parser.add_argument("--n-jobs", type=int, default=-1, help="Number of parallel jobs (-1 = all cores)")
    parser.add_argument("--no-summary", action="store_true", help="Don't generate summary report")
    parser.add_argument("--filter-output", action="store_true", help="Only output molecules passing all filters")

    args = parser.parse_args()

    # Load molecules
    df, mols = load_molecules(args.input, args.smiles_column)

    # Apply filters
    result_dfs = [df]

    # Rules
    if args.rules:
        rule_list = [r.strip() for r in args.rules.split(",")]
        df_rules = apply_rule_filters(mols, rule_list, args.n_jobs)
        result_dfs.append(df_rules)

    # Structural alerts
    if args.common_alerts:
        df_alerts = apply_structural_alerts(mols, "common", args.n_jobs)
        result_dfs.append(df_alerts)

    if args.nibr:
        df_nibr = apply_structural_alerts(mols, "nibr", args.n_jobs)
        result_dfs.append(df_nibr)

    if args.lilly:
        df_lilly = apply_structural_alerts(mols, "lilly", args.n_jobs)
        result_dfs.append(df_lilly)

    if args.pains:
        df_pains = apply_structural_alerts(mols, "pains", args.n_jobs)
        result_dfs.append(df_pains)

    # Complexity
    if args.complexity:
        df_complexity = apply_complexity_filter(
            mols, args.complexity_method, args.complexity_limit, args.n_jobs
        )
        result_dfs.append(df_complexity)

    # Property windows
    if any([args.mw_range, args.logp_range, args.tpsa_max,
            args.hbd_max, args.hba_max, args.rotatable_bonds_max]):
        df_props = apply_property_windows(mols, args)
        result_dfs.append(df_props)

    # Chemical groups
    if args.groups:
        group_list = [g.strip() for g in args.groups.split(",")]
        df_groups = apply_chemical_groups(mols, group_list)
        result_dfs.append(df_groups)

    # Combine results
    df_final = pd.concat(result_dfs, axis=1)

    # Filter output if requested
    if args.filter_output:
        pass_cols = [col for col in df_final.columns if col.startswith("passes_")]
        if pass_cols:
            df_final["passes_all"] = df_final[pass_cols].all(axis=1)
            df_final = df_final[df_final["passes_all"]]
            print(f"\nFiltered to {len(df_final)} molecules passing all filters")

    # Save results
    args.output.parent.mkdir(parents=True, exist_ok=True)
    df_final.to_csv(args.output, index=False)
    print(f"\nResults saved to: {args.output}")

    # Generate summary
    if not args.no_summary:
        generate_summary(df_final, args.output)

    print("\nDone!")


if __name__ == "__main__":
    main()