Blood-Brain Barrier Penetration Prediction

# Blood–Brain Barrier Penetration from SMILES **Task.** Binary classification: given a molecule's SMILES string, predict whether it penetrates the blood–brain barrier (`p_np` = 1) or not (`p_np` = 0). **Data.** [MoleculeNet BBBP](https://arxiv.org/abs/1703.00564) — 2,053 raw rows (`num, name, p_np, smiles`) pulled from the DeepChem S3 bucket, with the HuggingFace mirror `scikit-fingerprints/MoleculeNet_BBBP` as fallback. The label is heavily skewed (**≈76% penetrating / 24% non-penetrating**), so accuracy alone is misleading — AUROC and F1 are the headline metrics. Roughly 14 SMILES fail RDKit sanitization and are dropped, which is why published work reports ~2,039 usable molecules. **Protocol.** Clean → canonicalize → dedupe → **Bemis–Murcko scaffold split (80/10/10)** → define a *single* shared eval harness → **then** branch and train. Scaffold splitting means the test set contains molecular frameworks the model has never seen, which is far harder (and far more honest) than a random split. **Two modelling arms, forked in parallel:** | Branch | Approach | |---|---| | **A** | Fine-tune [`DeepChem/ChemBERTa-10M-MLM`](https://huggingface.co/DeepChem/ChemBERTa-10M-MLM) — a small RoBERTa SMILES encoder (hidden 384, 3 layers, vocab 600) — via `AutoModelForSequenceClassification` with class-weighted loss | | **B** | RDKit **Morgan/ECFP4** fingerprints (radius 2, 2048 bits) → **XGBoost** (`scale_pos_weight`) plus a regularized **logistic-regression** reference | Both arms are…

# === Setup: dependencies, seeds, versions ===
import subprocess, sys, importlib, os

# This sandbox ships PYTHONPATH="/pkg/:/root/" and HOME=/root, but /root is
# unreadable as uid 1000 -> pip's environment scan dies with PermissionError.
# Sanitize the env for any pip subprocess.
PIP_ENV = {**os.environ, "HOME": "/home/user", "PYTHONPATH": "/pkg"}

def ensure(pkg, import_name=None):
    """Install `pkg` only if it isn't importable. Returns True if installed."""
    name = import_name or pkg.split("==")[0].replace("-", "_")
    try:
        importlib.import_module(name)
        return False
    except ImportError:
        subprocess.check_call([sys.executable, "-m", "pip", "install", "-q", pkg], env=PIP_ENV)
        importlib.invalidate_caches()
        return True
…

Pinned environment
  python         3.11.13
  numpy          2.1.2
  pandas         2.2.3
  rdkit          2026.03.5
  scikit-learn   1.5.2
  xgboost        3.1.3
  torch          2.7.1+cu128
  transformers   5.13.1
  plotly         5.24.1

device = cuda   |   seed = 42

# === Load MoleculeNet BBBP ===
S3_URL = "https://deepchemdata.s3.us-west-1.amazonaws.com/datasets/BBBP.csv"
HF_MIRROR = "scikit-fingerprints/MoleculeNet_BBBP"

source_used = None
try:
    raw = pd.read_csv(S3_URL)
    source_used = f"DeepChem S3 ({S3_URL})"
except Exception as e:
    print(f"S3 fetch failed ({type(e).__name__}: {e}) — falling back to HF mirror")
    from datasets import load_dataset
    raw = load_dataset(HF_MIRROR, split="train").to_pandas()
    raw = raw.rename(columns={"SMILES": "smiles", "label": "p_np"})
    source_used = f"HuggingFace mirror ({HF_MIRROR})"

print(f"source      : {source_used}")
print(f"shape       : {raw.shape}")
print(f"columns     : {list(raw.columns)}")
…

source      : DeepChem S3 (https://deepchemdata.s3.us-west-1.amazonaws.com/datasets/BBBP.csv)
shape       : (2050, 4)
columns     : ['num', 'name', 'p_np', 'smiles']

dtypes:
num        int64
name      object
p_np       int64
smiles    object

class balance (raw):
  non-penetrant (0)      483    23.6%
  penetrant (1)         1567    76.4%
  imbalance ratio       3.24 : 1

unique canonical-input SMILES strings: 2050 of 2050

# === Clean: RDKit sanitize -> canonicalize -> dedupe (BEFORE any splitting) ===
from rdkit import Chem

audit = []          # row-by-row record of what we drop and why
n0 = len(raw)
audit.append(("raw rows from source", n0, ""))

df = raw.copy()
df["mol"] = df["smiles"].apply(Chem.MolFromSmiles)

bad_mask = df["mol"].isna()
bad_smiles = df.loc[bad_mask, ["num", "name", "smiles"]]
df = df.loc[~bad_mask].copy()
audit.append(("dropped: RDKit could not sanitize", -int(bad_mask.sum()), "invalid valence / unparseable"))

# canonical (isomeric) SMILES — the string both branches will actually consume
df["canonical_smiles"] = df["mol"].apply(lambda m: Chem.MolToSmiles(m, isomericSmiles=True))

…

step    n                          note
                             raw rows from source 2050                              
                dropped: RDKit could not sanitize  -11 invalid valence / unparseable
dropped: duplicate SMILES with CONFLICTING labels  -20       10 ambiguous structures
          dropped: duplicate SMILES, label agreed  -54         kept first occurrence
                dropped: fewer than 3 heavy atoms   -1      not a drug-like molecule
                            FINAL c…

# === Physicochemical descriptors (interpretable CNS-permeability axes) ===
from rdkit.Chem import Descriptors, Crippen, QED, rdMolDescriptors as rdmd

DESCRIPTORS = {
    "MolWt":     Descriptors.MolWt,
    "cLogP":     Crippen.MolLogP,
    "TPSA":      Descriptors.TPSA,
    "HBD":       Descriptors.NumHDonors,
    "HBA":       Descriptors.NumHAcceptors,
    "RotB":      Descriptors.NumRotatableBonds,
    "AromRings": rdmd.CalcNumAromaticRings,
    "HeavyAtoms": lambda m: m.GetNumHeavyAtoms(),
    "FracCSP3":  rdmd.CalcFractionCSP3,
    "QED":       QED.qed,
    "MolMR":     Crippen.MolMR,
    "FormalCharge": lambda m: Chem.GetFormalCharge(m),
}

…

computed 12 descriptors for 1964 molecules in 4.8s

Descriptor means by class, ranked by rank-biserial effect size:
  descriptor   mean_pen  mean_non  effect_size_rb    p_value
        TPSA      54.16     128.4         -0.6665 7.966e-105
         HBD      1.059     3.211         -0.6219  1.626e-98
         HBA      3.624      7.28         -0.5643  2.006e-77
         QED     0.6743    0.4542          0.5177  5.627e-64
       MolWt      317.7     441.3          -0.418  2.433e-42
       cLogP…

# === EDA: class balance, descriptor separation, correlation, CNS window ===
KEY = ["TPSA", "HBD", "HBA", "QED", "MolWt", "cLogP"]

fig = make_subplots(
    rows=3, cols=3,
    specs=[[{"type": "domain"}, {"type": "xy"}, {"type": "xy"}],
           [{"type": "xy"}, {"type": "xy"}, {"type": "xy"}],
           [{"type": "xy"}, {"type": "heatmap", "colspan": 2}, None]],
    subplot_titles=("Class balance", *[f"<b>{k}</b>" for k in KEY],
                    "TPSA vs cLogP density", "Descriptor correlation (Spearman)"),
    vertical_spacing=0.11, horizontal_spacing=0.09,
)

# (1) donut
vc = df["label_str"].value_counts()
fig.add_trace(go.Pie(labels=vc.index, values=vc.values, hole=0.58,
                     marker=dict(colors=[PALETTE[l] for l in vc.index],
                                 line=dict(color="white", width=2)),
…

# === 3D chemical space: ECFP4 PCA and the physicochemical CNS box ===
from rdkit.Chem import rdFingerprintGenerator
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler

FP_RADIUS, FP_BITS = 2, 2048
mfpgen = rdFingerprintGenerator.GetMorganGenerator(radius=FP_RADIUS, fpSize=FP_BITS)

def ecfp_matrix(mols):
    return np.array([mfpgen.GetFingerprintAsNumPy(m) for m in mols], dtype=np.uint8)

t0 = time.time()
FP_ALL = ecfp_matrix(df["mol"].tolist())
print(f"ECFP{2*FP_RADIUS} bit matrix {FP_ALL.shape} in {time.time()-t0:.1f}s "
      f"| mean bits set = {FP_ALL.sum(1).mean():.1f} | density = {FP_ALL.mean():.3%}")

pca = PCA(n_components=3, random_state=SEED).fit(FP_ALL.astype(np.float32))
P3 = pca.transform(FP_ALL.astype(np.float32))
…

ECFP4 bit matrix (1964, 2048) in 0.1s | mean bits set = 43.0 | density = 2.099%
PCA-3 explained variance: 5.1%, 3.0%, 2.7% (total 10.9%)