Run metadata-driven DESeq2 differential expression — differential

Fits a DESeq2 model from the tximport-derived counts of a study and returns one result table per comparison.

Usage

differential_deseq2(
  study,
  formula,
  term,
  metadata = NULL,
  reference = character(),
  contrast = character(),
  at = character(),
  test = c("wald", "lrt"),
  reduced = NULL,
  categorical = character(),
  numeric = character(),
  alpha = 0.05
)

Arguments

study: A rsemflow_study or study path.
formula: R formula or formula string.
term: Metadata variable or interaction term to report.
metadata: Optional replacement metadata path/data.frame.
reference: Character vector like condition::Control.
contrast: Character vector like condition::Drug::Control.
at: Optional conditioning specification like genotype::KO; currently supports one simple-effect condition with an explicit contrast.
test: wald or lrt.
reduced: Reduced formula required for LRT.
categorical: Metadata columns to force categorical.
numeric: Metadata columns to force numeric.
alpha: Significance threshold passed to DESeq2 results.

Value

A list with the fitted DESeqDataSet (dds), the metadata used, the formula(s), the design matrix, a contrast manifest, and the named result tables (results).

Details

Which tables are produced depends on term:

a categorical variable: every level versus the reference level;
a numeric variable: the per-unit coefficient;
an interaction such as genotype:treatment: every interaction coefficient.

Extra comparisons can be requested with contrast. With at, a single simple effect is returned instead (for example the treatment effect within one genotype in a ~ genotype * treatment model). With test = "lrt", one omnibus table is returned and log2FoldChange/lfcSE are set to NA because the likelihood-ratio test has no single direction.

Log fold changes are not shrunk.

Examples

# \donttest{
ex <- system.file("extdata", "example", package = "rsemflow")
study <- read_rsem_study(file.path(ex, "data"), file.path(ex, "metadata.tsv"))
de <- differential_deseq2(
  study,
  formula = ~ batch + condition,
  term = "condition",
  reference = "condition::Control"
)
head(de$results$condition_Drug_vs_Control)
#>             gene_id ensembl_gene_id      term        contrast coefficient test
#> 1 ENSG00000000001.2 ENSG00000000001 condition Drug vs Control        <NA> Wald
#> 2 ENSG00000000002.3 ENSG00000000002 condition Drug vs Control        <NA> Wald
#> 3 ENSG00000000003.1 ENSG00000000003 condition Drug vs Control        <NA> Wald
#> 4 ENSG00000000004.2 ENSG00000000004 condition Drug vs Control        <NA> Wald
#> 5 ENSG00000000005.3 ENSG00000000005 condition Drug vs Control        <NA> Wald
#> 6 ENSG00000000006.1 ENSG00000000006 condition Drug vs Control        <NA> Wald
#>   baseMean log2FoldChange     lfcSE     stat       pvalue         padj
#> 1 169.1780       2.397434 0.5764777 4.158763 3.199762e-05 0.0007999406
#> 2 208.8794       2.232664 0.6062340 3.682842 2.306482e-04 0.0038441366
#> 3 248.5807       2.127022 0.6689849 3.179476 1.475414e-03 0.0126338835
#> 4 288.2820       2.053341 0.7721432 2.659275 7.830909e-03 0.0522060612
#> 5 327.9833       1.998960 0.8977139 2.226723 2.596576e-02 0.1207709981
#> 6 367.6847       1.957152 1.0240589 1.911171 5.598259e-02 0.1696442033
# }