--- title: "Dominant Community State Types: From Normalization to a Gut Demonstration" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Dominant Community State Types: From Normalization to a Gut Demonstration} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{css, echo=FALSE} table td, table th { white-space: normal; overflow-wrap: anywhere; } ``` ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) library(linf) article_link <- function(name, title) { if (identical(Sys.getenv("IN_PKGDOWN"), "true")) { url <- paste0(name, ".html") } else if (isTRUE(getOption("linf.local.preview"))) { url <- paste0("../articles/", name, ".html") } else { url <- paste0("https://github.com/pgajer/linf/blob/main/vignettes/articles/", name, ".Rmd") title <- paste0(title, " (article source)") } sprintf("[%s](%s)", title, url) } ``` The `linf` package implements L-infinity normalization and Dominant Community State Types (dCSTs) for compositional data. This vignette demonstrates the complete workflow in two parts: first a quick-start with toy data to introduce the core functions, then a real-data analysis of gut microbiome samples from the American Gut Project showing how to inspect dominant and subdominant community structure. For entry points and policy contracts see the [function guide](function-guide.html); for dataset selection and safe joins see [example datasets](example-datasets.html). ## Part 1: Quick Start with Toy Data ### L-infinity normalization L-infinity normalization scales each row of a count matrix by its maximum value when that maximum exceeds `tol`. At the default `tol = 0`, positive rows have maximum one and exact-zero rows remain zero. With positive tolerance, rows at or below it remain unchanged, including small nonzero entries. ```{r toy-normalize} set.seed(1) S.counts <- matrix( rpois(30, lambda = 5), nrow = 10, ncol = 3, dimnames = list(paste0("s", 1:10), c("Taxon_A", "Taxon_B", "Taxon_C")) ) Z <- normalize.linf(S.counts) apply(Z, 1, max) # all 1 ``` ### Dominance sample sets Each sample is assigned to its dominant feature: the column achieving the within-sample maximum. Samples with the same dominant feature form a depth-1 *dominance sample set*. ```{r toy-dominant.features} dominant.features <- linf.dominant.features(Z) table(dominant.features$label) ``` ### Truncated dCSTs `linf.csts()` groups samples by dominant feature and collapses small groups (below threshold `n0`) into a `RARE_DOMINANT` bucket. This produces *Dominant Community State Types* (dCSTs). ```{r toy-csts} A <- matrix(c(5, 1, 0), nrow = 6, ncol = 3, byrow = TRUE) B <- matrix(c(1, 5, 0), nrow = 2, ncol = 3, byrow = TRUE) C <- matrix(c(1, 0, 5), nrow = 2, ncol = 3, byrow = TRUE) S <- rbind(A, B, C) S <- sweep(S, 1, rowSums(S), "/") colnames(S) <- c("Dom1", "Dom2", "Dom3") res <- linf.csts(S, n0 = 5) table(res$lineage.label, useNA = "ifany") ``` Only Dom1 (6 samples) meets the `n0 = 5` threshold. Dom2 and Dom3 are collapsed into `RARE_DOMINANT`. ## Part 2: Gut Microbiome dCST Demonstration ### Background Single-species dominance is common in some microbial ecosystems and less common in others. With default first-column tie handling, dCSTs provide a deterministic way to summarize the dominant feature and then refine large groups by subdominant features. The gut example below demonstrates those mechanics without treating the selected package subset as an epidemiologic sample. ### The data The `agp_gut` dataset bundled with this package contains `r nrow(agp_gut$counts)` gut microbiome samples from the American Gut Project (PRJEB11419), a citizen-science 16S rRNA V4 survey. The dataset is a dCST-stratified demonstration subset that includes all samples in four selected uncommon dCSTs and fills the remaining slots with a seed-42 simple random sample from the eligible background. Phenotypes are joined only after selection and do not influence membership. Because inclusion probabilities differ by dCST, it is not a probability sample of the underlying cohort; phenotype frequencies, effect sizes, and hypothesis tests must not be interpreted as population results. ```{r load-data} data(agp_gut) str(agp_gut, max.level = 1) dim(agp_gut$counts) # samples x taxa ``` Metadata includes self-reported health conditions parsed from the AGP questionnaire: ```{r meta-summary} disease_cols <- c("IBS", "IBD", "Obesity", "Cardiovascular_disease", "Autoimmune", "Acid_reflux") sapply(disease_cols, function(col) sum(agp_gut$meta[[col]], na.rm = TRUE)) ``` ### Step 1: Quality filtering `filter.asv()` removes low-depth samples and rare taxa in one call. We require at least 1,000 reads per sample and taxa present in at least 5% of samples. ```{r filter} filt <- filter.asv(agp_gut$counts, min.lib = 1000, prev.prop = 0.05, min.count = 2) dim(filt$counts) dim(filt$rel) ``` ### Step 2: L-infinity normalization ```{r normalize} M <- normalize.linf(filt$counts) # Verify: every row max is 1 stopifnot(all(abs(apply(M, 1, max) - 1) < 1e-10)) ``` ### Step 3: Depth-1 dCSTs ```{r dcst-depth1} csts <- linf.csts(M, n0 = 30) # dCST size distribution sort(table(csts$lineage.label), decreasing = TRUE) ``` The largest retained groups are *Bacteroides* and *Escherichia-Shigella*. The figure below labels all seven retained groups and the separate rare bucket. This tutorial uses the pure policy; the README figure uses absorption, so its group counts differ even though both start with the same filtered samples. The `RARE_DOMINANT` bucket collects samples dominated by taxa too infrequent to form their own dCST at this threshold. The *Escherichia-Shigella* values are processed 16S read counts. [McDonald et al. (2018)](https://doi.org/10.1128/mSystems.00031-18) describe microbial overgrowth during room-temperature shipping in American Gut, including *E. coli* blooms. That source does not establish the cause of the abundance in this bundled subset. The displayed dominance count alone cannot establish organismal abundance or a particular technical cause, including primer cross-reactivity. ### How sensitive is depth 1 to the support threshold? The following comparison uses the same `r nrow(M)` filtered samples and `r ncol(M)` features at every setting. Only the minimum support and the low-support policy change; no phenotype information enters these fits. ```{r support-sensitivity} sensitivity <- do.call(rbind, lapply(c(10L, 30L, 60L, 200L), function(support) { do.call(rbind, lapply(c("pure", "absorb"), function(policy) { fit <- linf.csts(M, n0 = support, low.freq.policy = policy, return.diagnostics = TRUE) counts <- summary(fit) stopifnot(counts$input.samples == nrow(M), counts$input.samples == counts$assigned.samples + counts$unassigned.samples) data.frame(n0 = support, Policy = policy, Groups = counts$retained.lineages, Rare = counts$rare.samples, Reassigned = if (policy == "absorb") sum(fit$diagnostics$reassigned) else 0L, Unassigned = counts$unassigned.samples) })) })) knitr::kable(sensitivity) ``` At `n0 = 30`, the pure view has `r sensitivity$Rare[sensitivity$n0 == 30 & sensitivity$Policy == "pure"]` samples in a synthetic rare category; absorption reassigns `r sensitivity$Reassigned[sensitivity$n0 == 30 & sensitivity$Policy == "absorb"]` samples to the retained groups. Raising the support threshold reduces the number of retained groups. At 200, no feature has enough provisional support: the pure view puts every sample in the rare category, while absorption leaves every sample unassigned. Rare samples are assigned to a synthetic category; "Unassigned" means an `NA` assignment and is counted separately. This is a descriptive sensitivity check, not a search for an optimal or clinically validated threshold. Stable group counts would not by themselves establish stable membership, biological validity or predictive accuracy. The stratified sampling limitation applies to every row of this comparison. ### Step 4: Depth-2 refinement `refine.linf.csts()` subdivides large dCSTs by their *subdominant* species. For instance, a Bacteroides-dominated sample might be further classified by whether its second-most-abundant taxon is *Faecalibacterium* or *Lachnospiraceae*. This captures co-dominance patterns that are especially important in the diverse gut environment. ```{r dcst-depth2} csts2 <- refine.linf.csts(M, csts, n0 = 30) # Show depth-2 dCSTs with >= 20 samples tab2 <- sort(table(csts2$lineage.label), decreasing = TRUE) tab2[tab2 >= 20] ``` ### Dominance strength A useful diagnostic: the relative abundance of the dominant species in each sample. In the gut, most samples have modest dominance (< 50%), consistent with the diverse nature of the ecosystem. ```{r dominance-dist, fig.cap = "Distribution of dominance strength across gut samples. Most samples show moderate dominance, with a tail of strongly dominated communities."} rel <- filt$rel dom_strength <- apply(rel, 1, max) hist(dom_strength, breaks = 50, col = "#2ecc71", border = "white", main = "Dominance Strength in Gut Microbiome", xlab = "Relative abundance of dominant species", ylab = "Number of samples") abline(v = 0.5, col = "red", lty = 2, lwd = 2) legend("topright", "50% dominance", col = "red", lty = 2, lwd = 2, bty = "n") ``` ### dCST size distribution ```{r dcst-barplot, fig.cap = "Pure-policy depth-1 assignments for 763 filtered AGP samples: seven retained groups and 159 samples in the synthetic rare bucket. Counts describe the stratified demonstration subset, not population frequencies."} tab1 <- sort(table(csts$lineage.label), decreasing = TRUE) # Use the last available taxonomy rank for display; keep full IDs in the fit. short_taxon <- function(id) { if (id == "RARE_DOMINANT") return("Rare dominant") ranks <- sub("^[a-z]__", "", strsplit(id, ";", fixed = TRUE)[[1L]]) ranks <- ranks[nzchar(ranks) & ranks != "__"] tail(ranks, 1L) } names(tab1) <- make.unique(vapply(names(tab1), short_taxon, character(1))) par(mar = c(5, 10, 3, 2)) bp <- barplot(rev(tab1), horiz = TRUE, col = ifelse(rev(names(tab1)) == "Rare dominant", "#e74c3c", "#3498db"), las = 1, cex.names = 0.8, xlim = c(0, max(tab1) * 1.18), xlab = "Number of samples", main = "Gut dCST sizes (pure policy, depth 1)") text(rev(tab1), bp, labels = rev(tab1), pos = 4, cex = 0.8) ``` ### Phenotype metadata The bundle retains selected self-reported phenotype fields so that users can inspect the object structure and practice data alignment. The following table is purely descriptive; dCST-stratified selection precludes population prevalence estimation or treating association calculations as cohort results. ```{r phenotype-summary} meta <- agp_gut$meta[match(rownames(M), agp_gut$meta$Run), ] stopifnot(all(meta$Run == rownames(M))) conditions <- c("IBS", "IBD", "Obesity", "Cardiovascular_disease", "Autoimmune", "Acid_reflux", "Lung_disease") phenotype_summary <- data.frame( Condition = conditions, Recorded_cases = vapply( conditions, function(x) sum(meta[[x]] == 1, na.rm = TRUE), integer(1) ), Non_missing = vapply( conditions, function(x) sum(!is.na(meta[[x]])), integer(1) ) ) knitr::kable( phenotype_summary, caption = "Recorded phenotypes in the selected demonstration subset" ) ``` ## Limitations This analysis has several important caveats: - **Self-reported metadata.** AGP health conditions are self-reported by citizen-science participants, with unknown accuracy. - **16S resolution.** The Escherichia-Shigella inflated abundance is a known V4 artifact. Shotgun metagenomics would resolve this. - **dCST-stratified demonstration subset.** Four uncommon dCSTs are included exhaustively and the eligible background is sampled with seed 42. Phenotypes do not affect membership, but differing inclusion probabilities mean the object must not be used for population prevalence estimates, effect-size estimation, or association testing. ## Summary This workflow summarizes dominant features and refines selected groups using remaining features. The support threshold, refinement selection, low-support policy and tie handling determine the stored hierarchy. Its descriptive patterns apply to this selected dataset; they do not establish population associations or clinical utility. ## Choose the next step | Your next question | Guide | |:--|:--| | Why do compositional arms crowd near a simplex corner? | `r article_link("simplex-corner-geometry", "Simplex corners and max normalization")` | | How do pure top-2/top-3 states, coverage and graph membership relate? | `r article_link("leading-taxa-state-graphs", "From leading taxa to a 2-CST graph")` | | How much of the chosen distance geometry survives a 3D fit? | [Embedding fidelity in grip](https://pgajer.github.io/grip/articles/embedding-fidelity.html) | These companion articles complement the quick-start workflow above. The state-graph article explains unordered leading-feature sets, which are separate from the ordered dCST hierarchy constructed here. Their optional visualization dependencies are not required for this introductory vignette.