Perform differential transcription testing using edgeR quasi-likelihood (QLT), edgeR likelihood-ratio (LR), limma-voom, limma-voom-with-quality-weights or DESeq2

test_differential_abundance() takes as input A `tbl` (with at least three columns for sample, feature and transcript abundance) or `SummarizedExperiment` (more convenient if abstracted to tibble with library(tidySummarizedExperiment)) and returns a consistent object (to the input) with additional columns for the statistics from the hypothesis test.

test_differential_abundance(
  .data,
  .formula,
  abundance = assayNames(.data)[1],
  contrasts = NULL,
  method = c("edgeR_quasi_likelihood", "edgeR_likelihood_ratio",
    "edger_robust_likelihood_ratio", "DESeq2", "limma_voom", "limma_voom_sample_weights",
    "glmmseq_lme4", "glmmseq_glmmtmb"),
  test_above_log2_fold_change = NULL,
  scaling_method = "TMM",
  omit_contrast_in_colnames = FALSE,
  prefix = "",
  formula_dispersion = NULL,
  ...,
  significance_threshold = NULL,
  fill_missing_values = NULL,
  .contrasts = NULL,
  .abundance = NULL
)

# S4 method for class 'SummarizedExperiment'
test_differential_abundance(
  .data,
  .formula,
  abundance = assayNames(.data)[1],
  contrasts = NULL,
  method = c("edgeR_quasi_likelihood", "edgeR_likelihood_ratio",
    "edger_robust_likelihood_ratio", "DESeq2", "limma_voom", "limma_voom_sample_weights",
    "glmmseq_lme4", "glmmseq_glmmtmb"),
  test_above_log2_fold_change = NULL,
  scaling_method = "TMM",
  omit_contrast_in_colnames = FALSE,
  prefix = "",
  formula_dispersion = NULL,
  ...,
  significance_threshold = NULL,
  fill_missing_values = NULL,
  .contrasts = NULL,
  .abundance = NULL
)

# S4 method for class 'RangedSummarizedExperiment'
test_differential_abundance(
  .data,
  .formula,
  abundance = assayNames(.data)[1],
  contrasts = NULL,
  method = c("edgeR_quasi_likelihood", "edgeR_likelihood_ratio",
    "edger_robust_likelihood_ratio", "DESeq2", "limma_voom", "limma_voom_sample_weights",
    "glmmseq_lme4", "glmmseq_glmmtmb"),
  test_above_log2_fold_change = NULL,
  scaling_method = "TMM",
  omit_contrast_in_colnames = FALSE,
  prefix = "",
  formula_dispersion = NULL,
  ...,
  significance_threshold = NULL,
  fill_missing_values = NULL,
  .contrasts = NULL,
  .abundance = NULL
)

Arguments

.data: A `tbl` (with at least three columns for sample, feature and transcript abundance) or `SummarizedExperiment` (more convenient if abstracted to tibble with library(tidySummarizedExperiment))
.formula: A formula representing the desired linear model. If there is more than one factor, they should be in the order factor of interest + additional factors.
abundance: The name of the transcript/gene abundance column (character, preferred)
contrasts: This parameter takes the format of the contrast parameter of the method of choice. For edgeR and limma-voom is a character vector. For DESeq2 is a list including a character vector of length three. The first covariate is the one the model is tested against (e.g., ~ factor_of_interest)
method: A character vector. Available methods are "edgeR_quasi_likelihood" (i.e., QLF), "edgeR_likelihood_ratio" (i.e., LRT), "edger_robust_likelihood_ratio", "DESeq2", "limma_voom", "limma_voom_sample_weights", "glmmseq_lme4", "glmmseq_glmmtmb". Only one method can be specified at a time. For glmmSeq, pass `formula_dispersion` (a fixed-effects formula for edgeR) to plug in tagwise dispersion, or omit it to let each gene estimate phi.
test_above_log2_fold_change: A positive real value. This works for edgeR and limma_voom methods. It uses the `treat` function, which tests that the difference in abundance is bigger than this threshold rather than zero https://pubmed.ncbi.nlm.nih.gov/19176553.
scaling_method: A character string. The scaling method passed to the back-end functions: edgeR and limma-voom (i.e., edgeR::normLibSizes; "TMM","TMMwsp","RLE","upperquartile"). Setting the parameter to \"none\" will skip the compensation for sequencing-depth for the method edgeR or limma_voom.
omit_contrast_in_colnames: If just one contrast is specified you can choose to omit the contrast label in the colnames.
prefix: A character string. The prefix you would like to add to the result columns. It is useful if you want to compare several methods.
formula_dispersion: A fixed-effects formula for edgeR tagwise dispersion (glmmSeq). If NULL, each gene estimates phi.
...: Further arguments passed to some of the internal experimental functions. For glmmSeq: `.scaling_factor` skips TMM and uses that column from [scale_abundance()].
significance_threshold: DEPRECATED - A real between 0 and 1 (usually 0.05).
fill_missing_values: DEPRECATED - A boolean. Whether to fill missing sample/transcript values with the median of the transcript. This is rarely needed.
.contrasts: DEPRECATED - This parameter takes the format of the contrast parameter of the method of choice. For edgeR and limma-voom is a character vector. For DESeq2 is a list including a character vector of length three. The first covariate is the one the model is tested against (e.g., ~ factor_of_interest)
.abundance: DEPRECATED. The name of the transcript/gene abundance column (symbolic, for backward compatibility)

Value

A consistent object (to the input) with additional columns for the statistics from the test (e.g., log fold change, p-value and false discovery rate).

A `SummarizedExperiment` object

Details

`r lifecycle::badge("maturing")`

This function provides the option to use edgeR https://doi.org/10.1093/bioinformatics/btp616, limma-voom https://doi.org/10.1186/gb-2014-15-2-r29, limma_voom_sample_weights https://doi.org/10.1093/nar/gkv412 or DESeq2 https://doi.org/10.1186/s13059-014-0550-8 to perform the testing. All methods use raw counts, irrespective of if scale_abundance or adjust_abundance have been calculated, therefore it is essential to add covariates such as batch effects (if applicable) in the formula.

Underlying method for edgeR framework:

.data |>

# Filter keep_abundant( factor_of_interest = !!(as.symbol(parse_formula(.formula)[1])), minimum_counts = minimum_counts, minimum_proportion = minimum_proportion ) |>

# Format select(!!.transcript,!!.sample,!!.abundance) |> spread(!!.sample,!!.abundance) |> as_matrix(rownames = !!.transcript) |>

# edgeR edgeR::DGEList(counts = .) |> edgeR::normLibSizes(method = scaling_method) |> edgeR::estimateDisp(design) |>

# Fit edgeR::glmQLFit(design) |> // or glmFit according to choice edgeR::glmQLFTest(coef = 2, contrast = my_contrasts) // or glmLRT according to choice

Underlying method for DESeq2 framework:

keep_abundant( factor_of_interest = !!as.symbol(parse_formula(.formula)[[1]]), minimum_counts = minimum_counts, minimum_proportion = minimum_proportion ) |>

# DESeq2 DESeq2::DESeqDataSet(design = .formula) |> DESeq2::DESeq() |> DESeq2::results()

Underlying method for glmmSeq framework:

counts = assay(.data, my_assay)

# Tagwise phi from edgeR if formula_dispersion is given.

glmmSeq( .formula, countdata = counts , metadata = metadata |> as.data.frame(), dispersion = dispersion, progress = TRUE, method = method |> str_remove("(?i)^glmmSeq_" ), )

References

Mangiola, S., Molania, R., Dong, R., Doyle, M. A., & Papenfuss, A. T. (2021). tidybulk: an R tidy framework for modular transcriptomic data analysis. Genome Biology, 22(1), 42. doi:10.1186/s13059-020-02233-7

McCarthy, D. J., Chen, Y., & Smyth, G. K. (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Research, 40(10), 4288-4297. doi:10.1093/nar/gks042

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. doi:10.1186/s13059-014-0550-8

Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15(2), R29. doi:10.1186/gb-2014-15-2-r29

Examples

## Load airway dataset for examples

  data('airway', package = 'airway')
  # Ensure a 'condition' column exists for examples expecting it

    SummarizedExperiment::colData(airway)$condition <- SummarizedExperiment::colData(airway)$dex



 # edgeR (default method)

 airway |>
 identify_abundant() |>
  test_differential_abundance( ~ condition, method = "edgeR_quasi_likelihood" )
#> Warning: All samples appear to belong to the same group.
#> =====================================
#> tidybulk says: All testing methods use raw counts, irrespective of if scale_abundance
#> or adjust_abundance have been calculated. Therefore, it is essential to add covariates
#> such as batch effects (if applicable) in the formula.
#> =====================================
#> This message is displayed once per session.
#> tidybulk says: The design column names are "(Intercept), conditionuntrt"
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$edgeR_quasi_likelihood_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$edgeR_quasi_likelihood_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=14224 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>        .feature  .sample | counts | SampleName cell  dex   albut Run   avgLength
#>        <chr>     <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1      ENSG0000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2      ENSG0000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3      ENSG0000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4      ENSG0000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5      ENSG0000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>        --------  ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 113788 ENSG0000… SRR103… | 34     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113789 ENSG0000… SRR103… | 82     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113790 ENSG0000… SRR103… | 17     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113791 ENSG0000… SRR103… | 12     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113792 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 21 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, logFC <chr>, logCPM <chr>, F <chr>, PValue <chr>,
#> #   FDR <chr>

 # You can also explicitly specify the method
 airway |>
 identify_abundant() |>
  test_differential_abundance( ~ condition, method = "edgeR_quasi_likelihood" )
#> Warning: All samples appear to belong to the same group.
#> tidybulk says: The design column names are "(Intercept), conditionuntrt"
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$edgeR_quasi_likelihood_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$edgeR_quasi_likelihood_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=14224 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>        .feature  .sample | counts | SampleName cell  dex   albut Run   avgLength
#>        <chr>     <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1      ENSG0000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2      ENSG0000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3      ENSG0000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4      ENSG0000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5      ENSG0000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>        --------  ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 113788 ENSG0000… SRR103… | 34     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113789 ENSG0000… SRR103… | 82     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113790 ENSG0000… SRR103… | 17     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113791 ENSG0000… SRR103… | 12     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113792 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 21 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, logFC <chr>, logCPM <chr>, F <chr>, PValue <chr>,
#> #   FDR <chr>

  # The function `test_differential_abundance` operates with contrasts too

 airway |>
 identify_abundant(factor_of_interest = condition) |>
 test_differential_abundance(
      ~ 0 + condition,
      contrasts = c( "conditiontrt - conditionuntrt"),
    method = "edgeR_quasi_likelihood"
 )
#> tidybulk says: The design column names are "conditiontrt, conditionuntrt"
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$edgeR_quasi_likelihood_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$edgeR_quasi_likelihood_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=15926 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>        .feature  .sample | counts | SampleName cell  dex   albut Run   avgLength
#>        <chr>     <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1      ENSG0000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2      ENSG0000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3      ENSG0000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4      ENSG0000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5      ENSG0000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>        --------  ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 127404 ENSG0000… SRR103… | 17     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127405 ENSG0000… SRR103… | 12     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127406 ENSG0000… SRR103… | 10     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127407 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127408 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 21 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, `logFC___conditiontrt - conditionuntrt` <chr>,
#> #   `logCPM___conditiontrt - conditionuntrt` <chr>,
#> #   `F___conditiontrt - conditionuntrt` <chr>, …

 # DESeq2 - equivalent for limma-voom

my_se_mini = airway
my_se_mini$condition  = factor(my_se_mini$condition)

# demontrating with `fitType` that you can access any arguments to DESeq()
my_se_mini  |>
   identify_abundant(factor_of_interest = condition) |>
       test_differential_abundance( ~ condition, method="deseq2", fitType="local")
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$deseq2_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$deseq2_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=15926 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>        .feature  .sample | counts | SampleName cell  dex   albut Run   avgLength
#>        <chr>     <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1      ENSG0000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2      ENSG0000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3      ENSG0000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4      ENSG0000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5      ENSG0000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>        --------  ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 127404 ENSG0000… SRR103… | 17     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127405 ENSG0000… SRR103… | 12     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127406 ENSG0000… SRR103… | 10     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127407 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 127408 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 22 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, baseMean <chr>, log2FoldChange <chr>, lfcSE <chr>,
#> #   stat <chr>, pvalue <chr>, padj <chr>

# testing above a log2 threshold, passes along value to lfcThreshold of results()
res <- my_se_mini  |>
   identify_abundant(factor_of_interest = condition) |>
        test_differential_abundance( ~ condition, method="deseq2",
            fitType="local",
            test_above_log2_fold_change=4 )
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$deseq2_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$deseq2_fit`

# Use random intercept and random effect models

 airway[1:50,] |>
  identify_abundant(factor_of_interest = condition) |>
  test_differential_abundance(
    ~ condition + (1 + condition | cell),
    method = "glmmseq_lme4",
    formula_dispersion = ~ condition + cell,
    cores = 1
  )
#> Registered S3 method overwritten by 'lme4':
#>   method           from
#>   na.action.merMod car 
#> Warning: package version mismatch: 
#> glmmTMB was built with TMB package version 1.9.21
#> Current TMB package version is 1.9.25
#> Please re-install glmmTMB from source or restore original ‘TMB’ package (see '?reinstalling' for more information)
#> tidybulk says: calculating tagwise (shrinked) and trended dispersion with edgeR::estimateDisp() using ~condition + cell.
#> tidybulk says: plugging rowData column dispersion_shrinked into glmmSeq.
#> 
#> n = 8 samples, 4 individuals
#> Time difference of 31.00323 secs
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$glmmseq_lme4_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$glmmseq_lme4_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=42 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>     .feature     .sample | counts | SampleName cell  dex   albut Run   avgLength
#>     <chr>        <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1   ENSG0000000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2   ENSG0000000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3   ENSG0000000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4   ENSG0000000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5   ENSG0000000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>     --------     ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 332 ENSG0000000… SRR103… | 1347   | GSM1275875 N061… trt   untrt SRR1… 98       
#> 333 ENSG0000000… SRR103… | 77     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 334 ENSG0000000… SRR103… | 1256   | GSM1275875 N061… trt   untrt SRR1… 98       
#> 335 ENSG0000000… SRR103… | 2809   | GSM1275875 N061… trt   untrt SRR1… 98       
#> 336 ENSG0000000… SRR103… | 921    | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 52 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, Dispersion <chr>, AIC <chr>, logLik <chr>, meanExp <chr>,
#> #   `(Intercept)` <chr>, conditionuntrt <chr>,
#> #   `N052611_cell__(Intercept)__lower` <chr>, …

# confirm that lfcThreshold was used
# Not run to keep the example fast
# res |>
#     mcols() |>
#     DESeq2::DESeqResults() |>
#     DESeq2::plotMA()


# The function `test_differential_abundance` operates with contrasts too

 my_se_mini |>
 identify_abundant() |>
 test_differential_abundance(
      ~ 0 + condition,
      contrasts = list(c("condition", "trt", "untrt")),
      method="deseq2",
     fitType="local"
 )
#> Warning: All samples appear to belong to the same group.
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the DE object do `metadata(.)$tidybulk$deseq2_object`
#> tidybulk says: to access the raw results (fitted GLM) do `metadata(.)$tidybulk$deseq2_fit`
#> # A SummarizedExperiment-tibble abstraction: Features=14224 | Samples=8 | 
#> #   Assays=counts
#> #                                 |----------------- COVARIATES ---------------|
#>        .feature  .sample | counts | SampleName cell  dex   albut Run   avgLength
#>        <chr>     <chr>   | <chr>  | <fct>      <fct> <fct> <fct> <fct> <chr>    
#> 1      ENSG0000… SRR103… | 679    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 2      ENSG0000… SRR103… | 467    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 3      ENSG0000… SRR103… | 260    | GSM1275862 N613… untrt untrt SRR1… 126      
#> 4      ENSG0000… SRR103… | 60     | GSM1275862 N613… untrt untrt SRR1… 126      
#> 5      ENSG0000… SRR103… | 3251   | GSM1275862 N613… untrt untrt SRR1… 126      
#>        --------  ------- - ------ - ---------- ----  ---   ----- ---   ---------
#> 113788 ENSG0000… SRR103… | 34     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113789 ENSG0000… SRR103… | 82     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113790 ENSG0000… SRR103… | 17     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113791 ENSG0000… SRR103… | 12     | GSM1275875 N061… trt   untrt SRR1… 98       
#> 113792 ENSG0000… SRR103… | 11     | GSM1275875 N061… trt   untrt SRR1… 98       
#> # ℹ 22 more variables: Experiment <fct>, Sample <fct>, BioSample <fct>,
#> #   condition <fct>, `|` <|>, gene_id <chr>, gene_name <chr>, entrezid <chr>,
#> #   gene_biotype <chr>, gene_seq_start <chr>, gene_seq_end <chr>,
#> #   seq_name <chr>, seq_strand <chr>, seq_coord_system <chr>, symbol <chr>,
#> #   .abundant <chr>, `baseMean___condition trt-untrt` <chr>,
#> #   `log2FoldChange___condition trt-untrt` <chr>,
#> #   `lfcSE___condition trt-untrt` <chr>, `stat___condition trt-untrt` <chr>, …