RNA-seq

Find the genes whose expression differs between conditions from raw RNA-seq read counts, with the negative-binomial model DESeq2 fits, edgeR's quasi-likelihood F-test or limma-voom, then ask which gene sets the differential genes over-represent, which sets sit at the top or bottom of the whole ranking, and which sets change more than the rest of the genes once the correlation between a set's genes is allowed for.

What it doesCopy Anchor Link

Works on a count table (genes in rows with the gene ID first, samples in columns, raw read counts rather than normalized values or TPM) and a sample table (one row per sample, its name as in the count table's header, then the variables the design names). For enrichment it also needs gene sets in GMT format, such as Gene Ontology, MSigDB or WikiPathways exports, whose gene IDs are the same kind as the count table's. Choose them from the project or load them from disk, or fetch them as they stand today: Gene Ontology terms from the GO Consortium's own annotation file for your organism (each term's direct annotations, NOT annotations left out, names from QuickGO), or Reactome pathways for the count table's genes, in the species Reactome finds for the name you give. The fetched sets carry their source, release and date; Reactome's hold only the genes it was sent, so an over-representation background of every annotated gene means the annotated genes among those. Then ask one of four questions. Over-representation and set enrichment read the last differential-expression run's significant, tested and ranked genes, or a gene list or ranked list you supply; the competitive test reads the ranked list too, or fits the count table itself. Its picture opens beside the panel: the volcano, MA plot or dispersion fit, the sample PCA or distance heatmap, the strongest gene sets, or a set's running enrichment score along the ranking.

What this computesCopy Anchor Link

Differential expression (DESeq2; Love, Huber and Anders 2014)

• size factors per sample: the median of ratios to the geometric mean over genes with no zero count, with zeros allowed (poscounts) for sparse data, or iterated: alternating a dispersion fit with the size factors that maximize the negative-binomial likelihood of all but the worst-fitting 5% of genes, for data where every gene has a zero
• normalization factors per gene and sample, when transcript-level quantifications (salmon, kallisto, RSEM through tximport) supply each gene's average transcript length in every sample: each gene's lengths over their geometric mean, scaled by median-of-ratios size factors of the length-corrected counts, so a change in which isoforms a sample expresses does not read as a change in expression; every later step (dispersions, the model, outliers, shrinkage, transformations) uses them in place of the size factors
• dispersion per gene: the gene-wise estimate by Cox–Reid adjusted profile likelihood, a fitted trend against mean count (parametric, asymptotic dispersion plus extra-Poisson over the mean; a local regression of log dispersion on log mean, locfit's local quadratic over the nearest 70% of genes weighted by mean count, which DESeq2 also substitutes on its own when the parametric curve does not fit, and the panel says so; or the mean of the gene-wise estimates), the prior variance of the gene-wise estimates around it (with three or fewer residual degrees of freedom, as in two against two, matched by simulation: DESeq2's seeded draws of the residual distribution at 200 candidate variances, scored by Kullback–Leibler divergence against the observed one and smoothed by loess), and the final maximum a posteriori dispersion shrunk toward the trend; gene-wise outliers keep their own estimate
• glmGamPoi dispersions (Ahlmann-Eltze and Huber 2021), DESeq2's choice for single-cell counts: each gene's dispersion by Cox–Reid adjusted maximum likelihood on a log scale, a local-median trend (a Gaussian-weighted median of the gene-wise estimates over windows of 10% of genes ordered by mean count, read at each gene's nearest neighbour), and quasi-likelihood shrinkage toward it: each gene's quasi-likelihood dispersion, the trend's, the shrunken one and the prior degrees of freedom the shrinkage estimated, all reported; with the likelihood-ratio design this becomes glmGamPoi's quasi-likelihood F test, the deviance gain of the full over the reduced design per removed degree of freedom over the shrunken dispersion. Outlier replacement is off and fitted means are floored at 1e-6, as DESeq2 sets them; observation weights and the variance stabilizing transformation do not combine with it (the regularized log does)
• observation weights per gene and sample, when a weight table comes with the counts (zinbwave's posterior probability that a single-cell count is not a dropout, for example): each gene's weights over their largest weigh its mean, its dispersion likelihoods (a sample at or below 0.01 leaves the Cox–Reid adjustment) and its model fit, and a gene whose weights leave the design degenerate is reported untested; with t-distributed Wald tests each gene's degrees of freedom are its summed weights less the number of coefficients, and a floor on fitted means (minmu) keeps very low single-cell expectations in range
• the negative-binomial generalized linear model of your design formula, its coefficients named by factor and level, with the reference level you name or, if you do not, the one that sorts first; a numeric column enters as a slope
• log2 fold change and its standard error for the last coefficient, a named coefficient, or a contrast of two levels of a factor
• the Wald test of each fold change, against zero or against a threshold with the alternative you choose (greater in absolute value, less in absolute value, greater, less), or the likelihood-ratio test of the full design against a reduced one, testing every removed term at once
• Cook's distances: genes with an outlying count lose their p-value; with enough replicates per group, the outlying counts are replaced and the gene refitted instead
• independent filtering: genes below a mean-count threshold, chosen to maximize discoveries at your target false discovery rate, are left out of the adjustment
• adjusted p-values (Benjamini–Hochberg or any other standard adjustment), and per gene why an adjusted p is missing: all counts zero, a Cook's outlier, below the filter, or observation weights that leave the model degenerate
• log2 fold-change shrinkage, for ranking and plotting genes by effect size: the fold changes pulled toward zero where the counts are too few or too noisy to support them, while the p-values stay those of the unshrunken test. apeglm (Zhu, Ibrahim and Love 2019) puts an adaptive heavy-tailed Cauchy prior on the fold change and reports each gene's posterior mode, its posterior standard deviation and its s-value, the false sign rate of genes at least as confident; ashr (Stephens 2017) fits a zero-centred mixture of normals prior to all genes at once and reports the posterior mean, posterior standard deviation, local false sign rate and s-value; normal is DESeq2's original zero-centred normal prior, its variance matched to the upper tail of the fold changes, with the model refitted under it (Wald test only). A contrast of two levels is shrunk by ashr, or by the normal prior on the model with one column per level, its prior averaged over the factor's levels; apeglm shrinks a coefficient only. With a fold-change threshold, apeglm and ashr report the s-value of a false sign or a fold change inside the threshold, and the normal prior re-tests the shrunken fit against the threshold. The MA plot then draws the shrunken fold changes. Instead of shrinking afterwards, the beta prior (DESeq2's default before version 1.16) fits the normal prior inside the model, one column per level with the prior averaged over each factor's levels, so the reported fold changes, standard errors and Wald test are all the shrunken (maximum a posteriori) ones, and the outlier cutoff counts every level's column
• the summary: genes up and down at your alpha, outliers, low counts and the filter's threshold, genes fitted by the fallback optimizer when the iterative fit did not converge
• the significant, up, down and tested genes, and the ranked list of tested genes by Wald statistic, which the two enrichment questions read

Sample-level quality checks (DESeq2's transformations, before trusting the test)

• a count transformation that makes the variance roughly independent of the mean: the variance stabilizing transformation (VST, from the fitted dispersion trend, either on a subset of 1,000 genes or on every gene), the regularized log (rlog, each gene's counts refitted with a Normal prior on the per-sample effects, with its prior variance), or log2 of normalized count + 1; blind to the design or using the fitted one
• principal component analysis on the most variable genes of the transformed counts (500 by default): every sample's PC1 and PC2 scores, each component's share of the variance, coloured by any sample-table variables, to spot batch effects and outlying replicates
• sample-to-sample Euclidean distances over every gene, and their complete-linkage hierarchical clustering, drawn as a heatmap with its dendrogram

Differential expression (edgeR quasi-likelihood; Lund et al. 2012, Chen et al. 2025)

• filterByExpr: genes kept only where the counts per million reach the level of your minimum count at the median library size in as many samples as the smallest group holds, with a minimum total count; the genes removed have no test
• TMM normalization (Robinson and Oshlack 2010): one normalization factor per library from the trimmed mean of log ratios against a reference sample, and the effective library size
• average log2 counts per million per gene, the abundance every trend is drawn against
• a common negative-binomial dispersion by Cox–Reid adjusted profile likelihood over the most abundant genes, the negative-binomial generalized linear model of your design, and the average quasi-likelihood dispersion
• per gene the quasi-likelihood dispersion from the bias-adjusted deviance and its adjusted residual degrees of freedom, squeezed by empirical Bayes toward an abundance trend: the prior (the trend's value), the posterior and the prior degrees of freedom
• log2 fold change, the quasi-likelihood F-statistic of the coefficient, its p-value on the total degrees of freedom, and the adjusted p
• the ranked list by the signed square root of F, a t-statistic on the same degrees of freedom

Differential expression (limma-voom; Law et al. 2014, Smyth 2004)

• the same filterByExpr and TMM steps, then log2 counts per million per observation
• the mean-variance trend: each gene's square-root residual standard deviation against its log2 count size, smoothed by lowess, which gives every observation a precision weight
• the weighted linear model of your design per gene, its residual variance, the empirical Bayes prior variance and prior degrees of freedom, and the moderated variance
• log2 fold change, average expression, the moderated t-statistic, its p-value, the adjusted p and the B-statistic, the log-odds that the gene is differentially expressed
• the ranked list by moderated t

Over-represented gene sets (one-sided hypergeometric test, or Wallenius' noncentral hypergeometric test when genes carry weights)

• per set: the count of list genes in it, the gene ratio and background ratio, the rich factor, fold enrichment, odds ratio and z-score, the p-value, the adjusted p and the q-value (Storey), with the hit genes
• the background: every gene the experiment tested, or every gene any set annotates; sets are sized within the background and kept between a minimum and maximum size
• which sets pass both the p and adjusted-p cutoff and the q-value cutoff, how many list genes no set annotates, and π₀, the estimated share of true nulls behind the q-values
• with gene weights (a selection bias such as a gene's mean count or length): each set's odds, its members' mean weight over the rest of the background's, which Wallenius' test (Fog 2008) uses in place of equal odds, so a set of well-counted genes is not called only for being easy to detect; how many background genes had no weight and took the median, and how many weights at or below zero were lifted to the smallest positive one; the weighted p-value is one minus the lower tail at a precision of 1e-7, as clusterProfiler computes it, so the strongest sets read p = 0 and the panel says how many

Gene set enrichment along the ranking (preranked GSEA; Subramanian et al. 2005, with the multilevel p-values of Korotkevich et al.)

• the enrichment score: the running sum's largest deviation from zero, each member weighted by its statistic to the power you set (0 is the unweighted Kolmogorov–Smirnov form)
• the normalized enrichment score (NES), against size-matched random sets, and two-tailed or one-tailed scoring
• p-values by adaptive multilevel sampling, accurate far below what plain permutation reaches, down to a floor you set, with the log2 error of each estimate; the adjusted p and q-value
• the rank at maximum, the leading edge: the members before the peak for a positive score or after it for a negative one, its share of the set and of the list, and the core genes
• a seed fixes the random sets, so a run is reproducible
• two implementations, chosen by engine: clusterProfiler's (the default, with q-values) or fgsea's own. They agree in method and differ in their random streams and arithmetic, so the same seed gives each package's own p-values: fgsea's adaptive multilevel estimate, or, with a permutation count, its plain permutation p-values (none below 1/(count + 1)); fgsea adjusts by Benjamini–Hochberg and reports no q-value
• clusterProfiler's two plain permutation tests, chosen by method: sample draws random gene sets of each set's size, permute relabels the whole ranked list; either a fixed number of draws, or adaptive, drawing in doubling batches while a set's p-value stays under a threshold and stopping at a ceiling. The p-value and NES are taken on the observed score's own side of zero, so no p-value falls below 1/(same-side draws + 1)

Gene sets against the rest, correlation-adjusted (camera; Wu and Smyth 2012)

• a competitive test: is a set more differentially expressed than the genes outside it? Genes in a set tend to move together, so their statistics are not independent; the test inflates the variance of the set's mean statistic by the variance inflation factor that the inter-gene correlation implies, instead of counting each gene as fresh evidence
• from the count table: the voom fit (filterByExpr, TMM, precision weights), each gene's own weighted least-squares fit, the moderated t after empirical Bayes squeezing of the residual variances, with its prior degrees of freedom, turned into z-scores
• from a ranked list: the supplied statistic as it stands
• the inter-gene correlation: fixed (0.01 by default) or estimated per set from the genes' residuals, floored at zero unless you allow a negative value
• per set: the number of genes tested, its direction (up or down against the rest), the p-value of a two-sample t on the z-scores, or of a correlation-adjusted rank-sum test, and the false discovery rate (Benjamini–Hochberg); a non-directional test asks only whether the set's genes change more in either direction
• the degrees of freedom of the test, and which sets pass the FDR cutoff

ControlsCopy Anchor Link

Count table, Sample table and Gene sets take the three documents; Fetch sets from (GO or Reactome), Organism (a taxon ID or a name) and Fetch gene sets fill the gene sets from the authority instead. Compute picks the question. For differential expression, Method picks DESeq2, edgeR's quasi-likelihood test or limma-voom, and the controls below follow it. Design is the formula over the sample table's columns with the variable of interest last, Reference names each factor's reference level (condition=untreated), Test picks Wald or likelihood ratio with its Reduced design, and Contrast or Name picks what to report (under edgeR and voom a contrast must compare a level with the factor's reference, and a name may be the model column, conditiontreated). For DESeq2, Transcript lengths (optional) takes tximport's average transcript length table, the same genes and samples as the counts, to normalize by gene and sample, and Observation weights (optional) a weight per gene and sample, with Use t for t-distributed Wald tests and minmu for the floor on fitted means. Advanced holds the size-factor and dispersion-trend types, the beta prior, alpha, the fold-change threshold and its alternative, the adjustment, independent filtering, Cook's filtering and its cutoff, and the replicate count from which outliers are replaced. Transform picks the transformation behind the sample checks (or none), Blind whether it ignores the design, PCA top how many of the most variable genes the PCA uses, and PCA group the sample-table variables that colour it (empty is the factor being compared). Shrink picks apeglm, ashr or normal shrinkage of the reported fold changes, or none; ashr and normal shrink a contrast too, apeglm a coefficient (name the comparison with Reference and Name), and all three honour the fold-change threshold. For edgeR and voom, Filter by expr turns filterByExpr on or off, with its Min count and Min total count. For over-representation, Gene list takes the significant, up or down genes of the last run or a list in the project, Background the run's tested genes or every annotated gene, and Gene weights none, the DESeq2 run's mean normalized counts or a 'gene weight' table in the project; set sizes and cutoffs follow. For set enrichment, Ranked list takes the run's test statistic or a ranked list in the project; the weight exponent, score type, sampling sizes, floor and seed follow, then Engine (clusterProfiler or fgsea), for clusterProfiler Method (multilevel, sample or permute) with its Permutations or Adaptive batches (Min perm, Max perm, Pval threshold), and for fgsea Nperm (0 keeps the multilevel estimate). For the competitive test, Statistic picks the ranked list (the run's test statistic or a file) or the count table's own voom fit, whose design, reference, coefficient and filter controls follow; Use ranks, Inter-gene correlation or Estimate correlation, Directional and the FDR cutoff set the test. Stop ends a run that is taking too long.

Where the numbers come fromCopy Anchor Link

Every gene's size factor, dispersion estimates, coefficients, standard errors, statistics, p-values, adjusted p-values and outlier calls is checked against DESeq2 1.52.0 running in R 4.6.1, on the pasilla knockdown (14,599 genes, seven libraries, the dataset of DESeq2's own vignette) under Wald and likelihood-ratio tests, each trend and size-factor type, thresholds and alternatives, and outlier replacement. Every gene's shrunken fold change, posterior standard deviation and s-value is checked against DESeq2's lfcShrink with apeglm 1.34.0, ashr 2.2-63 and its normal prior, on pasilla and on a 69-sample human lymphoblastoid study (Pickrell et al. 2010), and with contrasts and a threshold on pasilla, Pickrell and the airway dexamethasone study (Himes et al. 2014, four cell lines, eight libraries), whose design, like pasilla's two-against-two subsets, takes the simulated prior. The beta prior's fold changes, standard errors, tests and outlier replacement, and the iterated size factors, are checked gene by gene against DESeq2's own on pasilla, airway and Pickrell. Every transformed value (VST, rlog, log2 + 1, blind or not, either trend), every principal-component score, loading and share of variance, and every sample distance, merge height and leaf order of the clustering is checked against DESeq2's vst, rlog and plotPCA and R's prcomp, dist and hclust on pasilla, airway and Pickrell's 69 libraries. With transcript lengths, every normalization factor and every gene's dispersions, fold changes, tests, outlier replacement, shrunken fold changes and transformed values is checked against DESeq2 on tximport's summary of salmon quantifications of human macrophages, naive and after interferon gamma, from six donors (Alasoo et al. 2018). The local trend, every point of locfit's fitted curve and every gene's dispersions, tests and transformed values on it, is checked against DESeq2 and locfit 1.5-9.12 on pasilla, airway, Pickrell and the macrophages. With observation weights, every gene's mean, dispersions, coefficients, t degrees of freedom, p-values, outlier replacement, shrunken fold changes and transformed values is checked against DESeq2 on 130 single-cell libraries of the developing human cortex (Pollen et al. 2014) with zinbwave 1.34.0's weights, following zinbwave's own recipe; on those cells the parametric trend fails and the substituted local one is checked too. Every tested gene set of over-representation and of set enrichment, including the multilevel and the sample and permute p-values (fixed and adaptive) for a given seed, is checked against enrichit 0.2.6, the engine of clusterProfiler, and with the fgsea engine against fgsea 1.38's fgseaMultilevel and fgseaSimple (two-sided, one-sided, weighted, a single set), in R on pasilla's differential genes against every direct fly Gene Ontology annotation from the GO Consortium (the weighted test, and the Wallenius probabilities beneath it, against BiasedUrn 2.0.12), and the hand-off from differential expression to both enrichment questions is checked as one chain. Every gene's filter call, library size, normalization factor, average log2 counts per million, quasi-likelihood dispersions and their squeeze, fold change, F, p and adjusted p is checked against edgeR 4.10.6, and every gene's voom weight, coefficient, moderated t, p, adjusted p and B against limma 3.68.5, on pasilla with and without the filter and on the Pickrell study. Every gene set's size, correlation, direction, p-value and FDR from the competitive test is checked against limma 3.68.5's camera and cameraPR in R, on pasilla's voom fit with weights, without them, with the correlation estimated per set and as a rank test, and on DESeq2's Wald statistics, every direct fly Gene Ontology set that has a tested gene. Gene sets fetched from the GO Consortium are checked to reproduce those fly sets term for term, from the same annotation file; Reactome's are read as its AnalysisService returns them.