Let tidy_deseq_results() extract a coefficient by name

NEWS.md

@@ -1,5 +1,7 @@
 # clover 0.1.0.9000 (development version)
 
+* `tidy_deseq_results()` gains a `name` argument for extracting a single model coefficient, which is the only way to reach an interaction term. A design such as `~ genotype + charge_status + genotype:charge_status` carries differential charging in its interaction coefficient, and no `contrast` specification refers to it. Exactly one of `contrast` or `name` must be supplied.
+
 # clover 0.1.0
 
 * `read_bedmethyl()` reads per-position modification percentages from bedMethyl files produced by `modkit pileup`, enabling integration of nanopore modification data with the heatmap workflow.

R/deseq.R

@@ -207,16 +207,28 @@ run_deseq <- function(count_matrix, coldata, design, ...) {
 #' Extract results from a DESeq2 analysis and return a tidy tibble with
 #' tRNA identifiers and significance flags.
 #'
+#' Give either `contrast` or `name`, not both. A contrast compares two levels
+#' of one factor. A name addresses a single model coefficient, which is the
+#' only way to reach an interaction term: in a design such as
+#' `~ genotype + charge_status + genotype:charge_status`, the interaction
+#' coefficient is what carries differential charging, and no `contrast`
+#' specification refers to it. Use [DESeq2::resultsNames()] to list the
+#' available coefficients.
+#'
 #' @param dds A `DESeqDataSet` object (from [run_deseq()]).
 #' @param contrast A contrast specification: either a character vector of
 #'   length 3 (e.g., `c("condition", "mutant", "wildtype")`) or a list
 #'   for coefficient-based contrasts.
+#' @param name Name of a single model coefficient to extract, as returned by
+#'   [DESeq2::resultsNames()]. Use this for interaction terms.
 #' @param padj_cutoff Adjusted p-value threshold for significance.
 #'   Default `0.05`.
 #'
 #' @return A tibble with columns: `ref`, `log2FoldChange`, `lfcSE`,
 #'   `pvalue`, `padj`, and `significant` (logical).
 #'
+#' @seealso [run_deseq()], [charging_count_matrix()]
+#'
 #' @export
 #'
 #' @examples
@@ -229,11 +241,38 @@ run_deseq <- function(count_matrix, coldata, design, ...) {
 #' )
 #' dds <- run_deseq(counts, coldata, design = ~condition)
 #' tidy_deseq_results(dds, contrast = c("condition", "inf", "ctl"))
+#'
+#' # An interaction coefficient, by name
+#' DESeq2::resultsNames(dds)
+#' tidy_deseq_results(dds, name = "condition_inf_vs_ctl")
 #' }
-tidy_deseq_results <- function(dds, contrast, padj_cutoff = 0.05) {
+tidy_deseq_results <- function(
+  dds,
+  contrast = NULL,
+  name = NULL,
+  padj_cutoff = 0.05
+) {
   rlang::check_installed("DESeq2", reason = "to extract results.")
 
-  res <- DESeq2::results(dds, contrast = contrast)
+  if (is.null(contrast) == is.null(name)) {
+    cli_abort(
+      "Supply exactly one of {.arg contrast} or {.arg name}."
+    )
+  }
+
+  if (!is.null(name)) {
+    available <- DESeq2::resultsNames(dds)
+    if (!name %in% available) {
+      cli_abort(c(
+        "{.val {name}} is not a coefficient of this model.",
+        i = "Available: {.val {available}}."
+      ))
+    }
+    res <- DESeq2::results(dds, name = name)
+  } else {
+    res <- DESeq2::results(dds, contrast = contrast)
+  }
+
   res_df <- as.data.frame(res)
 
   tibble::tibble(

tests/testthat/test-deseq.R

@@ -126,3 +126,111 @@ test_that("run_deseq and tidy_deseq_results work end-to-end", {
   expect_type(res$significant, "logical")
   expect_equal(nrow(res), n_genes)
 })
+
+# A small fitted model to extract coefficients from. With `interaction = TRUE`
+# the columns are split into charged and uncharged halves, which is the shape of
+# a differential charging design.
+synthetic_dds <- function(interaction = FALSE, n_genes = 20) {
+  set.seed(42)
+  n_samples <- 6
+  lambdas <- sample(50:500, n_genes, replace = TRUE)
+
+  counts <- matrix(
+    as.integer(rpois(n_genes * n_samples, lambda = rep(lambdas, n_samples))),
+    nrow = n_genes,
+    ncol = n_samples
+  )
+  rownames(counts) <- paste0("tRNA-", seq_len(n_genes))
+  colnames(counts) <- c("wt1", "wt2", "wt3", "mut1", "mut2", "mut3")
+
+  coldata <- data.frame(
+    sample_id = colnames(counts),
+    condition = factor(rep(c("wt", "mut"), each = 3), levels = c("wt", "mut")),
+    row.names = colnames(counts)
+  )
+
+  if (!interaction) {
+    return(suppressWarnings(run_deseq(counts, coldata, design = ~condition)))
+  }
+
+  counts <- cbind(counts, counts)
+  colnames(counts) <- c(
+    paste0(colnames(coldata)[0], rownames(coldata), "_uncharged"),
+    paste0(rownames(coldata), "_charged")
+  )
+
+  coldata <- data.frame(
+    condition = factor(rep(coldata$condition, 2), levels = c("wt", "mut")),
+    charge_status = factor(
+      rep(c("uncharged", "charged"), each = n_samples),
+      levels = c("uncharged", "charged")
+    ),
+    row.names = colnames(counts)
+  )
+
+  suppressWarnings(
+    run_deseq(counts, coldata, design = ~ condition * charge_status)
+  )
+}
+
+test_that("tidy_deseq_results extracts a coefficient by name", {
+  skip_if_not_installed("DESeq2")
+
+  dds <- synthetic_dds()
+  coef <- setdiff(DESeq2::resultsNames(dds), "Intercept")[1]
+
+  by_name <- tidy_deseq_results(dds, name = coef)
+
+  expect_s3_class(by_name, "tbl_df")
+  expect_named(
+    by_name,
+    c("ref", "log2FoldChange", "lfcSE", "pvalue", "padj", "significant")
+  )
+  # The single coefficient of a two-level factor is that factor's contrast, so
+  # the two routes have to agree.
+  by_contrast <- tidy_deseq_results(dds, contrast = c("condition", "mut", "wt"))
+  expect_equal(by_name$log2FoldChange, by_contrast$log2FoldChange)
+  expect_equal(by_name$pvalue, by_contrast$pvalue)
+})
+
+test_that("tidy_deseq_results reaches an interaction coefficient", {
+  skip_if_not_installed("DESeq2")
+
+  # An interaction term is unreachable through `contrast`, which is the reason
+  # `name` exists: this is the shape of a differential charging model.
+  dds <- synthetic_dds(interaction = TRUE)
+  coef <- grep("\\.", DESeq2::resultsNames(dds), value = TRUE)
+  expect_length(coef, 1)
+
+  res <- tidy_deseq_results(dds, name = coef)
+
+  expect_s3_class(res, "tbl_df")
+  expect_equal(nrow(res), 20)
+})
+
+test_that("tidy_deseq_results requires exactly one of contrast and name", {
+  skip_if_not_installed("DESeq2")
+
+  dds <- synthetic_dds()
+
+  expect_error(tidy_deseq_results(dds), "exactly one")
+  expect_error(
+    tidy_deseq_results(
+      dds,
+      contrast = c("condition", "mut", "wt"),
+      name = "condition_mut_vs_wt"
+    ),
+    "exactly one"
+  )
+})
+
+test_that("tidy_deseq_results rejects an unknown coefficient name", {
+  skip_if_not_installed("DESeq2")
+
+  dds <- synthetic_dds()
+
+  expect_error(
+    tidy_deseq_results(dds, name = "no_such_coefficient"),
+    "not a coefficient"
+  )
+})