scverse · pankev-in · Jul 24, 2026 · Jul 24, 2026
@@ -381,6 +381,31 @@
 ds_Y_vs_X = DeseqStats(dds, contrast=["group", "Y", "X"], inference=inference)
 ds_Y_vs_X.summary()
 
+# %%
+# .. _lrt_ref:
+#
+# Likelihood ratio test
+# """""""""""""""""""""
+#
+# Instead of the Wald test, we may assess significance with a likelihood ratio test
+# (LRT), as in R DESeq2's ``DESeq(dds, test="LRT", reduced=...)``. The LRT compares the
+# full model to a nested ``reduced`` model that drops the term(s) being tested, and is
+# especially useful to test a factor with several levels jointly, or to compare nested
+# models in general. Here we test the effect of ``condition`` on top of ``group`` by
+# comparing the full ``~group + condition`` model to a reduced ``~group`` model.
+#
+# The p-value reflects the model comparison, while the reported ``log2FoldChange`` still
+# corresponds to the requested ``contrast`` of the full model.
+
+ds_lrt = DeseqStats(
+    dds,
+    contrast=["condition", "B", "A"],
+    test="LRT",
+    reduced="~group",
+    inference=inference,
+)
+ds_lrt.summary()
+
 # %%
 # LFC shrinkage (multifactor)
 # """""""""""""""""""""""""""

@@ -12,3 +12,18 @@ in `/datasets/synthetic/`, respectively using `~condition` and `~condition + gro
 - `r_test_size_factors.csv` contains DESeq2's `estimateSizeFactors` output,
 - `r_vst.csv` contains DESeq2's `varianceStabilizingTransformation` output with `blind=TRUE` and `fitType="parametric"`,
 - `r_vst_with_design.csv` contains DESeq2's `varianceStabilizingTransformation` output with `blind=FALSE` and `fitType="parametric"`.
+
+### Likelihood ratio test (LRT) reference data
+
+The `r_test_res_lrt*.csv` files contain the output of DESeq2's `results()` after running
+`DESeq(dds, test="LRT", reduced=...)`, and are used by the LRT tests. They were generated
+with DESeq2 1.46.0 by `generate_lrt_reference.R` (run `Rscript tests/data/generate_lrt_reference.R .`
+from the repository root to regenerate them). For each case, a `*_no_independent_filtering.csv`
+variant is also stored (i.e. `results(..., independentFiltering=FALSE)`), which isolates the raw
+LRT statistic and p-value from the p-value adjustment step.
+
+- `single_factor/r_test_res_lrt.csv`: full `~condition` vs reduced `~1` (`df=1`),
+- `multi_factor/r_test_res_lrt.csv`: full `~group + condition` vs reduced `~group` (`df=1`),
+- `multi_factor/r_test_res_lrt_reduced_intercept.csv`: full `~group + condition` vs reduced `~1` (`df=2`),
+- `multi_factor/r_test_res_lrt_outliers.csv`: same as above (`~group + condition` vs `~group`) but with
+  injected Cooks outliers that trigger outlier replacement and model refitting.
@@ -0,0 +1,66 @@
+## Generate the R DESeq2 reference outputs used by the PyDESeq2 likelihood ratio
+## test (LRT) suite.
+##
+## These files are the ground truth that `tests/test_pydeseq2.py` compares the
+## PyDESeq2 LRT implementation against, mirroring how the existing Wald-test
+## reference data was produced. They use DESeq2 defaults (parametric dispersion
+## trend, median-of-ratios size factors, Cooks refitting), matching the defaults
+## of `DeseqDataSet`/`DeseqStats`.
+##
+## Usage (from the repository root):
+##     Rscript tests/data/generate_lrt_reference.R .
+##
+## Requires: DESeq2 (reference data was generated with DESeq2 1.46.0).
+
+suppressMessages(library(DESeq2))
+
+repo <- commandArgs(trailingOnly = TRUE)[1]
+if (is.na(repo)) repo <- "."
+
+counts   <- read.csv(file.path(repo, "datasets/synthetic/test_counts.csv"),   row.names = 1)  # genes x samples
+metadata <- read.csv(file.path(repo, "datasets/synthetic/test_metadata.csv"), row.names = 1)  # samples x factors
+
+## Run an LRT and write both the default `results()` output (independent
+## filtering ON) and a version with independent filtering OFF (which isolates the
+## raw LRT statistic and p-value from the p-value adjustment step).
+run_lrt <- function(counts, metadata, full, reduced, name, outfile) {
+  dds <- DESeqDataSetFromMatrix(countData = counts, colData = metadata, design = full)
+  dds <- DESeq(dds, test = "LRT", reduced = reduced, quiet = TRUE)
+  write.csv(as.data.frame(results(dds, name = name)), outfile)
+  write.csv(
+    as.data.frame(results(dds, name = name, independentFiltering = FALSE)),
+    sub("\\.csv$", "_no_independent_filtering.csv", outfile)
+  )
+  invisible(dds)
+}
+
+## --- Clean data (no outliers) -------------------------------------------------
+meta <- metadata
+meta$condition <- factor(meta$condition, levels = c("A", "B"))
+meta$group     <- factor(meta$group,     levels = c("X", "Y"))
+
+# Case A: single factor, full = ~condition, reduced = ~1 (df = 1)
+run_lrt(counts, meta, ~condition, ~1, "condition_B_vs_A",
+        file.path(repo, "tests/data/single_factor/r_test_res_lrt.csv"))
+
+# Case B: multi factor, full = ~group + condition, reduced = ~group (df = 1)
+run_lrt(counts, meta, ~group + condition, ~group, "condition_B_vs_A",
+        file.path(repo, "tests/data/multi_factor/r_test_res_lrt.csv"))
+
+# Case C: multi factor, full = ~group + condition, reduced = ~1 (df = 2)
+run_lrt(counts, meta, ~group + condition, ~1, "condition_B_vs_A",
+        file.path(repo, "tests/data/multi_factor/r_test_res_lrt_reduced_intercept.csv"))
+
+## --- Data with Cooks outliers (triggers refitting) ---------------------------
+## Same mild outlier injection as the existing Wald outlier test.
+counts_out <- counts
+counts_out["gene1", "sample1"]  <- 2000L
+counts_out["gene7", "sample11"] <- 1000L
+
+dds_out <- DESeqDataSetFromMatrix(countData = counts_out, colData = meta,
+                                  design = ~group + condition)
+dds_out <- DESeq(dds_out, test = "LRT", reduced = ~group, quiet = TRUE)
+write.csv(as.data.frame(results(dds_out, name = "condition_B_vs_A")),
+          file.path(repo, "tests/data/multi_factor/r_test_res_lrt_outliers.csv"))
+
+cat("Done. n replaced (outlier case):", sum(mcols(dds_out)$replace, na.rm = TRUE), "\n")
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,6.31579426514429,0.0119667128198756,0.0299167820496891
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,12.6913815611018,0.000367344758633349,0.00122448252877783
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,5.38895826761268,0.0202645651531322,0.0405291303062643
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,15.8216712349465,6.96007574745871e-05,0.000619867084233537
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,14.7312636318204,0.000123973416846707,0.000619867084233537
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,0.0039664276938538,0.949782763590195,0.949782763590195
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,0.813353148799138,0.36713074461821,0.407923049575788
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,4.25158626412713,0.0392136864907029,0.0653561441511715
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,2.57942335791518,0.108261094491313,0.141449313668176
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,2.50953313686284,0.113159450934541,0.141449313668176
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,6.31579426514429,0.0119667128198756,0.0299167820496891
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,12.6913815611018,0.000367344758633349,0.00122448252877783
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,5.38895826761268,0.0202645651531322,0.0405291303062643
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,15.8216712349465,6.96007574745871e-05,0.000619867084233537
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,14.7312636318204,0.000123973416846707,0.000619867084233537
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,0.0039664276938538,0.949782763590195,0.949782763590195
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,0.813353148799138,0.36713074461821,0.407923049575788
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,4.25158626412713,0.0392136864907029,0.0653561441511715
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,2.57942335791518,0.108261094491313,0.141449313668176
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,2.50953313686284,0.113159450934541,0.141449313668176
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.55034480088269,0.760136902014561,0.285643269329525,6.8291022584383,0.00896843310671871,0.0224210827667968
+"gene2",21.3782230626323,0.53636303894464,0.152196255941874,12.3502395356451,0.000440929894170484,0.00146976631390161
+"gene3",5.02959691052961,-0.674061326257139,0.287293062831826,5.39886648506354,0.02014983406417,0.04029966812834
+"gene4",101.066398705215,-0.424994832954691,0.110699439617865,14.6717658138381,0.000127948581466862,0.000911728459446636
+"gene5",27.3066986731155,0.588386623184279,0.156762787257669,14.0047881763476,0.000182345691889327,0.000911728459446636
+"gene6",5.446988786503,-0.019685001214121,0.309645507173678,0.00399799057129258,0.949583622558906,0.949583622558906
+"gene7",28.4098070748334,0.13780930645167,0.152926774840609,0.809310627127047,0.36832414249836,0.409249047220401
+"gene8",40.5694778066866,-0.270864009796395,0.135240979433579,3.9986201979707,0.0455375283850031,0.0758958806416718
+"gene9",37.3644308033897,-0.212456964616852,0.136834847771096,2.40207465321987,0.121174453507509,0.151468066884386
+"gene10",11.6436906897954,0.386870211256403,0.245746587087819,2.45919232559606,0.116838846893463,0.151468066884386
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,11.650996821073,0.00295133281036829,0.00491888801728048
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,13.4478350080053,0.00120182084187058,0.00295217659418973
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,13.8294233735875,0.000993067731985363,0.00295217659418973
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,42.2876068754728,6.5669368585673e-10,6.5669368585673e-09
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,17.8105083122929,0.000135674197059866,0.000678370985299332
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,2.6977762507496,0.259528663617795,0.28836518179755
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,1.86785410620894,0.393007315373388,0.393007315373388
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,13.0367194654262,0.00147608829709487,0.00295217659418973
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,3.20467970692152,0.201424661980554,0.251780827475693
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,5.63505781320873,0.0597534166962626,0.0853620238518037
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.731468216113309,0.286302634331705,11.650996821073,0.00295133281036829,0.00491888801728048
+"gene2",21.2812387436366,0.535277661455529,0.149824057325277,13.4478350080053,0.00120182084187058,0.00295217659418973
+"gene3",5.01012348853472,-0.673741470948083,0.287396810353061,13.8294233735875,0.000993067731985363,0.00295217659418973
+"gene4",100.51796142035,-0.423470709852266,0.106218338222812,42.2876068754728,6.5669368585673e-10,6.5669368585673e-09
+"gene5",27.1424502740787,0.588017156967864,0.152768106977919,17.8105083122929,0.000135674197059866,0.000678370985299332
+"gene6",5.4130427476525,-0.0194903186060668,0.30784249287543,2.6977762507496,0.259528663617795,0.28836518179755
+"gene7",28.2940230404605,0.135114917108395,0.14954938570182,1.86785410620894,0.393007315373388,0.393007315373388
+"gene8",40.3583444203556,-0.271585373967762,0.131499503574464,13.0367194654262,0.00147608829709487,0.00295217659418973
+"gene9",37.1661826339853,-0.21401822402321,0.133007375492728,3.20467970692152,0.201424661980554,0.251780827475693
+"gene10",11.5893249023836,0.389531997473394,0.244928947229985,5.63505781320873,0.0597534166962626,0.0853620238518037
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.632812315254828,0.289103638926867,4.73285162924537,0.0295917841865267,0.0663276489248662
+"gene2",21.2812387436366,0.538551973774898,0.149962947101771,12.8392275422336,0.000339427708341099,0.0016971385417055
+"gene3",5.01012348853472,-0.632831858568865,0.295221453480919,4.53730195229139,0.0331638244624331,0.0663276489248662
+"gene4",100.51796142035,-0.412101731455792,0.118627877254952,12.0212585890596,0.000525971731194201,0.00175323910398067
+"gene5",27.1424502740787,0.582066403700141,0.15473043602656,14.0806522268775,0.000175136067078282,0.0016971385417055
+"gene6",5.4130427476525,0.00145730049165907,0.310309910140032,0.00819390042352097,0.927873869713322,0.927873869713322
+"gene7",28.2940230404605,0.134335684024563,0.149986853732723,0.805396982860088,0.369484643110571,0.410538492345078
+"gene8",40.3583444203556,-0.270655780074565,0.136401659050961,3.92714778918253,0.0475124627494821,0.0791874379158035
+"gene9",37.1661826339853,-0.212714657425418,0.133243837041424,2.54476500804708,0.110660374187379,0.144411939198395
+"gene10",11.5893249023836,0.386011089345352,0.244586104807745,2.47690282625604,0.115529551358716,0.144411939198395
@@ -0,0 +1,11 @@
+"","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj"
+"gene1",8.54131729397935,0.632812315254828,0.289103638926867,4.73285162924537,0.0295917841865267,0.0663276489248662
+"gene2",21.2812387436366,0.538551973774898,0.149962947101771,12.8392275422336,0.000339427708341099,0.0016971385417055
+"gene3",5.01012348853472,-0.632831858568865,0.295221453480919,4.53730195229139,0.0331638244624331,0.0663276489248662
+"gene4",100.51796142035,-0.412101731455792,0.118627877254952,12.0212585890596,0.000525971731194201,0.00175323910398067
+"gene5",27.1424502740787,0.582066403700141,0.15473043602656,14.0806522268775,0.000175136067078282,0.0016971385417055
+"gene6",5.4130427476525,0.00145730049165907,0.310309910140032,0.00819390042352097,0.927873869713322,0.927873869713322
+"gene7",28.2940230404605,0.134335684024563,0.149986853732723,0.805396982860088,0.369484643110571,0.410538492345078
+"gene8",40.3583444203556,-0.270655780074565,0.136401659050961,3.92714778918253,0.0475124627494821,0.0791874379158035
+"gene9",37.1661826339853,-0.212714657425418,0.133243837041424,2.54476500804708,0.110660374187379,0.144411939198395
+"gene10",11.5893249023836,0.386011089345352,0.244586104807745,2.47690282625604,0.115529551358716,0.144411939198395