Tweedieverse provides a unified differential analysis
workflow for omics measurements whose empirical variance structure is
well represented by the Tweedie family. The framework is intended for
heterogeneous high-throughput data, including microbiome abundances,
single-cell RNA-seq measurements, bulk RNA-seq summaries, metabolite
intensities, and other continuous or count-like assays.
The package uses Bioconductor containers as the primary interface. A
single omics assay is represented as a
SummarizedExperiment,
TreeSummarizedExperiment, or
SingleCellExperiment, and matched multi-omics assays are
represented as a MultiAssayExperiment. This keeps feature
assays, sample metadata, and experiment-level annotations together and
makes the workflow compatible with the broader Bioconductor
ecosystem.
This vignette demonstrates six Bioconductor-oriented workflows:
TreeSummarizedExperiment,SingleCellExperiment,SingleCellExperiment,muscat::aggregateData(),airway
dataset stored in a SummarizedExperiment,MultiAssayExperiment.Modern omics studies increasingly combine assays with different
measurement scales, sparsity patterns, and mean-variance relationships.
Microbiome profiles, single-cell RNA-seq measurements, bulk
transcriptomic summaries, and metabolomics intensities can all contain
many low-abundance or zero-valued features while still differing in the
variance model most appropriate for differential analysis.
Tweedieverse is designed for this setting by using the
Tweedie variance function as a common modeling framework while
preserving domain-specific choices such as normalization-derived
offsets.
Bioconductor containers provide a scientifically reproducible
representation of these data. Assays, sample-level covariates, feature
identifiers, and experiment metadata remain linked inside formal objects
such as SummarizedExperiment,
SingleCellExperiment,
TreeSummarizedExperiment, and
MultiAssayExperiment. This reduces ambiguity about sample
orientation, metadata matching, and multi-omics alignment, which are
common sources of analysis error in table-based workflows. It also
allows Tweedieverse results to fit naturally into existing
Bioconductor workflows for preprocessing, quality control,
visualization, and downstream biological interpretation.
After Bioconductor acceptance, the release version can be installed with:
if (!requireNamespace("BiocManager", quietly = TRUE)) {
install.packages("BiocManager")
}
BiocManager::install("Tweedieverse")The development version can be installed from GitHub:
if (!requireNamespace("remotes", quietly = TRUE)) {
install.packages("remotes")
}
remotes::install_github("himelmallick/Tweedieverse")
library(Tweedieverse)
library(SummarizedExperiment)
#> Loading required package: MatrixGenerics
#> Loading required package: matrixStats
#>
#> Attaching package: 'MatrixGenerics'
#> The following objects are masked from 'package:matrixStats':
#>
#> colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
#> colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
#> colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
#> colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
#> colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
#> colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
#> colWeightedMeans, colWeightedMedians, colWeightedSds,
#> colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
#> rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
#> rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
#> rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
#> rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
#> rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
#> rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
#> rowWeightedSds, rowWeightedVars
#> Loading required package: GenomicRanges
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#> Loading required package: generics
#>
#> Attaching package: 'generics'
#> The following objects are masked from 'package:base':
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following object is masked from 'package:utils':
#>
#> data
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget,
#> order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
#> rbind, Reduce, rownames, sapply, saveRDS, scale, sequence, table,
#> tapply, transform, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> expand.grid, I, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> Loading required package: Biobase
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#>
#> Attaching package: 'Biobase'
#> The following object is masked from 'package:MatrixGenerics':
#>
#> rowMedians
#> The following objects are masked from 'package:matrixStats':
#>
#> anyMissing, rowMedians
library(SingleCellExperiment)
library(TreeSummarizedExperiment)
#> Loading required package: Biostrings
#> Loading required package: XVector
#>
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#>
#> strsplit
library(MultiAssayExperiment)
library(muscat)
library(airway)The example containers used in this vignette are distributed with the
package under inst/extdata. This allows the vignette to run
during Bioconductor checks without network access.
For microbiome data, Tweedieverse() accepts a
TreeSummarizedExperiment. The assay contains taxa/features
in rows and samples in columns; colData contains sample
metadata.
The example table contains relative abundances, so this analysis
disables additional offset normalization by setting
normalization = "NONE". For raw microbiome counts, omit
that argument to use the microbiome default TSS size-factor
offset.
microbiome_tse
#> class: TreeSummarizedExperiment
#> dim: 6 12
#> metadata(0):
#> assays(1): counts
#> rownames(6): Bacteroides Faecalibacterium ... Roseburia Other_taxa
#> rowData names(0):
#> colnames(12): MB01 MB02 ... MB11 MB12
#> colData names(3): diagnosis antibiotics site
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
#> rowLinks: NULL
#> rowTree: NULL
#> colLinks: NULL
#> colTree: NULL
microbiome_fit <- Tweedieverse(
input_features = microbiome_tse,
output = NULL,
assay_name = "counts",
domain = "microbiome",
normalization = "NONE",
fixed_effects = c("diagnosis", "antibiotics"),
tweedie_p = NULL,
median_comparison = FALSE
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:28.506205 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:28.511712 INFO::Determining format of input files
#> 2026-09-28 22:22:28.516718 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:28.518405 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:28.518715 INFO::Total samples in data: 12
#> 2026-09-28 22:22:28.518952 INFO::Min samples required with min abundance for a feature not to be filtered: 1.200000
#> 2026-09-28 22:22:28.519594 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:28.519927 INFO::Filtered feature names:
#> 2026-09-28 22:22:28.520382 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:28.520657 INFO::Filtered feature names:
#> 2026-09-28 22:22:28.521144 INFO::Offset size factors computed with normalization: NONE
#> 2026-09-28 22:22:28.522654 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:28.523012 INFO::Formula for fixed effects: expr ~ diagnosis + antibiotics
#> 2026-09-28 22:22:28.523442 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:28.533665 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:28.555656 INFO::Fitting model to feature number 1, Bacteroides
#> 2026-09-28 22:22:28.766166 INFO::Fitting model to feature number 2, Faecalibacterium
#> 2026-09-28 22:22:28.805638 INFO::Fitting model to feature number 3, Escherichia
#> 2026-09-28 22:22:28.962818 INFO::Fitting model to feature number 4, Akkermansia
#> 2026-09-28 22:22:29.003861 INFO::Fitting model to feature number 5, Roseburia
#> 2026-09-28 22:22:29.03698 INFO::Fitting model to feature number 6, Other_taxa
#> 2026-09-28 22:22:29.904716 INFO::Counting prevalence for each feature
head(microbiome_fit)
#> feature metadata value coef stderr pval
#> 1 Bacteroides diagnosis nonIBD 0.6588056 0.05615203 9.327652e-07
#> 2 Escherichia diagnosis nonIBD -1.1261872 0.10529776 2.039466e-06
#> 3 Faecalibacterium diagnosis nonIBD -0.4054021 0.04498805 8.450622e-06
#> 4 Akkermansia diagnosis nonIBD 0.6221236 0.10651038 2.465199e-04
#> 5 Other_taxa diagnosis nonIBD 0.2642319 0.06716538 3.436856e-03
#> 6 Roseburia diagnosis nonIBD -0.1230465 0.03592271 7.564099e-03
#> qval base.model tweedie.index N N.not.zero percent.zero
#> 1 1.119318e-05 CPLM 1.99 12 12 0
#> 2 1.223679e-05 CPLM 1.01 12 12 0
#> 3 3.380249e-05 CPLM 1.01 12 12 0
#> 4 7.395596e-04 CPLM 1.01 12 12 0
#> 5 8.248454e-03 CPLM 1.99 12 12 0
#> 6 1.512820e-02 CPLM 1.01 12 12 0For single-cell RNA-seq data, Tweedieverse() accepts a
SingleCellExperiment. The UMI example uses integer molecule
counts.
UMI counts are count-like. Here we pin the variance power to
p = 1, the Poisson variance case, and use the single-cell
domain. normalization = "SCRAN" computes size factors for
the model offset without modifying the assay.
umi_sce
#> class: SingleCellExperiment
#> dim: 6 12
#> metadata(0):
#> assays(1): counts
#> rownames(6): GeneA GeneB ... GeneE GeneF
#> rowData names(0):
#> colnames(12): UMI01 UMI02 ... UMI11 UMI12
#> colData names(3): condition batch donor
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
umi_fit <- Tweedieverse(
input_features = umi_sce,
output = NULL,
assay_name = "counts",
domain = "single_cell",
normalization = "SCRAN",
fixed_effects = c("condition", "batch"),
tweedie_p = 1,
median_comparison = FALSE
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:29.971604 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:29.972442 INFO::Determining format of input files
#> 2026-09-28 22:22:29.972784 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:29.973866 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:29.97418 INFO::Total samples in data: 12
#> 2026-09-28 22:22:29.974433 INFO::Min samples required with min abundance for a feature not to be filtered: 1.200000
#> 2026-09-28 22:22:29.974829 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:29.975117 INFO::Filtered feature names:
#> 2026-09-28 22:22:29.975547 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:29.975807 INFO::Filtered feature names:
#> 2026-09-28 22:22:30.196584 INFO::Offset size factors computed with normalization: SCRAN
#> 2026-09-28 22:22:30.198139 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:30.19856 INFO::Formula for fixed effects: expr ~ condition + batch
#> 2026-09-28 22:22:30.19891 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:30.200914 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:30.21371 INFO::Fitting model to feature number 1, GeneA
#> 2026-09-28 22:22:30.217154 INFO::Fitting model to feature number 2, GeneB
#> 2026-09-28 22:22:30.219588 INFO::Fitting model to feature number 3, GeneC
#> 2026-09-28 22:22:30.222051 INFO::Fitting model to feature number 4, GeneD
#> 2026-09-28 22:22:30.224421 INFO::Fitting model to feature number 5, GeneE
#> 2026-09-28 22:22:30.226716 INFO::Fitting model to feature number 6, GeneF
#> 2026-09-28 22:22:30.233018 INFO::Counting prevalence for each feature
head(umi_fit)
#> feature metadata value coef stderr pval qval
#> 1 GeneB condition stimulated -0.5993041 0.04219361 1.810091e-07 2.172110e-06
#> 2 GeneA condition stimulated 0.5327817 0.04406900 7.226339e-07 4.335804e-06
#> 3 GeneE condition stimulated 1.0812511 0.12284247 1.023904e-05 4.095614e-05
#> 4 GeneF condition stimulated -0.1984687 0.03562332 3.467956e-04 1.040387e-03
#> 5 GeneC condition stimulated 1.9589759 0.40088608 8.634536e-04 2.072289e-03
#> base.model tweedie.index N N.not.zero percent.zero
#> 1 Tweedie GLM 1 12 12 0
#> 2 Tweedie GLM 1 12 12 0
#> 3 Tweedie GLM 1 12 12 0
#> 4 Tweedie GLM 1 12 12 0
#> 5 Tweedie GLM 1 12 9 25The non-UMI example uses read-count-like measurements. These can have
a more flexible mean-variance relationship than UMI counts, so we leave
tweedie_p = NULL and let Tweedieverse estimate the index
after filtering.
nonumi_sce
#> class: SingleCellExperiment
#> dim: 6 12
#> metadata(0):
#> assays(1): counts
#> rownames(6): GeneA GeneB ... GeneE GeneF
#> rowData names(0):
#> colnames(12): NONUMI01 NONUMI02 ... NONUMI11 NONUMI12
#> colData names(3): condition batch donor
#> reducedDimNames(0):
#> mainExpName: NULL
#> altExpNames(0):
nonumi_fit <- Tweedieverse(
input_features = nonumi_sce,
output = NULL,
assay_name = "counts",
domain = "single_cell",
normalization = "SCRAN",
fixed_effects = c("condition", "batch"),
tweedie_p = NULL,
median_comparison = FALSE
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:30.291321 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:30.292368 INFO::Determining format of input files
#> 2026-09-28 22:22:30.292799 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:30.29421 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:30.294565 INFO::Total samples in data: 12
#> 2026-09-28 22:22:30.294865 INFO::Min samples required with min abundance for a feature not to be filtered: 1.200000
#> 2026-09-28 22:22:30.295404 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:30.295766 INFO::Filtered feature names:
#> 2026-09-28 22:22:30.296299 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:30.296653 INFO::Filtered feature names:
#> 2026-09-28 22:22:30.361001 INFO::Offset size factors computed with normalization: SCRAN
#> 2026-09-28 22:22:30.362877 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:30.363382 INFO::Formula for fixed effects: expr ~ condition + batch
#> 2026-09-28 22:22:30.363734 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:30.366068 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:30.381649 INFO::Fitting model to feature number 1, GeneA
#> 2026-09-28 22:22:30.441818 INFO::Fitting model to feature number 2, GeneB
#> 2026-09-28 22:22:30.509487 INFO::Fitting model to feature number 3, GeneC
#> 2026-09-28 22:22:30.542335 INFO::Fitting model to feature number 4, GeneD
#> 2026-09-28 22:22:30.581379 INFO::Fitting model to feature number 5, GeneE
#> 2026-09-28 22:22:30.609962 INFO::Fitting model to feature number 6, GeneF
#> 2026-09-28 22:22:30.656899 INFO::Counting prevalence for each feature
head(nonumi_fit)
#> feature metadata value coef stderr pval qval
#> 1 GeneB condition stimulated -0.6126147 0.03385356 2.189998e-08 2.145964e-07
#> 2 GeneE condition stimulated 1.3174037 0.07698921 3.576607e-08 2.145964e-07
#> 3 GeneA condition stimulated 0.5157527 0.03600941 1.683966e-07 6.735866e-07
#> 4 GeneF condition stimulated -0.1540168 0.02942856 5.393497e-04 1.618049e-03
#> 5 GeneC condition stimulated 2.5139682 0.53577336 1.132512e-03 2.718029e-03
#> base.model tweedie.index N N.not.zero percent.zero
#> 1 CPLM 1.010 12 12 0
#> 2 CPLM 1.644 12 12 0
#> 3 CPLM 1.656 12 12 0
#> 4 CPLM 1.010 12 12 0
#> 5 CPLM 1.010 12 8 33Pseudobulk analysis aggregates single-cell counts to sample-level
counts before modeling. This is useful when the biological replicate,
rather than the cell, is the experimental unit. The muscat
package provides a Bioconductor-native workflow for this aggregation.
Its aggregateData() function summarizes a
SingleCellExperiment by variables in colData,
commonly cluster_id and sample_id.
This example follows the workflow used by the muscat
vignette: multi-subject single-cell data are first represented as a
prepared SingleCellExperiment, cells are aggregated into
cluster-by-sample pseudobulk counts, and the resulting sample-level
count matrix is analyzed with Tweedieverse(). The
package-provided example_sce object contains the required
cluster_id, sample_id, and
group_id fields and is a lightweight version of the Kang et
al. interferon stimulation example used throughout
muscat.
data("example_sce", package = "muscat")
pseudobulk_input <- example_sce[seq_len(100), ]
pseudobulk_sce <- muscat::aggregateData(
pseudobulk_input,
assay = "counts",
fun = "sum",
by = c("cluster_id", "sample_id"),
verbose = FALSE
)
SummarizedExperiment::assayNames(pseudobulk_sce)
#> [1] "B cells" "CD14+ Monocytes" "CD4 T cells"
#> [4] "CD8 T cells" "FCGR3A+ Monocytes"
pseudobulk_assay <- "CD4 T cells"
pseudobulk_counts <- SummarizedExperiment::assay(
pseudobulk_sce,
pseudobulk_assay
)
top_pseudobulk_genes <- order(
rowSums(pseudobulk_counts),
decreasing = TRUE
)[seq_len(20)]
pseudobulk_metadata <- unique(as.data.frame(
SummarizedExperiment::colData(pseudobulk_input)[, c("sample_id", "group_id")]
))
rownames(pseudobulk_metadata) <- pseudobulk_metadata$sample_id
pseudobulk_metadata <- pseudobulk_metadata[
colnames(pseudobulk_counts),
,
drop = FALSE
]
pseudobulk_metadata$condition <- pseudobulk_metadata$group_id
pseudobulk_se <- SummarizedExperiment::SummarizedExperiment(
assays = list(counts = pseudobulk_counts[top_pseudobulk_genes, , drop = FALSE]),
colData = pseudobulk_metadata
)
pseudobulk_fit <- Tweedieverse(
input_features = pseudobulk_se,
output = NULL,
assay_name = "counts",
domain = "bulk_rnaseq",
fixed_effects = c("condition"),
tweedie_p = 1,
median_comparison = FALSE,
max_significance = 1
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:30.960539 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:30.961448 INFO::Determining format of input files
#> 2026-09-28 22:22:30.961793 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:30.963142 INFO::Factor detected for categorial metadata 'condition'. Provide a reference argument or manually set factor ordering to change reference level.
#> 2026-09-28 22:22:30.963438 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:30.963685 INFO::Total samples in data: 4
#> 2026-09-28 22:22:30.963914 INFO::Min samples required with min abundance for a feature not to be filtered: 0.400000
#> 2026-09-28 22:22:30.964477 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:30.964769 INFO::Filtered feature names:
#> 2026-09-28 22:22:30.965368 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:30.96566 INFO::Filtered feature names:
#> calcNormFactors has been renamed to normLibSizes
#> 2026-09-28 22:22:30.969479 INFO::Offset size factors computed with normalization: TMM
#> 2026-09-28 22:22:30.970737 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:30.971115 INFO::Formula for fixed effects: expr ~ condition
#> 2026-09-28 22:22:30.97141 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:30.973247 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:30.985653 INFO::Fitting model to feature number 1, RPL11
#> 2026-09-28 22:22:30.98825 INFO::Fitting model to feature number 2, RPS27
#> 2026-09-28 22:22:30.990294 INFO::Fitting model to feature number 3, RPS8
#> 2026-09-28 22:22:30.992166 INFO::Fitting model to feature number 4, ISG15
#> 2026-09-28 22:22:30.994145 INFO::Fitting model to feature number 5, RPL5
#> 2026-09-28 22:22:30.996068 INFO::Fitting model to feature number 6, IFI6
#> 2026-09-28 22:22:30.998041 INFO::Fitting model to feature number 7, SH3BGRL3
#> 2026-09-28 22:22:30.999861 INFO::Fitting model to feature number 8, YBX1
#> 2026-09-28 22:22:31.001734 INFO::Fitting model to feature number 9, S100A11
#> 2026-09-28 22:22:31.003576 INFO::Fitting model to feature number 10, ENO1
#> 2026-09-28 22:22:31.005423 INFO::Fitting model to feature number 11, CD53
#> 2026-09-28 22:22:31.007287 INFO::Fitting model to feature number 12, S100A10
#> 2026-09-28 22:22:31.009071 INFO::Fitting model to feature number 13, S100A6
#> 2026-09-28 22:22:31.010867 INFO::Fitting model to feature number 14, TMEM59
#> 2026-09-28 22:22:31.012695 INFO::Fitting model to feature number 15, CD2
#> 2026-09-28 22:22:31.014555 INFO::Fitting model to feature number 16, NDUFS5
#> 2026-09-28 22:22:31.016346 INFO::Fitting model to feature number 17, CD52
#> 2026-09-28 22:22:31.018163 INFO::Fitting model to feature number 18, LAPTM5
#> 2026-09-28 22:22:31.01993 INFO::Fitting model to feature number 19, SYF2
#> 2026-09-28 22:22:31.021747 INFO::Fitting model to feature number 20, RPL22
#> 2026-09-28 22:22:31.028657 INFO::Counting prevalence for each feature
head(pseudobulk_fit)
#> feature metadata value coef stderr pval qval
#> 1 IFI6 condition stim 3.8839623 0.27615612 0.005017428 0.1003486
#> 2 ISG15 condition stim 3.8805689 0.57303845 0.021117751 0.2111775
#> 3 RPL5 condition stim -0.4307140 0.12650534 0.076498419 0.3059937
#> 4 ENO1 condition stim 0.4636025 0.13026340 0.070681909 0.3059937
#> 5 RPL22 condition stim -0.2870616 0.07274647 0.058629532 0.3059937
#> 6 SH3BGRL3 condition stim -0.1499796 0.07721882 0.191591787 0.4414440
#> base.model tweedie.index N N.not.zero percent.zero
#> 1 Tweedie GLM 1 4 4 0
#> 2 Tweedie GLM 1 4 4 0
#> 3 Tweedie GLM 1 4 4 0
#> 4 Tweedie GLM 1 4 4 0
#> 5 Tweedie GLM 1 4 4 0
#> 6 Tweedie GLM 1 4 4 0For bulk RNA-seq data, Tweedieverse() accepts a
SummarizedExperiment. This example uses the Bioconductor
airway dataset, which contains RNA-seq read counts from
airway smooth muscle cells treated with dexamethasone.
The full airway object contains tens of thousands of
genes. To keep the vignette fast, the example keeps a small set of
high-count genes while preserving the original airway
sample metadata. Bulk RNA-seq uses TMM as the default
normalization strategy in Tweedieverse; this computes
library-size factors for the model offset without modifying the count
assay.
data("airway", package = "airway")
airway_counts <- SummarizedExperiment::assay(airway, "counts")
top_airway_genes <- order(rowSums(airway_counts), decreasing = TRUE)[seq_len(20)]
bulk_rnaseq_se <- SummarizedExperiment::SummarizedExperiment(
assays = list(counts = airway_counts[top_airway_genes, , drop = FALSE]),
colData = SummarizedExperiment::colData(airway)
)
bulk_rnaseq_se
#> class: SummarizedExperiment
#> dim: 20 8
#> metadata(0):
#> assays(1): counts
#> rownames(20): ENSG00000115414 ENSG00000011465 ... ENSG00000137801
#> ENSG00000142173
#> rowData names(0):
#> colnames(8): SRR1039508 SRR1039509 ... SRR1039520 SRR1039521
#> colData names(9): SampleName cell ... Sample BioSample
bulk_rnaseq_fit <- Tweedieverse(
input_features = bulk_rnaseq_se,
output = NULL,
assay_name = "counts",
domain = "bulk_rnaseq",
fixed_effects = c("dex"),
tweedie_p = 1,
median_comparison = FALSE,
max_significance = 1
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:33.213149 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:33.213992 INFO::Determining format of input files
#> 2026-09-28 22:22:33.214346 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:33.215755 INFO::Factor detected for categorial metadata 'dex'. Provide a reference argument or manually set factor ordering to change reference level.
#> 2026-09-28 22:22:33.216029 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:33.216288 INFO::Total samples in data: 8
#> 2026-09-28 22:22:33.216511 INFO::Min samples required with min abundance for a feature not to be filtered: 0.800000
#> 2026-09-28 22:22:33.217014 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:33.217327 INFO::Filtered feature names:
#> 2026-09-28 22:22:33.217874 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:33.218156 INFO::Filtered feature names:
#> calcNormFactors has been renamed to normLibSizes
#> 2026-09-28 22:22:33.22085 INFO::Offset size factors computed with normalization: TMM
#> 2026-09-28 22:22:33.223284 INFO::Excluded metadata with entropy less or equal to 0: albut
#> 2026-09-28 22:22:33.223665 INFO::Formula for fixed effects: expr ~ dex
#> 2026-09-28 22:22:33.223948 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:33.225769 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:33.237626 INFO::Fitting model to feature number 1, ENSG00000115414
#> 2026-09-28 22:22:33.240154 INFO::Fitting model to feature number 2, ENSG00000011465
#> 2026-09-28 22:22:33.242089 INFO::Fitting model to feature number 3, ENSG00000198804
#> 2026-09-28 22:22:33.243942 INFO::Fitting model to feature number 4, ENSG00000156508
#> 2026-09-28 22:22:33.245768 INFO::Fitting model to feature number 5, ENSG00000164692
#> 2026-09-28 22:22:33.247589 INFO::Fitting model to feature number 6, ENSG00000116260
#> 2026-09-28 22:22:33.249372 INFO::Fitting model to feature number 7, ENSG00000087086
#> 2026-09-28 22:22:33.251164 INFO::Fitting model to feature number 8, ENSG00000115461
#> 2026-09-28 22:22:33.25297 INFO::Fitting model to feature number 9, ENSG00000077942
#> 2026-09-28 22:22:33.254708 INFO::Fitting model to feature number 10, ENSG00000166923
#> 2026-09-28 22:22:33.256556 INFO::Fitting model to feature number 11, ENSG00000087245
#> 2026-09-28 22:22:33.258349 INFO::Fitting model to feature number 12, ENSG00000091986
#> 2026-09-28 22:22:33.260095 INFO::Fitting model to feature number 13, ENSG00000103888
#> 2026-09-28 22:22:33.261901 INFO::Fitting model to feature number 14, ENSG00000210082
#> 2026-09-28 22:22:33.263633 INFO::Fitting model to feature number 15, ENSG00000124942
#> 2026-09-28 22:22:33.265377 INFO::Fitting model to feature number 16, ENSG00000108821
#> 2026-09-28 22:22:33.267184 INFO::Fitting model to feature number 17, ENSG00000196205
#> 2026-09-28 22:22:33.268941 INFO::Fitting model to feature number 18, ENSG00000150991
#> 2026-09-28 22:22:33.270748 INFO::Fitting model to feature number 19, ENSG00000137801
#> 2026-09-28 22:22:33.272574 INFO::Fitting model to feature number 20, ENSG00000142173
#> 2026-09-28 22:22:33.279066 INFO::Counting prevalence for each feature
head(bulk_rnaseq_fit)
#> feature metadata value coef stderr pval qval
#> 1 ENSG00000108821 dex untrt 0.8287049 0.08821054 8.260095e-05 0.001652019
#> 2 ENSG00000137801 dex untrt -1.3972630 0.23161755 9.374555e-04 0.009374555
#> 3 ENSG00000124942 dex untrt -0.2768458 0.07273457 8.902531e-03 0.059350204
#> 4 ENSG00000164692 dex untrt 0.2559281 0.08206528 2.062187e-02 0.095415880
#> 5 ENSG00000150991 dex untrt -0.2134059 0.07101706 2.385397e-02 0.095415880
#> 6 ENSG00000087086 dex untrt 0.3890967 0.14496888 3.634185e-02 0.121139513
#> base.model tweedie.index N N.not.zero percent.zero
#> 1 Tweedie GLM 1 8 8 0
#> 2 Tweedie GLM 1 8 8 0
#> 3 Tweedie GLM 1 8 8 0
#> 4 Tweedie GLM 1 8 8 0
#> 5 Tweedie GLM 1 8 8 0
#> 6 Tweedie GLM 1 8 8 0For multi-omics data, Tweedieverse() accepts a
MultiAssayExperiment. The packaged Franzosa et al. example
is stored as an .RData file in inst/extdata.
It contains matched microbiome species abundances and metabolite
profiles. Each omics layer is analyzed independently, and per-omics
parameters can be supplied as named vectors or lists.
franzosa_example_mae
#> A MultiAssayExperiment object of 2 listed
#> experiments with user-defined names and respective classes.
#> Containing an ExperimentList class object of length 2:
#> [1] taxaAbundanceProfile: SummarizedExperiment with 4 rows and 16 columns
#> [2] harmonizedMetaboliteProfile: SummarizedExperiment with 4 rows and 16 columns
#> Functionality:
#> experiments() - obtain the ExperimentList instance
#> colData() - the primary/phenotype DataFrame
#> sampleMap() - the sample coordination DataFrame
#> `$`, `[`, `[[` - extract colData columns, subset, or experiment
#> *Format() - convert into a long or wide DataFrame
#> assays() - convert ExperimentList to a SimpleList of matrices
#> exportClass() - save data to flat files
multi_fit <- Tweedieverse(
input_features = franzosa_example_mae,
output = NULL,
assay_name = "abundance",
domain = c(
taxaAbundanceProfile = "microbiome",
harmonizedMetaboliteProfile = "custom"
),
normalization = c(
taxaAbundanceProfile = "NONE",
harmonizedMetaboliteProfile = "MEDIAN"
),
fixed_effects = c("diseaseCat_clean", "Age"),
tweedie_p = list(
taxaAbundanceProfile = NULL,
harmonizedMetaboliteProfile = NULL
),
median_comparison = FALSE,
max_significance = 1
)
#> The specified assay has been extracted
#> 2026-09-28 22:22:33.314659 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:33.315502 INFO::Determining format of input files
#> 2026-09-28 22:22:33.315837 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:33.316867 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:33.31719 INFO::Total samples in data: 16
#> 2026-09-28 22:22:33.317451 INFO::Min samples required with min abundance for a feature not to be filtered: 1.600000
#> 2026-09-28 22:22:33.317852 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:33.318216 INFO::Filtered feature names:
#> 2026-09-28 22:22:33.318654 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:33.318928 INFO::Filtered feature names:
#> 2026-09-28 22:22:33.319233 INFO::Offset size factors computed with normalization: NONE
#> 2026-09-28 22:22:33.320294 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:33.320654 INFO::Formula for fixed effects: expr ~ diseaseCat_clean + Age
#> 2026-09-28 22:22:33.320926 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:33.324196 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:33.336991 INFO::Fitting model to feature number 1, Bacteroides_vulgatus
#> 2026-09-28 22:22:33.379011 INFO::Fitting model to feature number 2, Ruminococcus_gnavus
#> 2026-09-28 22:22:33.415004 INFO::Fitting model to feature number 3, Eubacterium_rectale
#> 2026-09-28 22:22:33.473737 INFO::Fitting model to feature number 4, Faecalibacterium_prausnitzii
#> 2026-09-28 22:22:33.523783 INFO::Counting prevalence for each feature
#> The specified assay has been extracted
#> 2026-09-28 22:22:33.530409 INFO::Writing function arguments to log file
#> 2026-09-28 22:22:33.531251 INFO::Determining format of input files
#> 2026-09-28 22:22:33.531597 INFO::Input format is data samples as columns and metadata samples as rows
#> 2026-09-28 22:22:33.532649 INFO::Filter data based on min abundance and min prevalence
#> 2026-09-28 22:22:33.53294 INFO::Total samples in data: 16
#> 2026-09-28 22:22:33.533204 INFO::Min samples required with min abundance for a feature not to be filtered: 1.600000
#> 2026-09-28 22:22:33.533671 INFO::Total filtered features with prevalence-abundance filtering: 0
#> 2026-09-28 22:22:33.533982 INFO::Filtered feature names:
#> 2026-09-28 22:22:33.534425 INFO::Total filtered features with variance filtering: 0
#> 2026-09-28 22:22:33.534702 INFO::Filtered feature names:
#> 2026-09-28 22:22:33.535511 INFO::Offset size factors computed with normalization: MEDIAN
#> 2026-09-28 22:22:33.536617 INFO::Excluded metadata with entropy less or equal to 0:
#> 2026-09-28 22:22:33.536986 INFO::Formula for fixed effects: expr ~ diseaseCat_clean + Age
#> 2026-09-28 22:22:33.53732 INFO::Applying z-score to standardize continuous metadata
#> 2026-09-28 22:22:33.539603 INFO::Running selected analysis method: CPLM
#> 2026-09-28 22:22:33.551031 INFO::Fitting model to feature number 1, C54:3 TAG
#> 2026-09-28 22:22:33.59975 INFO::Fitting model to feature number 2, 3-Carboxypropyl trimethylammonium
#> 2026-09-28 22:22:33.689693 INFO::Fitting model to feature number 3, C54:4 TAG
#> 2026-09-28 22:22:33.791343 INFO::Fitting model to feature number 4, C54:5 TAG
#> 2026-09-28 22:22:33.831235 INFO::Counting prevalence for each feature
names(multi_fit)
#> [1] "taxaAbundanceProfile" "harmonizedMetaboliteProfile"
head(multi_fit$taxaAbundanceProfile)
#> omics feature metadata value
#> 1 taxaAbundanceProfile Eubacterium_rectale Age Age
#> 2 taxaAbundanceProfile Eubacterium_rectale diseaseCat_clean Healthy
#> 3 taxaAbundanceProfile Ruminococcus_gnavus Age Age
#> 4 taxaAbundanceProfile Faecalibacterium_prausnitzii diseaseCat_clean Healthy
#> 5 taxaAbundanceProfile Bacteroides_vulgatus Age Age
#> 6 taxaAbundanceProfile Faecalibacterium_prausnitzii Age Age
#> coef stderr pval qval base.model tweedie.index N
#> 1 -2.7824200 0.6723543 0.001166404 0.009331229 CPLM 1.959 16
#> 2 -4.0777339 1.2742125 0.006964917 0.027859668 CPLM 1.959 16
#> 3 1.6620788 0.8757181 0.080128960 0.160257919 CPLM 1.990 16
#> 4 2.2109011 1.1622965 0.079534481 0.160257919 CPLM 1.869 16
#> 5 0.5819870 0.5659511 0.322546636 0.446136889 CPLM 1.951 16
#> 6 0.5975592 0.5963310 0.334602667 0.446136889 CPLM 1.869 16
#> N.not.zero percent.zero
#> 1 16 0
#> 2 16 0
#> 3 16 0
#> 4 16 0
#> 5 16 0
#> 6 16 0
head(multi_fit$harmonizedMetaboliteProfile)
#> omics feature
#> 1 harmonizedMetaboliteProfile C54:3 TAG
#> 2 harmonizedMetaboliteProfile 3-Carboxypropyl trimethylammonium
#> 3 harmonizedMetaboliteProfile C54:5 TAG
#> 4 harmonizedMetaboliteProfile C54:3 TAG
#> 5 harmonizedMetaboliteProfile C54:5 TAG
#> 6 harmonizedMetaboliteProfile 3-Carboxypropyl trimethylammonium
#> metadata value coef stderr pval qval base.model
#> 1 Age Age 0.4576428 0.2540086 0.09482193 0.2528585 CPLM
#> 2 diseaseCat_clean Healthy -3.6533834 1.9329772 0.08125656 0.2528585 CPLM
#> 3 Age Age -0.3291460 0.1750951 0.08273022 0.2528585 CPLM
#> 4 diseaseCat_clean Healthy 0.5723166 0.4908363 0.26456010 0.5291202 CPLM
#> 5 diseaseCat_clean Healthy -0.2784898 0.3378094 0.42458092 0.6793295 CPLM
#> 6 Age Age -0.4177465 1.0004811 0.68309312 0.6830931 CPLM
#> tweedie.index N N.not.zero percent.zero
#> 1 1.99 16 16 0
#> 2 1.99 16 16 0
#> 3 1.99 16 16 0
#> 4 1.99 16 16 0
#> 5 1.99 16 16 0
#> 6 1.99 16 16 0sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: Etc/UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats4 stats graphics grDevices utils datasets methods
#> [8] base
#>
#> other attached packages:
#> [1] airway_1.33.2 muscat_1.27.4
#> [3] MultiAssayExperiment_1.39.1 TreeSummarizedExperiment_2.21.0
#> [5] Biostrings_2.81.9 XVector_0.53.0
#> [7] SingleCellExperiment_1.35.2 SummarizedExperiment_1.43.0
#> [9] Biobase_2.73.2 GenomicRanges_1.65.4
#> [11] Seqinfo_1.3.2 IRanges_2.47.5
#> [13] S4Vectors_0.51.10 BiocGenerics_0.59.12
#> [15] generics_0.1.4 MatrixGenerics_1.25.0
#> [17] matrixStats_1.5.0 Tweedieverse_0.99.0
#> [19] BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] splines_4.6.1 bitops_1.1-0 tibble_3.3.1
#> [4] lifecycle_1.0.5 Rdpack_2.6.6 edgeR_4.99.6
#> [7] doParallel_1.0.17 lattice_0.23-1 MASS_7.3-66
#> [10] backports_1.5.1 magrittr_2.0.5 limma_3.99.0
#> [13] sass_0.4.10 rmarkdown_2.32 jquerylib_0.1.4
#> [16] yaml_2.3.12 metapod_1.21.0 otel_0.2.0
#> [19] DBI_1.3.0 buildtools_1.0.0 minqa_1.2.8
#> [22] RColorBrewer_1.1-3 abind_1.4-8 EnvStats_3.1.0
#> [25] cplm_0.7-12.1 glmmTMB_1.1.14 purrr_1.2.2
#> [28] yulab.utils_0.2.5 rappdirs_0.3.4 sandwich_3.1-3
#> [31] circlize_0.4.18 ggrepel_0.9.8 pbkrtest_0.5.5
#> [34] irlba_2.3.7 tidytree_0.4.8 maketools_1.3.2
#> [37] tweedie_3.1.0 dqrng_0.4.1 codetools_0.2-20
#> [40] DelayedArray_0.39.7 scuttle_1.23.2 tidyselect_1.2.1
#> [43] shape_1.4.6.1 farver_2.1.2 lme4_2.0-6
#> [46] ScaledMatrix_1.21.0 viridis_0.6.5 jsonlite_2.0.0
#> [49] GetoptLong_1.1.1 BiocNeighbors_2.7.3 scater_1.41.2
#> [52] iterators_1.0.14 foreach_1.5.2 tools_4.6.1
#> [55] progress_1.2.3 treeio_1.37.1 Rcpp_1.1.2
#> [58] blme_1.0-7 glue_1.8.1 gridExtra_2.3.1
#> [61] SparseArray_1.13.3 xfun_0.61 mgcv_1.9-4
#> [64] dplyr_1.2.1 withr_3.0.3 numDeriv_2016.8-1.1
#> [67] BiocManager_1.30.27 fastmap_1.2.0 boot_1.3-32
#> [70] bluster_1.23.1 caTools_1.18.4 digest_0.6.39
#> [73] rsvd_1.0.5 R6_2.6.1 colorspace_2.1-3
#> [76] gtools_3.9.5 RhpcBLASctl_0.23-42 tidyr_1.3.2
#> [79] variancePartition_1.43.1 corpcor_1.6.10 prettyunits_1.2.0
#> [82] S4Arrays_1.13.1 pkgconfig_2.0.3 gtable_0.3.6
#> [85] ComplexHeatmap_2.29.0 S7_0.2.2 sys_3.4.3
#> [88] remaCor_0.0.20 htmltools_0.5.9 TMB_1.9.25
#> [91] clue_0.3-68 scales_1.4.0 logging_0.10-111
#> [94] png_0.1-9 fANCOVA_0.6-1 reformulas_0.4.4
#> [97] scran_1.41.1 knitr_1.52 reshape2_1.4.5
#> [100] rjson_0.2.23 coda_0.19-4.1 nlme_3.1-171
#> [103] nloptr_2.2.1 cachem_1.1.0 zoo_1.9-1
#> [106] GlobalOptions_0.1.4 stringr_1.6.0 KernSmooth_2.23-27
#> [109] parallel_4.6.1 vipor_0.4.7 pillar_1.11.1
#> [112] grid_4.6.1 vctrs_0.7.3 gplots_3.3.0
#> [115] BiocSingular_1.29.1 beachmat_2.29.3 cluster_2.1.8.3
#> [118] beeswarm_0.4.0 evaluate_1.0.5 mvtnorm_1.4-2
#> [121] cli_3.6.6 locfit_1.5-9.12 compiler_4.6.1
#> [124] rlang_1.3.0 crayon_1.5.3 plyr_1.8.9
#> [127] fs_2.1.0 ggbeeswarm_0.7.3 stringi_1.8.9
#> [130] viridisLite_0.4.3 BiocParallel_1.47.0 lmerTest_3.2-1
#> [133] lazyeval_0.2.3 aod_1.3.3 scrapper_1.7.9
#> [136] Matrix_1.7-6 hms_1.1.4 ggplot2_4.0.3
#> [139] statmod_1.5.2 rbibutils_2.4.1 igraph_2.3.3
#> [142] broom_1.0.13 bslib_0.12.0 biglm_0.9-3
#> [145] ape_5.8-1