Differential expression from RNA-seq counts

What this involves

Differential expression analysis quantifies changes in transcript abundance between experimental conditions while accounting for biological and technical variability across sequencing libraries. To obtain defensible results, negative binomial generalised linear models are fitted directly to raw integer counts produced by tools such as featureCounts or pipelines such as nf-core/rnaseq. Quantifications normalised to TPM, FPKM, or CPM cannot be used, as the dispersion model requires the original integer scale.

Valid transcriptomic modelling requires exploratory checks before testing. Library sizes, sample-to-sample correlations, and principal component analysis on variance-stabilising transformed (VST) counts identify sample outliers and unmodelled batch structure. Technical covariates must be integrated into the design formula rather than pre-corrected. Empirical Bayes shrinkage algorithms such as apeglm or ashr are then applied to log fold changes, preventing genes with minimal read counts from dominating hit lists.

Gene identifier stability is maintained across every step by joining on versionless Ensembl gene identifiers or HGNC identifiers with MANE Select transcripts, rather than volatile gene symbols. The agent runs these workflows—including pydeseq2 analysis and nf-core/differentialabundance pipelines—on dedicated execution environments and returns structured summary tables containing gene identifiers, base means, shrunken fold changes, and Benjamini–Hochberg adjusted q-values.

How it goes

1. 1Inspect sample structure and qualityEvaluate library sizes, sample correlations, and PCA plots of variance-stabilising transformed (VST) counts to detect batch effects and verify replicate concordance before fitting models.
2. 2Specify experimental design in pydeseq2Instantiate a DeseqDataSet with raw integer counts and sample metadata, incorporating known batch factors and conditions directly into the design formula.
3. 3Fit dispersion parameters and calculate test statisticsEstimate size factors and gene-wise dispersions using pydeseq2, then extract condition contrasts with DeseqStats to generate Wald test statistics and unadjusted p-values.
4. 4Apply log fold change shrinkageShrink log fold changes using apeglm or ashr implementations to moderate noise-driven, inflated effect sizes in low-count transcripts.
5. 5Control false discovery rates and map identifiersCalculate Benjamini–Hochberg adjusted p-values (q-values) or test against explicit fold-change thresholds, joining output rows to Ensembl gene identifiers and HGNC symbols.

Start it here

Ask in plain language. The agent works on a Linux machine of its own — it installs and runs what the job needs, keeps every file, and reads the results back with the numbers explained. Nothing to set up, and a long run continues after you close the page.

What to watch for

• Never supply TPM, FPKM, CPM, or normalised values to the model, as negative binomial dispersion fitting requires raw integer counts.
• Do not pre-correct counts with algorithms like ComBat prior to testing; batch terms must be included in the design formula to avoid deflating p-values.
• If principal component analysis reveals that batch is the dominant driver of variance, stop and update the design formula before conducting hypothesis tests.
• Do not impute adjusted p-values as 1 or silently drop rows where padj is NA, as these represent informative independent filtering for low expression.
• Never apply post hoc filtering on both unadjusted p-values and fold changes simultaneously; test directly against a threshold using lfcThreshold to retain FDR control.
• Always join tables using Ensembl gene identifiers or HGNC identifiers rather than gene symbols, which drift and can be reallocated.

What it runs

Questions