R/methods.R, R/methods_SE.R
test_differential_abundance-methods.Rdtest_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,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)
# S4 method for class 'spec_tbl_df'
test_differential_abundance(
.data,
.formula,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)
# S4 method for class 'tbl_df'
test_differential_abundance(
.data,
.formula,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)
# S4 method for class 'tidybulk'
test_differential_abundance(
.data,
.formula,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)
# S4 method for class 'SummarizedExperiment'
test_differential_abundance(
.data,
.formula,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)
# S4 method for class 'RangedSummarizedExperiment'
test_differential_abundance(
.data,
.formula,
.sample = NULL,
.transcript = NULL,
.abundance = NULL,
contrasts = NULL,
method = "edgeR_quasi_likelihood",
test_above_log2_fold_change = NULL,
scaling_method = "TMM",
omit_contrast_in_colnames = FALSE,
prefix = "",
action = "add",
...,
significance_threshold = NULL,
fill_missing_values = NULL,
.contrasts = NULL
)A `tbl` (with at least three columns for sample, feature and transcript abundance) or `SummarizedExperiment` (more convenient if abstracted to tibble with library(tidySummarizedExperiment))
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.
The name of the sample column
The name of the transcript/gene column
The name of the transcript/gene abundance column
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)
A string character. Either "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"
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.
A character string. The scaling method passed to the back-end functions: edgeR and limma-voom (i.e., edgeR::calcNormFactors; "TMM","TMMwsp","RLE","upperquartile"). Setting the parameter to \"none\" will skip the compensation for sequencing-depth for the method edgeR or limma-voom.
If just one contrast is specified you can choose to omit the contrast label in the colnames.
A character string. The prefix you would like to add to the result columns. It is useful if you want to compare several methods.
A character string. Whether to join the new information to the input tbl (add), or just get the non-redundant tbl with the new information (get).
Further arguments passed to some of the internal experimental functions. For example for glmmSeq, it is possible to pass .dispersion, and .scaling_factor column tidyeval to skip the caluclation of dispersion and scaling and use precalculated values. This is helpful is you want to calculate those quantities on many genes and do DE testing on fewer genes. .scaling_factor is the TMM value that can be obtained with tidybulk::scale_abundance.
DEPRECATED - A real between 0 and 1 (usually 0.05).
DEPRECATED - A boolean. Whether to fill missing sample/transcript values with the median of the transcript. This is rarely needed.
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)
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 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 consistent object (to the input) with additional columns for the statistics from the hypothesis test (e.g., log fold change, p-value and false discovery rate).
A consistent object (to the input) with additional columns for the statistics from the hypothesis test (e.g., log fold change, p-value and false discovery rate).
A `SummarizedExperiment` object
A `SummarizedExperiment` object
`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::calcNormFactors(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 = .data assay(my_assay)
# Create design matrix for dispersion, removing random effects design = model.matrix( object = .formula |> lme4::nobars(), data = metadata )
dispersion = counts |> edgeR::estimateDisp(design = design)
glmmSeq( .formula, countdata = counts , metadata = metadata |> as.data.frame(), dispersion = dispersion, progress = TRUE, method = method |> str_remove("(?i)^glmmSeq_" ), )
# edgeR
tidybulk::se_mini |>
identify_abundant() |>
test_differential_abundance( ~ condition )
#> No group or design set. Assuming all samples belong to one 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), conditionTRUE"
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$edgeR_quasi_likelihood`
#> class: SummarizedExperiment
#> dim: 182 5
#> metadata(0):
#> assays(1): count
#> rownames(182): ACAP1 ACP5 ... ZNF286A ZNF324
#> rowData names(7): entrez .abundant ... PValue FDR
#> colnames(5): SRR1740034 SRR1740035 SRR1740043 SRR1740058 SRR1740067
#> colData names(5): Cell.type time condition days dead
# The function `test_differential_abundance` operates with contrasts too
tidybulk::se_mini |>
identify_abundant(factor_of_interest = condition) |>
test_differential_abundance(
~ 0 + condition,
contrasts = c( "conditionTRUE - conditionFALSE")
)
#> tidybulk says: The design column names are "conditionFALSE, conditionTRUE"
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$edgeR_quasi_likelihood`
#> class: SummarizedExperiment
#> dim: 394 5
#> metadata(0):
#> assays(1): count
#> rownames(394): ABCB4 ABCB9 ... ZNF324 ZNF442
#> rowData names(7): entrez .abundant ... PValue___conditionTRUE -
#> conditionFALSE FDR___conditionTRUE - conditionFALSE
#> colnames(5): SRR1740034 SRR1740035 SRR1740043 SRR1740058 SRR1740067
#> colData names(5): Cell.type time condition days dead
# DESeq2 - equivalent for limma-voom
my_se_mini = tidybulk::se_mini
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")
#> converting counts to integer mode
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$deseq2`
#> class: SummarizedExperiment
#> dim: 394 5
#> metadata(0):
#> assays(1): count
#> rownames(394): ABCB4 ABCB9 ... ZNF324 ZNF442
#> rowData names(8): entrez .abundant ... pvalue padj
#> colnames(5): SRR1740034 SRR1740035 SRR1740043 SRR1740058 SRR1740067
#> colData names(5): Cell.type time condition days dead
# 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 )
#> converting counts to integer mode
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$deseq2`
# Use random intercept and random effect models
se_mini[1:50,] |>
identify_abundant(factor_of_interest = condition) |>
test_differential_abundance(
~ condition + (1 + condition | time),
method = "glmmseq_lme4", cores = 1
)
#>
#> n = 5 samples, 2 individuals
#> Time difference of 30.54381 secs
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$glmmseq_lme4`
#> class: SummarizedExperiment
#> dim: 40 5
#> metadata(0):
#> assays(1): count
#> rownames(40): ABCB4 ABCB9 ... C5AR1 C5AR2
#> rowData names(26): entrez .abundant ... P_condition
#> P_condition_adjusted
#> colnames(5): SRR1740034 SRR1740035 SRR1740043 SRR1740058 SRR1740067
#> colData names(5): Cell.type time condition days dead
# confirm that lfcThreshold was used
if (FALSE) { # \dontrun{
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", "TRUE", "FALSE")),
method="deseq2",
fitType="local"
)
#> No group or design set. Assuming all samples belong to one group.
#> converting counts to integer mode
#> estimating size factors
#> estimating dispersions
#> gene-wise dispersion estimates
#> mean-dispersion relationship
#> final dispersion estimates
#> fitting model and testing
#> tidybulk says: to access the raw results (fitted GLM) do `attr(..., "internals")$deseq2`
#> class: SummarizedExperiment
#> dim: 182 5
#> metadata(0):
#> assays(1): count
#> rownames(182): ACAP1 ACP5 ... ZNF286A ZNF324
#> rowData names(8): entrez .abundant ... pvalue___condition TRUE-FALSE
#> padj___condition TRUE-FALSE
#> colnames(5): SRR1740034 SRR1740035 SRR1740043 SRR1740058 SRR1740067
#> colData names(5): Cell.type time condition days dead