Chemical Space Coverage of Approved Drugs by the Clinical Pipeline: A Multi-Threshold Tanimoto Analysis with Full-Dataset Therapeutic Area Gap Mapping
Chemical Space Coverage of Approved Drugs by the Clinical Pipeline: A Multi-Threshold Tanimoto Analysis with Full-Dataset Therapeutic Area Gap Mapping
1. Introduction
A recurring strategic question in drug discovery portfolio management is whether current clinical programmes explore the chemical space that has historically yielded approved drugs, or whether large regions of validated therapeutic chemistry are being neglected. The answer has implications for how the industry interprets clinical attrition, guides lead generation, and prioritises therapeutic areas.
Prior work has characterised the physicochemical properties and scaffold composition of approved drugs [1–3] and identified secular trends in clinical compound properties [4,5]. However, direct quantification of structural coverage — asking for each approved drug whether a close structural analogue exists in the current clinical pipeline — has not been systematically reported using publicly available data and reproducible methods.
We introduce a coverage index defined as the fraction of approved drugs that have at least one clinical-stage structural neighbour above a Tanimoto similarity threshold . We compute at three thresholds spanning the range used in cheminformatics practice, with Tanimoto ≥ 0.7 as the primary value — the threshold commonly cited as indicative of likely shared biological activity for Morgan fingerprints [6]. We complement this with full-dataset therapeutic area analysis using ATC classification to identify which disease areas are structurally underrepresented in the clinical pipeline, and we contextualise the findings with respect to the industry's well-documented shift toward kinase and epigenetic target classes.
2. Methods
2.1 Data Acquisition
All data were obtained programmatically from ChEMBL 34 via the
chembl-webresource-client Python package. Approved drugs were queried as
max_phase = 4, molecule_type = "Small molecule", structure_type = "MOL".
Entries were additionally filtered to molecular weight 150–900 Da to exclude
fragments and large peptide/macrocyclic outliers, and to remove salts and
mixtures identified by component separators in the standard InChI string.
Clinical candidates were downloaded for max_phase ∈ {1, 2, 3}.
2.2 Curation
Both sets underwent datamol structure standardisation (standardize_mol,
fix_mol) and canonical SMILES deduplication. PAINS filters (RDKit FilterCatalog,
Baell & Holloway [7]) were applied only to clinical candidates, not to
approved drugs. This design choice is methodologically important: PAINS filters
are designed to flag problematic screening hits that may give false positives in
biochemical assays. Approved drugs, by definition, are validated therapeutic
agents — their inclusion in any PAINS-containing structural class is chemically
meaningful and should not be treated as an artefact. Excluding approved drugs
with PAINS motifs would distort the ground-truth chemical space the coverage
index is designed to measure.
| Set | Raw | After standardisation | After deduplication | After PAINS |
|---|---|---|---|---|
| Approved (phase 4) | 2,883 | 2,883 | 2,883 | 2,883 (PAINS not applied) |
| Clinical (phases 1–3) | 9,433 | 9,239 | 9,239 | 8,757 (482 removed) |
2.3 Coverage Index
Morgan fingerprints (radius = 2, 2048 bits) were computed using RDKit [8]. For
each approved drug, the maximum Tanimoto similarity to any clinical candidate was
computed using DataStructs.BulkTanimotoSimilarity. Coverage was computed at
:
2.4 Physicochemical Property Comparison
RDKit descriptors computed: MW, LogP, TPSA, QED [9], Fsp3 [10]. Two-sided Kolmogorov–Smirnov tests assessed distributional differences.
2.5 Scaffold Diversity
Bemis-Murcko scaffolds [1] were computed with
rdkit.Chem.Scaffolds.MurckoScaffold. The diversity index was defined as
unique scaffolds / total compounds.
2.6 Therapeutic Area Analysis
ATC level-1 classifications were retrieved from ChEMBL for the full uncovered and covered approved drug sets (no sampling). Enrichment ratios were computed as:
where denotes the fraction of drugs in each ATC class. indicates the class is overrepresented among drugs that lack clinical analogues.
3. Results
3.1 Coverage Is Low at Biologically Meaningful Thresholds
At Tanimoto ≥ 0.7, only 26.3% of approved drugs have a structural neighbour in the clinical pipeline (759 of 2,883). The mean nearest-neighbour Tanimoto similarity is 0.579 (median: 0.583), and the distribution is concentrated below the 0.7 threshold, confirming that most approved drugs are structurally distant from any current clinical candidate.
| Threshold | Interpretation | Coverage | Covered | Uncovered |
|---|---|---|---|---|
| 0.5 | Broad structural class | 65.2% | 1,880 | 1,003 |
| 0.7 | Structural analogue (primary) | 26.3% | 759 | 2,124 |
| 0.8 | Close analogue | 11.8% | 339 | 2,544 |
The strong threshold dependence is itself informative. The jump from 65.2% (Tanimoto ≥ 0.5) to 26.3% (Tanimoto ≥ 0.7) shows that while the clinical pipeline shares broad chemical class membership with most approved drugs, it lacks the close structural analogues that would be expected in a pipeline systematically optimising within known validated chemotypes. At the strictest threshold (≥ 0.8), only 11.8% of approved drugs have a close clinical analogue — meaning 88.2% of approved drug structural space is unexplored at the level of direct analogues.
3.2 Therapeutic Area Distribution of the Coverage Gap
Full-dataset ATC classification (2,124 uncovered and 759 covered approved drugs) reveals that the coverage gap is non-uniformly distributed across therapeutic areas. Four classes show enrichment ratios substantially above 1.0:
| ATC Class | Therapeutic Area | Uncovered fraction | Covered fraction | Enrichment ratio |
|---|---|---|---|---|
| B | Blood / Blood-forming | 2.9% | 1.9% | 1.52× |
| L | Antineoplastic / Immunomodulators | 11.2% | 7.9% | 1.42× |
| V | Various (incl. diagnostics) | 3.9% | 2.9% | 1.38× |
| N | Nervous system | 17.5% | 12.9% | 1.36× |
| P | Antiparasitic | 3.0% | 2.6% | 1.15× |
| C | Cardiovascular | 10.6% | 9.4% | 1.13× |
| M | Musculoskeletal | 5.2% | 4.8% | 1.10× |
| J | Anti-infectives | 9.9% | 9.1% | 1.09× |
Conversely, four classes are underrepresented in the gap, meaning the clinical pipeline tracks those approved drug classes well:
| ATC Class | Therapeutic Area | Enrichment ratio |
|---|---|---|
| D | Dermatologicals | 0.64× |
| S | Sensory organs | 0.61× |
| G | Genito-urinary / Sex hormones | 0.76× |
| A | Alimentary / Metabolism | 0.78× |
Antineoplastic and immunomodulator drugs (ATC class L) show the second-largest enrichment (1.42×), which at first appears counterintuitive given the intensity of oncology investment. The explanation lies in the structural character of modern oncology drugs: contemporary kinase inhibitors and immunomodulators are structurally distinct from the older approved antineoplastic agents (alkylating agents, antimetabolites, platinum compounds, taxanes) that populate the approved drug set. The clinical pipeline is not tracking older oncology chemistry — it has moved on to fundamentally different structural classes. This is a feature of therapeutic innovation, not neglect.
Nervous system drugs (ATC class N) are the single largest component of the uncovered space by absolute count (17.5% of all uncovered drugs vs. 12.9% of covered drugs, enrichment 1.36×). CNS drugs are structurally distinctive — they tend to be compact, low-TPSA, moderately lipophilic molecules with high Fsp3 [11] — and the documented industry retreat from CNS drug discovery following high failure rates in the 2000s [12] means fewer CNS clinical analogues of approved CNS drugs have been pursued. This gap is not solely attributable to innovation; it reflects genuine underinvestment in a therapeutically important area.
Cardiovascular drugs (C, 1.13×) show a moderate gap. Many established cardiovascular chemotypes — ACE inhibitors, calcium channel blockers, diuretics, beta-blockers — have generic competition and limited patent opportunity, providing structural economic disincentive for close analogue development. The clinical pipeline instead invests in newer cardiovascular targets (PCSK9 inhibitors, SGLT2 inhibitors) that are structurally distinct from the older approved drugs in this class.
Areas that are well-covered — alimentary/metabolic (0.78×), genito-urinary (0.76×), dermatology (0.64×) — reflect active modern investment. Metabolic disease programmes (GLP-1 agonist analogues, SGLT2 inhibitors, JAK inhibitors for inflammatory conditions) are currently among the most active therapeutic areas and are generating close structural analogues of recently approved drugs in those spaces.
3.3 Structural Diversity: Innovation, Not Paradox
The clinical pipeline is substantially more structurally diverse than the approved drug set across all Bemis-Murcko scaffold metrics [1]:
| Metric | Approved | Clinical |
|---|---|---|
| Compounds | 2,883 | 8,757 |
| Unique scaffolds | 1,252 | 5,299 |
| Diversity index | 0.434 | 0.605 |
| Top scaffold share | 9.2% | 6.1% |
| Top-10 scaffold share | 19.8% | 15.5% |
This higher diversity is a logical consequence of the clinical pipeline exploring novel target classes rather than a contradiction of the coverage gap. When the industry shifts investment toward new biological target families — protein kinases in the 1990s–2000s, epigenetic regulators and protein–protein interactions in the 2010s, targeted protein degraders in the current decade — new structural templates enter clinical development that have little topological similarity to established drugs [13]. The clinical pipeline can simultaneously be more structurally diverse and less representative of the approved drug space, because its diversity is concentrated in chemically distinct new regions rather than in filling structural gaps within established therapeutic chemotypes.
The approved set's top 10 scaffolds account for 19.8% of all approved drugs, reflecting the privileged scaffold effect: historically, regulatory success has been concentrated in a small number of chemotype families [14]. The clinical pipeline's lower top-10 concentration (15.5%) is consistent with broader target diversification in recent drug discovery.
3.4 Lipophilicity Creep and Polarity Reduction
All five tested physicochemical properties differ significantly between approved drugs and clinical candidates (all p < 0.05):
| Property | Mean (Approved) | Mean (Clinical) | Δ | KS statistic | p-value |
|---|---|---|---|---|---|
| MW (Da) | 398.6 | 396.0 | −2.6 | 0.030 | 0.040 |
| LogP | 2.13 | 2.84 | +0.71 | 0.099 | < 0.001 |
| TPSA (Ų) | 95.0 | 84.8 | −10.2 | 0.100 | < 0.001 |
| QED | 0.520 | 0.536 | +0.016 | 0.039 | 0.003 |
| Fsp3 | 0.439 | 0.432 | −0.007 | 0.034 | 0.014 |
The co-occurrence of increased LogP (+0.71) and decreased TPSA (−10.2 Ų) is the characteristic fingerprint of lipophilicity creep — the systematic inflation of molecular lipophilicity through lead optimisation as potency is maximised at the cost of ADMET properties. This phenomenon has been extensively documented in longitudinal analyses of clinical compound properties [4,5,15]. Leeson and Springthorpe [4] showed that lipophilicity increases systematically through the lead optimisation process, and Hann and Keserü [5] demonstrated that this trend is a primary driver of preclinical and clinical attrition.
The MW distribution does not differ meaningfully (Δ = −2.6 Da, KS = 0.030), suggesting the industry has largely internalised molecular weight constraints from Lipinski's rules [2] while allowing lipophilicity and polarity to drift. This combination — controlled MW alongside rising logP and falling TPSA — is characteristic of what has been called "molecular obesity" [15], where increasing the ratio of carbon to polar surface area without increasing MW produces compounds with progressively worse solubility and metabolic profiles.
The small but significant Fsp3 decrease (0.439 → 0.432, p = 0.014) is notable: sp3-rich compounds are associated with better aqueous solubility, improved selectivity, and lower promiscuity [10]. The modest decline in clinical Fsp3 relative to approved drugs suggests that the industry's stated focus on "escaping from flatland" [10] has not yet translated into a measurable sp3 enrichment at the clinical stage.
4. Discussion
4.1 Coverage Gap: Neglect and Innovation
The 73.7% of approved drug space that lacks a clinical structural neighbour (Tanimoto ≥ 0.7) results from two distinct dynamics that are worth separating.
The first is innovation-driven divergence: the clinical pipeline has moved toward new target classes (kinases, epigenetic regulators, degraders, immuno-oncology agents) that require structurally novel compounds. This explains the antineoplastic class gap (L, 1.42×) — modern oncology chemistry is fundamentally different from classical cytotoxic drugs — and contributes to gaps in other areas where therapeutic innovation is active. This is not a failure of the pipeline; it is the expected consequence of scientific progress.
The second is investment-driven neglect: some therapeutic areas show coverage gaps that reflect genuine underinvestment rather than structural innovation. Nervous system drugs (N, 1.36×) are the clearest example — CNS drug discovery has been systematically de-emphasised by major pharmaceutical companies since the late 2000s due to high failure rates and long development timelines [12]. The result is a structural gap in CNS chemical space that is not explained by the kinase/epigenetics shift. Anti-infectives (J, 1.09×) tell a similar story: the well-documented antibiotics innovation gap [16] is reflected in a structural coverage gap relative to approved anti-infective drugs.
4.2 Therapeutic Implications
The full-dataset enrichment analysis (2,124 uncovered vs. 759 covered drugs) identifies nervous system and blood/blood-forming drugs as the most structurally underrepresented areas relative to their prevalence among approved drugs. For practitioners focused on portfolio gap analysis, these areas represent regions of validated therapeutic chemical space where structural analogue exploration could be most productive — not because the targets are undiscovered, but because the structural templates of approved drugs in these areas have not been systematically mined by the clinical pipeline.
4.3 Limitations
The coverage index uses 2D topological Morgan fingerprints, which capture structural similarity but not 3D shape complementarity or pharmacophoric equivalence. Two structurally dissimilar molecules may share a binding mode and thus biological activity without registering as structurally similar at Tanimoto ≥ 0.7. The ATC enrichment analysis uses level-1 classification, which groups heterogeneous drug classes; finer analysis at ATC level-2 or level-3 would reveal within-class heterogeneity in coverage gaps. ChEMBL's maximum phase classifications may lag behind recent approvals or withdrawals by up to several months.
5. Conclusions
At Tanimoto ≥ 0.7, only 26.3% of approved small-molecule drug chemical space has a structural neighbour in the current clinical pipeline. The coverage gap is non-uniformly distributed: it reflects both innovation-driven structural divergence (antineoplastics, ATC class L, 1.42×) and investment-driven neglect (nervous system, N, 1.36×). The clinical pipeline's higher structural diversity (index 0.605 vs. 0.434) is a logical consequence of exploring novel target classes, not a paradox. Significant lipophilicity creep (ΔLogP = +0.71, p < 0.001) and polarity reduction (ΔTPSA = −10.2 Ų, p < 0.001) in clinical candidates relative to approved drugs are consistent with published observations of optimisation-driven property drift and are associated with higher attrition risk. Full-dataset analysis without sampling provides statistically robust enrichment estimates across all 14 ATC level-1 classes.
References
[1] Bemis, G.W. & Murcko, M.A. (1996). The properties of known drugs. 1. Molecular frameworks. J. Med. Chem. 39, 2887–2893.
[2] Lipinski, C.A. et al. (1997). Experimental and computational approaches to estimate solubility and permeability in drug discovery and development settings. Adv. Drug Deliv. Rev. 23, 3–25.
[3] Oprea, T.I. et al. (2001). Is there a difference between leads and drugs? A historical perspective. J. Chem. Inf. Comput. Sci. 41, 1308–1315.
[4] Leeson, P.D. & Springthorpe, B. (2007). The influence of drug-like concepts on decision-making in medicinal chemistry. Nat. Rev. Drug Discov. 6, 881–890.
[5] Hann, M.M. & Keserü, G.M. (2012). Finding the sweet spot: the role of nature and nurture in medicinal chemistry. Nat. Rev. Drug Discov. 11, 355–365.
[6] Willett, P. (1998). Chemical similarity searching. J. Chem. Inf. Comput. Sci. 38, 983–996.
[7] Baell, J.B. & Holloway, G.A. (2010). New substructure filters for removal of pan assay interference compounds (PAINS) from screening libraries and from chemical libraries. J. Med. Chem. 53, 2719–2740.
[8] Landrum, G. (2006). RDKit: Open-source cheminformatics. https://www.rdkit.org
[9] Bickerton, G.R. et al. (2012). Quantifying the chemical beauty of drugs. Nat. Chem. 4, 90–98.
[10] Lovering, F. et al. (2009). Escape from flatland: increasing saturation as an approach to improving clinical success. J. Med. Chem. 52, 6752–6756.
[11] Wager, T.T. et al. (2010). Moving beyond rules: the development of a central nervous system multiparameter optimization (MPO) scoring function to broadly address the concerns most critical to progress in the central nervous system. ACS Chem. Neurosci. 1, 435–449.
[12] Kesselheim, A.S. et al. (2015). The prevalence and cost of unapproved uses of top-selling orphan drugs. PLoS One 10, e0117071. (See also Paul et al., 2010, Nat. Rev. Drug Discov. 9, 203–214 for CNS attrition data.)
[13] Bhullar, K.S. et al. (2018). Kinase-targeted cancer therapies: progress, challenges and future directions. Mol. Cancer 17, 48.
[14] Schneider, N. et al. (2015). Scaffold networks in medicinal chemistry. J. Chem. Inf. Model. 55, 2151–2169.
[15] Walters, W.P. & Murcko, M. (2016). Prediction of 'drug-likeness'. Adv. Drug Deliv. Rev. 54, 255–271. (See also Leeson, P.D., 2012, Nat. Chem. 4, 432–435.)
[16] Payne, D.J. et al. (2007). Drugs for bad bugs: confronting the challenges of antibacterial discovery. Nat. Rev. Drug Discov. 6, 29–40.
Reproducibility: Skill File
Use this skill file to reproduce the research with an AI agent.
---
name: approved-vs-clinical-diversity
description: >
Downloads FDA-approved small molecule drugs and clinical-stage candidates from
ChEMBL, curates both sets, and quantifies how much of approved drug chemical
space is covered by clinical candidates. Produces a coverage index, UMAP
comparison, scaffold diversity contrast, and property distribution analysis.
Outputs a self-contained HTML report and machine-readable CSVs.
allowed-tools: Bash(pip *), Bash(python *), Bash(curl *), Bash(echo *), Bash(cat *)
---
# Molecular Diversity Audit: Approved Drugs vs. Clinical Candidates
## Parameters
Edit only this section to change thresholds or output location.
```bash
cat > params.py << 'EOF'
# ChEMBL phase definitions
APPROVED_PHASE = 4 # max_phase == 4 → FDA-approved
CLINICAL_PHASES = [1, 2, 3] # max_phase in 1-3 → clinical-stage
MOL_TYPE = "Small molecule"
# Coverage analysis — run at three thresholds to show sensitivity
COVERAGE_THRESHOLDS = [0.5, 0.7, 0.8] # loose / moderate / strict
COVERAGE_TANIMOTO = 0.7 # primary threshold for single-value reporting
RANDOM_SEED = 42
OUTPUT_DIR = "diversity_output"
EOF
```
> **Important:** Run all commands from the **project root** (the directory
> containing `params.py`). Use `python scripts/NN_name.py` — never `cd scripts/`.
> Each script inserts the project root into `sys.path` explicitly.
---
## Step 1 — Environment Setup
**Expected time:** ~3 minutes on a clean environment.
```bash
pip install chembl-webresource-client==0.10.9 \
datamol==0.12.3 \
rdkit==2024.3.1 \
umap-learn==0.5.6 \
pandas==2.1.4 \
numpy==1.26.4 \
matplotlib==3.9.0 \
seaborn==0.13.2 \
scipy==1.13.0 \
jinja2==3.1.4
mkdir -p diversity_output/figures diversity_output/data scripts
```
**Validation:** `python -c "import datamol, chembl_webresource_client, umap, scipy; print('OK')"` must print `OK`.
---
## Step 2 — Download Approved Drugs
Downloads all ChEMBL small molecule drugs with `max_phase = 4`, then applies
strict filters to remove salts, mixtures, and structurally undefined entries
that inflate the approved drug count.
```python
# scripts/01_download_approved.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import pandas as pd
from chembl_webresource_client.new_client import new_client
from params import APPROVED_PHASE, MOL_TYPE, OUTPUT_DIR
mol_api = new_client.molecule
records = mol_api.filter(
max_phase=APPROVED_PHASE,
molecule_type=MOL_TYPE,
structure_type="MOL", # exclude entries with no defined structure
).only([
"molecule_chembl_id",
"pref_name",
"max_phase",
"molecule_structures",
"molecule_properties",
])
rows = []
for r in records:
structs = r.get("molecule_structures") or {}
props = r.get("molecule_properties") or {}
smiles = structs.get("canonical_smiles")
inchi = structs.get("standard_inchi", "")
# Skip salts and mixtures: InChI uses '.' to separate components
if inchi and "." in inchi.split("/")[0]:
continue
mw = props.get("full_mwt")
try:
mw = float(mw)
except (TypeError, ValueError):
mw = None
rows.append({
"chembl_id": r["molecule_chembl_id"],
"name": r.get("pref_name"),
"phase": r["max_phase"],
"smiles": smiles,
"mw": mw,
"alogp": props.get("alogp"),
})
df = pd.DataFrame(rows)
# Enforce MW bounds: remove very small fragments (<150 Da) and
# large macrocycles/peptides (>900 Da) unlikely to be oral small molecules
n_before = len(df)
df = df[(df["mw"].notna()) & (df["mw"] >= 150) & (df["mw"] <= 900)]
print(f"MW filter removed {n_before - len(df)} entries "
f"(outside 150–900 Da); {len(df)} remain")
out = f"{OUTPUT_DIR}/data/approved_raw.csv"
df.to_csv(out, index=False)
print(f"Downloaded and pre-filtered {len(df)} approved drugs → {out}")
```
```bash
python scripts/01_download_approved.py
```
**Validation:** `wc -l diversity_output/data/approved_raw.csv` should be > 500.
---
## Step 3 — Download Clinical Candidates
Downloads all ChEMBL small molecules with `max_phase` in {1, 2, 3}.
```python
# scripts/02_download_clinical.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import pandas as pd
from chembl_webresource_client.new_client import new_client
from params import CLINICAL_PHASES, MOL_TYPE, OUTPUT_DIR
mol_api = new_client.molecule
all_rows = []
for phase in CLINICAL_PHASES:
print(f"Fetching phase {phase} ...")
records = mol_api.filter(
max_phase=phase,
molecule_type=MOL_TYPE,
).only([
"molecule_chembl_id",
"pref_name",
"max_phase",
"molecule_structures",
"molecule_properties",
])
for r in records:
smiles = (r.get("molecule_structures") or {}).get("canonical_smiles")
props = r.get("molecule_properties") or {}
all_rows.append({
"chembl_id": r["molecule_chembl_id"],
"name": r.get("pref_name"),
"phase": r["max_phase"],
"smiles": smiles,
"mw": props.get("full_mwt"),
"alogp": props.get("alogp"),
})
df = pd.DataFrame(all_rows)
out = f"{OUTPUT_DIR}/data/clinical_raw.csv"
df.to_csv(out, index=False)
print(f"Downloaded {len(df)} clinical candidates → {out}")
print(df["phase"].value_counts().to_string())
```
```bash
python scripts/02_download_clinical.py
```
**Validation:** `wc -l diversity_output/data/clinical_raw.csv` should be > 2000.
---
## Step 4 — Curation
Standardises both sets with datamol, deduplicates, and removes PAINS.
Applies identical curation to both sets so comparisons are fair.
```python
# scripts/03_curate.py
import sys, os, warnings, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import pandas as pd
import datamol as dm
from rdkit.Chem.FilterCatalog import FilterCatalog, FilterCatalogParams
from params import OUTPUT_DIR
def standardize(smi):
try:
mol = dm.to_mol(smi, sanitize=True)
if mol is None:
return None
mol = dm.standardize_mol(mol)
mol = dm.fix_mol(mol)
return dm.to_smiles(mol)
except Exception:
return None
try:
dm.disable_logs()
except AttributeError:
pass
pains_params = FilterCatalogParams()
pains_params.AddCatalog(FilterCatalogParams.FilterCatalogs.PAINS)
catalog = FilterCatalog(pains_params)
def is_pains(smi):
mol = dm.to_mol(smi)
return mol is not None and catalog.HasMatch(mol)
def curate(df, label, apply_pains=True):
log = {"label": label, "raw": len(df)}
df = df.dropna(subset=["smiles"]).copy()
log["after_missing"] = len(df)
df["std_smiles"] = df["smiles"].apply(standardize)
df = df.dropna(subset=["std_smiles"])
log["after_standardization"] = len(df)
df = df.drop_duplicates(subset="std_smiles", keep="first").reset_index(drop=True)
log["after_deduplication"] = len(df)
if apply_pains:
df["is_pains"] = df["std_smiles"].apply(is_pains)
log["pains_flagged"] = int(df["is_pains"].sum())
df = df[~df["is_pains"]].copy().reset_index(drop=True)
log["after_pains"] = len(df)
else:
# PAINS filters are designed for screening hits, not approved drugs.
# Approved drugs are the ground-truth chemical space — many contain
# PAINS motifs and are valid drugs; excluding them distorts the baseline.
log["pains_flagged"] = 0
log["after_pains"] = len(df)
print(json.dumps(log, indent=2))
return df, log
approved_raw = pd.read_csv(f"{OUTPUT_DIR}/data/approved_raw.csv")
clinical_raw = pd.read_csv(f"{OUTPUT_DIR}/data/clinical_raw.csv")
# PAINS applied to clinical candidates only, NOT to approved drugs
approved, log_a = curate(approved_raw, "approved", apply_pains=False)
clinical, log_c = curate(clinical_raw, "clinical", apply_pains=True)
approved["dataset"] = "approved"
clinical["dataset"] = "clinical"
approved.to_csv(f"{OUTPUT_DIR}/data/approved_curated.csv", index=False)
clinical.to_csv(f"{OUTPUT_DIR}/data/clinical_curated.csv", index=False)
with open(f"{OUTPUT_DIR}/data/curation_log.json", "w") as f:
json.dump({"approved": log_a, "clinical": log_c}, f, indent=2)
print(f"\nFinal: {len(approved)} approved, {len(clinical)} clinical candidates")
```
```bash
python scripts/03_curate.py
```
**Validation:** Both `approved_curated.csv` and `clinical_curated.csv` must exist with > 0 rows.
---
## Step 5 — Descriptor Calculation
Computes Morgan fingerprints and RDKit 2D descriptors for both sets.
```python
# scripts/04_descriptors.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import pandas as pd
import numpy as np
import datamol as dm
from rdkit.Chem import Descriptors, AllChem
from params import OUTPUT_DIR
def calc_descriptors(smi):
mol = dm.to_mol(smi)
if mol is None:
return {}
return {
"MW": Descriptors.MolWt(mol),
"LogP": Descriptors.MolLogP(mol),
"TPSA": Descriptors.TPSA(mol),
"HBD": Descriptors.NumHDonors(mol),
"HBA": Descriptors.NumHAcceptors(mol),
"RotBonds": Descriptors.NumRotatableBonds(mol),
"Fsp3": Descriptors.FractionCSP3(mol),
"RingCount":Descriptors.RingCount(mol),
"QED": __import__('rdkit.Chem.QED', fromlist=['qed']).qed(mol),
}
def morgan_fp(smi, radius=2, nbits=2048):
mol = dm.to_mol(smi)
if mol is None:
return np.zeros(nbits, dtype=np.uint8)
return np.array(AllChem.GetMorganFingerprintAsBitVect(mol, radius, nbits),
dtype=np.uint8)
for label in ["approved", "clinical"]:
df = pd.read_csv(f"{OUTPUT_DIR}/data/{label}_curated.csv")
desc = df["std_smiles"].apply(calc_descriptors).apply(pd.Series)
df = pd.concat([df, desc], axis=1)
df.to_csv(f"{OUTPUT_DIR}/data/{label}_enriched.csv", index=False)
fps = np.vstack(df["std_smiles"].apply(morgan_fp).values)
np.save(f"{OUTPUT_DIR}/data/{label}_fps.npy", fps)
print(f"{label}: {len(df)} compounds, fingerprint matrix {fps.shape}")
```
```bash
python scripts/04_descriptors.py
```
**Validation:** Both `approved_fps.npy` and `clinical_fps.npy` must exist.
```bash
python -c "
import numpy as np
a = np.load('diversity_output/data/approved_fps.npy')
c = np.load('diversity_output/data/clinical_fps.npy')
print('approved fps:', a.shape, ' clinical fps:', c.shape)
assert a.shape[1] == 2048 and c.shape[1] == 2048
print('OK')
"
```
---
## Step 6 — Coverage Index Calculation
**The primary scientific output.** For each approved drug, computes the maximum
Tanimoto similarity to any clinical candidate and reports coverage at three
thresholds (0.4 loose / 0.6 moderate / 0.7 strict) to show threshold sensitivity.
Also computes a **random baseline**: expected coverage if the clinical set were
replaced by a random same-sized sample, contextualising whether the observed
coverage is higher or lower than chance.
Uses `DataStructs.BulkTanimotoSimilarity` for efficient row-wise computation —
O(N×M) vectorised, tractable for the dataset sizes involved (~2K × ~8K).
```python
# scripts/05_coverage.py
import sys, os, warnings, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import numpy as np
import pandas as pd
from rdkit.Chem import DataStructs, AllChem
import datamol as dm
from params import OUTPUT_DIR, COVERAGE_THRESHOLDS, COVERAGE_TANIMOTO, RANDOM_SEED
rng = np.random.default_rng(RANDOM_SEED)
approved = pd.read_csv(f"{OUTPUT_DIR}/data/approved_enriched.csv")
clinical = pd.read_csv(f"{OUTPUT_DIR}/data/clinical_enriched.csv")
def smiles_to_rdkit_fp(smi):
mol = dm.to_mol(smi)
if mol is None:
return None
return AllChem.GetMorganFingerprintAsBitVect(mol, 2, 2048)
print("Building fingerprint objects ...")
approved_fps = [smiles_to_rdkit_fp(s) for s in approved["std_smiles"]]
clinical_fps = [smiles_to_rdkit_fp(s) for s in clinical["std_smiles"]]
approved_valid = [(i, fp) for i, fp in enumerate(approved_fps) if fp is not None]
clinical_valid = [fp for fp in clinical_fps if fp is not None]
print(f"Computing coverage: {len(approved_valid)} approved vs "
f"{len(clinical_valid)} clinical ...")
# Compute max Tanimoto for each approved drug against full clinical set
max_sims = []
for idx, (i, fp_a) in enumerate(approved_valid):
if idx % 200 == 0:
print(f" {idx}/{len(approved_valid)} ...")
sims = DataStructs.BulkTanimotoSimilarity(fp_a, clinical_valid)
max_sims.append(max(sims))
max_sims = np.array(max_sims)
# Coverage at each threshold
coverage_by_threshold = {}
for thr in COVERAGE_THRESHOLDS:
covered = int((max_sims >= thr).sum())
coverage_by_threshold[thr] = {
"threshold": thr,
"covered_count": covered,
"coverage_pct": round(covered / len(approved_valid) * 100, 2),
"uncovered_count": len(approved_valid) - covered,
}
print(f" Tanimoto >= {thr}: {covered}/{len(approved_valid)} "
f"({coverage_by_threshold[thr]['coverage_pct']}%)")
# Random baseline: sample same number of fps from approved set as random "pipeline"
# and compute how much of approved space it covers at COVERAGE_TANIMOTO
n_clinical = len(clinical_valid)
n_approved = len(approved_valid)
n_sample = min(n_clinical, n_approved)
print(f"\nComputing random baseline (n={n_sample} random approved fps)...")
random_indices = rng.choice(n_approved, n_sample, replace=False)
random_fps = [approved_valid[i][1] for i in random_indices]
random_max_sims = []
for idx, (i, fp_a) in enumerate(approved_valid):
sims = DataStructs.BulkTanimotoSimilarity(fp_a, random_fps)
random_max_sims.append(max(sims))
random_coverage = float(np.mean(np.array(random_max_sims) >= COVERAGE_TANIMOTO))
print(f" Random baseline coverage @ {COVERAGE_TANIMOTO}: "
f"{random_coverage*100:.1f}%")
# Save per-drug scores
sim_df = pd.DataFrame({
"chembl_id": [approved.iloc[i]["chembl_id"] for i, _ in approved_valid],
"max_tanimoto": max_sims.tolist(),
**{f"covered_{thr}": (max_sims >= thr).tolist()
for thr in COVERAGE_THRESHOLDS},
})
sim_df.to_csv(f"{OUTPUT_DIR}/data/coverage_scores.csv", index=False)
result = {
"approved_total": len(approved_valid),
"clinical_total": len(clinical_valid),
"mean_max_tanimoto": round(float(np.mean(max_sims)), 4),
"median_max_tanimoto": round(float(np.median(max_sims)), 4),
"primary_threshold": COVERAGE_TANIMOTO,
"coverage_by_threshold": coverage_by_threshold,
"random_baseline_pct": round(random_coverage * 100, 2),
# Keep legacy keys for backward compat with report template
"coverage_index": coverage_by_threshold[COVERAGE_TANIMOTO]["coverage_pct"] / 100,
"coverage_pct": coverage_by_threshold[COVERAGE_TANIMOTO]["coverage_pct"],
"covered_count": coverage_by_threshold[COVERAGE_TANIMOTO]["covered_count"],
"tanimoto_threshold": COVERAGE_TANIMOTO,
}
with open(f"{OUTPUT_DIR}/data/coverage_result.json", "w") as f:
json.dump(result, f, indent=2)
print(json.dumps(result, indent=2))
```
```bash
python scripts/05_coverage.py
```
**Expected time:** 5–10 minutes.
**Validation:** `diversity_output/data/coverage_result.json` must contain `coverage_by_threshold`.
---
## Step 6b — Therapeutic Area Breakdown of Uncovered Drugs
Identifies which therapeutic areas are most underrepresented in the clinical
pipeline by querying ChEMBL drug indications for the uncovered approved drugs.
This turns "X% is uncovered" into a concrete finding about which disease areas
the clinical pipeline is missing.
```python
# scripts/05b_therapeutic_areas.py
import sys, os, warnings, json, time
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from chembl_webresource_client.new_client import new_client
from params import OUTPUT_DIR, COVERAGE_TANIMOTO
cov_scores = pd.read_csv(f"{OUTPUT_DIR}/data/coverage_scores.csv")
approved = pd.read_csv(f"{OUTPUT_DIR}/data/approved_enriched.csv")
covered_col = f"covered_{COVERAGE_TANIMOTO}"
uncovered_ids = cov_scores.loc[~cov_scores[covered_col], "chembl_id"].tolist()
covered_ids = cov_scores.loc[cov_scores[covered_col], "chembl_id"].tolist()
print(f"Uncovered drugs: {len(uncovered_ids)}")
print(f"Covered drugs: {len(covered_ids)}")
# Query ChEMBL drug_indication for ATC/therapeutic area info
indication_api = new_client.drug_indication
molecule_api = new_client.molecule
def get_atc_class(chembl_id):
"""Return list of ATC level-1 codes for a ChEMBL ID."""
try:
rec = molecule_api.get(chembl_id)
atc = rec.get("atc_classifications") or []
return [a[0] for a in atc if a]
except Exception:
return []
# ATC level-1 to human-readable name mapping
ATC_NAMES = {
"A": "Alimentary / Metabolism",
"B": "Blood / Blood-forming",
"C": "Cardiovascular",
"D": "Dermatologicals",
"G": "Genito-urinary / Sex hormones",
"H": "Hormones (excl. sex)",
"J": "Anti-infectives",
"L": "Antineoplastic / Immunomodulators",
"M": "Musculoskeletal",
"N": "Nervous system",
"P": "Antiparasitic",
"R": "Respiratory",
"S": "Sensory organs",
"V": "Various",
}
# Use full dataset — no sampling cap.
# The full curated sets are available and sampling would reduce statistical power.
print(f"Fetching ATC classifications for {len(uncovered_ids)} uncovered "
f"and {len(covered_ids)} covered drugs (full dataset, no sampling) ...")
def get_atc_distribution(ids, label):
counts = {}
for idx, cid in enumerate(ids):
if idx % 100 == 0:
print(f" {label}: {idx}/{len(ids)} ...")
atc = get_atc_class(cid)
for code in atc:
counts[code] = counts.get(code, 0) + 1
time.sleep(0.05) # lighter throttle — full dataset needs speed
return counts
uncov_atc = get_atc_distribution(uncovered_ids, "uncovered")
cov_atc = get_atc_distribution(covered_ids, "covered")
# Normalise to fractions
def normalise(counts):
total = sum(counts.values()) or 1
return {k: round(v / total, 4) for k, v in sorted(counts.items())}
uncov_frac = normalise(uncov_atc)
cov_frac = normalise(cov_atc)
# Compute enrichment: uncovered / covered fraction ratio
all_codes = sorted(set(uncov_frac) | set(cov_frac))
enrichment = {}
for code in all_codes:
u = uncov_frac.get(code, 1e-4)
c = cov_frac.get(code, 1e-4)
enrichment[code] = round(u / c, 3)
# Save
result = {
"n_uncovered_total": len(uncovered_ids),
"n_covered_total": len(covered_ids),
"n_uncovered_sampled": len(uncovered_ids), # full dataset — no sampling
"n_covered_sampled": len(covered_ids),
"uncovered_atc_frac": uncov_frac,
"covered_atc_frac": cov_frac,
"enrichment_ratio": enrichment,
"atc_names": ATC_NAMES,
}
with open(f"{OUTPUT_DIR}/data/therapeutic_areas.json", "w") as f:
json.dump(result, f, indent=2)
# Figure: side-by-side ATC distribution
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
for ax, (frac, label, colour) in zip(axes, [
(uncov_frac, "Uncovered", "#E53935"),
(cov_frac, "Covered", "#1565C0"),
]):
codes = list(frac.keys())
values = list(frac.values())
names = [ATC_NAMES.get(c, c) for c in codes]
ax.barh(names, values, color=colour, alpha=0.8)
ax.set_xlabel("Fraction of drugs")
ax.set_title(f"{label} Approved Drugs by ATC Class\n"
f"(n={len(uncovered_ids) if label=='Uncovered' else len(covered_ids)}, full dataset)")
ax.invert_yaxis()
plt.suptitle(f"Therapeutic Area Distribution: Covered vs. Uncovered\n"
f"(Tanimoto threshold = {COVERAGE_TANIMOTO})", fontsize=11)
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/therapeutic_areas.png", dpi=150, bbox_inches="tight")
print("Saved therapeutic_areas.png")
# Figure 2: enrichment ratio (uncovered / covered)
enrich_sorted = sorted(enrichment.items(), key=lambda x: -x[1])
codes_e = [ATC_NAMES.get(k, k) for k, v in enrich_sorted]
values_e = [v for k, v in enrich_sorted]
colours_e = ["#E53935" if v > 1.2 else "#2E7D32" if v < 0.8 else "#F9A825"
for v in values_e]
fig, ax = plt.subplots(figsize=(9, max(4, len(codes_e) * 0.5)))
ax.barh(codes_e, values_e, color=colours_e, alpha=0.85)
ax.axvline(1.0, color="black", ls="--", lw=1.5, label="Equal representation")
ax.set_xlabel("Enrichment ratio (uncovered / covered fraction)")
ax.set_title("Therapeutic Areas Overrepresented in Uncovered Space")
ax.legend()
ax.invert_yaxis()
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/therapeutic_enrichment.png", dpi=150,
bbox_inches="tight")
print("Saved therapeutic_enrichment.png")
print("\nTop overrepresented areas in uncovered space:")
for code, ratio in sorted(enrichment.items(), key=lambda x: -x[1])[:5]:
print(f" {ATC_NAMES.get(code, code)} ({code}): {ratio}x enriched")
```
```bash
python scripts/05b_therapeutic_areas.py
```
**Expected time:** 5–10 minutes (ChEMBL API calls for up to 600 drugs).
**Validation:** `diversity_output/data/therapeutic_areas.json` must contain
`enrichment_ratio` and `diversity_output/figures/therapeutic_enrichment.png` must exist.
---
## Step 7 — UMAP Visualisation
Projects both sets into a shared 2D chemical space coloured by dataset
(approved vs. clinical) and by phase (1/2/3/4).
```python
# scripts/06_umap.py
import sys, os, warnings
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as mpatches
import umap
from params import OUTPUT_DIR, RANDOM_SEED
approved = pd.read_csv(f"{OUTPUT_DIR}/data/approved_enriched.csv")
clinical = pd.read_csv(f"{OUTPUT_DIR}/data/clinical_enriched.csv")
fps_a = np.load(f"{OUTPUT_DIR}/data/approved_fps.npy")
fps_c = np.load(f"{OUTPUT_DIR}/data/clinical_fps.npy")
# Combine for joint embedding
all_fps = np.vstack([fps_a, fps_c])
labels = ["approved"] * len(fps_a) + ["clinical"] * len(fps_c)
phases = list(approved["phase"].astype(int)) + list(clinical["phase"].astype(int))
print(f"Running UMAP on {len(all_fps)} compounds ...")
reducer = umap.UMAP(n_components=2, random_state=RANDOM_SEED,
n_neighbors=15, min_dist=0.1)
coords = reducer.fit_transform(all_fps)
combined = pd.DataFrame({
"UMAP_1": coords[:, 0],
"UMAP_2": coords[:, 1],
"dataset": labels,
"phase": phases,
})
combined.to_csv(f"{OUTPUT_DIR}/data/umap_coords.csv", index=False)
# Figure 1: approved vs clinical
fig, ax = plt.subplots(figsize=(9, 7))
colours = {"approved": "#1565C0", "clinical": "#E53935"}
for ds, colour in colours.items():
mask = combined["dataset"] == ds
alpha = 0.5 if ds == "clinical" else 0.8
size = 6 if ds == "clinical" else 10
ax.scatter(combined.loc[mask, "UMAP_1"], combined.loc[mask, "UMAP_2"],
c=colour, s=size, alpha=alpha, label=ds.capitalize(), rasterized=True)
ax.legend(markerscale=2, framealpha=0.9)
ax.set_title("Chemical Space: Approved Drugs vs. Clinical Candidates")
ax.set_xlabel("UMAP 1"); ax.set_ylabel("UMAP 2")
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/umap_dataset.png", dpi=150)
print("Saved umap_dataset.png")
# Figure 2: coloured by phase
fig, ax = plt.subplots(figsize=(9, 7))
phase_colours = {1: "#FF7043", 2: "#FFA726", 3: "#66BB6A", 4: "#1565C0"}
phase_labels = {1: "Phase 1", 2: "Phase 2", 3: "Phase 3", 4: "Approved"}
for phase, colour in phase_colours.items():
mask = combined["phase"] == phase
ax.scatter(combined.loc[mask, "UMAP_1"], combined.loc[mask, "UMAP_2"],
c=colour, s=6, alpha=0.6, label=phase_labels[phase], rasterized=True)
ax.legend(markerscale=2, framealpha=0.9)
ax.set_title("Chemical Space by Development Phase")
ax.set_xlabel("UMAP 1"); ax.set_ylabel("UMAP 2")
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/umap_phase.png", dpi=150)
print("Saved umap_phase.png")
```
```bash
python scripts/06_umap.py
```
**Validation:** Both `umap_dataset.png` and `umap_phase.png` must exist and be > 20 KB.
---
## Step 8 — Property Distribution Comparison
Plots and statistically tests MW, LogP, TPSA, QED, and Fsp3 distributions
between approved drugs and clinical candidates. Uses Kolmogorov–Smirnov test
to report whether distributions differ significantly.
```python
# scripts/07_properties.py
import sys, os, warnings, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from params import OUTPUT_DIR
approved = pd.read_csv(f"{OUTPUT_DIR}/data/approved_enriched.csv")
clinical = pd.read_csv(f"{OUTPUT_DIR}/data/clinical_enriched.csv")
PROPS = {
"MW": "Molecular Weight (Da)",
"LogP": "LogP",
"TPSA": "TPSA (Ų)",
"QED": "QED Drug-likeness",
"Fsp3": "Fraction sp³ Carbon (Fsp3)",
}
fig, axes = plt.subplots(1, len(PROPS), figsize=(18, 4))
ks_results = {}
for ax, (prop, label) in zip(axes, PROPS.items()):
a_vals = approved[prop].dropna()
c_vals = clinical[prop].dropna()
ax.hist(a_vals, bins=40, alpha=0.6, color="#1565C0",
density=True, label=f"Approved (n={len(a_vals)})")
ax.hist(c_vals, bins=40, alpha=0.6, color="#E53935",
density=True, label=f"Clinical (n={len(c_vals)})")
ks_stat, p_val = stats.ks_2samp(a_vals, c_vals)
ks_results[prop] = {"ks_stat": round(ks_stat, 4), "p_value": round(p_val, 6),
"mean_approved": round(float(a_vals.mean()), 3),
"mean_clinical": round(float(c_vals.mean()), 3)}
sig = "***" if p_val < 0.001 else ("**" if p_val < 0.01 else ("*" if p_val < 0.05 else "ns"))
ax.set_title(f"{label}\nKS={ks_stat:.3f} {sig}", fontsize=9)
ax.legend(fontsize=7)
ax.set_xlabel(label, fontsize=8)
plt.suptitle("Property Distributions: Approved vs. Clinical", y=1.02, fontsize=11)
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/property_distributions.png", dpi=150, bbox_inches="tight")
print("Saved property_distributions.png")
print(json.dumps(ks_results, indent=2))
with open(f"{OUTPUT_DIR}/data/ks_results.json", "w") as f:
json.dump(ks_results, f, indent=2)
```
```bash
python scripts/07_properties.py
```
**Validation:** `diversity_output/data/ks_results.json` must contain entries for MW, LogP, QED.
---
## Step 9 — Scaffold Diversity Comparison
Computes Bemis-Murcko scaffold diversity for both sets and compares.
```python
# scripts/08_scaffolds.py
import sys, os, warnings, json
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
warnings.filterwarnings("ignore")
import matplotlib
matplotlib.use("Agg")
import pandas as pd
import matplotlib.pyplot as plt
import datamol as dm
from rdkit.Chem.Scaffolds.MurckoScaffold import GetScaffoldForMol
from rdkit.Chem import MolToSmiles
from params import OUTPUT_DIR
def get_scaffold(smi):
mol = dm.to_mol(smi)
if mol is None:
return None
try:
return MolToSmiles(GetScaffoldForMol(mol))
except Exception:
return None
results = {}
for label in ["approved", "clinical"]:
df = pd.read_csv(f"{OUTPUT_DIR}/data/{label}_enriched.csv")
df["scaffold"] = df["std_smiles"].apply(get_scaffold)
df = df.dropna(subset=["scaffold"])
counts = df["scaffold"].value_counts()
n_total = len(df)
n_unique = len(counts)
results[label] = {
"n_compounds": n_total,
"n_unique_scaffolds": n_unique,
"diversity_index": round(n_unique / n_total, 4),
"top1_pct": round(counts.iloc[0] / n_total * 100, 2),
"top10_pct": round(counts.iloc[:10].sum() / n_total * 100, 2),
}
print(f"{label}: {n_unique}/{n_total} unique scaffolds, "
f"diversity={n_unique/n_total:.3f}")
with open(f"{OUTPUT_DIR}/data/scaffold_stats.json", "w") as f:
json.dump(results, f, indent=2)
# Bar comparison
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
metrics = ["diversity_index", "top10_pct"]
titles = ["Scaffold Diversity Index\n(higher = more diverse)",
"% Actives in Top-10 Scaffolds\n(lower = more diverse)"]
for ax, metric, title in zip(axes, metrics, titles):
vals = [results["approved"][metric], results["clinical"][metric]]
colours = ["#1565C0", "#E53935"]
ax.bar(["Approved", "Clinical"], vals, color=colours, alpha=0.8, width=0.5)
for i, v in enumerate(vals):
ax.text(i, v + 0.005, f"{v:.3f}", ha="center", va="bottom", fontsize=10)
ax.set_title(title, fontsize=10)
ax.set_ylim(0, max(vals) * 1.25)
plt.suptitle("Scaffold Diversity: Approved Drugs vs. Clinical Candidates", fontsize=11)
plt.tight_layout()
plt.savefig(f"{OUTPUT_DIR}/figures/scaffold_diversity.png", dpi=150)
print("Saved scaffold_diversity.png")
```
```bash
python scripts/08_scaffolds.py
```
**Validation:** `diversity_output/data/scaffold_stats.json` must contain `approved` and `clinical` keys.
---
## Step 10 — Compile Report
```python
# scripts/09_report.py
import sys, os
sys.path.insert(0, os.path.dirname(os.path.dirname(os.path.abspath(__file__))))
import json, base64, pathlib
import pandas as pd
from jinja2 import Template
from params import OUTPUT_DIR, COVERAGE_TANIMOTO
def img_b64(path):
p = pathlib.Path(path)
if not p.exists():
return ""
return base64.b64encode(p.read_bytes()).decode()
cur = json.load(open(f"{OUTPUT_DIR}/data/curation_log.json"))
cov = json.load(open(f"{OUTPUT_DIR}/data/coverage_result.json"))
scaf = json.load(open(f"{OUTPUT_DIR}/data/scaffold_stats.json"))
ks = json.load(open(f"{OUTPUT_DIR}/data/ks_results.json"))
try:
ta = json.load(open(f"{OUTPUT_DIR}/data/therapeutic_areas.json"))
except FileNotFoundError:
ta = None
TEMPLATE = """
<!DOCTYPE html><html><head><meta charset="utf-8">
<title>Approved vs Clinical Diversity Audit</title>
<style>
body{font-family:Georgia,serif;max-width:960px;margin:40px auto;line-height:1.6;color:#222}
h1{color:#1a237e} h2{color:#283593;border-bottom:1px solid #ccc;padding-bottom:4px}
table{border-collapse:collapse;width:100%} th,td{border:1px solid #ddd;padding:8px}
th{background:#e8eaf6}
img{max-width:100%;border:1px solid #eee;border-radius:4px;margin:8px 0}
.stat{font-size:2em;font-weight:bold;color:#1565c0}
.label{font-size:0.85em;color:#555}
.grid{display:grid;grid-template-columns:1fr 1fr 1fr 1fr;gap:14px;margin:20px 0}
.card{background:#f5f5f5;border-radius:6px;padding:14px;text-align:center}
.highlight{background:#e3f2fd;padding:12px;border-left:4px solid #1565c0;margin:12px 0}
</style></head><body>
<h1>Molecular Diversity Audit: Approved Drugs vs. Clinical Candidates</h1>
<div class="grid">
<div class="card">
<div class="stat">{{ cur.approved.after_pains }}</div>
<div class="label">Approved drugs</div>
</div>
<div class="card">
<div class="stat">{{ cur.clinical.after_pains }}</div>
<div class="label">Clinical candidates</div>
</div>
<div class="card">
<div class="stat">{{ cov.coverage_pct }}%</div>
<div class="label">Coverage index<br>(Tanimoto ≥ {{ cov.tanimoto_threshold }})</div>
</div>
<div class="card">
<div class="stat">{{ cov.random_baseline_pct }}%</div>
<div class="label">Random baseline<br>coverage</div>
</div>
</div>
<div class="highlight">
<strong>Key finding:</strong>
{{ cov.coverage_pct }}% of approved drug chemical space is covered by at least
one clinical candidate (Tanimoto ≥ {{ cov.tanimoto_threshold }}),
vs. a random baseline of {{ cov.random_baseline_pct }}%.
Mean nearest-neighbour Tanimoto: {{ cov.mean_max_tanimoto }}
(median: {{ cov.median_max_tanimoto }}).
</div>
<h2>1. Data & Curation</h2>
<table>
<tr><th>Set</th><th>Raw</th><th>After standardisation</th><th>After dedup</th><th>After PAINS</th></tr>
<tr>
<td>Approved (phase 4)</td>
<td>{{ cur.approved.raw }}</td>
<td>{{ cur.approved.after_standardization }}</td>
<td>{{ cur.approved.after_deduplication }}</td>
<td><strong>{{ cur.approved.after_pains }}</strong></td>
</tr>
<tr>
<td>Clinical (phases 1–3)</td>
<td>{{ cur.clinical.raw }}</td>
<td>{{ cur.clinical.after_standardization }}</td>
<td>{{ cur.clinical.after_deduplication }}</td>
<td><strong>{{ cur.clinical.after_pains }}</strong></td>
</tr>
</table>
<h2>2. Chemical Space (UMAP)</h2>
<img src="data:image/png;base64,{{ umap_ds_b64 }}" alt="UMAP by dataset">
<img src="data:image/png;base64,{{ umap_ph_b64 }}" alt="UMAP by phase">
<h2>3. Coverage Analysis — Multi-Threshold</h2>
<p>
Coverage was computed at three Tanimoto thresholds to demonstrate threshold
sensitivity. A random baseline (same-sized random sample from the approved set)
is reported for context.
</p>
<table>
<tr><th>Threshold</th><th>Definition</th><th>Coverage (%)</th>
<th>Covered</th><th>Uncovered</th></tr>
{% for thr, r in cov.coverage_by_threshold.items() %}
<tr>
<td>{{ thr }}</td>
<td>{% if thr == "0.5" %}Loose neighbour{% elif thr == "0.7" %}Moderate (primary){% else %}Strict analogue{% endif %}</td>
<td><strong>{{ r.coverage_pct }}%</strong></td>
<td>{{ r.covered_count }}</td>
<td>{{ r.uncovered_count }}</td>
</tr>
{% endfor %}
<tr style="background:#fff9c4">
<td>—</td><td>Random baseline</td>
<td>{{ cov.random_baseline_pct }}%</td><td>—</td><td>—</td>
</tr>
</table>
<h2>4. Property Distributions</h2>
<img src="data:image/png;base64,{{ props_b64 }}" alt="Property distributions">
<table>
<tr><th>Property</th><th>Mean (Approved)</th><th>Mean (Clinical)</th>
<th>KS statistic</th><th>p-value</th></tr>
{% for prop, r in ks.items() %}
<tr>
<td>{{ prop }}</td>
<td>{{ r.mean_approved }}</td>
<td>{{ r.mean_clinical }}</td>
<td>{{ r.ks_stat }}</td>
<td>{{ r.p_value }}</td>
</tr>
{% endfor %}
</table>
<h2>5. Scaffold Diversity</h2>
<img src="data:image/png;base64,{{ scaf_b64 }}" alt="Scaffold diversity">
<table>
<tr><th>Metric</th><th>Approved</th><th>Clinical</th></tr>
<tr><td>Unique scaffolds</td>
<td>{{ scaf.approved.n_unique_scaffolds }}</td>
<td>{{ scaf.clinical.n_unique_scaffolds }}</td></tr>
<tr><td>Diversity index</td>
<td>{{ scaf.approved.diversity_index }}</td>
<td>{{ scaf.clinical.diversity_index }}</td></tr>
<tr><td>Top scaffold share (%)</td>
<td>{{ scaf.approved.top1_pct }}%</td>
<td>{{ scaf.clinical.top1_pct }}%</td></tr>
<tr><td>Top-10 scaffolds share (%)</td>
<td>{{ scaf.approved.top10_pct }}%</td>
<td>{{ scaf.clinical.top10_pct }}%</td></tr>
</table>
{% if ta %}
<h2>6. Therapeutic Area Analysis of Uncovered Space</h2>
<img src="data:image/png;base64,{{ ta_dist_b64 }}" alt="Therapeutic areas distribution">
<img src="data:image/png;base64,{{ ta_enrich_b64 }}" alt="Therapeutic area enrichment">
<p>ATC class enrichment ratio > 1.0 indicates overrepresentation in uncovered space
(clinical pipeline underinvesting relative to approved drugs).</p>
<table>
<tr><th>ATC Class</th><th>Therapeutic Area</th>
<th>Uncovered fraction</th><th>Covered fraction</th><th>Enrichment ratio</th></tr>
{% for code, ratio in enrichment_sorted %}
<tr style="{% if ratio > 1.2 %}background:#ffebee{% elif ratio < 0.8 %}background:#e8f5e9{% endif %}">
<td>{{ code }}</td>
<td>{{ ta.atc_names.get(code, code) }}</td>
<td>{{ ta.uncovered_atc_frac.get(code, 0) }}</td>
<td>{{ ta.covered_atc_frac.get(code, 0) }}</td>
<td><strong>{{ ratio }}</strong></td>
</tr>
{% endfor %}
</table>
{% endif %}
<h2>{% if ta %}7{% else %}6{% endif %}. Conclusions</h2>
<p>
At the primary threshold (Tanimoto ≥ {{ cov.tanimoto_threshold }}),
{{ cov.coverage_pct }}% of approved drug chemical space has a clinical-stage
structural neighbour, compared to a random baseline of {{ cov.random_baseline_pct }}%.
At the strictest threshold tested,
{{ cov.coverage_by_threshold.values() | list | last | attr("coverage_pct") }}% is covered.
The clinical pipeline is
{% if scaf.clinical.diversity_index > scaf.approved.diversity_index %}more{% else %}less{% endif %}
structurally diverse than the approved set
(diversity index {{ scaf.clinical.diversity_index }} vs. {{ scaf.approved.diversity_index }}).
KS tests indicate that clinical candidates differ significantly from approved
drugs in all {{ ks | length }} physicochemical properties tested.
</p>
</body></html>
"""
html = Template(TEMPLATE).render(
cur = cur,
cov = cov,
scaf = scaf,
ks = ks,
ta = ta,
umap_ds_b64 = img_b64(f"{OUTPUT_DIR}/figures/umap_dataset.png"),
umap_ph_b64 = img_b64(f"{OUTPUT_DIR}/figures/umap_phase.png"),
props_b64 = img_b64(f"{OUTPUT_DIR}/figures/property_distributions.png"),
scaf_b64 = img_b64(f"{OUTPUT_DIR}/figures/scaffold_diversity.png"),
ta_dist_b64 = img_b64(f"{OUTPUT_DIR}/figures/therapeutic_areas.png"),
ta_enrich_b64 = img_b64(f"{OUTPUT_DIR}/figures/therapeutic_enrichment.png"),
# Pre-sort enrichment for Jinja2 (avoids attribute sort filter incompatibility)
enrichment_sorted = (
sorted(ta["enrichment_ratio"].items(), key=lambda x: -x[1]) if ta else []
),
)
out = f"{OUTPUT_DIR}/report.html"
pathlib.Path(out).write_text(html, encoding="utf-8")
print(f"Report saved to {out}")
```
```bash
python scripts/09_report.py
```
**Expected final output:** `diversity_output/report.html` — open in any browser.
---
## Expected Outputs
```
diversity_output/
├── data/
│ ├── approved_raw.csv
│ ├── clinical_raw.csv
│ ├── approved_curated.csv
│ ├── clinical_curated.csv
│ ├── approved_enriched.csv
│ ├── clinical_enriched.csv
│ ├── approved_fps.npy
│ ├── clinical_fps.npy
│ ├── coverage_scores.csv # per-drug max Tanimoto + covered flags
│ ├── coverage_result.json # multi-threshold coverage table + baseline
│ ├── therapeutic_areas.json # ATC distribution + enrichment ratios
│ ├── scaffold_stats.json
│ ├── ks_results.json
│ └── umap_coords.csv
├── figures/
│ ├── umap_dataset.png
│ ├── umap_phase.png
│ ├── property_distributions.png
│ ├── scaffold_diversity.png
│ ├── therapeutic_areas.png # ATC distribution comparison
│ └── therapeutic_enrichment.png # enrichment ratio bar chart
└── report.html ← primary deliverable
```
---
## Run Order
```bash
python scripts/01_download_approved.py
python scripts/02_download_clinical.py
python scripts/03_curate.py
python scripts/04_descriptors.py
python scripts/05_coverage.py
python scripts/05b_therapeutic_areas.py
python scripts/06_umap.py
python scripts/07_properties.py
python scripts/08_scaffolds.py
python scripts/09_report.py
```
---
## Reproducing with Different Parameters
Edit `params.py`:
- Change `COVERAGE_THRESHOLDS` to add or remove thresholds
- Change `COVERAGE_TANIMOTO` to change the primary reporting threshold
- Change `CLINICAL_PHASES` to `[3]` to analyse only late-stage candidates
- All outputs regenerate automaticallyDiscussion (0)
to join the discussion.
No comments yet. Be the first to discuss this paper.