Version history

{
  "cells": [
    {
      "source": "# Dirichlet inflection-point analysis: DAM/MGnD expansion vs OPC depletion across Braak stages\n\n**Plan ref:** 166d17f0-b18c-46a1-8af6-7923191fb523  \n**Cohorts:** SEA-AD DLPFC (84 donors, CELLxGENE Census), Mathys 2023 Cell (GEO GSE174367, ~427 donors), Green 2024 Nature (ROSMAP community fractions)  \n**Lineages:** microglia (homeostatic MG-1, DAM/MGnD MG-2, IRM MG-3, LDA MG-4), oligodendrocyte + OPC  \n**Staging variable:** Braak stage 0–VI (ordinal)  \n**Covariates:** age, sex, PMI, RIN  \n**Model:** scCODA (Dirichlet regression), one model per lineage  \n\n## Kill criteria\n- CELLxGENE Census obs-level query returns no `Braak_stage` or `donor_id` field → BLOCKED, file data-request to Allen portal\n- GEO GSE174367 supplement lacks per-donor count matrix with Braak metadata → fall back to Mathys 2024 multiregion (DOI: 10.1038/s41586-024-07606-7) or BLOCK\n- scCODA convergence diagnostics (R-hat > 1.05 on any composition parameter) → report but do not interpret\n\n## Expected output artifacts\n1. `sea_ad_mg_composition_by_braak.csv` — per-donor microglia subtype proportions × Braak stage\n2. `sea_ad_oligo_composition_by_braak.csv` — per-donor oligo/OPC proportions × Braak stage\n3. `dirichlet_mg_results.csv` — scCODA posterior summary (effect size, 95% CI, expected FDR) per subtype\n4. `dirichlet_oligo_results.csv` — same for oligo/OPC lineage\n5. `inflection_point_summary.md` — Braak stage of inflection per subtype with uncertainty bounds\n6. UMAP colored by Braak stage (batch-corrected, Harmony/scVI)\n7. Dotplot: marker gene expression by subtype (P2RY12, TMEM119, SPP1, GPNMB, TREM2, LPL, PLIN2, IFIT1, PDGFRA, MBP, MOG)",
      "cell_id": "c-fe1b7bb4",
      "outputs": [],
      "cell_hash": "sha256:e5a9b39b4ae17ffb369fa15354ff696b0f5d97892f8eeea7fffdb0d2dde87d00",
      "cell_type": "markdown",
      "execution_count": null
    },
    {
      "source": "# Cell 1: CELLxGENE Census — SEA-AD DLPFC microglia and oligodendrocyte obs pull\n# KILL CHECK: confirm obs fields include donor_id and Braak_stage before proceeding\n\nimport cellxgene_census\nimport pandas as pd\n\nSEA_AD_DATASET_TITLE = \"Seattle Alzheimer's Disease Brain Cell Atlas\"\nTARGET_TISSUE = \"prefrontal cortex\"\nTARGET_CELL_CLASSES = [\"microglia\", \"oligodendrocyte\", \"oligodendrocyte precursor cell\"]\n\nwith cellxgene_census.open_soma() as census:\n    obs_df = cellxgene_census.get_obs(\n        census,\n        organism=\"Homo sapiens\",\n        value_filter=(\n            f\"tissue_general == '{TARGET_TISSUE}' and \"\n            f\"dataset_title.str.contains('SEA-AD', case=False, na=False) and \"\n            f\"cell_type.isin({TARGET_CELL_CLASSES})\"\n        ),\n        column_names=[\n            \"soma_joinid\", \"dataset_id\", \"donor_id\", \"tissue\", \"tissue_general\",\n            \"cell_type\", \"assay\", \"sex\", \"development_stage\",\n            # SEA-AD-specific obs fields (verify presence; kill if absent):\n            \"Braak_stage\", \"CERAD_score\", \"Thal_phase\",\n            \"PMI\", \"RIN\", \"age_at_death\",\n            # Subtype annotation field (may be 'subclass_label' or 'cluster_label'):\n            \"subclass_label\", \"cluster_label\", \"supertype_label\"\n        ]\n    )\n\nprint(\"obs shape:\", obs_df.shape)\nprint(\"columns:\", obs_df.columns.tolist())\nprint(\"donor_id nunique:\", obs_df['donor_id'].nunique() if 'donor_id' in obs_df.columns else 'MISSING — KILL')\nprint(\"Braak_stage nunique:\", obs_df['Braak_stage'].nunique() if 'Braak_stage' in obs_df.columns else 'MISSING — KILL')\nprint(\"cell_type value_counts:\\n\", obs_df['cell_type'].value_counts())\nprint(\"subclass_label value_counts:\\n\", obs_df.get('subclass_label', pd.Series()).value_counts().head(20))\n\n# Save for downstream scCODA cell\nobs_df.to_csv('sea_ad_obs_raw.csv', index=False)\nprint('Saved sea_ad_obs_raw.csv')",
      "cell_id": "c-1bdfe25c",
      "outputs": [],
      "cell_hash": "sha256:35608890c63050ee86a719b9b5a70329784261d92d4361e12091020627779f9e",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Cell 2: Compute per-donor subtype proportions and run scCODA Dirichlet regression\n# Requires sea_ad_obs_raw.csv from Cell 1\n# KILL CHECK: R-hat > 1.05 on any posterior parameter → flag, do not interpret\n\nimport pandas as pd\nimport numpy as np\nfrom pertpy.tools import Sccoda  # scCODA via pertpy >=0.7\nimport anndata as ad\n\nobs_df = pd.read_csv('sea_ad_obs_raw.csv')\n\n# --- KILL check ---\nfor required_col in ['donor_id', 'Braak_stage', 'subclass_label']:\n    if required_col not in obs_df.columns:\n        raise ValueError(f'KILL: required column {required_col!r} absent from Census obs. '\n                         'File data-request to Allen portal before proceeding.')\n\n# Braak stage as ordinal integer (0–6)\nbraak_map = {'Braak 0': 0, 'Braak I': 1, 'Braak II': 2,\n             'Braak III': 3, 'Braak IV': 4, 'Braak V': 5, 'Braak VI': 6}\nobs_df['braak_int'] = obs_df['Braak_stage'].map(braak_map)\nobs_df = obs_df.dropna(subset=['braak_int', 'donor_id', 'subclass_label'])\n\n# --- Microglia subset ---\nmg_df = obs_df[obs_df['cell_type'] == 'microglia'].copy()\n\n# Per-donor × subtype counts\nmg_counts = (\n    mg_df.groupby(['donor_id', 'subclass_label'])\n    .size().unstack(fill_value=0)\n)\n\n# Donor metadata (one row per donor)\ndonor_meta = (\n    mg_df.groupby('donor_id')\n    [['braak_int', 'sex', 'PMI', 'RIN', 'age_at_death']]\n    .first()\n)\nmg_counts = mg_counts.join(donor_meta)\nmg_counts.to_csv('sea_ad_mg_composition_by_braak.csv')\nprint('Microglia donor×subtype matrix shape:', mg_counts.shape)\nprint(mg_counts.head(3))\n\n# --- scCODA: microglia lineage ---\n# Build AnnData with composition as X and donor metadata in obs\ncell_cols = [c for c in mg_counts.columns if c not in ['braak_int','sex','PMI','RIN','age_at_death']]\nmg_adata = ad.AnnData(\n    X=mg_counts[cell_cols].values,\n    obs=mg_counts[['braak_int','sex','PMI','RIN','age_at_death']].reset_index(drop=False),\n    var=pd.DataFrame(index=cell_cols)\n)\nmg_adata.obs.columns = [c.replace(' ','_') for c in mg_adata.obs.columns]\n\nsccoda_mg = Sccoda()\nsccoda_mg.load(mg_adata, type='cell_level', generate_sample_level=False)\nsccoda_mg.prepare(sample_identifier='donor_id', covariate_formula='braak_int + sex + PMI + RIN')\nsccoda_mg.run_nuts(num_samples=10000, num_warmup=2000, rng_key=42)\nresults_mg = sccoda_mg.summary_prepare(credible_effects_threshold=0.95)\nresults_mg.to_csv('dirichlet_mg_results.csv')\nprint('scCODA microglia results:')\nprint(results_mg[['covariate','cell_type','effect','HDI_lower','HDI_upper','inclusion_probability']].to_string())\n\n# R-hat kill check\nif hasattr(sccoda_mg, 'mcmc') and hasattr(sccoda_mg.mcmc, 'diagnostics'):\n    diag = sccoda_mg.mcmc.diagnostics()\n    rhat_max = float(diag['r_hat'].max())\n    print(f'Max R-hat (microglia model): {rhat_max:.4f}')\n    if rhat_max > 1.05:\n        print('WARNING: R-hat > 1.05 — model not converged. Results flagged; do not interpret.')\n\nprint('Saved dirichlet_mg_results.csv and sea_ad_mg_composition_by_braak.csv')",
      "cell_id": "c-d9881f60",
      "outputs": [],
      "cell_hash": "sha256:479e82f7c6cb6b24be50cd6ded59ec2068c44db549f4f9b4fb9247dfae86ef86",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Cell 3: Oligodendrocyte + OPC lineage — per-donor proportion table and scCODA model\n# Requires sea_ad_obs_raw.csv from Cell 1\n# KILL CHECK: R-hat > 1.05 on any posterior parameter → flag, do not interpret\n\nimport pandas as pd\nimport numpy as np\nfrom pertpy.tools import Sccoda\nimport anndata as ad\n\nobs_df = pd.read_csv('sea_ad_obs_raw.csv')\n\n# Confirm required columns\nfor required_col in ['donor_id', 'Braak_stage', 'subclass_label', 'sex', 'age_at_death', 'PMI', 'RIN']:\n    if required_col not in obs_df.columns:\n        available = list(obs_df.columns)\n        raise ValueError(f'KILL: required column {required_col!r} absent. Available: {available}')\n\n# Braak stage as ordinal integer (0–6)\nbraak_map = {\n    'Braak 0': 0, 'Braak I': 1, 'Braak II': 2,\n    'Braak III': 3, 'Braak IV': 4, 'Braak V': 5, 'Braak VI': 6\n}\nobs_df['braak_int'] = obs_df['Braak_stage'].map(braak_map)\nif obs_df['braak_int'].isna().any():\n    unmapped = obs_df.loc[obs_df['braak_int'].isna(), 'Braak_stage'].unique().tolist()\n    raise ValueError(f'KILL: unmapped Braak_stage values: {unmapped}')\n\n# --- OLIGODENDROCYTE + OPC lineage ---\nOLIGO_LABELS = {\n    'oligodendrocyte precursor cell': 'OPC',\n    'oligodendrocyte': 'Oligo_mature'\n}\noligo_obs = obs_df[obs_df['cell_type'].isin(OLIGO_LABELS.keys())].copy()\noligo_obs['oligo_subtype'] = oligo_obs['cell_type'].map(OLIGO_LABELS)\n\n# Per-donor cell counts by oligo subtype\noligo_counts = (\n    oligo_obs.groupby(['donor_id', 'oligo_subtype'])\n    .size()\n    .unstack(fill_value=0)\n    .reset_index()\n)\n\n# Attach donor-level covariates\ndonor_meta = obs_df[['donor_id', 'braak_int', 'sex', 'age_at_death', 'PMI', 'RIN']]\\\n    .drop_duplicates('donor_id').set_index('donor_id')\noligo_counts = oligo_counts.set_index('donor_id').join(donor_meta, how='left').reset_index()\noligo_counts.to_csv('sea_ad_oligo_composition_by_braak.csv', index=False)\nprint(f'Oligo/OPC donor count table shape: {oligo_counts.shape}')\nprint(oligo_counts.head())\n\n# scCODA: oligodendrocyte lineage\nsubtype_cols_oligo = [c for c in oligo_counts.columns if c in ['OPC', 'Oligo_mature']]\nif len(subtype_cols_oligo) < 2:\n    raise ValueError(f'KILL: fewer than 2 oligo subtype columns detected: {subtype_cols_oligo}')\n\ncounts_matrix_oligo = oligo_counts[subtype_cols_oligo].values\ncov_df_oligo = oligo_counts[['braak_int', 'sex', 'age_at_death', 'PMI', 'RIN']].copy()\ncov_df_oligo['sex_bin'] = (cov_df_oligo['sex'] == 'female').astype(int)\n\nadata_oligo = ad.AnnData(\n    X=counts_matrix_oligo,\n    obs=cov_df_oligo.reset_index(drop=True)\n)\nadata_oligo.var_names = subtype_cols_oligo\n\nsccoda_oligo = Sccoda()\nsccoda_oligo.load(adata_oligo, type='cell_level',\n                  covariate_formula='braak_int + sex_bin + age_at_death + PMI + RIN',\n                  cell_type_identifier='cell_type',\n                  sample_identifier='donor_id')\nsccoda_oligo.fit()\n\n# R-hat kill check\nrhat_oligo = sccoda_oligo.get_rhat()\nif (rhat_oligo > 1.05).any():\n    print('WARNING: R-hat > 1.05 detected in oligo model — results flagged, do not interpret')\nelse:\n    print('Oligo model R-hat OK (all ≤ 1.05)')\n\nresults_oligo = sccoda_oligo.credible_effects()\nresults_oligo.to_csv('dirichlet_oligo_results.csv')\nprint('Oligo Dirichlet results:')\nprint(results_oligo)",
      "cell_id": "c-c181145b",
      "outputs": [],
      "cell_hash": "sha256:383cca98b846e4981cd8b01d9813ed8aa5597186975acf09edf5285d09ede31a",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Cell 4: Mathys 2023 replication — per-donor microglia and oligo proportion extraction\n# Source: Mathys et al. Cell 2023, DOI 10.1016/j.cell.2023.08.039\n# GEO accession: GSE174367 (TODO_VERIFY — check supplement for per-donor count table)\n# KILL CHECK: if per-donor Braak metadata absent from supplement → BLOCKED\n\nimport pandas as pd\nimport os\n\n# ---- Load Mathys 2023 per-donor cell-type proportion table ----\n# Expected: supplement Table S1 or equivalent with columns:\n#   donor_id, Braak_stage (or braak_int), cell_type_proportions...\n# If file is not present, raise KILL to trigger data-request\n\nMATHYS_PROP_FILE = 'mathys2023_per_donor_proportions.csv'\nif not os.path.exists(MATHYS_PROP_FILE):\n    raise FileNotFoundError(\n        f'KILL: Mathys 2023 replication file {MATHYS_PROP_FILE!r} not found. '\n        'Required: per-donor cell-type proportion table with Braak staging from '\n        'GEO GSE174367 supplement or paper Table S1. '\n        'File data-request to ROSMAP / GEO before proceeding with replication.'\n    )\n\nmathys_df = pd.read_csv(MATHYS_PROP_FILE)\nprint(f'Mathys 2023 table shape: {mathys_df.shape}')\nprint(f'Columns: {list(mathys_df.columns)}')\n\n# Confirm Braak and donor columns\nfor col in ['donor_id', 'Braak_stage']:\n    if col not in mathys_df.columns:\n        raise ValueError(f'KILL: column {col!r} missing from Mathys 2023 file. Columns present: {list(mathys_df.columns)}')\n\nprint(f'Braak stage distribution (Mathys 2023):')\nprint(mathys_df['Braak_stage'].value_counts().sort_index())\n\n# Expected microglia subtypes: homeostatic (Mic.1 / P2RY12+), DAM/MGnD (Mic.2 / SPP1+), etc.\n# Map to MG-1…MG-4 schema for cross-cohort comparison\nMG_LABEL_MAP = {\n    'Mic.1': 'MG-1_homeostatic',\n    'Mic.2': 'MG-2_DAM',\n    'Mic.3': 'MG-3_IRM',\n    'Mic.4': 'MG-4_LDA'\n}\nOLIGO_LABEL_MAP = {\n    'Oli.1': 'Oligo_mature',\n    'OPC': 'OPC'\n}\n\navailable_mg = [c for c in MG_LABEL_MAP if c in mathys_df.columns]\navailable_oligo = [c for c in OLIGO_LABEL_MAP if c in mathys_df.columns]\nprint(f'Mathys MG columns found: {available_mg}')\nprint(f'Mathys Oligo/OPC columns found: {available_oligo}')\n\nif not available_mg:\n    print('WARNING: no recognizable MG subtype columns — cross-cohort MG replication not possible without column remapping')\n\nmathys_df.to_csv('mathys2023_validated_proportions.csv', index=False)",
      "cell_id": "c-82c5cba1",
      "outputs": [],
      "cell_hash": "sha256:9069ce85b99a7a00aad6ed2c48f3fa5ac81df9670b5daeb44694bc88da76916d",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Cell 5: Green 2024 community-fraction validation + inflection-point detection\n# Source: Green et al. Nature 2024, DOI 10.1038/s41586-024-07871-6\n# Goal: extract microglia and oligodendrocyte community fractions by Braak stage,\n#        compare trajectory shape to SEA-AD posterior means from Cells 2-3,\n#        detect inflection point via logistic derivative on posterior means.\n# KILL CHECK: if Green 2024 community fractions do not map to microglia or oligo\n#             subtypes (i.e., community labels are unlabeled or spatial-only),\n#             flag BLOCKED and fall back to single-cohort SEA-AD result.\n\nimport pandas as pd\nimport numpy as np\nfrom scipy.optimize import curve_fit\nimport os\n\n# ---- Inflection-point detection on SEA-AD posterior means ----\n# Load posterior mean proportions output from Cells 2 (microglia) and 3 (oligo)\n# Expected files: sea_ad_microglia_composition_by_braak.csv,\n#                 sea_ad_oligo_composition_by_braak.csv\n# Columns: subtype, braak_int (0-6), posterior_mean, credible_interval_lo, credible_interval_hi\n\ndef logistic(x, L, k, x0):\n    \"\"\"Standard logistic curve for inflection-point fitting.\"\"\"\n    return L / (1 + np.exp(-k * (x - x0)))\n\ndef fit_inflection(braak_vals, proportion_vals, subtype_name):\n    \"\"\"Fit logistic to braak_int vs proportion; return inflection Braak stage.\"\"\"\n    try:\n        p0 = [proportion_vals.max(), 1.0, 3.0]\n        popt, _ = curve_fit(logistic, braak_vals, proportion_vals,\n                            p0=p0, maxfev=5000)\n        L, k, x0 = popt\n        inflection = x0  # Braak stage at maximum rate of change\n        return {'subtype': subtype_name, 'inflection_braak': round(inflection, 2),\n                'L': round(L, 4), 'k': round(k, 4)}\n    except RuntimeError:\n        return {'subtype': subtype_name, 'inflection_braak': float('nan'),\n                'L': float('nan'), 'k': float('nan'),\n                'note': 'curve_fit did not converge'}\n\nresults = []\nfor fname, lineage_label in [\n    ('sea_ad_microglia_composition_by_braak.csv', 'microglia'),\n    ('sea_ad_oligo_composition_by_braak.csv', 'oligo')\n]:\n    if not os.path.exists(fname):\n        print(f'SKIP: {fname} not yet present (upstream cell not executed)')\n        continue\n    df = pd.read_csv(fname)\n    for subtype in df['subtype'].unique():\n        sub = df[df['subtype'] == subtype].sort_values('braak_int')\n        if sub['posterior_mean'].isna().any() or len(sub) < 4:\n            print(f'SKIP inflection for {subtype}: insufficient data points')\n            continue\n        res = fit_inflection(sub['braak_int'].values,\n                             sub['posterior_mean'].values, subtype)\n        res['lineage'] = lineage_label\n        results.append(res)\n\nif results:\n    inflection_df = pd.DataFrame(results)\n    inflection_df.to_csv('sea_ad_inflection_points_by_subtype.csv', index=False)\n    print('Inflection-point table:')\n    print(inflection_df.to_string(index=False))\nelse:\n    print('INFO: No upstream proportion files found; inflection detection deferred to post-execution tick.')\n\n# ---- Green 2024 community fraction stub ----\n# Planned: load per-donor community fraction table from Green 2024 supplement\n# (Nature 2024 DOI 10.1038/s41586-024-07871-6, Supplementary Table S3 or equivalent)\n# Expected columns: donor_id, braak_int, community_id, community_label, fraction\n# Success criterion: at least one community maps to microglia or oligodendrocyte lineage\nGREEN_FILE = 'green2024_community_fractions.csv'\nif not os.path.exists(GREEN_FILE):\n    print(f'INFO: {GREEN_FILE} absent. Green 2024 replication arm DEFERRED.')\n    print('Action required: download Supplementary Table with per-donor community fractions from DOI 10.1038/s41586-024-07871-6.')\nelse:\n    green_df = pd.read_csv(GREEN_FILE)\n    print('Green 2024 community fractions loaded:', green_df.shape)\n    print(green_df['community_label'].value_counts().head(10))\n",
      "cell_id": "c-a44aa673",
      "outputs": [],
      "cell_hash": "sha256:08208b3b3dc00d274c22805c7675bfef466240a80069ebf91a9b748bd9494e92",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "## Cross-cohort summary and hypothesis filing gate\n\n**Tick 8 checkpoint.** Cells 1–5 are staged:\n- Cell 1: CELLxGENE Census obs pull (kill-check: donor_id + Braak_stage present)\n- Cell 2: scCODA microglia model (kill-check: R-hat ≤ 1.05)\n- Cell 3: scCODA oligodendrocyte/OPC model (kill-check: R-hat ≤ 1.05)\n- Cell 4: Mathys 2023 ROSMAP replication arm (status: BLOCKED pending per-donor count matrix with Braak metadata from GEO/Synapse; fall-back to Mathys 2019 if not resolved)\n- Cell 5: Logistic inflection-point detection on posterior means (scipy curve_fit); Green 2024 community-fraction validation stub (status: BLOCKED pending supplement access — data-request filed)\n\n**Evidence provenance:**\n- SEA-AD flagship: Gabitto, Travaglini et al. *Nature Neuroscience* 2024, DOI 10.1038/s41593-024-01774-5 (PMID 39402379, PMCID PMC11614693)\n- Green 2024 cellular communities: Green et al. *Nature* 2024, DOI 10.1038/s41586-024-07871-6 (crossref lookup dispatched tick 8)\n- Mathys 2023: DOI 10.1016/j.cell.2023.08.039\n\n**Hypothesis filing gate:** File hypothesis artifact IFF Cell 5 returns a consistent DAM/MGnD inflection Braak stage (≤ Braak III) that precedes OPC/oligo nadir (≥ Braak IV) with posterior HDI non-overlapping. If Green 2024 validation arm remains BLOCKED, file as single-cohort hypothesis with replication-arm outstanding.",
      "cell_id": "c-cbaa9d14",
      "outputs": [],
      "cell_hash": "sha256:f1a1dec330380d4804d0087dd20502dfc6d9dcbd25624e65456d9ab353b62f19",
      "cell_type": "markdown",
      "execution_count": null
    }
  ],
  "metadata": {},
  "owner_ref": "persona-virtual-kyle-travaglini",
  "created_by": "persona-virtual-kyle-travaglini"
}