""" drug_repurposing.py: Registry Forge drug repurposing module ============================================================ Assembles the data sets a drug repurposing study in an ALS registry starts from: who is in the motor neuron disease cohort and when they entered it, which medicines each person was exposed to under explicit criteria, what evidence each exposure rests on, and how long each person was followed to a defined data cut-off. The exposure criteria adapt the methodology of Reimer et al. (Lancet Digital Health, 2026): "Identification of drug repurposing candidates for amyotrophic lateral sclerosis using electronic health records: a retrospective cohort study." Cohort assembly is in scope here; propensity score matching, Cox models and other causal analyses are out of scope and belong in a prespecified analysis with dedicated tools. Method ------ 1. Cohort: every participant with a motor neuron disease code in problems[] (ICD-10-CM G12.20 to G12.29, ICD-9-CM 335.20 to 335.29, the SNOMED CT and Mondo codes in COHORT_CODES). Cohort entry is the earliest such code date. 2. Exposure (per medication): criterion 'any' (default) counts any record in the window from 6 months before to 12 months after cohort entry; the strict Reimer criteria A and B are available for dispense-event data. The evidence types behind each exposure (order, dispense, administration, reported use, C-CDA list entry) are kept, because a prescription row is not evidence that a medicine was taken. 3. Follow-up: from cohort entry to a defined data cut-off. A death on or before the cut-off ends follow-up as an event. Otherwise follow-up ends at the cut-off when the participant is confirmed alive at or after it, at the last confirmed alive date, or at the last dated record. A participant is never assumed alive at the cut-off. Deaths come from the Patient resource and from an optional vital_status.csv (for example National Death Index linkage results). 4. Summary: per medication, counts and Kaplan-Meier median survival measured from a landmark 12 months after cohort entry, among participants alive and under follow-up at the landmark. The landmark removes the immortal time the exposure window would otherwise create. These are crude descriptive figures, not effect estimates. 5. Panels: the ALS clinical-interest medicines (approved, withdrawn and of interest) are always shown. Medicines are grouped by ATC class through ATC_SEED. Outputs ------- drug_repurposing_cohort.csv one row per (participant, medication) for every medication meeting min_exposed; pseudonymized; keeps calendar dates, so it is restricted analysis data drug_repurposing_summary.csv one row per medication, counts and landmark medians, with small cells held back cure_id_review_worksheet.csv exposed (participant, medication) rows in the shape of a CURE ID case report, for clinical review; not a CURE ID submission drug_repurposing_report.html single-file report, counts only off_label_snippets.md representative source records for chart review Privacy: participant IDs are pseudonymized to PT-NNNN, ages are banded to the CURE ID age groups, and every count below k (default 5) is held back for cell size, together with any figure that would let a held-back count be worked out. This module is for hypothesis generation. It is not a substitute for a prespecified protocol, confounding control, or prospective trials. Registry Forge (Boyce D, Premasiri A, Sullivan S, Levine B, Vieira FG. JAMIA Open 2026;9(5):ooag200). MIT License. """ import os import csv import html import json import re import datetime as _dt from collections import defaultdict, Counter # ============================================================================ # CONFIG # ============================================================================ DEFAULT_BUNDLE_PATH = './dashboard_data.json' DEFAULT_OUT_DIR = './drug_repurposing_output' # Reimer 2026 thresholds MIN_EXPOSED_PATIENTS = 30 # medications below this are dropped CRITERION_A_PRE_MONTHS = 0 # criterion A: dispense within 12 mo AFTER dx CRITERION_A_POST_MONTHS = 12 CRITERION_A_MIN_DURATION_MONTHS = 6 # end must be at least 6 mo after dispense CRITERION_B_PRE_MONTHS = 6 # criterion B: 6 mo before to 12 mo after CRITERION_B_POST_MONTHS = 12 CRITERION_B_MIN_DISPENSES = 2 # Privacy: counts below k are held back for cell size, including for internal use DEFAULT_K_ANONYMITY = 5 # Survival summary: landmark (months after cohort entry) from which follow-up is measured LANDMARK_MONTHS = 12 # ============================================================================ # COHORT-DEFINING CODES: motor neuron disease spectrum # ============================================================================ COHORT_CODES = { # ICD-10-CM ('ICD-10-CM', 'G12.21'): 'Amyotrophic lateral sclerosis', ('ICD-10-CM', 'G12.20'): 'Motor neuron disease, unspecified', ('ICD-10-CM', 'G12.22'): 'Progressive bulbar palsy', ('ICD-10-CM', 'G12.23'): 'Primary lateral sclerosis (adult)', ('ICD-10-CM', 'G12.29'): 'Other motor neuron disease', # ICD-9-CM (legacy, present in older records) ('ICD-9-CM', '335.20'): 'Amyotrophic lateral sclerosis', ('ICD-9-CM', '335.21'): 'Progressive muscular atrophy', ('ICD-9-CM', '335.22'): 'Progressive bulbar palsy', ('ICD-9-CM', '335.23'): 'Pseudobulbar palsy', ('ICD-9-CM', '335.24'): 'Primary lateral sclerosis', ('ICD-9-CM', '335.29'): 'Other motor neuron disease', # SNOMED CT ('SNOMED-CT', '86044005'): 'Amyotrophic lateral sclerosis', ('SNOMED-CT', '230258005'): 'Motor neuron disease', # Mondo ('Mondo', 'MONDO:0004976'): 'Amyotrophic lateral sclerosis', ('Mondo', 'MONDO:0019056'): 'Motor neuron disease', } _SYSTEM_ALIASES = {'SNOMED CT': 'SNOMED-CT', 'SNOMED': 'SNOMED-CT', 'SNOMEDCT': 'SNOMED-CT', 'ICD10CM': 'ICD-10-CM', 'ICD-10': 'ICD-10-CM', 'ICD10': 'ICD-10-CM', 'ICD9CM': 'ICD-9-CM', 'ICD-9': 'ICD-9-CM', 'ICD9': 'ICD-9-CM', 'MONDO': 'Mondo'} def _norm_system(sysname): s = (sysname or '').strip() return _SYSTEM_ALIASES.get(s.upper().replace('_', ' '), _SYSTEM_ALIASES.get(s, s)) def _is_cohort_code(sysname, code): sysname = _norm_system(sysname) code = (code or '').strip() return (sysname, code.upper()) in COHORT_CODES or (sysname, code) in COHORT_CODES # ============================================================================ # RxNorm -> ATC class seed mapping # Seed covers the Reimer candidates, ALS medicines and medicines common in ALS care. # Extend this dict for the medicines in your registry. # ============================================================================ ATC_SEED = { # Statins — HMG-CoA reductase inhibitors (C10AA) 'simvastatin': ('C10AA01', 'HMG-CoA reductase inhibitors'), 'lovastatin': ('C10AA02', 'HMG-CoA reductase inhibitors'), 'pravastatin': ('C10AA03', 'HMG-CoA reductase inhibitors'), 'atorvastatin': ('C10AA05', 'HMG-CoA reductase inhibitors'), 'rosuvastatin': ('C10AA07', 'HMG-CoA reductase inhibitors'), 'pitavastatin': ('C10AA08', 'HMG-CoA reductase inhibitors'), 'fluvastatin': ('C10AA04', 'HMG-CoA reductase inhibitors'), # PDE5 inhibitors (G04BE) 'sildenafil': ('G04BE03', 'PDE5 inhibitors (drugs for erectile dysfunction)'), 'vardenafil': ('G04BE09', 'PDE5 inhibitors (drugs for erectile dysfunction)'), 'tadalafil': ('G04BE08', 'PDE5 inhibitors (drugs for erectile dysfunction)'), 'avanafil': ('G04BE10', 'PDE5 inhibitors (drugs for erectile dysfunction)'), # Alpha-adrenoreceptor antagonists (G04CA) 'tamsulosin': ('G04CA02', 'Alpha-adrenoreceptor antagonists'), 'terazosin': ('G04CA03', 'Alpha-adrenoreceptor antagonists'), 'alfuzosin': ('G04CA01', 'Alpha-adrenoreceptor antagonists'), 'silodosin': ('G04CA04', 'Alpha-adrenoreceptor antagonists'), 'doxazosin': ('C02CA04', 'Alpha-adrenoreceptor antagonists'), # cardiovasc indication # Centrally acting muscle relaxants (M03BX) — Reimer reported as class-level 'cyclobenzaprine':('M03BX08', 'Centrally acting muscle relaxants'), 'baclofen': ('M03BX01', 'Centrally acting muscle relaxants'), 'tizanidine': ('M03BX02', 'Centrally acting muscle relaxants'), # Medicines approved for ALS (ATC codes checked against the WHO ATC index, November 2026) 'riluzole': ('N07XX02', 'Other nervous system drugs (motor neuron disease)'), 'edaravone': ('N07XX14', 'Other nervous system drugs (motor neuron disease)'), 'tofersen': ('N07XX22', 'Other nervous system drugs (SOD1 ALS)'), 'sodium phenylbutyrate-taurursodiol': ('N07XX19', 'Other nervous system drugs (motor neuron disease; withdrawn in the US in 2024)'), # Dextromethorphan-quinidine combination (FDA-approved for pseudobulbar # affect in ALS; trade name Nuedexta). Reimer 2026 reported this combo # AND standalone dextromethorphan separately, so we track both. 'dextromethorphan-quinidine': ('N07XX59', 'Other nervous system drugs (pseudobulbar affect; brand name Nuedexta)'), 'nuedexta': ('N07XX59', 'Other nervous system drugs (pseudobulbar affect; brand name for dextromethorphan-quinidine)'), # Standalone dextromethorphan — Reimer reported this as a separate harm-direction # signal under R05DA (cough suppressants / antitussives). Even though it is the # active CNS ingredient in Nuedexta, when prescribed alone it is almost always # an antitussive and should be tracked under that class. 'dextromethorphan': ('R05DA09', 'Opium alkaloid antitussive (also active CNS ingredient in Nuedexta)'), # Standalone quinidine — Class IA antiarrhythmic. In Nuedexta it is present at # subtherapeutic dose to inhibit CYP2D6 metabolism of dextromethorphan; when # prescribed alone at therapeutic doses it is an antiarrhythmic. 'quinidine': ('C01BA01', 'Class IA antiarrhythmic (also CYP2D6-inhibitor ingredient in Nuedexta at subtherapeutic dose)'), 'meloxicam': ('M01AC06', 'Oxicam non-steroidal anti-inflammatory drug'), 'colecalciferol': ('A11CC05', 'Vitamin D and analogues'), 'cholecalciferol':('A11CC05', 'Vitamin D and analogues'), 'rivaroxaban': ('B01AF01', 'Direct factor Xa inhibitor (anticoagulant)'), 'hydrochlorothiazide':('C03AA03','Thiazide diuretic'), 'lisinopril': ('C09AA03', 'ACE inhibitor'), # Reimer harmful-direction candidates (for full disclosure) 'glycopyrrolate': ('A03AB02', 'Quaternary ammonium anticholinergic (for sialorrhea)'), 'ondansetron': ('A04AA01', '5HT3 antagonist (antiemetic)'), 'fentanyl': ('N02AB03', 'Opioid analgesic'), 'morphine': ('N02AA01', 'Opioid analgesic'), 'oxycodone': ('N02AA05', 'Opioid analgesic'), 'hydromorphone': ('N02AA03', 'Opioid analgesic'), 'methadone': ('N07BC02', 'Opioid analgesic / addiction treatment'), 'gabapentin': ('N02BF01', 'Gabapentinoid (neuropathic pain)'), 'pseudoephedrine':('R01BA02', 'Decongestant'), 'guaifenesin': ('R05CA03', 'Mucolytic'), 'cetirizine': ('R06AE07', 'H1 antihistamine'), 'salbutamol': ('R03AC02', 'Short-acting beta agonist'), 'albuterol': ('R03AC02', 'Short-acting beta agonist'), 'montelukast': ('R03DC03', 'Leukotriene receptor antagonist'), } # ============================================================================ # REGULATORY STATUS — used to filter the CURE ID off-label-use export # ============================================================================ # Medications with current FDA approval for an ALS indication. These are # on-label for ALS and therefore NOT off-label-use candidates for CURE ID. ALS_FDA_APPROVED = { 'riluzole', 'edaravone', 'tofersen', } # Medications previously FDA-approved for ALS but since withdrawn from the # US market. Excluded from the off-label-use export because they were # initially approved on-label; any continuing use is a separate regulatory # situation (compassionate use, ongoing supply, etc.), not a repurposing # candidate. Listed here for documentation and to keep the script's # decision-making transparent. ALS_FDA_WITHDRAWN = { # Brand name Relyvrio / Albrioza; withdrawn from market in 2024 after # the PHOENIX confirmatory trial did not meet its primary endpoint. 'sodium phenylbutyrate-taurursodiol', } # ============================================================================ # ALS CLINICAL-INTEREST MEDICATIONS (always shown in the dashboard panel) # ============================================================================ # Medications that always appear in the report's clinical-interest panel regardless of # whether they meet `min_exposed`. When fewer than k participants have a record, the # participant count, year range, record count and evidence types are all held back for # cell size, since a year range beside a small count can point to one person. ALS_CLINICAL_INTEREST = [ { 'med_key': 'riluzole', 'display': 'Riluzole', 'brands': 'Rilutek, Tiglutik, Exservan', 'atc': 'N07XX02', 'fda_status': 'approved', 'label': 'FDA-approved for ALS (1995)', }, { 'med_key': 'edaravone', 'display': 'Edaravone', 'brands': 'Radicava, Radicava ORS', 'atc': 'N07XX14', 'fda_status': 'approved', 'label': 'FDA-approved for ALS (2017)', }, { 'med_key': 'tofersen', 'display': 'Tofersen', 'brands': 'Qalsody', 'atc': 'N07XX22', 'fda_status': 'approved', 'label': 'FDA-approved for SOD1 ALS (2023, accelerated)', }, { 'med_key': 'dextromethorphan-quinidine', 'display': 'Dextromethorphan-quinidine', 'brands': 'Nuedexta', 'atc': 'N07XX59', 'fda_status': 'off_label_for_als', 'label': 'FDA-approved for PBA; off-label-of-interest in ALS', }, { 'med_key': 'sodium phenylbutyrate-taurursodiol', 'display': 'Sodium phenylbutyrate-taurursodiol', 'brands': 'Relyvrio, Albrioza', 'atc': 'N07XX19', 'fda_status': 'withdrawn', 'label': 'FDA-approved for ALS 2022; withdrawn 2024 after PHOENIX', }, ] # ============================================================================ # COMBINATION-DRUG DETECTION RULES # ============================================================================ # Each rule: (set of substrings that must ALL be present, case-insensitive, # in the raw lowercased medication name) -> normalized generic name to return. # # These rules run inside _norm_med_name BEFORE the dosage-stripping regex. # That order is critical: dosage stripping cuts the name at the first # " mg" match and drops everything after it, which silently loses the # second ingredient of any combination. For example: # # "dextromethorphan hydrobromide 20 MG / quinidine sulfate 10 MG Oral Capsule" # after dose strip -> "dextromethorphan hydrobromide" # first token -> "dextromethorphan" # ATC lookup -> R05DA09 (cough suppressant) ← WRONG # # With combination pre-detection both ingredients are seen first, and the # string resolves correctly to N07XX59 (Nuedexta / dextromethorphan-quinidine, # FDA-approved for pseudobulbar affect). # # Add new rules here as more combination drugs become relevant. The ALS- # relevant combinations currently covered: # - Nuedexta (dextromethorphan HBr 20 mg + quinidine sulfate 10 mg) # - Relyvrio / Albrioza (sodium phenylbutyrate + taurursodiol; # FDA-approved for ALS 2022, withdrawn 2024) COMBINATION_RULES = [ # Nuedexta — FDA-approved for pseudobulbar affect; off-label-of-interest in ALS (frozenset({'nuedexta'}), 'dextromethorphan-quinidine'), (frozenset({'dextromethorphan', 'quinidine'}), 'dextromethorphan-quinidine'), # Relyvrio / Albrioza — FDA-approved for ALS 2022, withdrawn 2024 (frozenset({'relyvrio'}), 'sodium phenylbutyrate-taurursodiol'), (frozenset({'albrioza'}), 'sodium phenylbutyrate-taurursodiol'), (frozenset({'phenylbutyrate', 'taurursodiol'}), 'sodium phenylbutyrate-taurursodiol'), (frozenset({'phenylbutyrate', 'tudca'}), 'sodium phenylbutyrate-taurursodiol'), ] # ============================================================================ # Brand-name → generic-name resolution # ---------------------------------------------------------------------------- # Real EHR data routinely names medications by brand (Lipitor, Flomax, # Nuedexta, Robinul) rather than generic. This dictionary is consulted by # _norm_med_name BEFORE the ATC_SEED lookup, so a record with display name # "Lipitor 20 mg" gets normalized to "atorvastatin" and matches the # atorvastatin entry in ATC_SEED. Brands that map to combination drugs # (e.g. Nuedexta → dextromethorphan-quinidine) are also handled here. # # This list focuses on (1) ALS-approved drugs, (2) drugs Reimer 2026 # identified as protective or harm-direction, and (3) drugs commonly # prescribed for ALS symptom management. Adopters extend this dict to # cover brand names common in their site's EHR. BRAND_TO_GENERIC = { # ALS FDA-approved 'rilutek': 'riluzole', 'tiglutik': 'riluzole', 'exservan': 'riluzole', 'radicava': 'edaravone', 'radicava ors': 'edaravone', 'qalsody': 'tofersen', 'nuedexta': 'dextromethorphan-quinidine', 'relyvrio': 'sodium phenylbutyrate-taurursodiol', 'albrioza': 'sodium phenylbutyrate-taurursodiol', # Statins (Reimer-protective) 'lipitor': 'atorvastatin', 'crestor': 'rosuvastatin', 'zocor': 'simvastatin', 'mevacor': 'lovastatin', 'altoprev': 'lovastatin', 'pravachol': 'pravastatin', 'livalo': 'pitavastatin', 'zypitamag': 'pitavastatin', 'lescol': 'fluvastatin', # PDE5 inhibitors (Reimer-protective) 'viagra': 'sildenafil', 'revatio': 'sildenafil', 'levitra': 'vardenafil', 'staxyn': 'vardenafil', 'cialis': 'tadalafil', 'adcirca': 'tadalafil', 'stendra': 'avanafil', # Alpha-adrenoreceptor antagonists (Reimer-protective) 'flomax': 'tamsulosin', 'hytrin': 'terazosin', 'uroxatral': 'alfuzosin', 'rapaflo': 'silodosin', 'cardura': 'doxazosin', # Centrally acting muscle relaxants (Reimer-protective) 'flexeril': 'cyclobenzaprine', 'amrix': 'cyclobenzaprine', 'fexmid': 'cyclobenzaprine', 'zanaflex': 'tizanidine', 'lioresal': 'baclofen', 'gablofen': 'baclofen', 'ozobax': 'baclofen', # Sialorrhea / ALS symptom management (Reimer harm-direction) 'robinul': 'glycopyrrolate', 'cuvposa': 'glycopyrrolate', 'glycate': 'glycopyrrolate', 'dartisla': 'glycopyrrolate', # Antiemetics (Reimer harm-direction) 'zofran': 'ondansetron', 'zuplenz': 'ondansetron', # Opioids (Reimer harm-direction; end-of-life ALS) 'sublimaze': 'fentanyl', 'duragesic': 'fentanyl', 'actiq': 'fentanyl', 'fentora': 'fentanyl', 'abstral': 'fentanyl', 'subsys': 'fentanyl', 'lazanda': 'fentanyl', # Gabapentinoids (Reimer harm-direction) 'neurontin': 'gabapentin', 'gralise': 'gabapentin', 'horizant': 'gabapentin', # OTC cough / cold (Reimer harm-direction) 'sudafed': 'pseudoephedrine', 'mucinex': 'guaifenesin', 'robitussin': 'guaifenesin', # Antihistamines (Reimer harm-direction) 'zyrtec': 'cetirizine', 'reactine': 'cetirizine', # Bronchodilators (Reimer harm-direction) 'ventolin': 'salbutamol', 'proventil': 'salbutamol', 'proair': 'salbutamol', 'accuneb': 'salbutamol', 'xopenex': 'levalbuterol', # Leukotriene receptor antagonists (Reimer harm-direction) 'singulair': 'montelukast', # ACE inhibitor brands (for the lisinopril seed entry) 'prinivil': 'lisinopril', 'zestril': 'lisinopril', 'qbrelis': 'lisinopril', # Anticoagulants 'xarelto': 'rivaroxaban', # Vitamin D — many alternate spellings 'cholecalciferol': 'colecalciferol', 'vitamin d3': 'colecalciferol', 'vitamin d': 'colecalciferol', 'd3': 'colecalciferol', } # ============================================================================ # CURE ID age band mapping (from the Cure ID intake CRF, Generic Other Disease) # ============================================================================ CURE_ID_AGE_BANDS = [ (0, 1/12, '<1 month'), (1/12, 3/12, '1 - 3 months'), (3/12, 6/12, '4 - 6 months'), (6/12, 11/12, '7 - 11 months'), (1, 5, '1 - 5 years'), (6, 10, '6 - 10 years'), (11, 15, '11 - 15 years'), (16, 20, '16 - 20 years'), (21, 30, '21 - 30 years'), (31, 40, '31 - 40 years'), (41, 50, '41 - 50 years'), (51, 60, '51 - 60 years'), (61, 70, '61 - 70 years'), (71, 80, '71 - 80 years'), (81, 89, '81 - 89 years'), (90, 150, '90+ years'), ] CURE_ID_RACE_MAP = { # Map race display strings from the bundle to CURE ID race options 'Asian': 'Asian', 'Black or African American': 'Black or African American', 'White': 'White', 'Native Hawaiian or Other Pacific Islander': 'Native Hawaiian or Other Pacific Islander', 'American Indian or Alaska Native': 'American Indian or Alaska Native', 'Other Race': 'Other', } # ============================================================================ # UTILITIES # ============================================================================ def log(msg): print(msg, flush=True) def _parse_iso_date(s): """Parse an ISO date / datetime string into a date. Returns None on failure.""" if not s: return None s = str(s).strip() if not s: return None # Try several common formats for fmt in ('%Y-%m-%d', '%Y-%m-%dT%H:%M:%S', '%Y-%m-%dT%H:%M:%SZ', '%Y-%m-%dT%H:%M:%S.%f', '%Y-%m-%dT%H:%M:%S.%fZ', '%Y/%m/%d', '%m/%d/%Y'): try: return _dt.datetime.strptime(s.split('+')[0], fmt).date() except ValueError: continue # Last resort: year-only m = re.match(r'^(\d{4})', s) if m: try: return _dt.date(int(m.group(1)), 1, 1) except ValueError: return None return None def _months_between(d1, d2): """Approximate months between two dates.""" if d1 is None or d2 is None: return None return (d2 - d1).days / 30.44 def _norm_med_name(name): """Normalize medication name for ATC lookup. Lowercases, strips dosage and formulation suffixes, then consults BRAND_TO_GENERIC to resolve brand names (Lipitor, Nuedexta, Flomax) to their generic equivalents so brand-named EHR records match the ATC_SEED keys (which use generic names). For combination drugs (e.g. dextromethorphan-quinidine), the slash and hyphen are normalized so 'dextromethorphan/quinidine' and 'dextromethorphan-quinidine' both resolve to the same string. For these multi-token combo names we keep the full hyphenated form; otherwise we take the first space-separated token. Combination drugs in RxNorm SCD form (e.g. "dextromethorphan hydrobromide 20 MG / quinidine sulfate 10 MG Oral Capsule [Nuedexta]") are detected by COMBINATION_RULES BEFORE the dosage-stripping step, which would otherwise truncate at "20 MG" and lose the second ingredient, causing Nuedexta to be misclassified as standalone dextromethorphan (R05DA09 cough suppressant) instead of the PBA combination (N07XX59). This was a real bug in earlier versions. """ if not name: return '' s = str(name).lower().strip() # === Combination-drug pre-detection. MUST happen before dosage stripping. for rule_words, generic_name in COMBINATION_RULES: if all(w in s for w in rule_words): return generic_name # Strip common suffixes (dosage, formulation) s = re.sub(r'\b\d+\s*(?:mg|mcg|µg|g|ml|meq|units?|iu)\b.*$', '', s, flags=re.I) s = re.sub(r'\b(?:tablet|capsule|injection|oral|solution|suspension|er|xr|sr|cr|extended[\s-]release).*$', '', s, flags=re.I) # Normalize slash to hyphen so 'dextromethorphan/quinidine' becomes # 'dextromethorphan-quinidine'. This must happen BEFORE the first-token # split or we'd lose the second ingredient. s = re.sub(r'\s*/\s*', '-', s) s = re.sub(r'\s+', ' ', s).strip() # Brand-name resolution: if the lowercased name matches a brand # (either as full string or as the first token), substitute the generic. # Try full-string match first to catch multi-word brands like "radicava ors". if s in BRAND_TO_GENERIC: return BRAND_TO_GENERIC[s] # Then try first-token match for "lipitor 20 mg-fragment-leftover" cases. first = s.split(' ')[0].strip() if ' ' in s else s if first in BRAND_TO_GENERIC: return BRAND_TO_GENERIC[first] # For known combination generic names (presence of a hyphen and both # halves are recognizable generics), keep the full hyphenated form so # 'dextromethorphan-quinidine' doesn't get truncated to 'dextromethorphan'. if '-' in first and first in ATC_SEED: return first # Otherwise return the first token (matches single-ingredient drugs). return first def _age_band(age_years): if age_years is None: return 'Unknown' for low, high, label in CURE_ID_AGE_BANDS: if low <= age_years < high + (1 if high < 1 else 0.99): return label return 'Unknown' def _to_pseudo_id(patient_id, mapping, prefix='PT'): """Convert a participant ID to PREFIX-NNNN form (stable within one run).""" if patient_id not in mapping: mapping[patient_id] = f'{prefix}-{len(mapping)+1:04d}' return mapping[patient_id] # ============================================================================ # COHORT ASSEMBLY # ============================================================================ def _coding_system_name(c): """Return the system name for a coding dict, trying every shape that run_pipeline.py / FHIR / synthetic test bundles can emit. run_pipeline.py emits codings as {'code', 'system_oid', 'system_name', 'display'}, with system_name being a human label like 'ICD-10-CM', 'SNOMED-CT', 'RxNorm'. Some synthetic fixtures or older bundles use plain 'system'. We fall through all three. """ return ((c.get('system_name') or c.get('system') or '').strip() or OID_TO_NAME.get(c.get('system_oid', ''), '')) # Minimal OID-to-name fallback for cases where only system_oid is populated OID_TO_NAME = { '2.16.840.1.113883.6.90': 'ICD-10-CM', '2.16.840.1.113883.6.103': 'ICD-9-CM', '2.16.840.1.113883.6.96': 'SNOMED-CT', '2.16.840.1.113883.6.88': 'RxNorm', '2.16.840.1.113883.6.1': 'LOINC', '2.16.840.1.113883.6.12': 'CPT-4', '2.16.840.1.113883.6.285':'HCPCS', } # ============================================================================ def identify_cohort(bundle, verbose=True): """Walk problems[] and return {patient_id: cohort_entry_date} for every participant with at least one motor neuron disease problem code (cohort entry = earliest). Uses a robust system-name reading (system_name / system / OID-derived) that matches the shape run_pipeline.py emits. Also accepts the top-level code_system field on problem rows for FHIR-derived records that ended up with an empty all_codings[] but have the prioritized code/system on the row itself. """ cohort = {} n_problems = 0 n_matched_codes = 0 n_unparseable_dates = 0 matched_systems = Counter() for prob in bundle.get('problems', []): n_problems += 1 codings = prob.get('all_codings') or [] # Also synthesize a coding from the prioritized top-level fields if # all_codings is empty (FHIR rows sometimes end up this way). if not codings and prob.get('code'): codings = [{ 'code': prob.get('code'), 'system_name': prob.get('code_system') or prob.get('system_name') or '', 'system_oid': prob.get('system_oid') or '', 'display': prob.get('display_name') or prob.get('display') or '', }] for c in codings: sysname = _coding_system_name(c) code = (c.get('code') or '').strip() if _is_cohort_code(sysname, code): n_matched_codes += 1 matched_systems[sysname] += 1 pid = prob.get('patient_id') if not pid: continue d = _parse_iso_date(prob.get('effective_date') or prob.get('onset_date') or prob.get('recorded_date')) if d is None: n_unparseable_dates += 1 continue if pid not in cohort or d < cohort[pid]: cohort[pid] = d break # one match per problem is enough if verbose: print(f' identify_cohort: scanned {n_problems} problem rows') print(f' identify_cohort: matched {n_matched_codes} motor neuron disease codes') if matched_systems: for sys_, n in matched_systems.most_common(): print(f' by system: {sys_}: {n}') if n_unparseable_dates: print(f' identify_cohort: dropped {n_unparseable_dates} matches with unparseable dates') print(f' identify_cohort: identified {len(cohort)} unique cohort patients') return cohort def last_evidence_dates(bundle): """{participant_id: latest dated record of any kind in the bundle}.""" out = {} for key, rows in bundle.items(): if not isinstance(rows, list) or key in ('patients',): continue for r in rows: if not isinstance(r, dict): continue pid = r.get('patient_id') d = _parse_iso_date(r.get('effective_date') or r.get('document_date') or r.get('end_date')) if pid and d and (pid not in out or d > out[pid]): out[pid] = d return out EVIDENCE_TYPE = {'MedicationRequest': 'order', 'MedicationDispense': 'dispense', 'MedicationAdministration': 'administration', 'MedicationStatement': 'reported_use'} def medication_evidence_types(bundle): """{participant_id: {normalized medication: set of evidence types}}. A C-CDA medication row is recorded as 'document_list' because C-CDA does not say whether it was taken.""" out = defaultdict(lambda: defaultdict(set)) for m in bundle.get('medications', []): pid = m.get('patient_id') name = _norm_med_name(m.get('display_name') or m.get('display') or '') if pid and name: out[pid][name].add(EVIDENCE_TYPE.get(m.get('fhir_resource_type') or m.get('medication_subtype') or '', 'document_list')) return out def read_vital_status(path): """vital_status.csv (templates/vital_status.csv): participant_id, death_date, last_confirmed_alive_date. Returns {participant_id: {'death': date, 'alive': date}}.""" if not path or not os.path.exists(path): return {} out = {} with open(path, encoding='utf-8-sig', newline='') as f: for r in csv.DictReader(f): pid = (r.get('participant_id') or '').strip() if pid: out[pid] = {'death': _parse_iso_date(r.get('death_date')), 'alive': _parse_iso_date(r.get('last_confirmed_alive_date'))} return out def follow_up_end(entry, death, alive, last_evidence, cutoff): """(end_date, event, basis, vital_status_at_cutoff). Never assumes a participant is alive at the cut-off, and never runs follow-up to the date the script runs.""" if death and death <= cutoff: return death, 1, 'death', 'deceased' if alive and alive >= cutoff: return cutoff, 0, 'administrative_cutoff', 'alive_at_cutoff' candidates = [d for d in (alive, last_evidence) if d and d <= cutoff] if alive and alive <= cutoff and (not last_evidence or alive >= last_evidence): return alive, 0, 'last_confirmed_alive', 'unknown_after_last_confirmation' if candidates: return max(candidates), 0, 'last_evidence_date', 'unknown' return entry, 0, 'no_follow_up', 'unknown' def kaplan_meier_median(times_events): """Median survival time from (time_days, event) pairs with right censoring. Returns (median_days or None, reached: bool). Reached is False when the curve stays above 0.5.""" data = sorted((t, e) for t, e in times_events if t is not None and t >= 0) if not data: return None, False s = 1.0 at_risk = len(data) i = 0 while i < len(data): t = data[i][0] deaths = sum(1 for tt, e in data if tt == t and e) removed = sum(1 for tt, _ in data if tt == t) if deaths: s *= (1 - deaths / at_risk) if s <= 0.5: return t, True at_risk -= removed while i < len(data) and data[i][0] == t: i += 1 return None, False def gather_patient_demographics(bundle): """Return {patient_id: {sex, race, ethnicity, dob, marital_status}}.""" demo = {} for p in bundle.get('patients', []): pid = p.get('patient_id') or p.get('id') if not pid: continue demo[pid] = { 'sex': (p.get('gender') or p.get('sex') or '').strip() or 'Unknown', 'race': (p.get('race') or '').strip() or 'Unknown', 'ethnicity': (p.get('ethnicity') or '').strip() or 'Unknown', 'dob': _parse_iso_date(p.get('birth_date') or p.get('dob')), 'marital_status': (p.get('marital_status') or {}).get('display', '') if isinstance(p.get('marital_status'), dict) else (p.get('marital_status') or '').strip(), 'deceased_date': _parse_iso_date(p.get('deceased_date_time') or p.get('death_date')), } return demo def collect_medication_events(bundle, verbose=True): """Return {patient_id: [(med_name_norm, raw_display, start_date, end_date), ...]}. Walks medications[] and normalizes each display name. Logs how many records resolved to a known ATC class (so empty cohort output explains itself). """ events = defaultdict(list) n_rows = 0 n_with_pid = 0 n_normed = 0 n_with_date = 0 n_with_end_date = 0 n_atc_matched = 0 unmatched_norms = Counter() for m in bundle.get('medications', []): n_rows += 1 pid = m.get('patient_id') if not pid: continue n_with_pid += 1 display_raw = m.get('display_name') or m.get('display') or '' med_name = _norm_med_name(display_raw) if not med_name: continue n_normed += 1 start = _parse_iso_date(m.get('effective_date') or m.get('start_date') or m.get('authored_on')) end = _parse_iso_date(m.get('end_date')) if start is not None: n_with_date += 1 if end is not None: n_with_end_date += 1 if med_name in ATC_SEED: n_atc_matched += 1 else: unmatched_norms[med_name] += 1 events[pid].append((med_name, display_raw, start, end)) if verbose: print(f' collect_medication_events: scanned {n_rows} medication rows') print(f' {n_with_pid} had patient_id') print(f' {n_normed} normalized to a non-empty name') print(f' {n_with_date} had a parseable start date') print(f' {n_with_end_date} had a parseable end_date') print(f' {n_atc_matched} matched an ATC_SEED entry') if n_with_end_date == 0 and n_rows > 0: print(f' NOTE: no medication rows have end_date. Reimer Criterion A ' f'(needs end_date) will never fire. Use criterion=\'any\' (default) ' f'for C-CDA / FHIR data.') if unmatched_norms: top_unmatched = unmatched_norms.most_common(15) print(f' Top 15 unmatched medication names (extend ATC_SEED / BRAND_TO_GENERIC to cover):') for name, n in top_unmatched: print(f' {n:5d} {name!r}') return events def apply_exposure_criteria(patient_meds, cohort_entry, criterion='any'): """Apply exposure criteria per medication. Parameters ---------- patient_meds : list of (med_name, raw_display, start_date, end_date) tuples cohort_entry : date criterion : str 'any' (default) Criterion C - any medication record with start_date in [-6, +12] months of cohort entry. Fits C-CDA / FHIR data where each record is one prescription event and end_date is typically missing. This is the right default for participant-contributed data and most patient-directed SMART on FHIR registries. 'reimer-a' Strict Reimer Criterion A only - start or end within 12mo of cohort entry AND end >= 6mo after start. Needs end_date so only fits VHA-style or dispense-event data. 'reimer-b' Strict Reimer Criterion B only - >=2 dispenses in window [-6, +12] months. Needs multiple rows per medication so only fits dispense-event data. 'reimer' Either Reimer A or B. Needs end_date or multiple rows. Returns {med_name: {'exposed', 'criterion_met', 'first_dispense', 'last_dispense', 'n_dispenses_window', 'raw_display'}}. """ by_med = defaultdict(list) for med, raw, start, end in patient_meds: by_med[med].append((raw, start, end)) out = {} for med, events in by_med.items(): # Sort by start date (None at end) events_sorted = sorted(events, key=lambda e: e[1] or _dt.date.max) starts = [e[1] for e in events_sorted if e[1] is not None] first_dispense = starts[0] if starts else None last_dispense = starts[-1] if starts else None criterion_met = None # Criterion C ('any'): at least one record with start in [-6,+12] mo of entry c_count = 0 for raw, start, end in events_sorted: if start is None: continue d = _months_between(cohort_entry, start) if d is None: continue if -CRITERION_B_PRE_MONTHS <= d <= CRITERION_B_POST_MONTHS: c_count += 1 c_met = c_count >= 1 # Criterion A (Reimer): dispense OR end within 12mo of cohort entry, # AND end at least 6mo after dispense a_met = False if criterion in ('reimer-a', 'reimer'): for raw, start, end in events_sorted: d_after_entry_start = _months_between(cohort_entry, start) if start else None d_after_entry_end = _months_between(cohort_entry, end) if end else None within = (d_after_entry_start is not None and 0 <= d_after_entry_start <= CRITERION_A_POST_MONTHS) or \ (d_after_entry_end is not None and 0 <= d_after_entry_end <= CRITERION_A_POST_MONTHS) duration_ok = (start is not None and end is not None and _months_between(start, end) is not None and _months_between(start, end) >= CRITERION_A_MIN_DURATION_MONTHS) if within and duration_ok: a_met = True break # Criterion B (Reimer): >=2 dispenses within window b_met = False if criterion in ('reimer-b', 'reimer'): b_met = c_count >= CRITERION_B_MIN_DISPENSES # Determine which criterion fired based on mode if criterion == 'any': if c_met: criterion_met = 'C' elif criterion == 'reimer-a': if a_met: criterion_met = 'A' elif criterion == 'reimer-b': if b_met: criterion_met = 'B' elif criterion == 'reimer': if a_met and b_met: criterion_met = 'A+B' elif a_met: criterion_met = 'A' elif b_met: criterion_met = 'B' out[med] = { 'exposed': bool(criterion_met), 'criterion_met': criterion_met, 'first_dispense': first_dispense, 'last_dispense': last_dispense, 'n_dispenses_window': c_count, 'raw_display': events_sorted[0][0] if events_sorted else '', } return out def assemble_cohort_table(cohort, demo, meds_per_patient, min_exposed, criterion='any', verbose=True, pseudonymize=True, cutoff=None, vital=None, last_evidence=None, groups=None, evidence_types=None): """Build the long-form per-(patient, medication) cohort table. When pseudonymize=True (default), patient_id values in the output rows are replaced with stable PT-NNNN pseudonyms. Raw patient_id values from EHR extracts (Epic patient IDs, FHIR resource IDs derived from Epic, MRNs, etc.) are HIPAA direct identifiers under Safe Harbor §164.514(b)(2)(i)(H) and should not appear in any output that may be shared beyond the IRB-approved circle. Setting pseudonymize=False emits the raw identifiers and is only appropriate for transient internal analyses where downstream chart-review linkage is needed; the resulting CSV must then be stored alongside PHI. """ # First pass: which medications meet min_exposed across the cohort med_exposed_count = Counter() per_patient_results = {} n_cohort_with_any_med = 0 for pid, entry_date in cohort.items(): meds = meds_per_patient.get(pid, []) if meds: n_cohort_with_any_med += 1 results = apply_exposure_criteria(meds, entry_date, criterion=criterion) per_patient_results[pid] = results for med, r in results.items(): if r['exposed']: med_exposed_count[med] += 1 keep_meds = {m for m, n in med_exposed_count.items() if n >= min_exposed} # Build the pseudonym mapping if requested. Stable per-run: the first # patient encountered gets PT-0001, the next gets PT-0002, etc. pseudo_map = {} if pseudonymize: for pid in sorted(cohort.keys()): _to_pseudo_id(pid, pseudo_map) if verbose: print(f' assemble_cohort_table: pseudonymize=True - {len(pseudo_map)} ' f'patient_id values will be replaced with PT-NNNN in cohort CSV') if verbose: print(f' assemble_cohort_table: criterion={criterion!r}') print(f' {n_cohort_with_any_med} of {len(cohort)} cohort patients had any medication record') print(f' {len(med_exposed_count)} unique medications had >=1 exposed cohort patient') if med_exposed_count: print(f' top exposed medications across cohort:') for med, n in med_exposed_count.most_common(10): marker = ' *KEPT*' if med in keep_meds else f' (dropped, < min_exposed={min_exposed})' shown = f'{n:4d}' if n >= DEFAULT_K_ANONYMITY else f' <{DEFAULT_K_ANONYMITY}' print(f' {shown} participants {med} {marker}') print(f' {len(keep_meds)} medications meet min_exposed={min_exposed}') if not keep_meds and med_exposed_count: print(f' HINT: no medication has {min_exposed} or more exposed participants. ' f'Lower --min-exposed only for exploratory work, and never below k.') rows = [] vital = vital or {} last_evidence = last_evidence or {} groups = groups or {} evidence_types = evidence_types or {} if cutoff is None: cutoff = max(last_evidence.values()) if last_evidence else _dt.date.today() if verbose: print(f' assemble_cohort_table: no --cutoff given; using the latest record date in the bundle ({cutoff})') for pid, entry_date in cohort.items(): d = demo.get(pid, {}) age_at_entry = None if d.get('dob') and entry_date: age_at_entry = (entry_date - d['dob']).days / 365.25 v = vital.get(pid, {}) deaths = [x for x in (v.get('death'), d.get('deceased_date')) if x] end, event, basis, status = follow_up_end(entry_date, min(deaths) if deaths else None, v.get('alive'), last_evidence.get(pid), cutoff) survival_days = (end - entry_date).days if end and entry_date and end >= entry_date else None results = per_patient_results.get(pid, {}) for med in keep_meds: r = results.get(med, {'exposed': False, 'criterion_met': None, 'first_dispense': None, 'last_dispense': None, 'n_dispenses_window': 0, 'raw_display': ''}) atc_code, atc_class = ATC_SEED.get(med, ('UNCLASSIFIED', 'Not in ATC seed mapping')) out_pid = pseudo_map.get(pid, pid) if pseudonymize else pid rows.append({ 'patient_id': out_pid, 'cohort_entry': entry_date.isoformat() if entry_date else '', 'age_at_entry': round(age_at_entry, 1) if age_at_entry else None, 'sex': d.get('sex', 'Unknown'), 'race': d.get('race', 'Unknown'), 'ethnicity': d.get('ethnicity', 'Unknown'), 'marital_status': d.get('marital_status', ''), 'medication': med, 'medication_display': r.get('raw_display', ''), 'atc_code': atc_code, 'atc_class': atc_class, 'exposed': 1 if r['exposed'] else 0, 'criterion': r.get('criterion_met') or '', 'n_dispenses_window': r.get('n_dispenses_window', 0), 'first_dispense': r['first_dispense'].isoformat() if r.get('first_dispense') else '', 'evidence_types': ';'.join(sorted(evidence_types.get(pid, {}).get(med, set()))), 'survival_days': survival_days, 'death_observed': event, 'follow_up_end': end.isoformat() if end else '', 'follow_up_end_basis': basis, 'vital_status_at_cutoff': status, 'data_cutoff': cutoff.isoformat(), }) return rows, keep_meds # ============================================================================ # SOURCE-RECORD SNIPPETS: traceback to the EHR records behind a medication # ============================================================================ # For each repurposing candidate, walk the bundle once more to pull a small # set of representative source records so a reviewer can verify findings and # prepare for chart review. Participant IDs are pseudonymized, dates are # reduced to year, free-text fields are truncated, resource UUIDs are not # written, and medicines with fewer than k participants are held back for # cell size. Medication display text can still hold free text typed by a # clinician, so the file is reviewed before it leaves the analysis environment. SNIPPET_DEFAULT_TARGETS = [ # Symptom-management mainstays for ALS 'baclofen', 'tizanidine', 'cyclobenzaprine', # antispasticity 'gabapentin', 'pregabalin', # neuropathic pain 'amitriptyline', 'nortriptyline', # sialorrhea, mood, pain 'ondansetron', # antiemetic (edaravone-related) 'glycopyrrolate', 'scopolamine', # sialorrhea 'morphine', 'oxycodone', 'hydromorphone', 'fentanyl', # palliative analgesia 'methadone', # Symptom-management adjuncts 'colecalciferol', 'cholecalciferol', 'vitamin d', # vitamin D supplementation 'albuterol', # respiratory # ALS-PBA (FDA-approved for PBA; off-label for ALS itself) 'dextromethorphan-quinidine', 'nuedexta', # Reimer-2026 protective-direction candidates worth checking 'atorvastatin', 'rosuvastatin', 'simvastatin', 'pravastatin', 'lovastatin', 'pitavastatin', 'fluvastatin', 'metformin', 'sertraline', 'fluoxetine', 'citalopram', 'escitalopram', ] def _snippet_year_only(date_str): """Reduce any date to year only for snippet output.""" if not date_str: return 'unknown' m = re.match(r'(\d{4})', str(date_str)) return m.group(1) if m else 'unknown' def _snippet_truncate(s, n=120): if not s: return '' s = str(s).strip() return s if len(s) <= n else s[:n-1] + '…' def _snippet_matches(display, code, target): """Substring match on display name or code, case-insensitive.""" t = target.lower() if display and t in display.lower(): return True if code and t in str(code).lower(): return True return False def inspect_off_label_records(bundle, targets=None, max_examples=5, min_records_to_show=DEFAULT_K_ANONYMITY): """Walk every medication record in the bundle and pull representative examples for each target drug. Returns {target_name: {'examples': [..], 'n_patients': N, 'n_records': M}} with privacy controls already applied. Parameters ---------- bundle : dict Parsed dashboard_data.json. targets : iterable[str] or None Drug names to search for. Defaults to SNIPPET_DEFAULT_TARGETS. max_examples : int Max example records to emit per drug. min_records_to_show : int Drugs with fewer than this many unique exposed patients are suppressed (k-anonymity). """ if targets is None: targets = SNIPPET_DEFAULT_TARGETS cache_pseudo = {} by_target = defaultdict(list) per_target_pt = defaultdict(set) meds = bundle.get('medications', []) if isinstance(meds, dict): flat = [] for pid, recs in meds.items(): for r in recs: if 'patient_id' not in r: r = {**r, 'patient_id': pid} flat.append(r) meds = flat for rec in meds: display = rec.get('display_name') or rec.get('display') or '' code = rec.get('code') or '' for target in targets: if _snippet_matches(display, code, target): pid = rec.get('patient_id') or rec.get('subject') or '' pseudo = _to_pseudo_id(pid, cache_pseudo, prefix='SN') example = { 'pseudo_patient_id': pseudo, 'display_name': _snippet_truncate(display), 'vocabulary': _snippet_truncate( rec.get('code_system') or rec.get('system_name') or rec.get('system') or '', 60), 'code': _snippet_truncate(code, 40), 'year': _snippet_year_only( rec.get('effective_date') or rec.get('authored_on') or rec.get('date')), 'source_format': rec.get('source') or 'unknown', 'record_type': EVIDENCE_TYPE.get(rec.get('fhir_resource_type') or '', 'document_list'), 'status': rec.get('status') or '', } by_target[target].append(example) per_target_pt[target].add(pseudo) break # one record matches one target; don't double-count out = {} for target, examples in by_target.items(): if len(per_target_pt[target]) >= min_records_to_show: seen_pts = set() kept = [] for ex in examples: if ex['pseudo_patient_id'] not in seen_pts: kept.append(ex) seen_pts.add(ex['pseudo_patient_id']) if len(kept) >= max_examples: break out[target] = { 'examples': kept, 'n_patients': len(per_target_pt[target]), 'n_records': len(examples), } return out def render_off_label_snippets_markdown(snippets, cohort_label='this cohort'): """Render the snippet dict as a markdown report ready to drop into a manuscript or docs page.""" out = [] out.append('# Source-record examples for repurposing candidates') out.append('') out.append(f'Representative medication records for cohort participants in **{cohort_label}**, ' f'drawn from the parsed bundle. Participant identifiers are pseudonymized (SN-NNNN, ' f'assigned independently of the cohort CSV), dates are reduced to year, display text ' f'is truncated, and medicines recorded for fewer than {DEFAULT_K_ANONYMITY} participants ' f'are held back for cell size. Display text can still hold words typed by a clinician, ' f'so review this file before it leaves the analysis environment. It is for source ' f'verification and chart review preparation, not for analysis.') out.append('') if not snippets: out.append('*No qualifying records found.*') return '\n'.join(out) for target, info in sorted(snippets.items(), key=lambda kv: -kv[1]['n_patients']): out.append(f'## {target}') out.append('') out.append(f'{info["n_patients"]} participants, {info["n_records"]} records, ' f'up to 5 examples shown') out.append('') out.append('| Participant | Display name in source | Vocabulary | Code | Year | Source format | Record type |') out.append('|---|---|---|---|---|---|---|') for ex in info['examples']: row = ('| ' + ' | '.join([ ex['pseudo_patient_id'], ex['display_name'] or '*(no display)*', ex['vocabulary'] or 'not recorded', ex['code'] or 'not recorded', ex['year'], ex['source_format'], ex['record_type'], ]) + ' |') out.append(row) out.append('') return '\n'.join(out) # ============================================================================ # CURE ID REVIEW WORKSHEET # ============================================================================ CRITERION_LABEL = {'C': 'any record from 6 months before to 12 months after cohort entry', 'A': 'Reimer criterion A', 'B': 'Reimer criterion B', 'A+B': 'Reimer criteria A and B'} def cure_id_intake_rows(cohort_rows, exposed_only=True, off_label_only=True, k=1): """Local CURE ID review worksheet: one row per exposed (participant, medication), in the shape of a CURE ID case report so a clinical reviewer can decide what, if anything, is reportable. It is not a CURE ID submission. With off_label_only (default), medicines approved for ALS (ALS_FDA_APPROVED) and those approved and later withdrawn (ALS_FDA_WITHDRAWN) are left out; every other medicine is marked for indication review. Every row starts as needs_review. """ pseudo = {} rows = [] med_counts = Counter(r['medication'] for r in cohort_rows if r['exposed']) for r in cohort_rows: if exposed_only and not r['exposed']: continue if k > 1 and med_counts[r['medication']] < k: continue med = (r['medication'] or '').strip().lower() if off_label_only and (med in ALS_FDA_APPROVED or med in ALS_FDA_WITHDRAWN): continue pid = r['patient_id'] # cohort rows are already pseudonymized (PT-NNNN) sex_norm = (r.get('sex') or 'Unknown').strip().title() if sex_norm not in ('Female', 'Male'): sex_norm = 'Unknown' races = [CURE_ID_RACE_MAP[x.strip()] for x in (r.get('race') or '').split(' / ') if x.strip() in CURE_ID_RACE_MAP] or ['Unknown'] days = r.get('survival_days') if r['death_observed'] and days is not None: outcome = ('Died within 12 months of cohort entry' if days <= 365 else 'Died within 24 months of cohort entry' if days <= 730 else 'Died more than 24 months after cohort entry') elif r.get('vital_status_at_cutoff') == 'alive_at_cutoff': outcome = 'Alive at the data cut-off' else: outcome = 'Vital status after last contact unknown' rows.append({ 'pseudo_patient_id': pid, 'age_group': _age_band(r.get('age_at_entry')), 'sex': sex_norm, 'races': ', '.join(races), 'medication_name': r['medication'], 'medication_atc_class': r['atc_class'], 'medication_atc_code': r['atc_code'], 'evidence_types': r.get('evidence_types', ''), 'exposure_criterion': CRITERION_LABEL.get(r['criterion'], r['criterion'] or ''), 'off_label_review': 'Not approved for ALS: review the indication and actual use', 'outcome_at_cutoff': outcome, 'indication': '', 'clinical_reviewer': '', 'review_status': 'needs_review', 'worksheet_status': 'Local review worksheet; not a CURE ID submission', }) return rows # ============================================================================ # SUMMARY + REPORT # ============================================================================ def build_summary(cohort_rows, keep_meds, k=DEFAULT_K_ANONYMITY, landmark_months=LANDMARK_MONTHS): """Per-medication summary with complementary suppression and landmark Kaplan-Meier medians. - Medians are Kaplan-Meier medians, which account for censoring, or "not reached". - Follow-up is measured from a landmark (default 12 months after cohort entry) among participants alive and under follow-up then, because exposure is assessed through month 12; measuring from entry would credit exposed participants with immortal time. - If either group is below k, every figure for that medication is held back; death counts between 1 and k-1 are held back too, along with the count that would let them be worked out. """ held = 'held back for cell size' landmark_days = round(landmark_months * 30.44) by_med = defaultdict(lambda: {'exposed': [], 'unexposed': []}) for r in cohort_rows: by_med[r['medication']]['exposed' if r['exposed'] else 'unexposed'].append(r) def landmark_set(rows): return [((r['survival_days'] - landmark_days), r['death_observed']) for r in rows if r['survival_days'] is not None and r['survival_days'] > landmark_days] summary = [] for med in sorted(keep_meds): e, u = by_med[med]['exposed'], by_med[med]['unexposed'] atc_code, atc_class = ATC_SEED.get(med, ('UNCLASSIFIED', 'Not in ATC seed mapping')) small = len(e) < k or len(u) < k le, lu = landmark_set(e), landmark_set(u) def deaths(rows): n = sum(1 for r in rows if r['death_observed']) return held if 0 < n < k else n def km(rows_l): if len(rows_l) < k: return held m, reached = kaplan_meier_median(rows_l) return round(m / 30.44, 1) if reached else 'not reached' if small: row = dict.fromkeys(['n_exposed', 'n_unexposed', 'n_deaths_exposed', 'n_deaths_unexposed', 'n_at_landmark_exposed', 'n_at_landmark_unexposed', 'km_median_months_from_landmark_exposed', 'km_median_months_from_landmark_unexposed'], held) else: row = {'n_exposed': len(e), 'n_unexposed': len(u), 'n_deaths_exposed': deaths(e), 'n_deaths_unexposed': deaths(u), 'n_at_landmark_exposed': len(le) if len(le) >= k or not le else held, 'n_at_landmark_unexposed': len(lu) if len(lu) >= k or not lu else held, 'km_median_months_from_landmark_exposed': km(le), 'km_median_months_from_landmark_unexposed': km(lu)} # Complementary suppression: exposed + unexposed always add up to a cohort total # (deaths, participants at the landmark), so one held-back cell would be # recoverable from the other. Hold both. for a, b in (('n_deaths_exposed', 'n_deaths_unexposed'), ('n_at_landmark_exposed', 'n_at_landmark_unexposed')): if held in (row[a], row[b]): row[a] = row[b] = held # Participants lost before the landmark (group size minus participants at the # landmark) are a count too; hold the landmark counts when that difference is small. if any(0 < n - m < k for n, m in ((len(e), len(le)), (len(u), len(lu)))): row['n_at_landmark_exposed'] = row['n_at_landmark_unexposed'] = held summary.append({'medication': med, 'atc_code': atc_code, 'atc_class': atc_class, 'als_status': ('approved for ALS' if med in ALS_FDA_APPROVED else 'withdrawn' if med in ALS_FDA_WITHDRAWN else ''), **row, 'landmark_months': landmark_months}) return summary HTML_TEMPLATE = """ {title}

Registry Forge: drug repurposing cohort report

{cohort_name} · data cut-off {cutoff} · generated {generated}
Privacy: participant IDs are pseudonymized and every count below {k} is shown as “held back for cell size”, along with any figure that would let a held-back count be worked out. The companion file drug_repurposing_cohort.csv keeps calendar dates, so it is restricted analysis data. Review this report in full before sharing it.

Cohort overview

{n_cohort}Participants in the motor neuron disease cohort {n_deaths}Deaths on or before the cut-off {n_alive}Confirmed alive at the cut-off {n_unknown}Vital status unknown after last contact {n_meds}Medicines with at least {min_exposed} exposed participants

Cohort: any motor neuron disease code in the problem list (ICD-10-CM G12.20 to G12.29, ICD-9-CM 335.20 to 335.29, plus the SNOMED CT and Mondo codes listed in the module). Cohort entry is the earliest such code. Exposure criterion: {criterion}. Follow-up runs from cohort entry to the data cut-off, ending earlier at death, at the last confirmed alive date or at the last dated record; nobody is assumed alive at the cut-off.

{group_table_html}
{candidate_panel_html}

Medication summary by ATC class

Kaplan-Meier medians are measured from a landmark {landmark} months after cohort entry, among participants alive and under follow-up at the landmark, because exposure is assessed through month {landmark}. “Not reached” means fewer than half of that group died during follow-up. These are crude descriptive figures with no adjustment for site of onset, disease stage, age or any other difference between the groups, so a difference between the two columns says nothing about the effect of the medicine.

{atc_groups_html}

Methods and limitations

This module assembles the cohort and assigns exposure following the approach of Reimer et al. (Lancet Digital Health, 2026). It does not run propensity score matching or survival regression; those steps need a prespecified analysis plan, dedicated tools and covariates such as site of onset, time from symptom onset, ALSFRS-R and FVC. drug_repurposing_cohort.csv is the input to that work.

Limitations include small groups, confounding by indication (a medicine prescribed for a comorbidity marks a different group of participants), incomplete medication histories in the records a registry receives, evidence that a medicine was ordered rather than taken, deaths missed when vital status has not been linked, and coded diagnoses that have not been checked against the clinical record. Treat every result as a question for a prespecified study.

cure_id_review_worksheet.csv lists each exposed participant and medicine in the shape of a CURE ID case report so a clinical reviewer can decide whether any use is reportable. It is a local review worksheet, not a CURE ID submission. Medicines approved for ALS, and the withdrawn phenylbutyrate-taurursodiol, are not in it.

""" def compute_candidate_summary(meds_per_patient, evidence_types=None, k=DEFAULT_K_ANONYMITY): """Presence, participant count, year range and evidence types for each medicine in ALS_CLINICAL_INTEREST, within the cohort. Counts, year ranges, record counts and evidence types are all held back when fewer than k participants have a record.""" held = 'held back for cell size' evidence_types = evidence_types or {} keys = {m['med_key'] for m in ALS_CLINICAL_INTEREST} alias = {'cholecalciferol': 'colecalciferol'} pts, recs, years, ev = defaultdict(set), Counter(), defaultdict(list), defaultdict(set) for pid, events in meds_per_patient.items(): for med_name, raw_display, start, end in events: key = alias.get(med_name, med_name) if key not in keys: continue pts[key].add(pid) recs[key] += 1 years[key] += [d.year for d in (start, end) if d is not None] ev[key] |= evidence_types.get(pid, {}).get(med_name, set()) rows = [] for entry in ALS_CLINICAL_INTEREST: key = entry['med_key'] n = len(pts.get(key, ())) if n == 0: out = {'present': False, 'n_participants': '0', 'year_range': '', 'n_records': '0', 'evidence': ''} elif n < k: out = {'present': True, 'n_participants': held, 'year_range': held, 'n_records': held, 'evidence': held} else: y = years.get(key, []) out = {'present': True, 'n_participants': str(n), 'year_range': f'{min(y)} to {max(y)}' if y else 'no dates recorded', 'n_records': f'{recs[key]:,}', 'evidence': ', '.join(sorted(ev.get(key, ()))) or 'not recorded'} rows.append({**entry, **out}) return rows def _cell(v): if isinstance(v, str) and v.startswith('held back'): return 'held back for cell size' if isinstance(v, int): return f'{v:,}' return html.escape(str(v)) if v not in (None, '') else '–' def render_candidate_panel(rows, k): if not rows: return '' trs = [] for r in rows: trs.append( f'{html.escape(r["display"])}
{html.escape(r["brands"])} · ATC {r["atc"]}' f'{html.escape(r["label"])}' f'{"Recorded in cohort" if r["present"] else "Not recorded"}' f'{_cell(r["n_participants"])}{_cell(r["year_range"])}' f'{_cell(r["n_records"])}{_cell(r["evidence"])}') return f"""

ALS clinical-interest medicines

Shown whether or not they meet the minimum exposure count. The status column describes the regulatory history only. Evidence types say what each record shows: an order, a dispense, an administration, a participant-reported use, or a C-CDA medication list entry that does not say whether the medicine was taken. Counts below {k} are held back for cell size, together with the year range, record count and evidence types for that medicine.

{''.join(trs)}
MedicineStatusIn cohortParticipantsYears with recordsRecordsEvidence types
""" def suppress_with_total(counts, k): """{label: n} -> {label: n or 'held back for cell size'} for categories that add up to a published total. Counts from 1 to k-1 are held back; when that leaves a single held-back cell, the next smallest is held back too so the first cannot be worked out by subtraction.""" held = 'held back for cell size' out = {g: (held if 0 < n < k else n) for g, n in counts.items()} if sum(1 for v in out.values() if v == held) == 1: rest = sorted((n, g) for g, n in counts.items() if out[g] != held and n > 0) if rest: out[rest[0][1]] = held return out def render_html(cohort_rows, summary, cohort_name, k, min_exposed=MIN_EXPOSED_PATIENTS, candidate_rows=None, cohort_overview=None, cutoff=None, criterion='any', landmark_months=LANDMARK_MONTHS): ov = cohort_overview or {} status = suppress_with_total({'n_deaths': ov.get('n_deaths', 0), 'n_alive': ov.get('n_alive', 0), 'n_unknown': ov.get('n_unknown', 0)}, k) def held_or(n): return n if n == 0 or n >= k else 'held back for cell size' by_class = defaultdict(list) for r in summary: by_class[r['atc_class']].append(r) parts = [] for atc_class in sorted(by_class): trs = [] for r in sorted(by_class[atc_class], key=lambda x: x['medication']): flag = f' {html.escape(r["als_status"])}' if r['als_status'] else '' km_e = r['km_median_months_from_landmark_exposed'] km_u = r['km_median_months_from_landmark_unexposed'] trs.append( f'{html.escape(r["medication"])}{flag}
{r["atc_code"]}' f'{_cell(r["n_exposed"])}{_cell(r["n_unexposed"])}' f'{_cell(r["n_deaths_exposed"])}{_cell(r["n_deaths_unexposed"])}' f'{_cell(r["n_at_landmark_exposed"])}{_cell(r["n_at_landmark_unexposed"])}' f'{_cell(km_e)}{_cell(km_u)}') parts.append(f"""

{html.escape(atc_class)}

{''.join(trs)}
MedicineExposedUnexposedDeaths, exposedDeaths, unexposedAt landmark, exposedAt landmark, unexposedKM median months, exposedKM median months, unexposed
""") return HTML_TEMPLATE.format( title=f'Drug repurposing cohort report: {cohort_name}', cohort_name=html.escape(cohort_name), cutoff=cutoff.isoformat() if cutoff else 'not set', generated=_dt.datetime.now(_dt.timezone.utc).strftime('%Y-%m-%d %H:%M UTC'), k=k, min_exposed=min_exposed, criterion=html.escape(criterion), landmark=landmark_months, n_cohort=_cell(held_or(ov.get('n_cohort', 0))), n_deaths=_cell(status['n_deaths']), n_alive=_cell(status['n_alive']), n_unknown=_cell(status['n_unknown']), n_meds=len(summary), group_table_html='', candidate_panel_html=render_candidate_panel(candidate_rows, k) if candidate_rows else '', atc_groups_html='\n'.join(parts) or '

No medicine met the minimum exposure count in this cohort.

', ) def _write_csv(path, rows): with open(path, 'w', newline='', encoding='utf-8') as f: if rows: w = csv.DictWriter(f, fieldnames=list(rows[0].keys()), lineterminator='\n') w.writeheader() for r in rows: w.writerow({key: (f"'{v}" if isinstance(v, str) and v[:1] in ('=', '+', '-', '@') else v) for key, v in r.items()}) # ============================================================================ # MAIN # ============================================================================ def main(bundle_path=DEFAULT_BUNDLE_PATH, out_dir=DEFAULT_OUT_DIR, min_exposed=5, k=DEFAULT_K_ANONYMITY, cohort_name='ALS registry cohort', criterion='any', verbose=True, pseudonymize=True, cutoff=None, vital_status_path=None, landmark_months=LANDMARK_MONTHS): """Run the drug repurposing cohort assembly against a Registry Forge bundle. Parameters ---------- bundle_path : path to dashboard_data.json produced by run_pipeline.py. out_dir : directory for the output files. min_exposed : minimum exposed participants to keep a medicine. Reimer used 30 in a cohort of about 11,000; smaller registries need a lower value, never below k. k : counts below k are held back for cell size everywhere they appear. Default 5. criterion : 'any' (default; fits C-CDA and FHIR registry data), or 'reimer-a', 'reimer-b', 'reimer' for dispense-event data. cutoff : data cut-off date (datetime.date or 'YYYY-MM-DD'). Follow-up ends here. When omitted, the latest dated record in the bundle is used and a warning is printed. vital_status_path : optional vital_status.csv with death and last-confirmed-alive dates. pseudonymize : True (default) writes PT-NNNN pseudonyms. False also writes drug_repurposing_cohort_identified.csv beside the pseudonymized file; treat it as PHI. """ if k < 1: raise ValueError('k must be at least 1') min_exposed = max(min_exposed, k) if isinstance(cutoff, str): cutoff = _parse_iso_date(cutoff) os.makedirs(out_dir, exist_ok=True) log('=' * 72) log('Registry Forge: drug repurposing cohort assembly') log('=' * 72) log(f'Bundle: {bundle_path}') log(f'Output: {out_dir}') log(f'Settings: min_exposed={min_exposed} k={k} criterion={criterion!r} landmark={landmark_months} months') with open(bundle_path, encoding='utf-8') as f: bundle = json.load(f) log(f' participants: {len(bundle.get("patients", [])):,}') log(f' problems: {len(bundle.get("problems", [])):,}') log(f' medications: {len(bundle.get("medications", [])):,}') cohort = identify_cohort(bundle, verbose=verbose) log(f'Cohort: {len(cohort)} participants with a motor neuron disease code') if not cohort: log(' No cohort participants identified. Common causes:') log(' - problems[] is empty or missing in the bundle') log(' - no problem coding uses a code listed in COHORT_CODES (ICD-10-CM G12.2x, ICD-9-CM 335.2x, SNOMED CT or Mondo)') log(' - every matching problem lacks a parseable effective_date or onset_date') demo = gather_patient_demographics(bundle) meds = collect_medication_events(bundle, verbose=verbose) vital = read_vital_status(vital_status_path) if vital_status_path and not vital: log(f' WARNING: no rows read from {vital_status_path}') last_ev = last_evidence_dates(bundle) ev_types = medication_evidence_types(bundle) rows_raw, keep_meds = assemble_cohort_table(cohort, demo, meds, min_exposed, criterion=criterion, verbose=verbose, pseudonymize=False, cutoff=cutoff, vital=vital, last_evidence=last_ev, evidence_types=ev_types) cutoff = _parse_iso_date(rows_raw[0]['data_cutoff']) if rows_raw else (cutoff or (max(last_ev.values()) if last_ev else _dt.date.today())) pseudo = {} for pid in sorted(cohort): _to_pseudo_id(pid, pseudo) cohort_rows = [{**r, 'patient_id': pseudo.get(r['patient_id'], r['patient_id'])} for r in rows_raw] log(f' medicines meeting min_exposed={min_exposed}: {len(keep_meds)}') # Per-participant overview (one row per participant, independent of medication) per_pid = {} for pid, entry in cohort.items(): d = demo.get(pid, {}) v = vital.get(pid, {}) deaths = [x for x in (v.get('death'), d.get('deceased_date')) if x] per_pid[pid] = follow_up_end(entry, min(deaths) if deaths else None, v.get('alive'), last_ev.get(pid), cutoff) overview = {'n_cohort': len(cohort), 'n_deaths': sum(1 for x in per_pid.values() if x[1] == 1), 'n_alive': sum(1 for x in per_pid.values() if x[3] == 'alive_at_cutoff'), 'n_unknown': sum(1 for x in per_pid.values() if x[3].startswith('unknown'))} log(f' deaths on or before {cutoff}: {overview["n_deaths"] if overview["n_deaths"] >= k or overview["n_deaths"] == 0 else f"<{k}"}') summary = build_summary(cohort_rows, keep_meds, k, landmark_months=landmark_months) worksheet = cure_id_intake_rows(cohort_rows, exposed_only=True, k=1) _write_csv(os.path.join(out_dir, 'drug_repurposing_cohort.csv'), cohort_rows) if not pseudonymize: _write_csv(os.path.join(out_dir, 'drug_repurposing_cohort_identified.csv'), rows_raw) log(' WARNING: drug_repurposing_cohort_identified.csv holds registry participant IDs. Treat it as PHI.') _write_csv(os.path.join(out_dir, 'drug_repurposing_summary.csv'), summary) _write_csv(os.path.join(out_dir, 'cure_id_review_worksheet.csv'), worksheet) log(f' wrote drug_repurposing_cohort.csv ({len(cohort_rows):,} rows), drug_repurposing_summary.csv ' f'({len(summary):,} rows), cure_id_review_worksheet.csv ({len(worksheet):,} rows)') cohort_meds = {pid: ev for pid, ev in meds.items() if pid in cohort} candidate_rows = compute_candidate_summary(cohort_meds, ev_types, k=k) if verbose: log('ALS clinical-interest medicines in the cohort (always shown):') for r in candidate_rows: log(f' {r["display"]:<22s} {r["n_participants"]:>24s} participants') report = render_html(cohort_rows, summary, cohort_name, k, min_exposed=min_exposed, candidate_rows=candidate_rows, cohort_overview=overview, cutoff=cutoff, criterion=criterion, landmark_months=landmark_months) with open(os.path.join(out_dir, 'drug_repurposing_report.html'), 'w', encoding='utf-8') as f: f.write(report) cohort_bundle = {**bundle, 'medications': [m for m in bundle.get('medications', []) if m.get('patient_id') in cohort]} snippets = inspect_off_label_records(cohort_bundle, min_records_to_show=k) with open(os.path.join(out_dir, 'off_label_snippets.md'), 'w', encoding='utf-8') as f: f.write(render_off_label_snippets_markdown(snippets, cohort_label=cohort_name)) log(f' wrote drug_repurposing_report.html and off_label_snippets.md ({len(snippets)} medicines shown)') log('Done.') return cohort_rows def _cli(argv=None): import argparse ap = argparse.ArgumentParser(description='Registry Forge: drug repurposing cohort assembly') ap.add_argument('--bundle', default=DEFAULT_BUNDLE_PATH, help='dashboard_data.json from run_pipeline.py') ap.add_argument('--out', default=DEFAULT_OUT_DIR, help='output folder') ap.add_argument('--cutoff', help='data cut-off date, YYYY-MM-DD (strongly recommended)') ap.add_argument('--vital-status', help='vital_status.csv with death and last-confirmed-alive dates') ap.add_argument('--min-exposed', type=int, default=5) ap.add_argument('--k', type=int, default=DEFAULT_K_ANONYMITY) ap.add_argument('--landmark-months', type=int, default=LANDMARK_MONTHS) ap.add_argument('--criterion', default='any', choices=['any', 'reimer-a', 'reimer-b', 'reimer']) ap.add_argument('--cohort-name', default='ALS registry cohort') ap.add_argument('--also-write-identified', action='store_true', help='also write a cohort CSV with registry participant IDs (PHI)') a = ap.parse_args(argv) main(bundle_path=a.bundle, out_dir=a.out, min_exposed=a.min_exposed, k=a.k, cohort_name=a.cohort_name, criterion=a.criterion, pseudonymize=not a.also_write_identified, cutoff=a.cutoff, vital_status_path=a.vital_status, landmark_months=a.landmark_months) if __name__ == '__main__': _cli()