# How to Create a DESeq2 Object from a Count Matrix: A Template By the end of this guide you will be able to take a raw count matrix and its sample metadata, build a DESeqDataSet (the DESeq2 object), verify that the object holds the right data in the right shape, filter out uninformative genes, and confirm that the design formula you intend to fit is actually attached to the object. You will also get a reusable template you can paste into any project, plus a table of the arguments that matter and a list of the errors you are most likely to hit. To follow along you need R (version 4.0 or newer is realistic for current Bioconductor releases), the DESeq2 package from Bioconductor, and two inputs: a count matrix with genes as rows and samples as columns, and a sample table (colData) with one row per sample. Everything in this walkthrough runs in an R console or an R Markdown notebook. If you are building your skills in R for bioinformatics, this single step, converting a plain table of numbers into a formal statistical object, is the hinge between upstream counting pipelines and every downstream analysis you will ever run on bulk RNA-seq data. ## What a Count Matrix Must Look Like A count matrix is a table of integers. Each row is a feature (usually a gene, sometimes a transcript, exon, or contig). Each column is one sequencing library, one sample. The cell value is the number of sequencing reads or fragments assigned to that gene in that sample. In bulk RNA-seq, each gene in each individual is summarized by a single number of this kind, a fragment count, which is exactly what count-based differential expression tools consume [1]. RNA-seq produces this digital form of data, aligned sequence read counts, as its primary measurement, and the comparative [statistical methods](/blog/guides/statistical-methods) are built on that foundation [2]. Before any object creation, four format rules must hold. Violating any of them is the source of most support questions on the Bioconductor support site, so it pays to check them explicitly. 1. **The values are raw, unnormalized counts.** Pass raw integer counts, not TPM, not FPKM, not RPKM, and not counts averaged across replicates. Sequencing assays do not measure absolute transcript abundance. The data are constrained to an arbitrary total, the library size, imposed by the [sequencing depth](/knowledge/diagnostics/molecular/how-to-calculate-sequencing-depth-and-coverage-for-your-ngs-run) of each sample, which is why normalization happens inside the statistical method rather than before it [3]. DESeq2 estimates its own size factors from the raw counts, so pre-normalized input destroys the estimate. 2. **Values are non-negative integers stored as integers or whole-number numerics.** Floating point counts with fractional parts produce a warning at object creation and get rounded. If your counts came out of a transcript-level estimator or a normalization step, check for fractional values first. 3. **Row names are unique gene or feature identifiers.** If a matrix has duplicate row names, the object cannot be built. Collapse or disambiguate them before you start. 4. **Column names exactly match the row names of your sample table.** Not just in content, in order and spelling. This is checked automatically, and a mismatch is the most common creation error of all. One practical caution comes from the upstream side. The counting step you choose, including how you handle reads that map to more than one gene (multireads), feeds directly into the counts you put in this matrix. Dropping multireads underestimates counts, and how you resolve them affects the reliability of downstream differential expression calls [4]. The number of genes called significant can vary widely depending on which combination of alignment, quantification, and counting tools produced your matrix, an effect documented across 495 distinct workflows on real human data [5]. The DESeq2 object does not fix upstream choices, it inherits them. Automated end-to-end pipelines make this explicit by treating read counting as a distinct stage before differential analysis [6]. ## Step 0: Install and Load DESeq2 DESeq2 is distributed through Bioconductor, the open-source software project for genomics analysis in R. If you are following a Bioconductor tutorial for the first time, note that Bioconductor packages are installed through BiocManager, not through plain install.packages(). ```r if (!requireNamespace("BiocManager", quietly = TRUE)) { install.packages("BiocManager") } BiocManager::install("DESeq2") library(DESeq2) packageVersion("DESeq2") ``` Expected output of the last line looks like this (your version may be newer): ``` [7] '1.42.0' ``` If the install completes but library(DESeq2) fails, restart R so the freshly installed package is on the search path. Do not load DESeq2 and other count-based tools such as edgeR in the same session and assume identical behavior, their vocabularies differ even though the underlying statistical framework of modeling counts with gene-specific variability is shared [7]. ## Step 1: Prepare the Count Matrix Here is a self-contained example you can run right now. It simulates a small negative-binomial count matrix for 1,000 genes across 6 samples in two conditions, which is enough to exercise every argument of the object constructor. ```r set.seed(42) genes <- paste0("gene_", sprintf("%04d", 1:1000)) samples <- paste0("sample_", 1:6) cts <- matrix( rnbinom(1000 * 6, mu = 100, size = 10), nrow = 1000, dimnames = list(genes, samples) ) class(cts) dim(cts) cts[1:3, 1:3] ``` Expected output: ``` [7] "matrix" "array" [7] 1000 6 sample_1 sample_2 sample_3 gene_01 107 104 106 gene_02 93 88 103 gene_03 118 131 104 ``` In a real project your counts usually arrive as a CSV or TSV file exported from featureCounts, HTSeq-count, Salmon, or a cloud platform. Read it with the gene identifier in the first column and convert immediately to a matrix: ```r cts <- read.csv("counts.csv", row.names = 1) cts <- as.matrix(cts) mode(cts) <- "integer" stopifnot(!anyNA(cts), all(cts >= 0)) ``` Three checks worth running at this point, because they catch the errors that surface later as confusing failures: ```r anyDuplicated(rownames(cts)) # must return 0 all(cts == round(cts)) # must return TRUE, no fractional counts str(cts[1:5, 1:2]) # integer, not character ``` If the third check reports a character matrix, some gene name contained a comma or the CSV was read with the wrong separator. Fix the read call before proceeding, because a character matrix will be rejected at object creation. ## Step 2: Prepare colData colData is a data frame of sample metadata with one row per sample. Every column you plan to use in your model must live here. Build it from a metadata file or define it directly, and set row names to the sample identifiers: ```r coldata <- data.frame( condition = factor(rep(c("control", "treated"), each = 3)), batch = factor(c("A", "A", "A", "B", "B", "B")), row.names = samples ) coldata ``` Expected output: ``` condition batch sample_1 control A sample_2 control A sample_3 control A sample_4 treated B sample_5 treated B sample_6 treated B ``` Two discipline rules here. First, every design variable must be a factor with a defined reference level, not a character vector. DESeq2 will stop you at creation if you pass characters, and the conversion is a one-liner you control. Second, the row names of colData must match the column names of the count matrix exactly, and the order must match too, because the constructor pairs them positionally. Verify the pairing before you build the object: ```r identical(colnames(cts), rownames(coldata)) ``` This must return TRUE. If it returns FALSE because the orders differ, reorder rather than rename: ```r coldata <- coldata[match(colnames(cts), rownames(coldata)), , drop = FALSE] stopifnot(identical(colnames(cts), rownames(coldata))) ``` A word on what goes into the design. The point of the model formula is to estimate biological variation between samples separately from the treatment effect, so any known source of unwanted variation that is balanced across conditions, such as a processing batch, belongs in the model as a blocking variable [7]. Repeated measurements on the same experimental unit over time introduce within-unit correlation that standard count models do not handle, so if your design involves multiple time points per subject, plan the analysis approach accordingly rather than folding the time points in as if they were independent samples [8]. ## Step 3: Create the DESeq2 Object This is the template. One function call does the work: ```r dds <- DESeqDataSetFromMatrix( countData = cts, colData = coldata, design = ~ batch + condition ) dds ``` Expected output: ``` class: DESeqDataSet dim: 1000 6 metadata(1): version assays(1): counts rownames(1000): gene_0001 gene_0001 ... gene_0999 gene_1000 rowData names(0): colnames(6): sample_1 sample_2 ... sample_5 sample_6 colData names(2): condition batch ``` (The rownames line above is abbreviated for readability. On a freshly simulated matrix every gene name is distinct.) The printed object is a DESeqDataSet, an S4 class that extends the SummarizedExperiment class. It bundles three things: the count assay, the sample metadata in colData, and the design formula. From here on, every subsetting, transformation, and statistical test operates on all three simultaneously, which is what prevents the classic failure mode of a model fitted against the wrong sample labels. ### Each Argument Explained ![Comparison of three DESeq2 design choices with what each models, when to use it, and what results reports.](/article-images/illustrations/create-deseq2-object-design-formulas.webp) | Argument | What it does | Typical value | Notes | |--|--|--|--| | countData | The raw count matrix or data frame, genes by samples | cts | Values must be non-negative integers. Data frames are converted internally. | | colData | Data frame of sample metadata, samples by variables | coldata | Row names must equal colnames(countData). Design variables must be factors. | | design | Model formula built from colData columns | ~ batch + condition | Determines the comparisons results() will report. Variables must exist in colData. | | tidy | Treat countData as a table with feature IDs in the first column | FALSE (default) | Set TRUE to read a wide CSV directly, for example from read.csv with tidy input. | | rowRanges | Optional genomic coordinates for each feature | NULL or a GRanges object | Useful when counts came from transcript quantification tools that supply gene coordinates. | | ignoreRank | Internal check on the design matrix rank | FALSE (default) | Leave at the default unless you understand the rank check. | The design formula deserves its own paragraph, since it is the argument with the longest downstream shadow. It must reference columns of colData by name, it defines which coefficients appear when you call results(), and it should include known nuisance factors before the variable of interest. The convention is that the last term in the formula is the comparison of interest, so ~ batch + condition places condition as the effect you extract. For a two-group comparison, either ~ condition or the no-intercept shorthand ~ 0 + condition works, but the intercept form keeps results() behavior predictable. The rationale for modeling blocking variables alongside the treatment effect, so that biological variation is estimated separately from technical structure, comes from the generalized linear model framework for multifactor count experiments [7]. One more detail on factor encoding. If you do not set levels explicitly, R sorts them alphabetically and the default level becomes the reference. To force a biological reference: ```r dds$condition <- relevel(dds$condition, ref = "control") levels(dds$condition) ``` Expected output: ``` [7] "control" "treated" ``` Log2 fold changes reported later will then be treated relative to control, which is what readers expect. ## Step 4: Verify the Object Never trust that the object is right because the call did not error. Run four checks, each of which has caught real mistakes in real projects. First, confirm the counts inside the object match your original matrix: ```r identical(counts(dds), cts) # [7] TRUE ``` If this returns FALSE, the sample order inside the object differs from your input order, most often because colData was reordered to match. Confirm the pairing is still correct with head(colData(dds)) against head(cts). Second, confirm the metadata landed: ```r colData(dds) ``` Expected output: ``` DataFrame with 6 rows and 2 columns condition batch sample_1 control A sample_2 control A sample_3 control A sample_4 treated B sample_5 treated B sample_6 treated B ``` Third, confirm the design is attached and its coefficient names are what you expect: ```r design(dds) # ~batch + condition ``` Fourth, sanity check the counts themselves before any modeling: ```r summary(colSums(counts(dds))) ``` Library sizes in the same experiment typically differ, and that difference is normal. What matters is that DESeq2 corrects for it through size factors rather than you pre-normalizing, because scaling normalization of raw counts demonstrably improves differential expression inference compared with naive comparisons [9]. Watch also for one sample whose total is wildly out of line with the rest. A library total a hundredfold below its peers usually signals failed sequencing, not biology, and no object construction can rescue it. ## Step 5: Pre-filter Low-Count Genes Filtering is not required, but it saves memory and fitting time and it modestly sharpens multiple-testing correction. The standard practice is a simple total-count threshold applied after object creation, so that the object and the filter travel together: ```r keep <- rowSums(counts(dds)) >= 10 dds <- dds[keep, ] dim(dds) ``` Expected output depends on your data. On the simulated matrix above, most genes have mean counts near 100, so nearly all 1,000 survive. On real data, the drop is often severe, sometimes keeping fewer than half the rows in low-depth experiments. Two cautions on filtering. A fixed threshold like 10 total counts works across many datasets, but the threshold should be chosen with your library sizes in mind, not copied blindly. More fundamentally, filtering decisions are analysis choices that feed the final result, and diagnostic inspection of count distributions before and after filtering is part of a quality-aware analysis rather than an optional garnish [10]. Keep the filter simple, document it, and apply it identically to all samples. ## Step 6: What Happens Next Inside the Object Once the DESeq2 object exists, the standard pipeline is one function call, dds <- DESeq(dds), which runs three internal steps in sequence: size factor estimation, dispersion estimation per gene, and fitting of the negative binomial generalized linear model per gene. The size factor step is where the raw-count requirement pays off. DESeq2 computes a scaling factor per sample from the ratios of each count to a per-gene geometric mean, so samples with larger library sizes get scaled down to a common footing before any test is run [9]. You rarely need to run those steps manually. After dds <- DESeq(dds), check the model with resultsNames(dds), then extract results with results(dds, contrast = c("condition", "treated", "control")). The object you built in Step 3 is what makes that call produce the right comparison, against the right reference level, with the right covariates adjusted. A note for teams working in Python. The same statistical workflow exists as PyDESeq2, an open-source Python implementation of the DESeq2 pipeline that accepts count data through its own API and is designed for integration with Python data science tooling [11]. The object you create there is not the S4 class described here, but the inputs, raw integer counts plus a properly encoded design, are the same, and everything in this article about count matrix format applies unchanged. ## Parameters That Matter, At a Glance The table in Step 3 covers the constructor. Three downstream parameters also affect what your object ultimately yields, and setting them consciously avoids silent surprises. | Decision | Where it is set | Common choice | Effect | |--|--|--|--| | Reference level | relevel or factor(levels) before or after creation | The biologically baseline condition | Flips the sign of log2 fold changes in results() | | Contrast | results(dds, contrast = ...) | c("condition", "treated", "control") | Selects the comparison reported, independent of level ordering | | Filter threshold | rowSums(counts(dds)) >= threshold | 10 total counts across all samples | Removes genes with too little data to inform dispersion | ## Common Errors and Their Fixes These are the failures you will actually encounter, with the exact text where DESeq2 or base R supplies it. **Error: "some variables in design formula are characters, convert them to factors first: condition".** Your design column is a character vector. Fix it before creation: ```r coldata$condition <- factor(coldata$condition) ``` **Error: "counts matrix should be numeric mode, non-negative values".** Your matrix contains negative numbers, missing values encoded as text, or character strings where numbers belong. Check with anyNA(cts), any(cts < 0), and str(cts[1:5, 1:2]). A count of zero is valid, NA is not. **Error: "object of type 'closure' is not subsettable".** You typed counts[1:5, 1:2] instead of counts(dds)[1:5, 1:2]. counts is a function, not a variable, and base R's error message reflects that. **An error stating that row names of colData and column names of the count matrix do not match.** This is the SummarizedExperiment validity check doing its job. The fix is the match() reordering shown in Step 2. Do not fix it by renaming columns, fix it by reordering rows of colData. **A warning that some rows have non-integer values and will be rounded.** Your counts contain fractions, usually because an earlier step estimated or normalized them. Trace back to the counting step and read in raw counts instead. Rounding at this stage hides a data provenance problem rather than solving it. **An error stating that the model matrix is not full rank.** Your design formula asks for more coefficients than your data can support. The classic cause is a variable that is perfectly confounded with the condition, for example a batch column where every treated sample is in batch B and every control is in batch A. No statistical method can separate those effects, so the design must change, typically by dropping the confounded variable. **Error: duplicate row names.** The constructor refuses duplicate gene identifiers. If duplicates are genuinely the same gene measured twice, collapse them deliberately, for example by summing: ```r cts_collapsed <- rowsum(cts, group = rownames(cts)) ``` If they are different features that share a name, make them unique with make.unique(rownames(cts)) and document that choice. ## Common Mistakes and Limitations **Feeding normalized or averaged values into the constructor.** TPM, FPKM, counts divided by library size, or counts averaged across replicate lanes all break the size factor estimation that Step 6 depends on. Raw counts only. If your collaborator sends a normalized matrix, ask for the raw table. **Forgetting that the choice of count matrix already made decisions for you.** Benchmarking work comparing many combinations of alignment, quantification, and testing tools found wide variation in the number of significant genes called, measured as differences in precision and recall against known truths [5]. Creating a DESeq2 object is a small step in code but a large commitment in methodology. **Building the object before checking sample metadata.** A swap of two column labels in a spreadsheet silently produces a perfectly valid object with the wrong pairing of conditions to libraries. The verification checks in Step 4 exist because of this exact failure. Library size summaries and, where possible, a PCA of variance-stabilized counts, catch most label swaps quickly. **Treating sample size as a formality.** The number of biological replicates has a larger effect on differential expression power than library size does, with the notable exception of low-expressed genes where both matter comparably [12]. Three replicates per group can run without errors and still underpower the analysis badly. Plan replicate numbers before sequencing, not after object creation. **Assuming reproducibility from any single analysis.** Work on tumor transcriptomics found that poor reproducibility of differentially expressed gene lists occurs even at relatively large sample sizes, driven largely by high variation of gene expression between samples of the same condition, including the influence of outlier counts [13]. A clean object and a clean run are necessary, not sufficient, for results that generalize. **Skipping basic quality control before modeling.** Diagnostic plots of count distributions, sample distances, and normalization behavior should inform preprocessing decisions, and quality-aware tools exist precisely to make this inspection routine rather than optional [10]. Check the data, then create the object, then run the model. **Rounding fractional counts without asking why.** A warning about non-integer rows almost always means your matrix is not raw counts. The honest fix is upstream, not a silent rounding. ## Next Steps With a verified DESeqDataSet in hand, the canonical sequence is dds <- DESeq(dds), then resultsNames(dds) to confirm coefficients, then res <- results(dds, contrast = c("condition", "treated", "control")) with an adjusted p-value cutoff chosen for your experiment, then visualization with plotMA and PCA of transformed counts. From there, annotate gene identifiers, check whether top hits make biological sense, and consider whether pathway-level summaries serve your question better than long gene lists. For a multi-condition, multi-factor experiment, revisit the design formula against the blocking principles discussed above before adding interaction terms. ## Frequently Asked Questions ### Can I use TPM or FPKM values instead of raw counts? No. DESeq2 requires raw integer counts because it estimates per-sample size factors and per-gene dispersions from the raw data, and pre-normalized values break that estimation. If you only have TPM values, go back to the counting step and export raw counts. ### Should I filter genes before or after creating the DESeq2 object? Either order produces the same statistical result, but the standard workflow filters after creation using rowSums(counts(dds)) >= 10 or a similar threshold. Filtering after creation keeps the object and its gene filter in one place and is easier to document. ### What should the design formula be for paired samples? Include the pairing variable before the condition, for example ~ patient + condition. This models the subject-to-subject variation as a block, which is the point of blocking variables in multifactor count models. Failing to include the pairing term inflates the residual variability and loses power. ### How do I handle duplicate gene names in my count matrix? Collapse duplicates that are truly the same gene by summing their counts with rowsum(), or make identifiers unique with make.unique() if the rows represent distinct features. Do not leave duplicates in place, because object creation will fail. ### Can I create a DESeq2 object from a data frame instead of a matrix? Yes. DESeqDataSetFromMatrix accepts a data frame for countData and converts it to a matrix internally. Make sure the first column is not accidentally treated as a data column, which is what the tidy argument is for when gene IDs sit in column one. ### What does the tidy argument do in DESeqDataSetFromMatrix? tidy = TRUE tells the constructor that countData is a table with feature identifiers in the first column, which should be moved into row names. The default is tidy = FALSE, which expects a matrix with row names already set. ### Can I run DESeq2 in Python? Yes. PyDESeq2 is an open-source Python implementation of the DESeq2 workflow for bulk [RNA-seq differential expression](/knowledge/bioinformatics/rna-seq-differential-expression-deseq2-edger), released under the MIT license. The count matrix and design requirements are the same as in R. ### Do I need the same number of replicates in each group? No, unbalanced designs are fine, but replicate number drives [statistical power](/blog/guides/statistical-power-what-it-is-and-why-it-matters-in-research) more than library size does, so plan for adequate replicates in every group you intend to compare. Very small groups will produce wide [confidence intervals](/blog/data-analysis/confidence-interval-formula-examples) regardless of sequencing depth. ## Sources 1. [Individual Level Differential Expression Analysis for Single Cell RNA-seq data](https://doi.org/10.1101/2021.05.10.443350) 2. [NPEBseq: nonparametric empirical bayesian-based procedure for differential expression analysis of RNA-seq data](https://pubmed.ncbi.nlm.nih.gov/23981227/) 3. [Benchmarking differential expression analysis tools for RNA-Seq: normalization-based vs. log-ratio transformation-based methods](https://pubmed.ncbi.nlm.nih.gov/30021534/) 4. [A fuzzy method for RNA-Seq differential expression analysis in presence of multireads](https://pubmed.ncbi.nlm.nih.gov/28185579/) 5. [Empirical assessment of analysis workflows for differential expression analysis of human samples using RNA-Seq](https://pubmed.ncbi.nlm.nih.gov/28095772/) 6. [SeqExpressionAnalyser: An R Package for Automated End-to-End RNA-Seq Analysis From Reads to Differential Expression](https://pubmed.ncbi.nlm.nih.gov/41602550/) 7. [Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation](https://pubmed.ncbi.nlm.nih.gov/22287627/) 8. [rmRNAseq: differential expression analysis for repeated-measures RNA-seq data](https://pubmed.ncbi.nlm.nih.gov/32449749/) 9. [A scaling normalization method for differential expression analysis of RNA-seq data](https://pubmed.ncbi.nlm.nih.gov/20196867/) 10. [Data quality aware analysis of differential expression in RNA-seq with NOISeq R/Bioc package](https://pubmed.ncbi.nlm.nih.gov/26184878/) 11. [PyDESeq2: a python package for bulk RNA-seq differential expression analysis](https://pubmed.ncbi.nlm.nih.gov/37669147/) 12. [Optimization of an RNA-Seq Differential Gene Expression Analysis Depending on Biological Replicate Number and Library Size](https://pubmed.ncbi.nlm.nih.gov/29491871/) 13. [High heterogeneity undermines generalization of differential expression results in RNA-Seq analysis](https://pubmed.ncbi.nlm.nih.gov/33509298/)