---
title: "Structural-absence and conditional-abundance testing with SACAT"
author:
  - name: Yiqian Zhang
package: SACAT
output:
  BiocStyle::html_document:
    toc: true
    toc_depth: 2
vignette: >
  %\VignetteIndexEntry{Structural-absence and conditional-abundance testing with SACAT}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

# Overview

SACAT analyzes microbiome counts to separate evidence for differences in
structural absence (SA) from differences in relative abundance conditional
on presence (RA). Sequencing depth enters the model because a zero can
reflect either absence or nondetection. RA effects are standardized over
covariates and centered against a reference from other eligible taxa.
The primary omnibus p-value is `min(1, 2 * min(p_SA, p_RA))`.
Benjamini–Hochberg (BH) adjustment is applied separately to each result family.

# Installation

Install SACAT from GitHub using R 4.6.0 or later and a C++ compiler.

```{r installation, eval=FALSE}
install.packages(c("BiocManager", "remotes"))
BiocManager::install("Statwright/SACAT")
```

```{r packages}
library(SACAT)
library(SummarizedExperiment)
```

# The GEMS age comparison

The Global Enteric Multicenter Study collected stool samples from children
with diarrhea and controls without diarrhea. Here we use controls from
Bangladesh, The Gambia, Kenya, and Mali to compare 191 infants aged 1–11
months with 137 children aged 24–58 months.

```{r example-data}
data("gems_microbiome", package = "SACAT")
dim(gems_microbiome)
table(gems_microbiome$group)
```

The `counts` assay contains raw 16S counts for the 41 named bacterial genera
detected in at least 20% of the 328 children. The `reads` column in
`colData(gems_microbiome)` gives each sample's library size from the full
count table, before genus selection. Use this original depth when fitting
a subset of taxa.

# Fit and inspect the results

The model compares older children with infants while adjusting for country.
`sacat_se()` reads the counts and sample annotations from the
`SummarizedExperiment`. Use `sacat()` for a count matrix and metadata frame.

```{r fit-model}
fit <- sacat_se(
    gems_microbiome,
    formula = ~ country + group,
    group = "group",
    reference = "Infant",
    library_size = "reads",
    store_plot_data = TRUE
)
fit
```

We can inspect the adjusted p-values for three genera after fitting all 41.
Each column uses BH adjustment across the full set of genera.

```{r primary-results}
genera <- c("Haemophilus", "Lactobacillus", "Faecalibacterium")
columns <- c("p_adj_structural_absence", "p_adj_relative_abundance",
             "p_adj_omnibus")
rows <- match(genera, fit$results$taxon)
p_values <- fit$results[rows, columns]
p_values[] <- lapply(p_values, function(p) sprintf("%.3g", p))
knitr::kable(
    data.frame(Genus = genera, p_values),
    col.names = c("Genus", "SA", "RA", "Omnibus"),
    row.names = FALSE
)
```

In older children, *Haemophilus* showed evidence of lower structural absence
(SA only). *Lactobacillus* showed lower relative abundance conditional on
presence (RA only). For *Faecalibacterium*, older children had both lower
structural absence and higher relative abundance conditional on presence.
The omnibus test detected all three associations. For other genera,
`fit$results` provides the effect estimates and p-values, and
`fit$diagnostics` reports test availability and numerical warnings.

# Plot the association profile

```{r association-profile, dev="svg", fig.width=11.3, fig.height=4.8, out.width="100%"}
plot(fit, features = genera, group_labels = c("1-11 mo", "24-58 mo"))
```

**A. Component evidence.** Filled points mark BH-adjusted p-values at or
below 0.05. Positive signed statistics indicate higher structural absence
or higher reference-adjusted conditional abundance in older children.
Stars beside the genus names summarize the adjusted omnibus p-values.

**B. Detection and zero sources.** The blue portion gives the fraction of
children in whom the genus was detected. The other two portions divide the
observed zeros between estimated structural absence (SA) and nondetection
despite presence (ND), using fitted probabilities.

**C. Abundance.** Each point is a positive observed count divided by that
sample's full library size, expressed as a percentage. Diamonds show fitted
abundance conditional on presence, obtained by averaging fitted log abundance
over the same country profiles for each age group and exponentiating. The
log ratio of the diamonds gives the raw abundance contrast. The RA test
centers this contrast against the reference estimated from other genera.

# Data source and session information

Pop M, Walker AW, Paulson J, et al. (2014). Diarrhea in young children from
low-income countries leads to large-scale alterations in intestinal microbiota
composition. *Genome Biology* **15**, R76. [Study][gems-paper]

The counts and sample annotations come from the Bioconductor experiment data
package [msd16s][gems-data]. See `?gems_microbiome` for the selection criteria.

[gems-paper]: https://doi.org/10.1186/gb-2014-15-6-r76
[gems-data]: https://bioconductor.org/packages/msd16s/

```{r session-info}
session <- sessionInfo()
for (field in intersect(c("BLAS", "LAPACK"), names(session))) {
    session[[field]] <- basename(session[[field]])
}
print(session)
```