---
title: "Tweedieverse"
author: "Nalin Arora, Ziyu Liu, Ali Rahnavard,Himel Mallick"
date: "`r BiocStyle::doc_date()`"
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Tweedieverse}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
```

## Overview

`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:

- microbiome relative abundances stored in a `TreeSummarizedExperiment`,
- single-cell RNA-seq UMI counts stored in a `SingleCellExperiment`,
- single-cell RNA-seq non-UMI/read counts stored in a `SingleCellExperiment`,
- pseudobulk RNA-seq counts generated with `muscat::aggregateData()`,
- bulk RNA-seq counts from the Bioconductor `airway` dataset stored in a
  `SummarizedExperiment`,
- matched microbiome and metabolomics assays stored in a
  `MultiAssayExperiment`.

## Bioconductor Motivation

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.

## Installation

After Bioconductor acceptance, the release version can be installed with:

```{r install-bioc, eval=FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

BiocManager::install("Tweedieverse")
```

The development version can be installed from GitHub:

```{r install-github, eval=FALSE}
if (!requireNamespace("remotes", quietly = TRUE)) {
  install.packages("remotes")
}

remotes::install_github("himelmallick/Tweedieverse")
```

```{r packages}

library(Tweedieverse)
library(SummarizedExperiment)
library(SingleCellExperiment)
library(TreeSummarizedExperiment)
library(MultiAssayExperiment)
library(muscat)
library(airway)
```

## Loading Packaged Example Data

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.

```{r load-example-data}
load(system.file(
  "extdata",
  "tweedieverse_data.RData",
  package = "Tweedieverse"
))
```

## Microbiome Relative Abundance Example

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.

```{r microbiome-fit}
microbiome_tse

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
)

head(microbiome_fit)
```

## Single-Cell RNA-seq UMI Count Example

For 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.

```{r umi-fit}
umi_sce

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
)

head(umi_fit)
```

## Single-Cell RNA-seq Non-UMI / Read Count Example

The 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.

```{r nonumi-fit}
nonumi_sce

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
)

head(nonumi_fit)
```

## Pseudobulk Single-Cell RNA-seq Example

Pseudobulk 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`.

```{r pseudobulk-fit}
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)

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
)

head(pseudobulk_fit)
```

## Bulk RNA-seq Count Example

For 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.

```{r bulk-rnaseq-fit}
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

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
)

head(bulk_rnaseq_fit)
```

## Multi-omics Container

For 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.

```{r multiassay}
franzosa_example_mae

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
)

names(multi_fit)
head(multi_fit$taxaAbundanceProfile)
head(multi_fit$harmonizedMetaboliteProfile)
```

## Session Info

```{r session-info}
sessionInfo()
```