Extracts, filters, annotates, and exports differential expression results from DESeq2 or edgeR with proper handling of padj=NA (independent filtering, Cook's outliers, all-zero), multiple-testing correction choice (BH vs Storey q-value vs IHW vs lfsr), TREAT vs post-hoc fold-change filtering, p-value histogram diagnostics, gene annotation via org.db/biomaRt/mygene, GSEA preranked input, ORA background construction, replication reality (Schurch 2016 small-n result), and SABV/sex-stratified reporting. Use when extracting and interpreting DE results, troubleshooting padj=NA, choosing FDR method, preparing ranked lists for pathway analysis, annotating gene IDs, or comparing DESeq2 vs edgeR outputs.
Standardmäßig ist der Prompt ausgewählt, der zuerst die Quelle prüft. Sie können zu einem direkten Befehl wechseln oder eine lokale Kopie herunterladen.
Quelldateien prüfen
Lesen Sie SKILL.md und alle von SkillsMP angezeigten Begleitdateien, bevor Sie sich für eine Installation entscheiden.
Mit Codex oder Claude installieren Kopieren Sie diesen Prompt, fügen Sie ihn in Codex, Claude oder einen anderen Assistant ein und lassen Sie die Skill-Seite prüfen und installieren.
Ein direkter Befehl überspringt den Prüf-Prompt. Prüfen Sie die Quelle, bevor Sie ihn ausführen.
Extracts, filters, annotates, and exports differential expression results from DESeq2 or edgeR with proper handling of padj=NA (independent filtering, Cook's outliers, all-zero), multiple-testing correction choice (BH vs Storey q-value vs IHW vs lfsr), TREAT vs post-hoc fold-change filtering, p-value histogram diagnostics, gene annotation via org.db/biomaRt/mygene, GSEA preranked input, ORA background construction, replication reality (Schurch 2016 small-n result), and SABV/sex-stratified reporting. Use when extracting and interpreting DE results, troubleshooting padj=NA, choosing FDR method, preparing ranked lists for pathway analysis, annotating gene IDs, or comparing DESeq2 vs edgeR outputs.
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, NApadj, baseMean below auto threshold
Removed before BH adjustment to maximize rejections at alpha
results(dds, independentFiltering = FALSE) OR filterFun = ihw
Cook's distance outlier
NApvalue, NApadj, 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
NApvalue 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.
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.
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)
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.
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.
ashr lfsr (local false sign rate -- probability the estimated direction is wrong):
res_ashr <- lfcShrink(dds, coef ='condition_treated_vs_control',
type ='ashr', svalue =TRUE)
res_ashr$svalue # FDR-like, based on lfsr; requires svalue=TRUE
svalue=TRUE is required to populate the svalue column; the default returns the standard pvalue/padj columns only. lfsr and padj are NOT interchangeable. lfsr asks "P(sign wrong)"; padj asks "expected fraction of false discoveries". When reporting, state which.
P-value Histogram Diagnostics
Goal: Diagnose model misspecification, hidden batch effects, or over-correction by inspecting the raw p-value distribution.
Approach: Plot raw p-values; under a correctly specified null, the histogram is uniform with an upward spike near zero (the true DE genes).
library(ggplot2)
ggplot(res_df, aes(x = pvalue))+
geom_histogram(bins =50, fill ='steelblue', color ='white')+
labs(x ='P-value', y ='Frequency', title ='P-value distribution')+
theme_bw()
Inspect PCA for batch; add covariate; check plotDispEsts
Conservative (depleted near 0, spike near 1)
Over-correction; too many covariates; wrong dispersion
Simplify model; check dispersion plot for excess shrinkage
Spike only at p = 1
Discrete artifact from very-low-count genes
Pre-filter more aggressively
Bimodal with spike at 0.5
Unusual; suggests a discrete categorical test masquerading
Investigate
The histogram is one of the cheapest sanity checks in a DE pipeline; always plot it before believing the gene list.
Filtering and Ordering
Goal: Subset to significant genes and rank by p-value, fold change, or expression level for downstream use.
Approach: dplyr-style filter + arrange; handle NA padj explicitly per the three-meanings table at the top.
sig <- res_df %>%
filter(!is.na(padj), padj <0.05,abs(log2FoldChange)>1, baseMean >10)%>%
arrange(padj)# Up- vs down-regulated
up <- sig %>% filter(log2FoldChange >0)
down <- sig %>% filter(log2FoldChange <0)# Summary
n_tested <-sum(!is.na(res$padj))
n_sig <-sum(res$padj <0.05, na.rm =TRUE)
cat(sprintf('Tested: %d Significant (padj<0.05): %d Up: %d Down: %d\n',
n_tested, n_sig,sum(sig$log2FoldChange >0),sum(sig$log2FoldChange <0)))
Gene Annotation
Goal: Map gene IDs to symbols, descriptions, and cross-database identifiers for human-readable results.
Approach: Prefer AnnotationDbi::mapIds with org.db (fast, local, version-pinned); fall back to biomaRt or mygene for symbols/aliases not in org.db; for prokaryotes, use Prokka/Bakta GFF.
The sub('\\..*', '', ...) strips the Ensembl version. CAUTION: this regex destroys the _PAR_Y suffix in GENCODE 25-43 PAR genes -- use sub('\\.[0-9]+(_PAR_Y)?$', '\\1', ...) to preserve. See expression-matrix/gene-id-mapping for full details.
For HGNC symbols changed since 2020 (SEPT1 -> SEPTIN1, MARCH1 -> MARCHF1, MARC1 -> MTARC1, DEC1 -> DELEC1) old symbol-keyed downstream tools silently drop genes. Always join on stable Ensembl or Entrez IDs; use symbols as display labels only.
Goal: Produce a ranked list of all genes (no significance filter) for fgsea / clusterProfiler GSEA.
Approach: Rank by Wald statistic (DESeq2 stat) or shrunken LFC. NEVER use a filtered set as GSEA input -- GSEA's permutation null requires the full background.
stat (Wald Z) is preferred over raw LFC for GSEA because it combines effect and precision in one number. Unshrunken LFC is dominated by low-count noise.
ORA Input
Goal: Run over-representation analysis (enrichGO, enrichKEGG) on a significant gene list with the correct background.
Approach: Subset to padj<0.05; background = ALL TESTED genes (post-independent-filtering), NOT the genome.
Common mistake: omitting universe= lets clusterProfiler default to "all annotated genes for this organism" -- which includes thousands of genes never tested. The resulting enrichment p-values are wrong (too small). The background MUST be the tested set.
Concordance >70% at the top 500: robust. <60%: suspect filtering, normalization, or design difference -- not a tool difference. Run both pipelines with the same filtering and design to isolate.
Per-Method Failure Modes
Dropped a key gene by removing NAs
Trigger: Pipeline does res_df <- na.omit(res_df); downstream gene of interest is missing from results.
Mechanism: Gene was flagged by independent filtering OR Cook's distance; padj is NA but the biology is real.
Symptom: A gene with clear differential expression in the count matrix is absent from the results table.
Fix: Diagnose which filter fired (independent filtering vs Cook's vs all-zero); rerun results() with the appropriate filter off (independentFiltering = FALSE or cooksCutoff = FALSE).
Reported FDR on a magnitude-filtered gene set
Trigger: Methods section says "genes with padj < 0.05 and abs(LFC) > 1 (FDR < 5%)".
Mechanism: BH controls FDR for the |LFC| > 0 hypothesis, not the |LFC| > 1 hypothesis. The post-hoc filter adds no FDR control.
Symptom: Reviewer challenges the FDR claim; replication studies show many of the filtered genes are not the magnitude expected.
Fix: Use TREAT (glmTreat) or lfcThreshold= to test the magnitude hypothesis with proper FDR control. Re-do the methods sentence to match what was actually computed.
ORA universe wrong
Trigger: ORA p-values look implausibly small for a small significant gene set.
Mechanism:universe= argument omitted; clusterProfiler defaulted to all annotated genes in the organism, including thousands never in the tested set.
Symptom: Many enriched pathways at strict thresholds; results don't replicate; reviewer questions the background.
Fix: Explicitly pass universe = bg_entrez where bg_entrez is the set of tested gene IDs (i.e., those with non-NA padj).
"Significant" gene list in n=3 study doesn't replicate
Trigger: Small RNA-seq study finds 200 DE genes; validation in independent cohort recovers 60.
Mechanism: Schurch 2016 RNA 22:839: at n=3/group, all tools miss 20-40% of true positives compared to n=30. Variability of the gene list itself is high.
Symptom: 30-50% replication of the gene list across independent runs of the SAME data.
Fix: Frame the small-n DE list as hypothesis-generating, not as a stable set of facts. Validate top hits orthogonally before drawing conclusions. Use TREAT for biologically meaningful thresholds to require larger effects.
Sex-confounded design gives spurious chrX/chrY signal
Trigger: Mixed-sex cohort; sex not in the design; many chrY genes call as DE.
Mechanism: Sex distribution differs across the experimental groups; the "treatment effect" partially captures sex.
Symptom: chrY genes (DDX3Y, RPS4Y1, UTY) and XIST dominate the top DE list.
Fix: Include sex in the design (~ sex + condition); rerun. For chrX/chrY-specific analyses, sex MUST be in the model or the analysis is uninterpretable. Mauvais-Jarvis et al. 2020 Lancet 396:565 reviews the SABV requirement.
Common errors
Error / symptom
Cause
Fix
$FDR not found on DESeq2 result
DESeq2 uses padj; edgeR uses FDR
Check tool, use correct column
summary(res) shows different cutoff than results(alpha=)
summary(res, alpha=) defaults to 0.1
Pass alpha explicitly to summary()
All padj NA
All genes filtered (rare; usually a data problem)
Check independentFilteringResults(res); inspect baseMean distribution
Direction of LFC reversed
Reference level not set; alphabetical default
relevel() BEFORE DESeq()
Gene symbol mapping rate <50%
Mixed Ensembl versions; recent HGNC renames
Verify Ensembl release, check for SEPT/MARCH/MARC renames
enrichGO reports thousands of pathways
Wrong universe=
Pass universe = bg_entrez (tested set, not genome)
References
Benjamini Y, Hochberg Y. 1995. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B 57(1):289-300. doi:10.1111/j.2517-6161.1995.tb02031.x
Storey JD. 2003. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat 31(6):2013-2035. doi:10.1214/aos/1074290335
Ignatiadis N, Klaus B, Zaugg JB, Huber W. 2016. Data-driven hypothesis weighting increases detection power in genome-scale multiple testing. Nat Methods 13(7):577-580. doi:10.1038/nmeth.3885
Bourgon R, Gentleman R, Huber W. 2010. Independent filtering increases detection power for high-throughput experiments. PNAS 107(21):9546-9551. doi:10.1073/pnas.0914005107
Stephens M. 2017. False discovery rates: a new deal. Biostatistics 18(2):275-294. doi:10.1093/biostatistics/kxw041
McCarthy DJ, Smyth GK. 2009. Testing significance relative to a fold-change threshold is a TREAT. Bioinformatics 25(6):765-771. doi:10.1093/bioinformatics/btp053
Schurch NJ et al. 2016. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA 22(6):839-851. doi:10.1261/rna.053959.115
Love MI, Huber W, Anders S. 2014. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol 15(12):550. doi:10.1186/s13059-014-0550-8
Mauvais-Jarvis F et al. 2020. Sex and gender: modifiers of health, disease, and medicine. Lancet 396(10250):565-582. doi:10.1016/S0140-6736(20)31561-0
Bruford EA et al. 2020. Guidelines for human gene nomenclature. Nat Genet 52:754-758. doi:10.1038/s41588-020-0669-3
Ziemann M, Eren Y, El-Osta A. 2016. Gene name errors are widespread in the scientific literature. Genome Biol 17:177. doi:10.1186/s13059-016-1044-7
Wu T et al. 2021. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. The Innovation 2(3):100141. doi:10.1016/j.xinn.2021.100141
