#' Differential analysis of multi-omics data using Tweedie GLMs
#'
#' Fit a per-feature Tweedie generalized linear model to omics features.
#' @param 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.
#' @param 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.
#' @param output The output folder to write results.
#' @param 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.
#' @param 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}.
#' @param 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}.
#' @param 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).
#' @param 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).
#' @param base_model The per-feature base model. Default is "CPLM". Must be one of "CPLM", "ZICP", "ZSCP", or "ZACP".
#' @param link A specification of the GLM link function. Default is "log". Must be one of "log", "identity", "sqrt", or "inverse".
#' @param fixed_effects Metadata variable(s) describing the fixed effects coefficients.
#' @param random_effects Metadata variable(s) describing the random effects part of the model.
#' @param 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.
#' @param 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.
#' @param adjust_offset If TRUE (default), an offset term will be included as the logarithm of \code{scale_factor}.
#' @param 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}.
#' @param max_significance The q-value threshold for significance. Default is 0.05.
#' @param correction The correction method for computing the q-value (see \code{\link[stats]{p.adjust}} for options, default is 'BH').
#' @param standardize Should continuous metadata be standardized? Default is TRUE. Bypassed for categorical variables.
#' @param cores An integer that indicates the number of R processes to run in parallel. Default is 1.
#' @param 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}}).
#' @param na.action How to handle missing values? See \code{\link{na.action}}. Default is \code{\link{na.exclude}}.
#' @param plot_heatmap Logical. If TRUE (default is FALSE), generate a heatmap of the (top \code{heatmap_first_n}) significant associations.
#' @param plot_scatter Logical. If TRUE (default is FALSE), generate scatter/box plots of individual associations.
#' @param heatmap_first_n In heatmap, plot top N features with significant associations (default is 50).
#' @param 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).
#'
#' @importFrom grDevices colorRampPalette dev.off jpeg pdf
#' @importFrom stats coef fitted as.formula na.exclude p.adjust plogis sd update
#' @importFrom utils capture.output read.table write.table
#' @importFrom dplyr %>% everything
#' @importFrom parallel clusterExport
#' @return A data frame containing coefficient estimates, p-values, and q-values (multiplicity-adjusted p-values) are returned.
#'
#' @author Himel Mallick, \email{him4004@@med.cornell.edu}
#'
#' @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'))
#'
#' }
#' @keywords microbiome, metagenomics, multiomics, scRNASeq, tweedie, singlecell
#' @export
Tweedieverse <- function(input_features,
input_metadata = NULL,
output,
assay_name = "counts",
abd_threshold = 0.0,
prev_threshold = 0.1,
var_threshold = 0.0,
entropy_threshold = 0.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) {
#################################
# Specify all available options #
#################################
model_choices <- c("CPLM", "ZICP", "ZACP", "ZSCP")
link_choices <- c("log", "identity", "sqrt", "inverse")
criteria_ZACP_choices <- c("AIC", "BIC")
correction_choices <-
c("BH", "holm", "hochberg", "hommel", "bonferroni", "BY")
optimizer_choices <- c("nlminb", "bobyqa", "L-BFGS-B")
#######################################################################
#=====================================================================#
# Read in the data and metadata, create output folder, initialize log #
#=====================================================================#
#######################################################################
#########################################
# Support multiple Bioconductor classes #
#########################################
valid_classes <- c("SummarizedExperiment", "RangedSummarizedExperiment",
"SingleCellExperiment", "TreeSummarizedExperiment")
##############################################################
# Extract features and metadata based on user-provided input #
##############################################################
input_class<-class(input_features)
if (input_class %in% valid_classes) {
data <- extractAssay(input_features, assay_name)
if(is.null(input_metadata)) {
metadata <- data.frame(colData(input_features))
} else{
metadata<-input_metadata
}
} else if (!(is.character(input_features)) & !(is.data.frame(input_features))) {
stop(cat(paste('Input data of class <', class(input_features), '> not supported. Please use SummarizedExperiment, SingleCellExperiment, RangedSummarizedExperiment, TreeSummarizedExperiment, or data.frame')))
} else{
# if a character string then this is a file name, else it
# is a data frame
if (is.character(input_features)) {
data <-
data.frame(
read.delim::fread(input_features, header = TRUE, sep = '\t'),
header = TRUE,
fill = T,
comment.char = "" ,
check.names = F,
row.names = 1
)
if (nrow(input_features) == 1) {
# read again to get row name
data <- read.delim(
input_features,
header = TRUE,
fill = T,
comment.char = "" ,
check.names = F,
row.names = 1
)
}
} else {
data <- input_features
}
if (is.character(input_metadata)) {
input_metadata <-
data.frame(
read.delim::fread(input_metadata, header = TRUE, sep = '\t'),
header = TRUE,
fill = T,
comment.char = "" ,
check.names = F,
row.names = 1
)
if (nrow(input_metadata) == 1) {
input_metadata <- read.delim(
input_metadata,
header = TRUE,
fill = T,
comment.char = "" ,
check.names = F,
row.names = 1
)
}
} else {
metadata <- input_metadata
}
}
# create an output folder and figures folder if it does not exist
if (!file.exists(output)) {
print("Creating output folder")
dir.create(output)
}
#if (plot_heatmap || plot_scatter) {
figures_folder <- file.path(output, "figures")
if (!file.exists(figures_folder)) {
print("Creating output figures folder")
dir.create(figures_folder)
}
#}
# Create log file (write info to stdout and debug level to log file)
# Set level to finest so all log levels are reviewed
log_file <- file.path(output, "Tweedieverse.log")
# Remove log file if already exists (to avoid append)
if (file.exists(log_file)) {
print(paste("Warning: Deleting existing log file:", log_file))
unlink(log_file)
}
logging::basicConfig(level = 'FINEST')
logging::addHandler(logging::writeToFile,
file = log_file, level = "DEBUG")
logging::setLevel(20, logging::getHandler('basic.stdout'))
#####################
# Log the arguments #
#####################
logging::loginfo("Writing function arguments to log file")
logging::logdebug("Function arguments")
if (is.character(input_features)) {
logging::logdebug("Input data file: %s", input_features)
}
if (is.character(input_metadata)) {
logging::logdebug("Input metadata file: %s", input_metadata)
}
logging::logdebug("Output folder: %s", output)
logging::logdebug("Abundance threshold: %f", abd_threshold)
logging::logdebug("Prevalence threshold: %f", prev_threshold)
logging::logdebug("Variance threshold: %f", var_threshold)
logging::logdebug("Base model: %s", base_model)
logging::logdebug("Link function: %s", link)
logging::logdebug("Fixed effects: %s", fixed_effects)
logging::logdebug("Random effects: %s", random_effects)
logging::logdebug("ZSCP cutoff: %f", cutoff_ZSCP)
logging::logdebug("ZACP criteria: %s", criteria_ZACP)
logging::logdebug("Offset adjustment: %s", adjust_offset)
logging::logdebug("Scale factor: %s", scale_factor)
logging::logdebug("Max significance: %f", max_significance)
logging::logdebug("Correction method: %s", correction)
logging::logdebug("Standardize: %s", standardize)
logging::logdebug("Cores: %d", cores)
logging::logdebug("Optimization routine: %s", optimizer)
#######################################
# Check if valid options are selected #
#######################################
# Check if the selected link is valid
if (!link %in% link_choices) {
option_not_valid_error("Please select a link from the list of available options",
toString(link_choices))
}
# Check if the selected base_model is valid
if (!base_model %in% model_choices) {
option_not_valid_error(
paste(
"Please select an analysis method",
"from the list of available options"
),
toString(model_choices)
)
}
# Check if the selected criteria_ZACP is valid
if (!criteria_ZACP %in% criteria_ZACP_choices) {
option_not_valid_error(
paste(
"Please select a criteria",
"from the list of available options"
),
toString(criteria_ZACP_choices)
)
}
# Check if the selected correction is valid
if (!correction %in% correction_choices) {
option_not_valid_error(
paste(
"Please select a correction method",
"from the list of available options"
),
toString(correction_choices)
)
}
# Check if the selected optimizer is valid
if (!optimizer %in% optimizer_choices) {
option_not_valid_error(
paste(
"Please select an optimizer method",
"from the list of available options"
),
toString(correction_choices)
)
}
############################################################
# Check if the selected numerical options are within range #
############################################################
prop_options <- c(prev_threshold, cutoff_ZSCP, max_significance)
if (any(prop_options < 0) || any(prop_options > 1)) {
stop(
paste(
"One of the following is outside [0, 1]:",
"prev_threshold, cutoff_ZSCP, max_significance"
)
)
}
###############################################################
# Determine orientation of data in input and reorder to match #
###############################################################
logging::loginfo("Determining format of input files")
samples_row_row <- intersect(rownames(data), rownames(metadata))
if (length(samples_row_row) > 0) {
# this is the expected formatting so do not modify data frames
logging::loginfo(paste(
"Input format is data samples",
"as rows and metadata samples as rows"
))
} else {
samples_column_row <- intersect(colnames(data), rownames(metadata))
if (length(samples_column_row) > 0) {
logging::loginfo(paste(
"Input format is data samples",
"as columns and metadata samples as rows"
))
# transpose data frame so samples are rows
data <- type.convert(as.data.frame(t(data)))
logging::logdebug("linked data so samples are rows")
} else {
samples_column_column <-
intersect(colnames(data), colnames(metadata))
if (length(samples_column_column) > 0) {
logging::loginfo(
paste(
"Input format is data samples",
"as columns and metadata samples as columns"
)
)
data <- type.convert(as.data.frame(t(data)))
metadata <- type.convert(as.data.frame(t(metadata)))
logging::logdebug("linked data and metadata so samples are rows")
} else {
samples_row_column <-
intersect(rownames(data), colnames(metadata))
if (length(samples_row_column) > 0) {
logging::loginfo(
paste(
"Input format is data samples",
"as rows and metadata samples as columns"
)
)
metadata <- type.convert(as.data.frame(t(metadata)))
logging::logdebug("linked metadata so samples are rows")
} else {
logging::logerror(
paste(
"Unable to find samples in data and",
"metadata files.",
"Rows/columns do not match."
)
)
logging::logdebug("input_features rows: %s",
paste(rownames(data), collapse = ","))
logging::logdebug("input_features columns: %s",
paste(colnames(data), collapse = ","))
logging::logdebug("Metadata rows: %s",
paste(rownames(metadata), collapse = ","))
logging::logdebug("Metadata columns: %s",
paste(colnames(data), collapse = ","))
stop()
}
}
}
}
# Replace unexpected characters in feature names
# colnames(data) <- make.names(colnames(data))
# Check for samples without metadata
extra_feature_samples <-
setdiff(rownames(data), rownames(metadata))
if (length(extra_feature_samples) > 0)
logging::logdebug(
paste(
"The following samples were found",
"to have features but no metadata.",
"They will be removed. %s"
),
paste(extra_feature_samples, collapse = ",")
)
# Check for metadata samples without features
extra_metadata_samples <-
setdiff(rownames(metadata), rownames(data))
if (length(extra_metadata_samples) > 0)
logging::logdebug(
paste(
"The following samples were found",
"to have metadata but no features.",
"They will be removed. %s"
),
paste(extra_metadata_samples, collapse = ",")
)
# Get a set of the samples with both metadata and features
intersect_samples <- intersect(rownames(data), rownames(metadata))
logging::logdebug(
"A total of %s samples were found in both the data and metadata",
length(intersect_samples)
)
# Now order both data and metadata with the same sample ordering
logging::logdebug("Reordering data/metadata to use same sample ordering")
data <- data[intersect_samples, , drop = FALSE]
metadata <- metadata[intersect_samples, , drop = FALSE]
########################################################################
# Assign reference values to categorical metadata (fixed effects only) #
########################################################################
if (is.null(reference)) {
reference <- ","
}
split_reference <- unlist(strsplit(reference, "[,;]"))
# for each fixed effect, check that a reference level has been set if necessary: number of levels > 2 and metadata isn't already an ordered factor
for (i in fixed_effects) {
# don't check for or require reference levels for numeric metadata
if (is.numeric(metadata[,i])) {
next
}
# respect ordering if a factor is explicitly passed in with no reference set
if (is.factor(metadata[,i]) && !(i %in% split_reference)) {
logging::loginfo(paste("Factor detected for categorial metadata '",
i, "'. Provide a reference argument or manually set factor ordering to change reference level.", sep=""))
next
}
# set metadata as a factor (ordered alphabetically)
metadata[,i] <- as.factor(metadata[,i])
mlevels <- levels(metadata[,i])
# get reference level for variable being considered, returns NA if not found
ref <- split_reference[match(i, split_reference)+1]
# if metadata has 2 levels, allow but don't require setting reference level, otherwise require it
if ((length(mlevels) == 2)) {
if(!is.na(ref)) {
metadata[,i] = relevel(metadata[,i], ref = ref)
}
} else if (length(mlevels) > 2) {
if (!is.na(ref)) {
metadata[,i] = relevel(metadata[,i], ref = ref)
} else {
stop(paste("Please provide the reference for the variable '",
i, "' which includes more than 2 levels: ",
paste(as.character(mlevels), collapse=", "), ".", sep=""))
}
} else {
stop("Provided categorical metadata has fewer than 2 unique, non-NA values.")
}
}
#########################################################
# Non-specific filtering based on user-provided options #
#########################################################
unfiltered_data <- data
unfiltered_metadata <- metadata
# require at least total samples * min prevalence values
# for each feature to be greater than min abundance
logging::loginfo("Filter data based on min abundance and min prevalence")
total_samples <- nrow(unfiltered_data)
logging::loginfo("Total samples in data: %d", total_samples)
min_samples <- total_samples * prev_threshold
logging::loginfo(
paste(
"Min samples required with min abundance",
"for a feature not to be filtered: %f"
),
min_samples
)
# Filter by abundance using zero as value for NAs
data_zeros <- unfiltered_data
data_zeros[is.na(data_zeros)] <- 0
filtered_data <-
unfiltered_data[,
colSums(data_zeros > abd_threshold) > min_samples,
drop = FALSE]
total_filtered_features <-
ncol(unfiltered_data) - ncol(filtered_data)
logging::loginfo(
"Total filtered features with prevalence-abundance filtering: %d",
total_filtered_features
)
filtered_feature_names <-
setdiff(names(unfiltered_data), names(filtered_data))
logging::loginfo("Filtered feature names: %s",
toString(filtered_feature_names))
#################################
# Filter data based on variance #
#################################
sds <- apply(filtered_data, 2, na.rm = T, sd)
final_features <-
filtered_data[, which(sds > var_threshold), drop = FALSE]
total_filtered_features_var <-
ncol(filtered_data) - ncol(final_features)
logging::loginfo("Total filtered features with variance filtering: %d",
total_filtered_features_var)
filtered_feature_names_var <-
setdiff(names(filtered_data), names(final_features))
logging::loginfo("Filtered feature names: %s",
toString(filtered_feature_names_var))
########################################################################
# Set the scale factor to rowsum if not provided and create a modified #
# metadata table without the scale factor variable (if present) ########
########################################################################
if (is.null(scale_factor)) {
offset <- rowSums(final_features)
} else {
if (!scale_factor %in% colnames(unfiltered_metadata)) {
stop(
paste(
"The specified scale_factor variable is not present in the metadata table:\n",
scale_factor
)
)
} else {
offset <- unfiltered_metadata[, scale_factor]
unfiltered_metadata <-
dplyr::select(unfiltered_metadata,-scale_factor)
}
}
####################################
# Filter metadata based on entropy #
####################################
# Reduce metadata to only include those pass entropy threshold
temp_filtered_metadata <- unfiltered_metadata[, apply(unfiltered_metadata, 2, entropy) > entropy_threshold, drop=F]
excluded_metadata <- setdiff(colnames(unfiltered_metadata), colnames(temp_filtered_metadata))
logging::loginfo(
paste(
"Excluded metadata with",
"entropy less or equal to %s: %s"
),
entropy_threshold, paste(excluded_metadata, collapse = ",")
)
filtered_metadata <- temp_filtered_metadata
###############################################
# Compute the formula based on the user input #
###############################################
#####################
# Determine formula #
#####################
random_effects_formula <- NULL
# Use all metadata if no fixed effects are provided
if (is.null(fixed_effects)) {
fixed_effects <- colnames(filtered_metadata)
} else {
fixed_effects <- unlist(strsplit(fixed_effects, ",", fixed = TRUE))
# remove any fixed effects not found in metadata names
to_remove <- setdiff(fixed_effects, colnames(filtered_metadata))
if (length(to_remove) > 0)
logging::logwarn(
paste(
"Feature name not found in metadata",
"so not applied to formula as fixed effect: %s"
),
paste(to_remove, collapse = " , ")
)
fixed_effects <- setdiff(fixed_effects, to_remove)
if (length(fixed_effects) == 0) {
logging::logerror("No fixed effects included in formula.")
stop()
}
}
if (!is.null(random_effects)) {
random_effects <-
unlist(strsplit(random_effects, ",", fixed = TRUE))
# subtract random effects from fixed effects
fixed_effects <- setdiff(fixed_effects, random_effects)
# remove any random effects not found in metadata
to_remove <-
setdiff(random_effects, colnames(filtered_metadata))
if (length(to_remove) > 0)
logging::logwarn(
paste(
"Feature name not found in metadata",
"so not applied to formula as random effect: %s"
),
paste(to_remove, collapse = " , ")
)
random_effects <- setdiff(random_effects, to_remove)
# create formula
if (length(random_effects) > 0) {
random_effects_formula_text <-
paste("expr ~ (1 | ",
paste(
random_effects,
")",
sep = '',
collapse = " + (1 | "
),
sep = '')
logging::loginfo("Formula for random effects: %s",
random_effects_formula_text)
random_effects_formula <-
tryCatch(
as.formula(random_effects_formula_text),
error = function(e)
stop(
paste(
"Invalid formula for random effects: ",
random_effects_formula_text
)
)
)
}
}
# Reduce metadata to only include fixed/random effects in formula
effects_names <- union(fixed_effects, random_effects)
filtered_metadata <-
filtered_metadata[, effects_names, drop = FALSE]
# Create the fixed effects formula text
formula_text <-
paste("expr ~ ", paste(fixed_effects, collapse = " + "))
logging::loginfo("Formula for fixed effects: %s", formula_text)
formula <-
tryCatch(
as.formula(formula_text),
error = function(e)
stop(
paste(
"Invalid formula.",
"Please provide a different formula: ",
formula_text
)
)
)
#############################################################
# Standardize metadata (excpet the offset variable), if set #
#############################################################
if (standardize) {
logging::loginfo("Applying z-score to standardize continuous metadata")
filtered_metadata <-
filtered_metadata %>% dplyr::mutate_if(is.numeric, scale)
} else {
logging::loginfo("Bypass z-score application to metadata")
}
##################################
# Merge metadata and offset back #
##################################
final_metadata <- as.data.frame(cbind.data.frame(filtered_metadata, offset))
##############################################################
# Apply the base model to the filtered data with user inputs #
##############################################################
logging::loginfo("Running selected analysis method: %s", base_model)
fit_data <- fit.Tweedieverse(
features = final_features,
metadata = final_metadata,
base_model = base_model,
link = link,
formula = formula,
random_effects_formula = random_effects_formula,
cutoff_ZSCP = cutoff_ZSCP,
criteria_ZACP = criteria_ZACP,
adjust_offset = adjust_offset,
correction = correction,
cores = cores,
optimizer = optimizer,
na.action = na.action
)
###################################################
# Count the N and Zero-inflation for each feature #
###################################################
logging::loginfo("Counting prevalence for each feature")
try(fit_data$results$N <-
apply(
fit_data$results,
1,
FUN = function(x)
length(final_features[, x[1]])
))
try(fit_data$results$N.not.zero <-
apply(
fit_data$results,
1,
FUN = function(x)
length(which(final_features[, x[1]] > 0))
))
try(fit_data$results$percent.zero <-
apply(
fit_data$results,
1,
FUN = function(x)
round(mean(final_features[, x[1]] == 0, na.rm = TRUE), 2) *
100
))
#########################
# Write out the results #
#########################
results_file <- file.path(output, "all_results.tsv")
logging::loginfo("Writing all results to file (ordered by increasing q-values): %s",
results_file)
ordered_results <-
fit_data$results[order(fit_data$results$qval),]
ordered_results <-
ordered_results[!is.na(ordered_results$qval),] # Remove NA's
ordered_results <-
dplyr::select(
ordered_results,
c(
'feature',
'metadata',
'value',
"coef",
"stderr",
"pval",
"qval"
),
everything()
)
write.table(
ordered_results,
file = results_file,
sep = "\t",
quote = FALSE,
row.names = FALSE
)
# Write results passing threshold to file
# (removing any that are NA for the q-value)
significant_results <-
ordered_results[ordered_results$qval <= max_significance,]
significant_results_file <-
file.path(output, "significant_results.tsv")
logging::loginfo(
paste(
"Writing the significant results",
"(those which are less than or equal to the threshold",
"of %f ) to file (ordered by increasing q-values): %s"
),
max_significance,
significant_results_file
)
write.table(
significant_results,
file = significant_results_file,
sep = "\t",
quote = FALSE,
row.names = FALSE
)
#######################################################
# Create visualizations for results passing threshold #
#######################################################
logging::loginfo("Writing Tweedie inxed plot to file: %s",
output)
tryCatch({
tweedie_index_plot(ordered_results, figures_folder)
}, error = function(err) {
logging::logerror("Unable to do make a Tweedie inxed plot of results!!!")
logging::logerror(err)
# dev.off()
})
if (plot_heatmap & length(unique(significant_results_file["metdata"])) > 1) {
heatmap_file <- file.path(output, "Tweedieverse_Heatmap.pdf")
logging::loginfo("Writing heatmap of significant results to file: %s",
heatmap_file)
tryCatch({
save_heatmap(significant_results_file,
heatmap_file,
figures_folder,
first_n = heatmap_first_n)
}, error = function(err) {
logging::logerror("Unable to do make a hetamp of results!!!")
logging::logerror(err)
# dev.off()
})
}
if (plot_scatter) {
logging::loginfo(
paste(
"Writing association plots",
"(one for each significant association)",
"to output folder: %s"
),
output
)
association_plots(
metadata,
final_features/offset,
significant_results_file,
output,
figures_folder
)
}
return(significant_results)
}
option_not_valid_error <- function(message, valid_options) {
logging::logerror(paste(message, ": %s"), toString(valid_options))
stop("Option not valid", call. = FALSE)
}
## Quiets concerns of R CMD check
utils::globalVariables(c("name"))
Datasets
Models
