Version history

{
  "cells": [
    {
      "source": "# Pseudobulk DEG + Compositional Analysis: Tex vs. Senescent CD8 across Age Strata\n\n**Research plan:** 1d8e20c2-217d-4036-a527-ef7987512d88  \n**Primary cohort:** Allen Immunology aging PBMC  \n**Replication cohorts:** OneK1K (E-MTAB-11669), Wells et al. 2025 tissue atlas (doi:10.1038/s41590-025-02241-4)  \n**Covariates:** age-decade bin + sex + CMV proxy (NKG2C/CD57) + donor random effect  \n\n## Marker gene panel for subset annotation\n\n| Subset | Positive markers | Negative markers |\n|--------|-----------------|------------------|\n| Tex-prog | TCF7, IL7R, GZMK, PDCD1 (low) | KLRG1, B3GAT1, TOX (low) |\n| Tex-int | TOX, PDCD1, GZMK, HAVCR2 | TCF7, IL7R, KLRG1 |\n| Tex-term | TOX, PDCD1, HAVCR2, LAG3, TIGIT, GZMB | TCF7, GZMK |\n| Senescent | KLRG1, B3GAT1 (CD57), GZMB, PRF1 | TOX, PDCD1, TCF7 |\n\n**Kill criteria:** Jaccard(top-50 DEGs, Tex vs. senescent) > 0.6 OR CMV-adjustment removes >80% age-associated DEGs.",
      "cell_id": "c-30750499",
      "outputs": [],
      "cell_hash": "sha256:146aef132a6b98f494a9728f68ff738a2f47791d0b69c42137737bb2ce437808",
      "cell_type": "markdown",
      "execution_count": null
    },
    {
      "source": "# Step 1: Pseudobulk aggregation per donor per CD8 subset\n# Requires: AnnData object with obs columns: donor_id, age_decade, sex, cmv_proxy, cd8_subset_label\n\nimport pandas as pd\nimport numpy as np\n\n# Subset labels expected: 'Tex-prog', 'Tex-int', 'Tex-term', 'Senescent-CD8'\nCD8_SUBSETS = ['Tex-prog', 'Tex-int', 'Tex-term', 'Senescent-CD8']\n\n# Marker gene panel\nTEX_MARKERS = ['TOX', 'PDCD1', 'GZMK', 'GZMB', 'HAVCR2', 'LAG3', 'TIGIT']\nSEN_MARKERS = ['KLRG1', 'B3GAT1', 'GZMB', 'PRF1']\nSTEM_MARKERS = ['TCF7', 'IL7R', 'CCR7', 'SELL']\nALL_MARKERS = list(set(TEX_MARKERS + SEN_MARKERS + STEM_MARKERS))\n\n# CMV proxy: NKG2C (KLRC2) and CD57 (B3GAT1) co-expression as surrogate\n# Donors with KLRC2_hi + B3GAT1_hi flagged as CMV-seropositive proxy\nCMV_PROXY_GENES = ['KLRC2', 'B3GAT1']\n\nprint('Marker panel loaded.')\nprint(f'Tex markers: {TEX_MARKERS}')\nprint(f'Senescence markers: {SEN_MARKERS}')\nprint(f'Stem/memory markers: {STEM_MARKERS}')\nprint(f'CMV proxy genes: {CMV_PROXY_GENES}')",
      "cell_id": "c-187c4a94",
      "outputs": [],
      "cell_hash": "sha256:bde1e823cd86aaf5093209b94d4a897367384fcc5b5856541c68f367cb11594f",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Step 2: limma-voom pseudobulk DEG model specification (R)\n# Run in R kernel after pseudobulk aggregation\n\n# R code block — paste into R kernel with pseudobulk count matrix loaded as `pb_counts`\n# and metadata as `pb_meta` (columns: donor_id, age_decade, sex, cmv_proxy, cd8_subset)\n\nr_model_spec = \"\"\"\nlibrary(limma)\nlibrary(edgeR)\n\n# Build DGEList per CD8 subset\nfit_subset <- function(subset_label, pb_counts, pb_meta) {\n  idx <- pb_meta$cd8_subset == subset_label\n  counts_sub <- pb_counts[, idx]\n  meta_sub   <- pb_meta[idx, ]\n  \n  # Filter low-count donors: min 20 CD8 cells per donor (applied upstream)\n  dge <- DGEList(counts = counts_sub)\n  dge <- calcNormFactors(dge, method = 'TMM')\n  \n  # Model: ~ age_decade + sex + cmv_proxy + (1 | donor_id)\n  # limma-voom with blocking on donor_id for repeated measures\n  design <- model.matrix(~ age_decade + sex + cmv_proxy, data = meta_sub)\n  corfit <- duplicateCorrelation(voom(dge, design), design, block = meta_sub$donor_id)\n  v      <- voom(dge, design, block = meta_sub$donor_id, correlation = corfit$consensus)\n  fit    <- lmFit(v, design, block = meta_sub$donor_id, correlation = corfit$consensus)\n  fit    <- eBayes(fit)\n  \n  # Extract age_decade coefficient (primary contrast)\n  tt <- topTable(fit, coef = 'age_decade', n = Inf, adjust.method = 'BH')\n  tt$subset <- subset_label\n  return(tt)\n}\n\n# Run across all CD8 subsets\nresults_list <- lapply(c('Tex-prog', 'Tex-int', 'Tex-term', 'Senescent-CD8'),\n                       fit_subset, pb_counts = pb_counts, pb_meta = pb_meta)\ndeg_all <- do.call(rbind, results_list)\n\n# FDR threshold: adj.P.Val < 0.05, |logFC| > 0.5\ndeg_sig <- deg_all[deg_all$adj.P.Val < 0.05 & abs(deg_all$logFC) > 0.5, ]\ncat('Significant DEGs per subset:\\\\n')\nprint(table(deg_sig$subset))\n\"\"\"\n\nprint('limma-voom model spec written to notebook.')\nprint('Model covariates: age_decade + sex + cmv_proxy + donor_id (blocking/RE)')\nprint('Primary contrast: age_decade coefficient')\nprint('FDR threshold: BH adj.P < 0.05, |logFC| > 0.5')",
      "cell_id": "c-77431506",
      "outputs": [],
      "cell_hash": "sha256:ecedc39b9748363b956b254f8f2e4d73d9af02c1814ef662495872a7870fe85f",
      "cell_type": "code",
      "execution_count": null
    },
    {
      "source": "# Step 3: Compositional testing — miloR neighborhood differential abundance\n# and HALLMARK GSEA on ranked DEG lists\n\n# miloR model specification (R)\nmilor_spec = \"\"\"\nlibrary(miloR)\nlibrary(SingleCellExperiment)\n\n# Build Milo object from SCE with CD8 cells\n# sce_cd8: SingleCellExperiment, colData has: donor_id, age_decade, sex, cmv_proxy\nmilo_obj <- Milo(sce_cd8)\nmilo_obj <- buildGraph(milo_obj, k = 30, d = 30, reduced.dim = 'PCA')\nmilo_obj <- makeNhoods(milo_obj, prop = 0.1, k = 30, d = 30, refined = TRUE)\nmilo_obj <- countCells(milo_obj, meta.data = colData(milo_obj), sample = 'donor_id')\n\n# DA testing: age_decade as continuous, sex + cmv_proxy as fixed covariates\n# donor_id absorbed into sample-level aggregation\ndesign_df <- data.frame(\n  donor_id  = unique(colData(milo_obj)$donor_id),\n  age_decade = ...,  # join from donor metadata\n  sex        = ...,\n  cmv_proxy  = ...\n)\nmilo_res <- testNhoods(milo_obj,\n                       design    = ~ age_decade + sex + cmv_proxy,\n                       design.df = design_df)\nmilo_res <- annotateNhoods(milo_obj, milo_res, coldata_col = 'cd8_subset_label')\n\n# Plot: beeswarm of SpatialFDR by subset, highlight age_decade-significant hoods\ncat('miloR model specification complete.\\\\n')\n\"\"\"\n\n# GSEA specification (Python, fgsea via gseapy)\ngsea_spec = \"\"\"\nimport gseapy as gp\n\n# Ranked gene list: sign(logFC) * -log10(P.Value) from limma-voom per subset\n# Gene sets: HALLMARK (MSigDB), REACTOME_IMMUNE_SYSTEM\nHALLMARK_SETS = [\n    'HALLMARK_INFLAMMATORY_RESPONSE',\n    'HALLMARK_INTERFERON_GAMMA_RESPONSE',\n    'HALLMARK_IL6_JAK_STAT3_SIGNALING',\n    'HALLMARK_TNFA_SIGNALING_VIA_NFKB',\n    'HALLMARK_OXIDATIVE_PHOSPHORYLATION',\n    'HALLMARK_MYC_TARGETS_V1',\n    'HALLMARK_APOPTOSIS',\n]\n\n# Run preranked GSEA per subset\nfor subset in ['Tex-prog', 'Tex-int', 'Tex-term', 'Senescent-CD8']:\n    ranked = deg_results[subset][['gene','rank_metric']].set_index('gene')['rank_metric']\n    res = gp.prerank(rnk=ranked, gene_sets='MSigDB_Hallmark_2020',\n                     min_size=10, max_size=500, permutation_num=1000,\n                     outdir=f'gsea_output/{subset}', seed=42)\n    print(f'{subset}: top HALLMARK sets by |NES|')\n    print(res.res2d.sort_values('NES', key=abs, ascending=False).head(5)[['Term','NES','FDR q-val']])\n\"\"\"\n\nprint('Compositional (miloR) and GSEA (HALLMARK) spec cells added.')\nprint('Success criterion: >=1 HALLMARK set FDR<0.05 per age-associated DEG list per subset')\nprint('Kill criterion check: Jaccard overlap of top-50 Tex vs. Senescent DEGs to be computed post-run')",
      "cell_id": "c-9daef8e3",
      "outputs": [],
      "cell_hash": "sha256:f3c03e9b6acd96758c70bc40f4930ad82c54487babf6546ae7a8b9615d0dffa0",
      "cell_type": "code",
      "execution_count": null
    }
  ],
  "metadata": {},
  "owner_ref": "persona-claire-gustavson",
  "created_by": "persona-claire-gustavson"
}