Skip to content
Bench to Insights
Blog

DESeq2 and pyDESeq2

DESeq2 is the most widely used package for differential expression analysis of RNA-seq data. It uses a negative binomial model to test whether gene expression differs between experimental conditions. pyDESeq2 is a Python reimplementation of the same algorithm, from Muzellec and colleagues at Owkin, and it exists because a pipeline already written in Python should not have to cross into R for one step.

Both tabs below run the same analysis on the airway dataset, in the same container, from the same count matrix. That makes the comparison at the end of the page a measurement rather than a claim.

The dataset, and getting it into Python

Section titled “The dataset, and getting it into Python”

The airway dataset comes from Himes et al. 2014. The experiment treated human airway smooth muscle cells with dexamethasone, a corticosteroid used to treat asthma. Four cell lines were each split into treated and untreated conditions, giving eight samples.

It ships as a Bioconductor package, which is the first problem a Python user hits. The R tab therefore writes the count matrix and the sample table to CSV, and the Python tab reads those files. pyDESeq2 wants samples as rows and genes as columns, the transpose of the R layout, so the export transposes on the way out.

library(airway)
library(DESeq2)


data("airway")
airway


counts_matrix <- assay(airway, "counts")
dim(counts_matrix)


# Hand the same matrix to Python. pyDESeq2 wants samples as rows, so transpose.
dir.create("outputs", showWarnings = FALSE, recursive = TRUE)
write.csv(t(counts_matrix), "outputs/airway_counts.csv")
write.csv(
  as.data.frame(colData(airway))[, c("cell", "dex")],
  "outputs/airway_metadata.csv"
)
class: RangedSummarizedExperiment
dim: 63677 8
metadata(1): ''
assays(1): counts
rownames(63677): ENSG00000000003 ENSG00000000005 ... ENSG00000273492
  ENSG00000273493
rowData names(10): gene_id gene_name ... seq_coord_system symbol
colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
colData names(9): SampleName cell ... Sample BioSample
[1] 63677     8

Building the object and filtering

Section titled “Building the object and filtering”

Both libraries hold the counts, the sample table and the design in one object. The design formula is written the same way in both: ~ cell + dex models expression as a function of the donor cell line and the treatment, so the coefficient of interest is the treatment effect after the donor is accounted for.

The filter keeps genes with at least ten reads in total. It removes over 41,000 genes that carry no information about a difference between conditions, which speeds up the fit and lightens the multiple testing correction.

dds <- DESeqDataSetFromMatrix(
  countData = counts_matrix,
  colData = colData(airway),
  design = ~ cell + dex
)


keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
nrow(dds)
[1] 22369

Running the model

Section titled “Running the model”

One call does the whole pipeline in both languages: size factor estimation to correct for sequencing depth, dispersion estimation to model the overdispersion that count data always has, a negative binomial fit per gene, and a Wald test on the coefficient.

dds <- DESeq(dds)
res <- results(dds, contrast = c("dex", "trt", "untrt"))


res_df <- as.data.frame(res)
res_df <- res_df[order(res_df$padj), ]
round(head(res_df[, c("baseMean", "log2FoldChange", "stat", "padj")], 5), 4)
                  baseMean log2FoldChange    stat padj
ENSG00000152583   997.4447         4.5750 24.8314    0
ENSG00000165995   495.0957         3.2911 24.7353    0
ENSG00000120129  3409.0384         2.9478 24.1872    0
ENSG00000101347 12703.4128         3.7670 24.1488    0
ENSG00000189221  2341.7807         3.3537 23.5872    0

Counting what changed

Section titled “Counting what changed”

The adjusted p-value is the column to filter on, because every gene was tested and a raw p-value of 0.01 is expected by chance in hundreds of them.

tested <- res_df[!is.na(res_df$padj), ]
significant <- tested[tested$padj < 0.05, ]


# Tested, significant, up, down.
nrow(tested)
nrow(significant)
sum(significant$log2FoldChange > 0)
sum(significant$log2FoldChange < 0)


# Keep the R numbers for the comparison at the end of the page.
write.csv(
  data.frame(
    engine = "R DESeq2",
    tested = nrow(tested),
    significant = nrow(significant),
    up = sum(significant$log2FoldChange > 0),
    down = sum(significant$log2FoldChange < 0),
    top_gene = rownames(res_df)[1]
  ),
  "outputs/r_summary.csv",
  row.names = FALSE
)


# The per-gene R estimates, so the Python tab can compare the two directly.
write.csv(
  res_df[, c("baseMean", "log2FoldChange", "padj")],
  "outputs/r_results.csv"
)
[1] 18032
[1] 4000
[1] 2193
[1] 1807

Where the two agree, and where they do not

Section titled “Where the two agree, and where they do not”

The comparison is the reason for running both. The table is built in the Python tab from the R summary written above and from the Python results in memory, so every cell in it comes from this page’s own run.

r_summary = pd.read_csv("outputs/r_summary.csv")


python_row = {
    "engine": "pyDESeq2",
    "tested": len(tested),
    "significant": len(significant),
    "up": int((significant["log2FoldChange"] > 0).sum()),
    "down": int((significant["log2FoldChange"] < 0).sum()),
    "top_gene": results_df.index[0],
}
comparison = pd.concat(
    [r_summary, pd.DataFrame([python_row])], ignore_index=True
)
print(comparison.to_string(index=False))


# How far apart are the two fold change estimates for the genes both tested?
r_res = pd.read_csv("outputs/r_results.csv", index_col=0)
shared = r_res.join(results_df, how="inner", lsuffix="_r", rsuffix="_py")
delta = (shared["log2FoldChange_r"] - shared["log2FoldChange_py"]).abs()
spearman = shared["log2FoldChange_r"].corr(
    shared["log2FoldChange_py"], method="spearman"
)


print(f"shared genes: {len(shared)}")
print(f"max abs log2FC difference: {delta.max():.4f}")
print(f"median abs log2FC difference: {delta.median():.6f}")
print(f"Spearman correlation: {spearman:.6f}")
  engine  tested  significant   up  down        top_gene
R DESeq2   18032         4000 2193  1807 ENSG00000152583
pyDESeq2   18032         3922 2156  1766 ENSG00000152583
shared genes: 22369
max abs log2FC difference: 0.7969
median abs log2FC difference: 0.000208
Spearman correlation: 0.999778

The two implementations do not return identical numbers, and they are not supposed to. They differ in the numerical details of dispersion estimation and in how the Wald p-values are computed, which moves genes near the significance boundary across it in one engine and not the other. What matters is whether the difference changes a conclusion, and the fold change agreement above answers that.

Use the R implementation when the rest of your analysis is Bioconductor, because the downstream annotation and gene set packages expect its objects. Use pyDESeq2 when the rest of the pipeline is Python and the alternative is exporting to R and back for one step.