RNA-seq count matrix analysis
Summary
Construct a DESeq2 data object from RNA-seq count matrices, estimate dispersions via negative binomial GLM, and extract differential expression results with automated independent filtering. This skill bridges raw count data to statistically rigorous fold-change and p-value inference.
When to use
You have a count matrix (genes × samples) from HTSeq, featureCounts, or transcript abundance quantification (Salmon, kallisto), a sample metadata table with experimental design, and you need to test for differential expression while controlling false discovery rate via independent filtering on mean normalized counts.
When NOT to use
- Input is already a feature table, normalized expression matrix, or continuous abundance (e.g., TPM, FPKM) — DESeq2 requires raw counts.
- Sample size is very small (n < 3 per group) — dispersion estimation becomes unreliable.
- Data are from a single sample or have no experimental grouping — differential expression testing requires at least two conditions.
Inputs
- count matrix (genes × samples, integer values from HTSeq, featureCounts, or similar)
- sample metadata table with condition/batch covariates
- experimental design formula (e.g., ~cell + dex)
- optional: transcript quantification output (Salmon/kallisto) and tx2gene mapping
Outputs
- DESeqDataSet object with estimated dispersions and fitted negative binomial coefficients
- results table with log2 fold change, standard error, p-value, padj, and baseMean per gene
- summary of independent filtering: number of genes passing filter, rejection threshold
- optional: shrunken log fold changes (via lfcShrink with apeglm type)
How to apply
First, construct a DESeqDataSet object using DESeqDataSetFromMatrix() (for count matrices) or DESeqDataSetFromTximport() (for quantification outputs). Optionally pre-filter low-count genes (e.g., keeping rows with ≥10 reads in ≥smallestGroupSize samples) to reduce memory and computation time. Run DESeq() to estimate dispersions and fit negative binomial GLMs across your design formula. Extract results using results(dds) with default alpha=0.1; the function automatically performs independent filtering based on the mean of normalized counts for each gene, rejecting null hypotheses only for genes above an optimal count threshold. Order results by adjusted p-value (padj) and compare filtered gene counts to document filtering impact. Optionally apply IHW (Independent Hypothesis Weighting) via results(dds, filterFun=ihw) for improved power when filter statistics and test statistics are correlated.
Related tools
- DESeq2 (Core package for constructing DESeqDataSet, estimating dispersions, fitting negative binomial GLMs, and extracting differential expression results with independent filtering) — https://github.com/thelovelab/DESeq2
- tximport (Import transcript quantification output (Salmon, kallisto, RSEM) and convert to gene-level counts)
- tximeta (Import quantification with automatic metadata and produce SummarizedExperiment for DESeqDataSet construction)
- IHW (Alternative to DESeq2's default independent filtering; applies Independent Hypothesis Weighting for improved power)
- HTSeq (Generate count matrices from aligned BAM files (alternative input source))
- featureCounts (Rsubread) (Generate count matrices from aligned BAM files (alternative, faster input source))
- ashr (Adaptive Shrinkage for stabilizing fold-change estimates via empirical Bayes shrinkage) — https://github.com/stephens999/ashr
Examples
library(DESeq2); dds <- DESeqDataSetFromMatrix(countData=cts, colData=coldata, design=~condition); dds <- DESeq(dds); res <- results(dds, alpha=0.1); res_ordered <- res[order(res$pvalue),]
Evaluation signals
- DESeqDataSet object successfully loads with counts(dds) returning integer matrix, design matches formula, and colData matches sample metadata.
- DESeq() completes without error; dispersions are estimated for all genes (rowData(dds)$dispersion is not NA) and show reasonable shrinkage pattern (observed vs. fitted).
- results() table has nrow equal to filtered gene count (< nrow(dds)), with non-null baseMean, log2FoldChange, lfcSE, pvalue, and padj columns; no NAs in padj except for filtered genes.
- Independent filtering threshold is reported (e.g., 'independent filtering removed X genes'); compare retained gene count before/after filtering and verify baseMean distribution is shifted upward post-filter.
- Adjusted p-values (padj) show expected multiple-testing correction (padj ≥ pvalue for all tests); genes with padj < 0.1 (default alpha) represent controlled FDR rejection set.
Limitations
- Dispersion estimation requires multiple replicates per condition; very small sample sizes lead to unreliable estimation and inflated Type I error.
- Independent filtering is automatic but data-driven; threshold varies by dataset and may not align with biological cutoffs (e.g., baseMean threshold is unrelated to fold-change magnitude).
- DESeq2 assumes negative binomial distribution; violations (e.g., zero-inflation, technical confounders) can invalidate inference unless addressed by pre-filtering or design adjustment.
- Log fold-change estimates are unshrunken by default; small-sample noise inflates effect sizes; use lfcShrink() with apeglm type for stability.
- Design formula must account for known confounders (batch, cell type); unmeasured confounding is not detected and biases results.
Evidence
- [other] Core DESeq2 workflow and independent filtering rationale: "The package DESeq2 provides methods to test for differential expression by use of negative binomial generalized linear models"
- [other] Independent filtering mechanism in results() function: "the results function automatically performs independent filtering based on the mean of normalized counts for each gene"
- [other] Pre-filtering low-count genes reduces memory and computation: "by removing rows in which there are very few reads, we reduce the memory size of the
dds data object, and we increase the speed"
- [other] DESeqDataSet construction from count matrix: "dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)"
- [other] Results with log fold changes and p-values: "res <- results(dds, name="condition_trt_vs_untrt")"
- [other] tximport for quantification-to-counts conversion: "txi <- tximport(files, type="salmon", tx2gene=tx2gene)"
- [other] DESeqDataSet from tximport output: "ddsTxi <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)"
1---2name: rnaseq-count-matrix-analysis3description: Use when you have a count matrix (genes × samples) from HTSeq, featureCounts, or transcript abundance quantification (Salmon, kallisto), a sample metadata table with experimental design, and you need to test for differential expression while controlling false discovery rate via independent.4license: CC-BY-4.05---67# RNA-seq count matrix analysis89## Summary1011Construct a DESeq2 data object from RNA-seq count matrices, estimate dispersions via negative binomial GLM, and extract differential expression results with automated independent filtering. This skill bridges raw count data to statistically rigorous fold-change and p-value inference.1213## When to use1415You have a count matrix (genes × samples) from HTSeq, featureCounts, or transcript abundance quantification (Salmon, kallisto), a sample metadata table with experimental design, and you need to test for differential expression while controlling false discovery rate via independent filtering on mean normalized counts.1617## When NOT to use1819- Input is already a feature table, normalized expression matrix, or continuous abundance (e.g., TPM, FPKM) — DESeq2 requires raw counts.20- Sample size is very small (n < 3 per group) — dispersion estimation becomes unreliable.21- Data are from a single sample or have no experimental grouping — differential expression testing requires at least two conditions.2223## Inputs2425- count matrix (genes × samples, integer values from HTSeq, featureCounts, or similar)26- sample metadata table with condition/batch covariates27- experimental design formula (e.g., ~cell + dex)28- optional: transcript quantification output (Salmon/kallisto) and tx2gene mapping2930## Outputs3132- DESeqDataSet object with estimated dispersions and fitted negative binomial coefficients33- results table with log2 fold change, standard error, p-value, padj, and baseMean per gene34- summary of independent filtering: number of genes passing filter, rejection threshold35- optional: shrunken log fold changes (via lfcShrink with apeglm type)3637## How to apply3839First, construct a DESeqDataSet object using DESeqDataSetFromMatrix() (for count matrices) or DESeqDataSetFromTximport() (for quantification outputs). Optionally pre-filter low-count genes (e.g., keeping rows with ≥10 reads in ≥smallestGroupSize samples) to reduce memory and computation time. Run DESeq() to estimate dispersions and fit negative binomial GLMs across your design formula. Extract results using results(dds) with default alpha=0.1; the function automatically performs independent filtering based on the mean of normalized counts for each gene, rejecting null hypotheses only for genes above an optimal count threshold. Order results by adjusted p-value (padj) and compare filtered gene counts to document filtering impact. Optionally apply IHW (Independent Hypothesis Weighting) via results(dds, filterFun=ihw) for improved power when filter statistics and test statistics are correlated.4041## Related tools4243- **DESeq2** (Core package for constructing DESeqDataSet, estimating dispersions, fitting negative binomial GLMs, and extracting differential expression results with independent filtering) — https://github.com/thelovelab/DESeq244- **tximport** (Import transcript quantification output (Salmon, kallisto, RSEM) and convert to gene-level counts)45- **tximeta** (Import quantification with automatic metadata and produce SummarizedExperiment for DESeqDataSet construction)46- **IHW** (Alternative to DESeq2's default independent filtering; applies Independent Hypothesis Weighting for improved power)47- **HTSeq** (Generate count matrices from aligned BAM files (alternative input source))48- **featureCounts (Rsubread)** (Generate count matrices from aligned BAM files (alternative, faster input source))49- **ashr** (Adaptive Shrinkage for stabilizing fold-change estimates via empirical Bayes shrinkage) — https://github.com/stephens999/ashr5051## Examples5253```54library(DESeq2); dds <- DESeqDataSetFromMatrix(countData=cts, colData=coldata, design=~condition); dds <- DESeq(dds); res <- results(dds, alpha=0.1); res_ordered <- res[order(res$pvalue),]55```5657## Evaluation signals5859- DESeqDataSet object successfully loads with counts(dds) returning integer matrix, design matches formula, and colData matches sample metadata.60- DESeq() completes without error; dispersions are estimated for all genes (rowData(dds)$dispersion is not NA) and show reasonable shrinkage pattern (observed vs. fitted).61- results() table has nrow equal to filtered gene count (< nrow(dds)), with non-null baseMean, log2FoldChange, lfcSE, pvalue, and padj columns; no NAs in padj except for filtered genes.62- Independent filtering threshold is reported (e.g., 'independent filtering removed X genes'); compare retained gene count before/after filtering and verify baseMean distribution is shifted upward post-filter.63- Adjusted p-values (padj) show expected multiple-testing correction (padj ≥ pvalue for all tests); genes with padj < 0.1 (default alpha) represent controlled FDR rejection set.6465## Limitations6667- Dispersion estimation requires multiple replicates per condition; very small sample sizes lead to unreliable estimation and inflated Type I error.68- Independent filtering is automatic but data-driven; threshold varies by dataset and may not align with biological cutoffs (e.g., baseMean threshold is unrelated to fold-change magnitude).69- DESeq2 assumes negative binomial distribution; violations (e.g., zero-inflation, technical confounders) can invalidate inference unless addressed by pre-filtering or design adjustment.70- Log fold-change estimates are unshrunken by default; small-sample noise inflates effect sizes; use lfcShrink() with apeglm type for stability.71- Design formula must account for known confounders (batch, cell type); unmeasured confounding is not detected and biases results.7273## Evidence7475- [other] Core DESeq2 workflow and independent filtering rationale: "The package DESeq2 provides methods to test for differential expression by use of negative binomial generalized linear models"76- [other] Independent filtering mechanism in results() function: "the *results* function automatically performs independent filtering based on the mean of normalized counts for each gene"77- [other] Pre-filtering low-count genes reduces memory and computation: "by removing rows in which there are very few reads, we reduce the memory size of the `dds` data object, and we increase the speed"78- [other] DESeqDataSet construction from count matrix: "dds <- DESeqDataSetFromMatrix(countData = cts, colData = coldata, design= ~ batch + condition)"79- [other] Results with log fold changes and p-values: "res <- results(dds, name="condition_trt_vs_untrt")"80- [other] tximport for quantification-to-counts conversion: "txi <- tximport(files, type="salmon", tx2gene=tx2gene)"81- [other] DESeqDataSet from tximport output: "ddsTxi <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)"
HolobiomicsLab