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.

install.packages(c("BiocManager", "remotes"))
BiocManager::install("Statwright/SACAT")

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.

data("gems_microbiome", package = "SACAT")
dim(gems_microbiome)
#> [1]  41 328
table(gems_microbiome$group)
#> 
#> Infant  Older 
#>    191    137

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.

fit <- sacat_se(
    gems_microbiome,
    formula = ~ country + group,
    group = "group",
    reference = "Infant",
    library_size = "reads",
    store_plot_data = TRUE
)
#> Warning: Structural-absence inference returned with numerical warnings for 2
#> taxon/taxa; inspect diagnostics$warning_structural_absence.
fit
#> SACAT fit
#>   Samples: 328; taxa: 41; retained taxa: 41
#>   Contrast: Older - Infant
#>   Structural results returned: 37/41
#>   Regular structural score tests: 37/41
#>   Nonregular structural results: 0/41
#>   Structural-absence discoveries (BH-adjusted p <= 0.05): 11
#>   Relative-abundance tests formed: 41/41
#>   Relative-abundance discoveries (BH-adjusted p <= 0.05): 17
#>   Omnibus discoveries (BH-adjusted p <= 0.05): 21
#>   Results: x$results; formation details: x$diagnostics

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.

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
)

Genus	SA	RA	Omnibus
Haemophilus	0.022	0.158	0.0164
Lactobacillus	0.073	0.000291	0.000388
Faecalibacterium	3.1e-07	9.54e-05	3.2e-07

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

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

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

session <- sessionInfo()
for (field in intersect(c("BLAS", "LAPACK"), names(session))) {
    session[[field]] <- basename(session[[field]])
}
print(session)
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   libblas.so.3 
#> LAPACK: 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] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [3] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [5] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [7] BiocGenerics_0.59.13        generics_0.1.4             
#>  [9] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [11] SACAT_0.99.2                BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] Matrix_1.7-6        jsonlite_2.0.0      compiler_4.6.2     
#>  [4] BiocManager_1.30.27 Rcpp_1.1.2          jquerylib_0.1.4    
#>  [7] statmod_1.5.2       yaml_2.3.12         fastmap_1.2.0      
#> [10] lattice_0.23-1      XVector_0.53.0      R6_2.6.1           
#> [13] S4Arrays_1.13.2     knitr_1.52          DelayedArray_0.39.8
#> [16] maketools_1.3.2     bslib_0.12.0        rlang_1.3.0        
#> [19] cachem_1.1.0        xfun_0.62           sass_0.4.10        
#> [22] sys_3.4.3           otel_0.3.0          SparseArray_1.13.5 
#> [25] cli_3.6.6           digest_0.6.39       grid_4.6.2         
#> [28] lifecycle_1.0.5     evaluate_1.0.5      buildtools_1.0.0   
#> [31] abind_1.4-8         rmarkdown_2.32      tools_4.6.2        
#> [34] htmltools_0.5.9