--- a
+++ b/man/Tweedieverse.Rd
@@ -0,0 +1,324 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/Tweedieverse.R
+\name{Tweedieverse}
+\alias{Tweedieverse}
+\title{Differential analysis of multi-omics data using Tweedie GLMs}
+\usage{
+Tweedieverse(
+  input_features,
+  input_metadata = NULL,
+  output,
+  assay_name = "counts",
+  abd_threshold = 0,
+  prev_threshold = 0.1,
+  var_threshold = 0,
+  entropy_threshold = 0,
+  base_model = "CPLM",
+  link = "log",
+  fixed_effects = NULL,
+  random_effects = NULL,
+  cutoff_ZSCP = 0.3,
+  criteria_ZACP = "BIC",
+  adjust_offset = TRUE,
+  scale_factor = NULL,
+  max_significance = 0.05,
+  correction = "BH",
+  standardize = TRUE,
+  cores = 1,
+  optimizer = "nlminb",
+  na.action = na.exclude,
+  plot_heatmap = FALSE,
+  plot_scatter = FALSE,
+  heatmap_first_n = 50,
+  reference = NULL
+)
+}
+\arguments{
+\item{input_features}{A tab-delimited input file or an R data frame of features (can be in rows/columns)
+and samples (or cells). Samples are expected to have matching names with \code{input_metadata}. 
+\code{input_features} can also be an object of class \code{SummarizedExperiment} or \code{SingleCellExperiment} 
+that contains the expression or abundance matrix and other metadata; the \code{assays} 
+slot contains the expression or abundance matrix and is named \code{"counts"}.  
+This matrix should have one row for each feature and one sample for each column.  
+The \code{colData} slot should contain a data frame with one row per 
+sample and columns that contain metadata for each sample. 
+Additional information about the experiment can be contained in the
+\code{metadata} slot as a list.}
+
+\item{input_metadata}{A tab-delimited input file or an R data frame of metadata (rows/columns).
+Samples are expected to have matching sample names with \code{input_features}. 
+This file is ignored when \code{input_features} is a \code{SummarizedExperiment} 
+or a \code{SingleCellExperiment} object with \code{colData} containing the same information.}
+
+\item{output}{The output folder to write results.}
+
+\item{assay_name}{If the input is provided as one of the accepted Bioconductor objects 
+(e.g., SummarizedExperiment, RangedSummarizedExperiment, SingleCellExperiment, or TreeSummarizedExperiment), 
+this argument selects the name of the assay slot in the input object that contains the omics measurements.}
+
+\item{abd_threshold}{If prevalence-abundance filtering is desired, only features that are present (or detected)
+in at least \code{prev_threshold} percent of samples at \code{abd_threshold} minimum abundance (read count or proportion)
+are retained. Default value for \code{abd_threshold} is \code{0.0}.
+To disable prevalence-abundance filtering, set \code{abd_threshold = -Inf}.}
+
+\item{prev_threshold}{If prevalence-abundance filtering is desired, only features that are present (or detected)
+in at least \code{prev_threshold} percent of samples at \code{abd_threshold} minimum abundance (read count or proportion)
+are retained. Default value for \code{prev_threshold} is \code{0.1}.}
+
+\item{var_threshold}{If variance filtering is desired, only features that have variances greater than
+\code{var_threshold} are retained. This step is done after the prevalence-abundance filtering.
+Default value for \code{var_threshold} is \code{0.0} (i.e. no variance filtering).}
+
+\item{entropy_threshold}{If entropy-based filtering is desired for metadata, only features that have entropy greater than
+\code{entropy_threshold} are retained. Default value for \code{entropy_threshold} is \code{0.0} (i.e. no entropy filtering).}
+
+\item{base_model}{The per-feature base model. Default is "CPLM". Must be one of "CPLM", "ZICP", "ZSCP", or "ZACP".}
+
+\item{link}{A specification of the GLM link function. Default is "log". Must be one of "log", "identity", "sqrt", or "inverse".}
+
+\item{fixed_effects}{Metadata variable(s) describing the fixed effects coefficients.}
+
+\item{random_effects}{Metadata variable(s) describing the random effects part of the model.}
+
+\item{cutoff_ZSCP}{For \code{base_model = "ZSCP"}, the cutoff to stratify features for
+adaptive ZI modeling based on sparsity (zero-inflation proportion). Default is 0.3. Must be between 0 and 1.}
+
+\item{criteria_ZACP}{For \code{base_model = "ZACP"}, the criteria to select the
+best fitting model per feature.  The possible options are 'AIC' and BIC' (default).
+More criteria will be supported in a future release.}
+
+\item{adjust_offset}{If TRUE (default), an offset term will be included as the logarithm of \code{scale_factor}.}
+
+\item{scale_factor}{Name of the numerical variable containing library size (for non-normalized data) or scale factor
+(for normalized data) across samples to be included as an offset in the base model (when \code{adjust_offset = TRUE}).
+If not found in metadata, defaults to the sample-wise total sums, unless \code{adjust_offset = FALSE}.}
+
+\item{max_significance}{The q-value threshold for significance. Default is 0.05.}
+
+\item{correction}{The correction method for computing the q-value (see \code{\link[stats]{p.adjust}} for options, default is 'BH').}
+
+\item{standardize}{Should continuous metadata be standardized? Default is TRUE. Bypassed for categorical variables.}
+
+\item{cores}{An integer that indicates the number of R processes to run in parallel. Default is 1.}
+
+\item{optimizer}{The optimization routine to be used for estimating the parameters of the Tweedie model.
+Possible choices are \code{"nlminb"} (the default, see \code{\link[stats]{nlminb}}),
+\code{"bobyqa"} (\code{\link[minqa]{bobyqa}}), and \code{"L-BFGS-B"} (\code{\link[stats]{optim}}).
+Ignored for random effects modeling which uses an alternative Template Model Builder (TMB) approach (\code{\link[glmmTMB]{glmmTMB}}).}
+
+\item{na.action}{How to handle missing values? See \code{\link{na.action}}. Default is \code{\link{na.exclude}}.}
+
+\item{plot_heatmap}{Logical. If TRUE (default is FALSE), generate a heatmap of the (top \code{heatmap_first_n}) significant associations.}
+
+\item{plot_scatter}{Logical. If TRUE (default is FALSE), generate scatter/box plots of individual associations.}
+
+\item{heatmap_first_n}{In heatmap, plot top N features with significant associations (default is 50).}
+
+\item{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).}
+}
+\value{
+A data frame containing coefficient estimates, p-values, and q-values (multiplicity-adjusted p-values) are returned.
+}
+\description{
+Fit a per-feature Tweedie generalized linear model to omics features.
+}
+\examples{
+
+\dontrun{
+
+##############################################################################
+# Example 1 - Differential Abundance Analysis of Synthetic Microbiome Counts #
+##############################################################################
+
+#######################################
+# Install and Load Required Libraries #
+#######################################
+
+library(devtools)
+devtools::install_github('biobakery/sparseDOSSA@varyLibSize')
+library(sparseDOSSA)
+library(stringi)
+
+######################
+# Specify Parameters #
+######################
+
+n.microbes <- 200 # Number of Features
+n.samples <- 100 # Number of Samples
+spike.perc <- 0.02 # Percentage of Spiked-in Bugs
+spikeStrength<-"20" # Effect Size
+
+###########################
+# Specify Binary Metadata #
+###########################
+
+n.metadata <- 1
+UserMetadata<-as.matrix(rep(c(0,1), each=n.samples/2))
+UserMetadata<-t(UserMetadata) # Transpose
+
+###################################################
+# Spiked-in Metadata (Which Metadata to Spike-in) #
+###################################################
+
+Metadatafrozenidx<-1
+spikeCount<-as.character(length(Metadatafrozenidx))
+significant_metadata<-paste('Metadata', Metadatafrozenidx, sep='')
+
+#############################################
+# Generate SparseDOSSA Synthetic Abundances #
+#############################################
+
+DD<-sparseDOSSA::sparseDOSSA(number_features = n.microbes,
+number_samples = n.samples,
+UserMetadata=UserMetadata,
+Metadatafrozenidx=Metadatafrozenidx,
+datasetCount = 1,
+spikeCount = spikeCount,
+spikeStrength = spikeStrength,
+noZeroInflate=TRUE,
+percent_spiked=spike.perc,
+seed = 1234)
+
+##############################
+# Gather SparseDOSSA Outputs #
+##############################
+
+sparsedossa_results <- as.data.frame(DD$OTU_count)
+rownames(sparsedossa_results)<-sparsedossa_results$X1
+sparsedossa_results<-sparsedossa_results[-1,-1]
+colnames(sparsedossa_results)<-paste('Sample', 1:ncol(sparsedossa_results), sep='')
+data<-as.matrix(sparsedossa_results[-c((n.metadata+1):(2*n.microbes+n.metadata)),])
+data<-data.matrix(data)
+class(data) <- "numeric"
+truth<-c(unlist(DD$truth))
+truth<-truth[!stri_detect_fixed(truth,":")]
+truth<-truth[(5+n.metadata):length(truth)]
+truth<-as.data.frame(truth)
+significant_features<-truth[seq(1,
+(as.numeric(spikeCount)+1)*(n.microbes*spike.perc), (as.numeric(spikeCount)+1)),]
+significant_features<-as.vector(significant_features)
+
+####################
+# Extract Features #
+####################
+
+features<-as.data.frame(t(data[-c(1:n.metadata),]))
+
+####################
+# Extract Metadata #
+####################
+
+metadata<-as.data.frame(data[1,])
+colnames(metadata)<-rownames(data)[1]
+
+###############################
+# Mark True Positive Features #
+###############################
+
+wh.TP = colnames(features) \%in\% significant_features
+colnames(features)<-paste("Feature", 1:n.microbes, sep = "")
+newname = paste0(colnames(features)[wh.TP], "_TP")
+colnames(features)[wh.TP] <- newname;
+colnames(features)[grep('TP', colnames(features))]
+
+####################
+# Run Tweedieverse #
+###################
+
+###################
+# Default options #
+###################
+
+CPLM <-Tweedieverse(
+features,
+metadata,
+output = './demo_output/CPLM') # Assuming demo_output exists
+
+###############################################
+# User-defined prevalence-abundance filtering #
+###############################################
+
+ZICP<-Tweedieverse(
+features,
+metadata,
+output = './demo_output/ZICP', # Assuming demo_output exists
+base_model = 'ZICP',
+abd_threshold = 0.0,
+prev_threshold = 0.2)
+
+####################################
+# User-defined variance filtering  #
+####################################
+
+sds<-apply(features, 2, sd)
+var_threshold = median(sds)/2
+ZSCP<-Tweedieverse(
+features,
+metadata,
+output = './demo_output/ZSCP', # Assuming demo_output exists
+base_model = 'ZSCP',
+var_threshold = var_threshold)
+
+##################
+# Multiple cores #
+##################
+
+ZACP<-Tweedieverse(
+features,
+metadata,
+output = './demo_output/ZACP', # Assuming demo_output exists
+base_model = 'ZACP',
+cores = 4)
+
+##########################################################################
+# Example 2 - Multivariable Association on HMP2 Longitudinal Microbiomes #
+##########################################################################
+
+######################
+# HMP2 input_features Analysis #
+######################
+
+#############
+# Load input_features #
+#############
+
+library(data.table)
+input_features <- fread("https://raw.githubusercontent.com/biobakery/Maaslin2/master/inst/extdata/HMP2_taxonomy.tsv", sep ="\t")
+input_metadata <-fread("https://raw.githubusercontent.com/biobakery/Maaslin2/master/inst/extdata/HMP2_metadata.tsv", sep ="\t")
+
+###############
+# Format data #
+###############
+
+library(tibble)
+features<- column_to_rownames(input_features, 'ID')
+metadata<- column_to_rownames(input_metadata, 'ID')
+
+#############
+# Fit Model #
+#############
+
+library(Tweedieverse)
+HMP2 <- Tweedieverse(
+features,
+metadata,
+output = './demo_output/HMP2', # Assuming demo_output exists
+fixed_effects = c('diagnosis', 'dysbiosisnonIBD','dysbiosisUC','dysbiosisCD', 'antibiotics', 'age'),
+random_effects = c('site', 'subject'),
+base_model = 'CPLM',
+adjust_offset = FALSE, # No offset as the values are relative abundances
+cores = 8, # Make sure your computer has the capability
+standardize = FALSE,
+reference = c('diagnosis,nonIBD'))
+
+}
+}
+\author{
+Himel Mallick, \email{him4004@med.cornell.edu}
+}
+\keyword{metagenomics,}
+\keyword{microbiome,}
+\keyword{multiomics,}
+\keyword{scRNASeq,}
+\keyword{singlecell}
+\keyword{tweedie,}