# === BLOCK 3: limma-voom differential expression model === # Per-subset pseudobulk limma-voom in R (pydeseq2 as Python fallback) # # R model formula (primary): # ~ age_continuous + sex + CMV_status + batch # Random effect: donor nested in cohort (edgeR::voomWithQualityWeights if multi-cohort) # # Contrast of interest: # 1. age_continuous coefficient: DE genes tracking age within each CD8 subset # 2. Interaction: age × CMV_status — genes where aging effect differs by CMV # 3. Cross-subset comparison: Tex-term vs TEMRA age-DE overlap # # R pseudocode: # library(limma); library(edgeR) # for subset in ['Tex_prog','Tex_int','Tex_term','TEMRA','Tcm','Tnaive']: # pb <- read_h5ad(paste0('pseudobulk_cd8/pb_',subset,'.h5ad')) # via anndata2ri # dge <- DGEList(counts=t(pb$X), samples=pb$obs) # dge <- filterByExpr(dge, min.count=10, min.total.count=15) # dge <- calcNormFactors(dge) # design <- model.matrix(~ age + sex + CMV_status + batch, data=dge$samples) # v <- voom(dge, design, plot=FALSE) # fit <- lmFit(v, design) # fit <- eBayes(fit) # # Extract age coefficient # res <- topTable(fit, coef='age', number=Inf, adjust='BH') # res$subset <- subset # write.csv(res, paste0('limma_results/de_age_', subset, '.csv')) # # Python DESeq2 fallback (if R unavailable): # from pydeseq2.dds import DeseqDataSet # from pydeseq2.ds import DeseqStats # for subset, pb in pb_frames.items(): # dds = DeseqDataSet(counts=..., metadata=pb.obs, # design_factors=['age_decade','sex','CMV_status']) # dds.deseq2() # stat_res = DeseqStats(dds, contrast=['age_decade','old','young']) # stat_res.summary() # # Output schema per subset CSV: # gene, logFC, AveExpr, t, P.Value, adj.P.Val (BH), B, subset # Significance filter: |logFC| > 0.5, adj.P.Val < 0.05 # Key genes to annotate: TOX, PDCD1, HAVCR2, LAG3, TIGIT, TCF7, GZMK, GZMB, # PRF1, KLRG1, B3GAT1 (CD57), IL7R, CCR7, SELL, MKI67, PCNA

Discussion