---
title: "dCSTs for Vaginal Microbiome Data"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{dCSTs for Vaginal Microbiome Data}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

## Introduction

Community state types (CSTs) are the standard way to describe the structure of
vaginal microbial communities. They were originally defined by hierarchical
clustering of 16S rRNA data (Ravel et al. 2011) and later standardized by the
VALENCIA nearest-centroid classifier (France et al. 2020).

**Dominant Community State Types (dCSTs)** group samples by their largest
features and retain groups meeting a sample-support threshold. Fitting uses
abundances rather than supplied CST labels; the retained states still depend
on the fitting sample set. L-infinity normalization divides positive rows by
their maximum at default tolerance, preserving their within-row rankings.

This vignette compares dCSTs with supplied Valencia labels on source training
samples. It is a descriptive comparison, not held-out validation. See the
[function guide](function-guide.html) for fitting versus frozen transfer and the
[dataset guide](example-datasets.html) for provenance and alignment.

**Reference:** Gajer, P. & Ravel, J. (2025). A New Approach to Compositional
Data Analysis using $L^\infty$-normalization with Applications to Vaginal
Microbiome. arXiv:2503.21543.

For the geometry behind this normalization, see the companion article
`r article_link("simplex-corner-geometry", "Simplex corners and max normalization")`,
which derives the departure-angle formulas and tests them in full abundance space
using this same bundled dataset.

For unordered top-2/top-3 memberships, coverage, and the distinction between
taxon graphs and pair-state graphs, see
`r article_link("leading-taxa-state-graphs", "From leading taxa to a 2-CST graph")`.

## The Valencia 2k dataset

The `valencia2k` dataset is a stratified subsample of 2,000 vaginal samples from
the Valencia CST-classifier training set (France et al. 2020). It contains
species-level relative abundances for 178 taxa along with the original Valencia
CST and sub-CST assignments.

```{r load-data}
data(valencia2k)

## Compositional matrix: 2000 samples x 178 taxa
dim(valencia2k$rel)

## Valencia CST assignments
table(valencia2k$cst$Val_CST)

## Top 10 most abundant taxa (by mean relative abundance)
means <- sort(colMeans(valencia2k$rel), decreasing = TRUE)
head(round(means, 4), 10)
```

The dominant taxa are familiar vaginal species: *Lactobacillus iners*,
*L. crispatus*, *Gardnerella vaginalis*, and others. This is exactly the kind
of data where dominant-feature assignments are most intuitive: most samples
are dominated by a single species, and the identity of that species is
biologically meaningful.

## Reconstructing counts

The bundled matrix is compositional (rows sum to 1). For workflows that require
count-like input (e.g., `filter.asv()`), we can reconstruct approximate counts
using the per-sample read depths. These reconstructed and rounded values are
not original observed counts:

```{r reconstruct-counts}
count_mat <- sweep(valencia2k$rel, 1, valencia2k$reads, "*")
count_mat <- round(count_mat)
storage.mode(count_mat) <- "integer"

## Check: library sizes
summary(rowSums(count_mat))
```

The bundled Valencia matrix was already filtered upstream. The subset-building
script takes that prepared matrix as input; we skip `filter.asv()` here and
work directly with the compositional matrix.

## $L^\infty$ normalization

$L^\infty$-normalization divides each row by its maximum value. After
normalization at default tolerance, every positive row has at least one
coordinate equal to 1. Exact-zero rows remain zero, and rows at or below a
positive `tol` remain unchanged. All rows in this example are positive.

```{r linf-normalize}
Z <- normalize.linf(valencia2k$rel)

## Every row max should be 1
summary(apply(Z, 1, max))
```

## Dominance sample sets: dominant species assignment

The dominant species is the column that achieves the maximum in each row.
Samples with the same dominant species form a depth-1 dominance sample set.
With the default first-column tie rule, dominant-feature assignment is
deterministic and sample-independent: it depends only on
the sample's own abundance profile, never on other samples.

```{r linf-dominant.features}
dominant.features <- linf.dominant.features(Z)

## How many distinct dominant taxa?
cat("Distinct dominant taxa:", length(dominant.features$observed.levels), "\n")

## Dominance sample-set size distribution (top 15)
dominance_tab <- sort(table(dominant.features$label), decreasing = TRUE)
head(dominance_tab, 15)
```

We see the expected vaginal community structure: a few well-supported
dominance sample sets (the *Lactobacillus* species, *Gardnerella*, BVAB1) and
a long tail of rarer dominant taxa.

## Depth-1 dCSTs

In practice, many provisional dominance sample sets contain very few samples.
**Dominant Community State Types (dCSTs)** retain states at or above a minimum
support threshold ($n_0$) and handle samples from low-support sets according
to the selected policy.

The package supports two policies for handling low-support provisional sets:

- **"pure"**: named dCSTs stay pure, and samples from low-support sets are
  labeled `RARE_DOMINANT`
- **"absorb"**: samples from low-support sets are reassigned to a retained
  state by restricted argmax

```{r depth1-dcsts}
dcst1 <- linf.csts(Z, n0 = 30, low.freq.policy = "pure")

## Retained depth-1 states
dcst1$retained.feature.labels

## dCST assignment (pure policy)
dcst1_tab <- sort(table(dcst1$lineage.label), decreasing = TRUE)
dcst1_tab
```

### Concordance with Valencia CSTs

The following cross-tabulation describes all 2,000 bundled rows, which were
subsampled from the Valencia classifier training source. The dCST hierarchy
was fitted on those same rows, without using their supplied CST labels.
This overlap precludes interpreting agreement as held-out validation:

```{r concordance-depth1}
## Build concordance table: rows = dCSTs, columns = Valencia CSTs
dcst_labels <- dcst1$lineage.label.absorb  # use absorb view for clean comparison
metadata.order <- match(rownames(Z), valencia2k$cst$sample_id)
stopifnot(!anyNA(metadata.order), !anyDuplicated(valencia2k$cst$sample_id))
val_cst     <- valencia2k$cst$Val_CST[metadata.order]

concordance <- table(dCST = dcst_labels, Valencia = val_cst)

## Show as percentage: what fraction of each dCST falls in each Valencia CST?
pct <- round(100 * prop.table(concordance, margin = 1), 1)
pct
```

Several dCSTs map predominantly to a single supplied Valencia CST in this
source subset:

- *Lactobacillus crispatus* dCST $\to$ CST I (dominated by *L. crispatus*)
- *Lactobacillus iners* dCST $\to$ CST III (dominated by *L. iners*)
- *Gardnerella vaginalis* dCST $\to$ CST IV-B
- BVAB1 dCST $\to$ CST IV-A
- *Lactobacillus gasseri* dCST $\to$ CST II
- *Lactobacillus jensenii* dCST $\to$ CST V

```{r concordance-reverse}
## Reverse view: what fraction of each Valencia CST falls in each dCST?
pct_rev <- round(100 * prop.table(concordance, margin = 2), 1)
pct_rev
```

The two tables describe overlap between different grouping rules on shared
source samples. They do not estimate performance on an independent cohort.

## Depth-2 refinement

Large dCSTs can be refined by dropping the dominant species and repeating the
procedure on the remaining taxa. This produces a two-level hierarchy that
reveals sub-community structure.

```{r depth2-refine}
dcst2 <- refine.linf.csts(Z, dcst1, n0 = 15, refinement.factor = 2,
                           low.freq.policy = "pure", verbose = FALSE)

## Depth-2 dCST table
dcst2_tab <- sort(table(dcst2$lineage.label), decreasing = TRUE)
dcst2_tab
```

### Concordance with Valencia sub-CSTs

```{r concordance-depth2}
dcst2_labels <- dcst2$lineage.label.absorb
val_subcst   <- valencia2k$cst$Val_subCST[metadata.order]

concordance2 <- table(dCST = dcst2_labels, Valencia = val_subcst)

## Percentage by dCST (row-wise)
pct2 <- round(100 * prop.table(concordance2, margin = 1), 1)

## Show only dCSTs with at least 20 samples for readability
dcst2_sizes <- rowSums(concordance2)
pct2[dcst2_sizes >= 20, , drop = FALSE]
```

The depth-2 refinement recovers sub-CST structure. For example, the
*Lactobacillus crispatus* dCST splits based on the second-most-abundant taxon,
and these sub-dCSTs align with Valencia sub-CSTs I-A and I-B.

## The full pipeline

For convenience, `linf.dcst.landmark.pipeline()` runs the entire workflow in
one call: normalization, depth-1 dCSTs, depth-2 refinement, and landmark
computation.

```{r full-pipeline}
out <- linf.dcst.landmark.pipeline(
  valencia2k$rel,
  feature.ids    = colnames(valencia2k$rel),
  feature.labels = colnames(valencia2k$rel),
  n0.depth1      = 30,
  n0.depth2      = 15,
  refinement.factor = 2,
  low.freq.policy   = "pure",
  landmark.view     = "absorb",
  verbose           = FALSE
)

names(out)
```

## Landmark points

All four landmark methods select observed rows by the lineage's leaf-feature
value: its maximum, minimum, or the observation nearest its mean or median.
They do not construct averaged community profiles. The pipeline supplies
L-infinity-normalized values, so these measurements are relative to the sample
maximum, not proportions of total community abundance.

```{r landmarks}
## Depth-1 landmarks
lm1 <- out$landmarks.depth1$landmarks

## Show landmarks for a few key dCST dominance-lineages
key_lineages <- c("Lactobacillus_crispatus", "Lactobacillus_iners",
                  "Gardnerella_vaginalis")
key_lm <- lm1[lm1$lineage.id %in% key_lineages, ]

## For each dCST, show the endpoint.max landmark's target value
## (target-feature abundance relative to the row maximum)
ep_max <- key_lm[key_lm$landmark.type == "endpoint.max",
                 c("lineage.label", "point.name", "target.value")]
ep_max

## And the endpoint.min (lowest target value among assigned rows)
ep_min <- key_lm[key_lm$landmark.type == "endpoint.min",
                 c("lineage.label", "point.name", "target.value")]
ep_min
```

A depth-1 pure dominant feature has normalized value one in every positive
member, so endpoint selection can consist entirely of ties. Absorbed samples
may have a smaller target value because their assigned feature is not their
original maximum. Neither endpoint establishes near-total or marginal
unit-sum dominance. Inspect unit-sum relative abundances for that question.

## Summary

This workflow fits supported dominance groups, refines selected parents and
selects observed landmark rows. Its concordance tables describe overlap with
Valencia labels on source training samples. To assign separate samples while
preserving the fitted states, use `transfer.dcsts()` and inspect unassigned
(`NA`) depths; fitting on new samples would learn a different hierarchy.
The [function guide](function-guide.html) explains these contracts and the
[dataset guide](example-datasets.html#reference-and-query) gives a disjoint
reference/query demonstration.

## References

- France, M. T., Ma, B., Gajer, P., et al. (2020). VALENCIA: a nearest
  centroid classification method for vaginal microbial communities based on
  composition. *Microbiome*, 8(1), 166.
- Gajer, P. & Ravel, J. (2025). A New Approach to Compositional Data Analysis
  using $L^\infty$-normalization with Applications to Vaginal Microbiome.
  arXiv:2503.21543.
- Ravel, J., Gajer, P., Abdo, Z., et al. (2011). Vaginal microbiome of
  reproductive-age women. *PNAS*, 108(Suppl 1), 4680--4687.