Inspiration
Traditional medicine encodes centuries of careful pharmacological observation, but it indexes that knowledge by plant, not by molecule or protein target. Drug discovery works the other way round. The two bodies of knowledge do not share an index, so connecting them is mostly done by hand, one ethnobotanical lead at a time.
The bridge already exists and is barely used: LOTUS, via Wikidata, holds hundreds of thousands of referenced organism–compound occurrence statements. Joined to a target model, a plant becomes a set of testable molecular hypotheses.
But the thing that actually motivated the build was a second observation. Virtual screening results are extremely easy to make look excellent and genuinely hard to make trustworthy. A random train/test split, a missing leakage control, and an enrichment factor quoted on a publication-biased dataset will together produce a spectacular number that means nothing. We wanted to build the version that reports what it actually knows — including where it fails.
What it does
PhytoTarget runs a complete pipeline and ships an interactive dashboard:
- Curates bioactivity for a disease target from ChEMBL — median pChEMBL, not maximum, and the ambiguous 10 µM–316 nM band is dropped rather than relabelled.
- Builds a natural-product library live from LOTUS/Wikidata by SPARQL: 102,470 referenced occurrence pairs, 21,377 source taxa, collapsing to 23,765 unique standardised structures.
- Trains and audits a model — scaffold split as the headline, random split reported beside it, y-scrambling as a leakage control, AVE bias, isotonic calibration.
- Screens the whole library, attaching an applicability-domain flag, a PAINS count and a known-active marker to every single row.
- Answers the inverse question: which source organisms are statistically over-represented among the top novel predictions? Hypergeometric test per genus and species, Benjamini–Hochberg corrected across 12,203 taxa tested.
- Builds a knowledge graph — organism → compound → target → disease — where every edge carries its evidence, exportable to Cytoscape/Gephi.
Step 5 is the part we think is new. A ranked list of 500 molecules is hard to act on. "These 62 genera contain significantly more predicted actives than chance would give you" is a field collection plan.
How we built it
Python, RDKit, LightGBM, scikit-learn, Streamlit, Plotly, NetworkX, ReportLab.
- Features: ECFP4 (Morgan r=2, 2048 bits) + 12 interpretable descriptors.
- One standardisation path for both sources — largest organic fragment, neutralised, canonical InChIKey — so a compound arriving from ChEMBL and from Wikidata collapses to one row instead of being counted twice.
- Model: LightGBM with isotonic calibration, because the output is a ranked list a human will threshold; a compound shown at P = 0.9 should be active ~90% of the time.
- Enrichment: hypergeometric over-representation with BH-FDR, implemented inline.
- Reproducibility:
python run_all.py, config-driven targets, 13 tests covering the invariants the claims rest on, and a pitch deck whose every number is read from pipeline output at build time so it cannot drift from the code.
Results
| Acetylcholinesterase | Monoamine oxidase B | |
|---|---|---|
| ROC-AUC, scaffold split (headline) | 0.943 | 0.922 |
| ROC-AUC, random split (optimistic) | 0.974 | 0.973 |
| AUC inflation from random splitting | +0.031 | +0.050 |
| y-scrambled ROC-AUC (must be ≈0.50) | 0.501 ± 0.022 | 0.505 ± 0.022 |
| AVE split bias | +0.502 | +0.376 |
| Brier / ECE | 0.109 / 0.073 | 0.130 / 0.112 |
| EF@1%, decoy-augmented | 54.9 (ceiling 55) | 66.7 (ceiling 94) |
| Inside applicability domain | 15.7% | 7.8% |
| Enriched taxa at FDR < 0.05 | 62 of 12,203 | 46 of 12,203 |
Retrospective recovery, with the model never told which rows these were:
| Compound | Source organism | Rank of 23,765 | Percentile |
|---|---|---|---|
| galantamine (licensed AChE drug) | Galanthus nivalis | 626 | top 2.6% |
| harmine | Peganum harmala | 927 | top 3.9% |
| huperzine A | Huperzia serrata | 3,730 | top 15.7% |
| harmaline | Peganum harmala | 5,358 | top 22.6% |
| physostigmine | Physostigma venenosum | 8,287 | top 34.9% |
Top novel in-domain hits (AChE): terreulactone D (P = 0.843) and territrem E (P = 0.818), both Aspergillus terreus meroterpenoids. The territrem family is independently documented as exceptionally potent AChE inhibitors, and they were not among the model's training actives — a checkable result rather than a claim.
Challenges we ran into
Four were substantive, and all four changed the result rather than just costing time.
1. Cross-source compound matching silently failed. Known-active flagging matched
0 of 2,216 ChEMBL actives against the library. The cause took a while to see:
Wikidata's P233 SMILES are frequently drawn flat — galantamine's record has no
stereocentres at all — while ChEMBL stores stereo-specified structures, so their full
InChIKeys never agree (...-UHFFFAOYSA-N vs ...-VTHUDJRQSA-N). Matching now uses the
InChIKey skeleton (connectivity layer). That deliberately merges stereoisomers, and
the looser match is the correct one for a flag whose only job is to stop a rediscovery
being presented as a novel hit.
2. Enrichment factor was mathematically incapable of being informative. ChEMBL is publication-biased — people publish what works — so our AChE set is 64% actives. EF@1% is capped at 1/active_fraction, so it cannot exceed 1.56 no matter how good the model is. Our first fix made it worse: downsampling actives to 1% prevalence left two actives in the test set and a standard deviation larger than the mean. The real fix was to invert it — keep every held-out active and hide them among the real 23,765-compound library, which puts prevalence at 1–2% with the active count still in the hundreds.
3. WDQS kept timing out. The transitive plant filter (wdt:P171* wd:Q25314) 504s,
as does resolving compound names through a SPARQL label scan. We pulled occurrence
pairs flat and resolved names through the much cheaper entity-search REST endpoint
instead. One consequence we could not work around is listed under "what we'd fix".
4. A full disk silently killed training, twice. Runs died mid-fit with exit code
0 and no traceback, then with OpenBLAS error: Memory allocation still failed. The
root cause was the C: drive at ~10 MB free, so Windows could not grow the pagefile,
while CalibratedClassifierCV(cv=5) over LightGBM with n_jobs=-1 started five fits
each spawning one thread per core. Capping the thread pools to 1 fixed it. A related
lesson: piping to tail masks the real exit code.
We also found statsmodels.stats missing from the environment — a truncated install,
almost certainly from the same full disk. Rather than depend on a package for twelve
lines of sorting, we implemented Benjamini–Hochberg inline, which removed the
dependency entirely.
Accomplishments that we're proud of
- The y-scrambling control. Ten refits on permuted labels collapse to 0.501 ± 0.022 — tight on chance. That single number is what makes the 0.943 worth reading, and it is the control most screening write-ups skip.
- Reporting the inflation rather than hiding it. We compute the random-split number specifically so we can show it is +0.031 higher than the honest one.
- Surfacing territrems. The model independently ranked a documented potent AChE-inhibitor class at the very top without them being training actives.
- Publishing the numbers that argue against us (below) instead of tuning them away.
What we learned
The most useful thing was how much work it takes to stop a model flattering itself, and how many of those failures are invisible rather than loud. Three of our four major bugs produced plausible output: a known-active count of zero, an enrichment factor of 1.56, and a process that exited 0. None of them raised an error. The only reason we caught them was building the validation harness before trusting the results, and treating a suspiciously clean number as a bug report.
We also learned that the honest version of a result is usually the more interesting one. "Galantamine ranks in the top 2.6%" is a fine headline. "Galantamine ranks in the top 2.6% but physostigmine ranks at 35%, because physostigmine inhibits by covalently carbamylating the catalytic serine and a 2D fingerprint has no representation of a covalent mechanism" is a far better one — it tells you exactly which class of lead this tool will miss.
What's next
What we would fix first. The library is not plants-only. We intended a plant filter, but the transitive query times WDQS out, so the library is all LOTUS source organisms — fungi, algae and cyanobacteria included. This is visible in the output (Lyngbya majuscula is a cyanobacterium; Cladosporium is a fungus) and we documented it rather than quietly rebranding. It cuts both ways: it also happens to be how the strongest AChE hits were found. A pre-resolved taxonomic kingdom table would fix it without touching WDQS at query time.
Beyond that:
- Reduce the applicability-domain problem, which is the real ceiling here: only 15.7% of the library is inside the domain, because natural-product scaffolds are mostly unlike ChEMBL's synthetic chemistry. Conformal prediction or an NP-pretrained representation would help more than a bigger model would.
- Add covalent and mechanism-aware features, so carbamates like physostigmine stop being systematic misses.
- Structure-based rescoring of the top in-domain hits for orthogonal evidence.
- More targets — the config makes adding one a six-line block, no code change.
- Correct the AVE bias honestly rather than by re-splitting, e.g. by reporting performance stratified by nearest-neighbour distance.
Built with
python · rdkit · lightgbm · scikit-learn · pandas · numpy · scipy ·
streamlit · plotly · networkx · reportlab · chembl-api · wikidata ·
sparql · lotus · cheminformatics · machine-learning · drug-discovery
Try it out
pip install -r requirements.txt
python run_all.py # or --skip-fetch to reuse committed data
streamlit run app/streamlit_app.py
Committed outputs/ means the ranked hits, enrichment tables and knowledge graphs can
be inspected without rerunning the ~25-minute pipeline.
Ethics, licensing and scope
Traditional knowledge referenced here originates with the communities who developed and maintained it over generations; this tool points back toward that knowledge rather than claiming it. Ethnobotanical notes shown in the dashboard are literature context for human interpretation and are never model input, so they cannot leak into the statistics they sit beside.
- ChEMBL (EMBL-EBI) — bioactivity, CC BY-SA 3.0
- LOTUS via Wikidata — occurrences, CC0
- RDKit — BSD-3
- Code: MIT
A predicted binder is not a drug. Nothing here addresses absorption, metabolism, toxicity, tissue exposure, or whether engaging the target helps a patient. This is a research prototype and not clinical advice.
