bio-differential-expression-de-results

What this skill does

## Version Compatibility

Reference examples tested with: DESeq2 1.42+, edgeR 4.0+, IHW 1.34+, qvalue 2.34+, ashr 2.2+, AnnotationDbi 1.66+, org.Hs.eg.db 3.18+, biomaRt 2.58+, mygene 1.38+ (Python), dplyr 1.1+, openxlsx 4.2+

Before using code patterns, verify installed versions match. If versions differ:
- R: `packageVersion('<pkg>')` then `?function_name` to verify parameters
- Python: `pip show <package>` then `help(module.function)` to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.

# DE Results

**"What are my significant genes?"** -> Extract DE estimates and p-values from the fitted model, handle missing padj correctly, apply FDR control appropriate to the design, and produce the table or ranked list the downstream tool actually needs.

## The Single Most Important Modern Insight -- `padj = NA` has three distinct meanings

A `NA` in the `padj` column is not a missing value; it is a flag indicating which filter excluded the gene. The three causes -- independent filtering, Cook's distance outlier, and all-zero in a group -- have completely different remediations. Dropping all NA rows blindly silently discards real signal, most often from low-count master regulators (transcription factors expressed at ~10 counts) that pass biology but fail the data-driven baseMean threshold.

| `padj = NA` cause | DESeq2 detection | What it means | Fix if undesired |
|-------------------|------------------|---------------|-------------------|
| Independent filtering | finite `pvalue`, `NA` `padj`, baseMean below auto threshold | Removed before BH adjustment to maximize rejections at `alpha` | `results(dds, independentFiltering = FALSE)` OR `filterFun = ihw` |
| Cook's distance outlier | `NA` `pvalue`, `NA` `padj`, baseMean > 0, group has >=3 reps | One sample has Cook's > `qf(0.99, p, m-p)` | `results(dds, cooksCutoff = FALSE)` |
| All-zero or near-zero in a group | `NA` `pvalue` AND baseMean very low | Insufficient information to test | Filter at preprocess time; or accept |

Independent filtering (Bourgon, Gentleman, Huber 2010 *PNAS* 107:9546) chooses the baseMean threshold to maximize rejections. The filter MUST be independent of the test statistic under the null -- this is why baseMean (the across-sample mean) is the canonical choice. Using "min count in treatment group" as a filter VIOLATES the independence requirement and inflates type-I error. Most pipelines unknowingly do this; do not.

A second axis: at n>=7 per group, `DESeq()` also REPLACES outlier counts via `replaceOutliers()` and refits (default `minReplicatesForReplace = 7`). Cook's filtering is NOT computed for continuous covariates -- a continuous-covariate analysis has effectively no outlier filtering.

## Algorithmic Taxonomy

| Method | What it computes | When to use | Failure mode |
|--------|------------------|-------------|--------------|
| BH (`p.adjust(method='BH')`, DESeq2 default `pAdjustMethod='BH'`) | FDR at fixed alpha; Benjamini-Hochberg 1995 | Default for most RNA-seq DE | Assumes independence or PRDS; many overlapping tests violate |
| Storey q-value (`qvalue::qvalue`) | q-value using estimated pi_0 | Genome-scale with many true nulls | Pi_0 estimation can fail at small test counts |
| IHW (`results(filterFun=ihw)`, Ignatiadis 2016 *Nat Methods* 13:577) | Weighted BH with covariate-informed weights | Modern default for DESeq2; +5-20% discoveries at same FDR | Covariate MUST be independent under null |
| ashr local false sign rate (`lfsr` / `svalue`, with `svalue=TRUE`) | P(sign of estimate is wrong) | When effect-direction certainty is what matters | Conservative lower bound on FDR; not interchangeable with padj |
| BY (`p.adjust(method='BY')`) | Benjamini-Yekutieli; arbitrary dependence | Strongly correlated tests | Uniformly conservative; rarely needed for DE |
| Holm / Bonferroni | FWER | Small confirmatory test sets | Far too conservative for genome-scale |
| TREAT / `lfcThreshold=` | FDR for "|LFC| > tau" hypothesis | Pre-specified biologically meaningful threshold | tau must be set BEFORE looking at data |

## Decision Tree by Scenario

| Scenario | Recommended approach | Why |
|----------|---------------------|-----|
| Standard bulk DE, two groups | DESeq2 results with default BH; report padj < 0.05 | Default works |
| Want more power at same FDR | `results(dds, filterFun = ihw)` | IHW typically gains 5-20% |
| Pre-specified biological fold-change matters | `results(dds, lfcThreshold = log2(1.5), altHypothesis = 'greaterAbs')` OR `glmTreat(fit, lfc = log2(1.5))` | Post-hoc `padj<0.05 & abs(LFC)>1` does NOT control FDR for the magnitude claim |
| Ranking for GSEA preranked | `stat` (Wald Z) for DESeq2 OR shrunken LFC | Never use unshrunken LFC -- low-count noise dominates |
| ORA input | Subset by padj<0.05; background = ALL TESTED genes (post-independent-filtering) | Background = "all genes in genome" is wrong; pre-filtering already excluded many |
| Many NA padj including biologically interesting genes | Diagnose: independent filtering vs Cook's vs all-zero; turn off the offending filter only for that gene set | Blanket `na.omit` discards signal |
| Multi-condition design | LRT for "any change" first; pairwise per-level Wald for effect sizes | LRT padj is omnibus; LRT LFC is one specific coefficient |
| Small n (<=3/group) | Report as exploratory, top hits only | Schurch 2016: tools miss 20-40% of true positives at n=3 |
| Human / mouse with mixed sexes | Include sex as covariate; run sex-stratified sensitivity | SABV mandate; sex effect is real and chromosomal |
| Prokaryotic | Use Prokka/Bakta GFF, KEGG strain code | Ensembl/org.db are eukaryote-only |

## Extracting Results

**Goal:** Pull DE estimates and p-values from a fitted DESeq2 or edgeR object into a usable data frame with explicit contrast naming.

**Approach:** `results()` (DESeq2) or `topTags()` (edgeR) with explicit `name=` / `coef=`; convert to data.frame; preserve row order if planning to join with annotation.

```r
library(DESeq2)
library(dplyr)

resultsNames(dds)
res <- results(dds, name = 'condition_treated_vs_control', alpha = 0.05)
res_shrunk <- lfcShrink(dds, coef = 'condition_treated_vs_control', type = 'apeglm')

res_df <- as.data.frame(res)
res_df$gene <- rownames(res_df)
```

```r
library(edgeR)
tt <- topTags(qlf, n = Inf, sort.by = 'none')$table
tt$gene <- rownames(tt)
```

`sort.by = 'none'` in `topTags` preserves the original gene order -- critical when joining with an annotation table by row index. Default is sort by p-value.

Column name reminder (a recurring cross-tool bug):

| Tool | LFC column | Adjusted p-value column |
|------|------------|------------------------|
| DESeq2 | `log2FoldChange` | `padj` |
| edgeR | `logFC` | `FDR` |
| limma `topTable` | `logFC` | `adj.P.Val` |
| limma `topTreat` | `logFC` | `adj.P.Val` (post-TREAT) |

## TREAT vs Post-hoc LFC Filtering

**Goal:** Make a defensible FDR claim about "biologically meaningful fold change" genes.

**Approach:** Use TREAT or `lfcThreshold=` to test a magnitude hypothesis with proper FDR control. Post-hoc filtering of `padj<0.05 & abs(LFC)>tau` does NOT control FDR for the magnitude claim.

```r
res_treat <- results(dds, lfcThreshold = log2(1.5), altHypothesis = 'greaterAbs', alpha = 0.05)

# edgeR equivalent
tr <- glmTreat(fit, coef = 2, lfc = log2(1.5))
```

What a reviewer is really probing with a "200 genes >2x changed at FDR 5%" claim: is the FDR for the change>2x claim or for the change-non-zero claim? Post-hoc filtering controls FDR only for the latter. TREAT (or `lfcThreshold=`) controls FDR for the former. McCarthy & Smyth 2009 *Bioinformatics* 25:765 is the canonical citation.

## IHW for Better Power

**Goal:** Gain 5-20% more discoveries at the same FDR by weighting p-values with a covariate (typically baseMean) that informs power but is independent of the null.

**Approach:** `results(d