Scientific Informatics: Live Integrations and Code

#!/usr/bin/env python3
"""
ChEMBL → Scientific Database ETL Pipeline
Fetches bioactivity data for Kymera target proteins,
validates and loads into research data store.

Pattern mirrors CDD Vault / D360 integration workflows.
"""

import numpy as np
import pandas as pd
import requests
import logging
from dataclasses import asdict, dataclass, field
from typing import Iterator

logging.basicConfig(level=logging.INFO, format="%(asctime)s %(levelname)s %(message)s")
log = logging.getLogger(__name__)

CHEMBL_BASE = "https://www.ebi.ac.uk/chembl/api/data"

# Kymera priority targets
TARGETS = {
    "STAT6":  "CHEMBL5401",
    "IRAK4":  "CHEMBL3778",
    "IRF5":   "CHEMBL5365",
}

@dataclass
class ActivityRecord:
    target: str
    molecule_id: str
    smiles: str | None
    assay_type: str
    standard_type: str           # IC50, Ki, Kd, etc.
    value_nm: float | None
    pchembl: float | None
    assay_description: str
    journal: str | None
    year: int | None
    flags: list[str] = field(default_factory=list)

def fetch_activities(target_name: str, chembl_id: str) -> Iterator[dict]:
    """Paginate ChEMBL activity endpoint for a given target."""
    offset, limit = 0, 100
    while True:
        url = (
            f"{CHEMBL_BASE}/activity.json"
            f"?target_chembl_id={chembl_id}"
            f"&assay_type=B&pchembl_value__isnull=false"
            f"&limit={limit}&offset={offset}"
        )
        r = requests.get(url, timeout=30)
        r.raise_for_status()
        data = r.json()
        rows = data.get("activities", [])
        log.info(f"{target_name}: fetched {len(rows)} records (offset={offset})")
        yield from rows
        if not data["page_meta"]["next"]:
            break
        offset += limit

def parse_record(target: str, raw: dict) -> ActivityRecord:
    """Parse and validate a raw ChEMBL activity dict."""
    try:
        value = float(raw["standard_value"]) if raw.get("standard_value") else None
    except (ValueError, TypeError):
        value = None

    try:
        pchembl = float(raw["pchembl_value"]) if raw.get("pchembl_value") else None
    except (ValueError, TypeError):
        pchembl = None

    flags = []
    if raw.get("potential_duplicate"):
        flags.append("POTENTIAL_DUPLICATE")
    if value is not None and value < 0:
        flags.append("NEGATIVE_VALUE")
    if raw.get("data_validity_comment"):
        flags.append(f"DVC:{raw['data_validity_comment']}")

    return ActivityRecord(
        target=target,
        molecule_id=raw.get("molecule_chembl_id", ""),
        smiles=raw.get("canonical_smiles"),
        assay_type=raw.get("assay_type", ""),
        standard_type=raw.get("standard_type", ""),
        value_nm=value,
        pchembl=pchembl,
        assay_description=raw.get("assay_description", ""),
        journal=raw.get("document_journal"),
        year=raw.get("document_year"),
        flags=flags,
    )

def validate_dataframe(df: pd.DataFrame) -> pd.DataFrame:
    """
    Curate and validate bioactivity data:
    - Remove duplicates and flagged records
    - Filter to accepted activity types
    - Enforce value range (1 pM – 100 µM = 0.001–100000 nM)
    """
    initial = len(df)
    df = df[~df["flags"].apply(lambda f: "POTENTIAL_DUPLICATE" in f)]
    df = df[df["standard_type"].isin(["IC50", "Ki", "Kd", "EC50", "GI50"])]
    df = df[df["value_nm"].between(0.001, 100_000)]
    df = df.dropna(subset=["smiles", "value_nm"])
    log.info(f"Validation: {initial} → {len(df)} records ({initial - len(df)} removed)")
    return df

def run_pipeline(output_path: str = "kymera_bioactivity.parquet") -> pd.DataFrame:
    """
    Full ETL: Extract → Transform → Validate → Load
    Output: curated Parquet file ready for downstream analysis
    or CDD Vault / D360 import.
    """
    all_records = []
    for target_name, chembl_id in TARGETS.items():
        for raw in fetch_activities(target_name, chembl_id):
            all_records.append(asdict(parse_record(target_name, raw)))

    df = pd.DataFrame(all_records)
    df = validate_dataframe(df)
    df["pIC50"] = df["pchembl"].fillna(9 - np.log10(df["value_nm"]))
    df.to_parquet(output_path, index=False)
    log.info(f"Loaded {len(df)} curated records → {output_path}")
    return df

if __name__ == "__main__":
    df = run_pipeline()
    print(df.groupby("target")[["value_nm", "pchembl"]].describe().round(2))