dCSTs for Vaginal Microbiome Data

knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5)
library(linf)

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) offer an alternative that requires no clustering and no training data. A dCST is defined by $L^\infty$-normalization: divide each sample's abundance vector by its maximum value, then group samples by which taxon achieves that maximum. This is a purely geometric operation on compositional data --- it identifies the dominant species in each sample and groups samples accordingly.

This vignette demonstrates the dCST workflow on a real vaginal 16S dataset and shows that dCSTs recover the established CST structure with high concordance.

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.

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.

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:

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))

Since the Valencia 13k data was already filtered (min 3,000 reads, taxa present in $\ge$ 1% of samples), 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, every sample has at least one coordinate equal to 1, and all other coordinates lie in $[0, 1]$.

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. The assignment is deterministic and sample-independent: it depends only on the sample's own abundance profile, never on other samples.

cells <- linf.cells(Z)

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

## Dominance sample-set size distribution (top 15)
dominance_tab <- sort(table(cells$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 only states 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:

The legacy policy name "rare" is still accepted as a deprecated alias for "pure".

dcst1 <- linf.csts(Z, n0 = 30, low.freq.policy = "pure")

## Retained depth-1 states
dcst1$kept.cells.lbl

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

Concordance with Valencia CSTs

The key test: do dCSTs correspond to the established Valencia CSTs? Let's build a cross-tabulation:

## Build concordance table: rows = dCSTs, columns = Valencia CSTs
dcst_labels <- dcst1$cell.label.absorb  # use absorb view for clean comparison
val_cst     <- valencia2k$cst$Val_CST

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

The concordance is striking. Each dCST maps predominantly to a single Valencia CST:

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

This two-way concordance confirms that dCSTs are an alternative construct of CSTs: they recover the same groupings without any clustering algorithm or training data.

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.

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$cell.label), decreasing = TRUE)
dcst2_tab

Concordance with Valencia sub-CSTs

dcst2_labels <- dcst2$cell.label.absorb
val_subcst   <- valencia2k$cst$Val_subCST

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.

out <- linf.dcst.landmark.pipeline(
  count_mat,
  feature.ids    = colnames(count_mat),
  feature.labels = colnames(count_mat),
  n0.depth1      = 30,
  n0.depth2      = 15,
  refinement.factor = 2,
  low.freq.policy   = "pure",
  landmark.view     = "absorb",
  verbose           = FALSE
)

names(out)

Landmark points

Each dCST dominance-lineage has characteristic landmark points: the most extreme sample (highest value of the dominant taxon), the least extreme (lowest value while still dominant), and mean/median representatives. These help characterize what a "typical" or "extreme" member of each community type looks like.

## 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$cell.id %in% key_lineages, ]

## For each dCST, show the endpoint.max landmark's target value
## (how dominant is the dominant species in the most extreme sample?)
ep_max <- key_lm[key_lm$landmark.type == "endpoint.max",
                 c("cell.label", "point.name", "target.value")]
ep_max

## And the endpoint.min (least dominant while still assigned to this lineage)
ep_min <- key_lm[key_lm$landmark.type == "endpoint.min",
                 c("cell.label", "point.name", "target.value")]
ep_min

The endpoint.max landmarks show samples where the dominant species comprises nearly the entire community, while endpoint.min landmarks show the boundary cases --- samples where dominance is marginal.

Summary

This vignette demonstrated the dCST workflow on real vaginal 16S data:

  1. $L^\infty$-normalization maps compositional data to the $L^\infty$-simplex
  2. Dominant-feature assignment assigns each sample to its dominant taxon
  3. Truncated dCSTs retain well-supported states and handle rare dominance patterns
  4. Depth-2 refinement reveals sub-community structure
  5. Landmark points characterize extreme and typical community members

The concordance with established Valencia CSTs is high, validating dCSTs as a principled, deterministic alternative to clustering-based community typing. Key advantages:

References



Try the linf package in your browser

Any scripts or data that you put into this service are public.

linf documentation built on Aug. 5, 2026, 9:08 a.m.