| Title: | Differential Analysis of Omics Data Based on the Tweedie Distribution |
|---|---|
| Description: | A toolkit for differential analysis of omics data. Implements a CPLM-based statistical workflow built on the Tweedie distribution. Unlike traditional single-omics tools, Tweedieverse is technology-agnostic and can be applied to both count and continuous measurements arising from diverse high-throughput technologies (e.g. transcript abundances from bulk and single-cell RNA-Seq studies in the form of UMI counts or non-UMI counts, microbiome taxonomic and functional profiles in the form of counts or relative abundances, and compound abundance levels or peak intensities from metabolomics and other mass spectrometry-based experiments, among others). The software includes multiple customization options such as the inclusion of random effects and multiple covariates along with several data exploration capabilities and visualization modules in a unified estimation umbrella. |
| Authors: | Nalin Arora [aut, ctb, cre] (ORCID: <https://orcid.org/0009-0009-1340-688X>), Ziyu Liu [ctb], Ali Rahnavard [aut, ctb] (ORCID: <https://orcid.org/0000-0002-9710-0248>), Himel Mallick [aut, ctb] (ORCID: <https://orcid.org/0000-0003-4956-2429>) |
| Maintainer: | Nalin Arora <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 0.99.0 |
| Built: | 2026-09-28 22:22:35 UTC |
| Source: | https://github.com/BiocStaging/Tweedieverse |
Compute Shannon entropy for a vector.
entropy(target)entropy(target)
target |
A vector of values. |
A numeric entropy value.
entropy(c("A", "A", "B", "B", "C"))entropy(c("A", "A", "B", "B", "C"))
Adjust Tweedieverse(or any other differential analysis methods) coefficient estimates and p-values by testing each taxon against the *median* effect for the same metadata variable - a simple post-hoc strategy to curb false discoveries driven by the compositional nature of microbiome count data (after the approach adopted in **MaAsLin 3**).
median_comparison_tweedie( df, p_cutoff = 0.95, subtract_median = FALSE, n_sims = 10000, median_threshold = 0 )median_comparison_tweedie( df, p_cutoff = 0.95, subtract_median = FALSE, n_sims = 10000, median_threshold = 0 )
df |
A 'data.frame' returned by **Tweedieverse** or any other differential analysis methods containing (at minimum) the columns 'taxon', 'metadata', 'effect_size', 'pval', 'stderr', and 'qval'. |
p_cutoff |
Numeric. Upper bound on the original p-value to include an effect in the median calculation (default '0.95' = all non-missing). |
subtract_median |
Logical. If 'TRUE', subtracts the group median from every coefficient before returning it in 'coef_median'; otherwise the original 'effect_size' is copied unchanged (default 'FALSE'). |
n_sims |
Integer. Number of Monte-Carlo simulations used to estimate the covariance between each coefficient and the group median (default '10 000'). |
median_threshold |
Numeric. Absolute difference below which a coefficient is considered effectively equal to the median and assigned 'pval_median = 1' (default '0'). |
For every distinct value in 'metadata' the algorithm
keeps coefficients with 'pval < p_cutoff' and computes their median;
simulates 'n_sims' draws of coefficients using a normal approximation ('N(effect_size, stderr^2)') and records the empirical distribution of the simulated medians;
derives a variance-inflated *offset* that accounts for the covariance between each coefficient and the group median;
performs a two-sided Z-test of 'H0 : beta_i = offset_i', returning the resulting p-value in 'pval_median'.
The input 'df' with two new columns:
Median-centred coefficient (or the original 'effect_size' if 'subtract_median = FALSE').
Two-sided p-value from the median-comparison test.
The procedure is heuristic and relies on normal approximations. Results may be unstable for very small sample sizes or when 'stderr' values are zero or missing.
[MaAsLin 3 GitHub](https://github.com/biobakery/maaslin3)
toy <- data.frame( taxon = c("tax1", "tax2"), metadata = c("grp", "grp"), effect_size = c(0.4, 0.2), pval = c(0.01, 0.2), stderr = c(0.1, 0.1), qval = c(0.02, 0.25) ) median_comparison_tweedie(toy, n_sims = 100)toy <- data.frame( taxon = c("tax1", "tax2"), metadata = c("grp", "grp"), effect_size = c(0.4, 0.2), pval = c(0.01, 0.2), stderr = c(0.1, 0.1), qval = c(0.02, 0.25) ) median_comparison_tweedie(toy, n_sims = 100)
A compact ggplot2 theme used by Tweedieverse visualization functions.
theme_omicsEye()theme_omicsEye()
A list of ggplot2 theme elements.
theme_omicsEye()theme_omicsEye()
Fit a per-feature Tweedie generalized linear model to omics features.
Tweedieverse( input_features, input_metadata = NULL, output = NULL, assay_name = "counts", abd_threshold = 0, prev_threshold = 0.1, var_threshold = 0, entropy_threshold = 0, base_model = "CPLM", link = "log", tweedie_p = NULL, p0_transform = "NONE", p0_transform_pseudocount = 1, fixed_effects = NULL, random_effects = NULL, domain = "microbiome", normalization = NULL, adjust_offset = TRUE, scale_factor = NULL, max_significance = 0.05, correction = "BH", Maaslin2_run = TRUE, median_comparison = FALSE, median_subtraction = FALSE, run_presence_absence_model = FALSE, method_args = NULL, method.args = NULL, standardize = TRUE, BPPARAM = BiocParallel::SerialParam(), optimizer = "nlminb", na.action = na.exclude, plot_heatmap = FALSE, plot_scatter = FALSE, heatmap_first_n = 50, reference = NULL )Tweedieverse( input_features, input_metadata = NULL, output = NULL, assay_name = "counts", abd_threshold = 0, prev_threshold = 0.1, var_threshold = 0, entropy_threshold = 0, base_model = "CPLM", link = "log", tweedie_p = NULL, p0_transform = "NONE", p0_transform_pseudocount = 1, fixed_effects = NULL, random_effects = NULL, domain = "microbiome", normalization = NULL, adjust_offset = TRUE, scale_factor = NULL, max_significance = 0.05, correction = "BH", Maaslin2_run = TRUE, median_comparison = FALSE, median_subtraction = FALSE, run_presence_absence_model = FALSE, method_args = NULL, method.args = NULL, standardize = TRUE, BPPARAM = BiocParallel::SerialParam(), optimizer = "nlminb", na.action = na.exclude, plot_heatmap = FALSE, plot_scatter = FALSE, heatmap_first_n = 50, reference = NULL )
input_features |
A domain-appropriate Bioconductor container, such as
|
input_metadata |
Optional R data frame of metadata. Samples are expected to have
matching sample names with |
output |
The output folder to write results. |
assay_name |
If the input is provided as one of the accepted Bioconductor objects,
this argument selects the name of the assay slot in the input object that contains the omics measurements.
For |
abd_threshold |
If prevalence-abundance filtering is desired, only features that are present (or detected)
in at least |
prev_threshold |
If prevalence-abundance filtering is desired, only features that are present (or detected)
in at least |
var_threshold |
If variance filtering is desired, only features that have variances greater than
|
entropy_threshold |
If entropy-based filtering is desired for metadata, only features that have entropy greater than
|
base_model |
The per-feature base model. Only "CPLM" is supported. |
link |
A specification of the GLM link function. Default is "log". Must be one of "log", "identity", "sqrt", or "inverse". |
tweedie_p |
Numeric Tweedie variance power. Default is |
p0_transform |
Transformation to apply to a copy of the filtered feature table when the
resolved Tweedie index is |
p0_transform_pseudocount |
Numeric pseudocount used for |
fixed_effects |
Metadata variable(s) describing the fixed effects coefficients. |
random_effects |
Metadata variable(s) describing the random effects part of the model. |
domain |
Domain used to choose the default size-factor normalization. Must be one of
|
normalization |
Size-factor strategy used only to compute the model offset. The input feature
table is never normalized or transformed. If |
adjust_offset |
If TRUE (default), an offset term will be included as the logarithm of the
size factor estimated by |
scale_factor |
Name of a numerical metadata variable containing user-supplied sample size
factors. When supplied, this overrides |
max_significance |
The q-value threshold for significance. Default is 0.05. |
correction |
The correction method for computing the q-value (see |
Maaslin2_run |
Logical. For |
median_comparison |
If TRUE, coefficients will be tested against a null value corresponding to the median coefficient for a covariate in the |
median_subtraction |
If TRUE, coefficients minus median will be used for compositionality adjustment. |
run_presence_absence_model |
If TRUE, also fit a DAssemble-style presence-absence logistic regression model
alongside the selected abundance model, including the MaAsLin2 path used by default for |
method_args |
Optional named list of method-specific arguments. For |
method.args |
Alias for |
standardize |
Should continuous metadata be standardized? Default is TRUE. Bypassed for categorical variables. |
BPPARAM |
|
optimizer |
The optimization routine to be used for estimating the parameters of the Tweedie model.
Possible choices are |
na.action |
How to handle missing values? See |
plot_heatmap |
Logical. If TRUE (default is FALSE), generate a heatmap of the (top |
plot_scatter |
Logical. If TRUE (default is FALSE), generate scatter/box plots of individual associations. |
heatmap_first_n |
In heatmap, plot top N features with significant associations (default is 50). |
reference |
The factor to use as a reference for a variable with more than two levels provided as a string of 'variable,reference' semi-colon delimited for multiple variables (default is NULL). |
For single-omics input, a data frame containing coefficient estimates, p-values,
and q-values (multiplicity-adjusted p-values) is returned. For MultiAssayExperiment
input, a named list of omics-specific result data frames is returned.
Himel Mallick, [email protected]
set.seed(123) features <- as.data.frame(matrix(rpois(24, lambda = 5), nrow = 12, ncol = 2)) colnames(features) <- c("feature1", "feature2") rownames(features) <- paste0("sample", seq_len(nrow(features))) metadata <- data.frame( group = rep(c("A", "B"), each = 6), row.names = rownames(features) ) experiment <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = t(as.matrix(features))), colData = metadata ) fit <- Tweedieverse( input_features = experiment, output = NULL, fixed_effects = "group", median_comparison = FALSE ) head(fit) maaslin2_median_fit <- Tweedieverse( input_features = experiment, output = NULL, fixed_effects = "group", tweedie_p = 0, Maaslin2_run = TRUE, median_comparison = TRUE, median_subtraction = TRUE ) # For Bioconductor-container workflow examples, see: # vignette("Tweedieverse", package = "Tweedieverse")set.seed(123) features <- as.data.frame(matrix(rpois(24, lambda = 5), nrow = 12, ncol = 2)) colnames(features) <- c("feature1", "feature2") rownames(features) <- paste0("sample", seq_len(nrow(features))) metadata <- data.frame( group = rep(c("A", "B"), each = 6), row.names = rownames(features) ) experiment <- SummarizedExperiment::SummarizedExperiment( assays = list(counts = t(as.matrix(features))), colData = metadata ) fit <- Tweedieverse( input_features = experiment, output = NULL, fixed_effects = "group", median_comparison = FALSE ) head(fit) maaslin2_median_fit <- Tweedieverse( input_features = experiment, output = NULL, fixed_effects = "group", tweedie_p = 0, Maaslin2_run = TRUE, median_comparison = TRUE, median_subtraction = TRUE ) # For Bioconductor-container workflow examples, see: # vignette("Tweedieverse", package = "Tweedieverse")
Create a heatmap summarizing association statistics from Tweedieverse results.
Tweedieverse_heatmap( output_results, title = NA, cell_value = "qval", data_label = "data", metadata_label = "metadata", border_color = "grey93", color = colorRampPalette(c("darkblue", "grey90", "darkred")), col_rotate = 90, first_n = 50, write_to = NA )Tweedieverse_heatmap( output_results, title = NA, cell_value = "qval", data_label = "data", metadata_label = "metadata", border_color = "grey93", color = colorRampPalette(c("darkblue", "grey90", "darkred")), col_rotate = 90, first_n = 50, write_to = NA )
output_results |
Tweedieverse result data frame or path to a TSV file. |
title |
Optional plot title. |
cell_value |
Which value to display ('"qval"', '"pval"', or '"coef"'). |
data_label |
Label for feature axis. |
metadata_label |
Label for metadata axis. |
border_color |
Heatmap border color. |
color |
Color palette function. |
col_rotate |
Rotation angle for column labels. |
first_n |
Number of top features to include. |
write_to |
Optional output directory for serialized plot objects. |
A 'pheatmap' object or 'NULL' when insufficient data are available.
toy <- data.frame( feature = c("tax1", "tax2", "tax1", "tax2"), metadata = c("grp", "grp", "time", "time"), value = c("A", "B", "A", "B"), pval = c(0.01, 0.05, 0.02, 0.03), qval = c(0.02, 0.08, 0.03, 0.04), coef = c(1.2, -0.7, 0.9, -0.4) ) Tweedieverse_heatmap(toy, first_n = 2)toy <- data.frame( feature = c("tax1", "tax2", "tax1", "tax2"), metadata = c("grp", "grp", "time", "time"), value = c("A", "B", "A", "B"), pval = c(0.01, 0.05, 0.02, 0.03), qval = c(0.02, 0.08, 0.03, 0.04), coef = c(1.2, -0.7, 0.9, -0.4) ) Tweedieverse_heatmap(toy, first_n = 2)