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

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

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

## ----linf-normalize-----------------------------------------------------------
Z <- normalize.linf(valencia2k$rel)

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

## ----linf-cells---------------------------------------------------------------
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)

## ----depth1-dcsts-------------------------------------------------------------
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-depth1-------------------------------------------------------
## 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

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

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

## ----concordance-depth2-------------------------------------------------------
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]

## ----full-pipeline------------------------------------------------------------
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)

## ----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$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