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

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

## ----toy-dominant.features----------------------------------------------------
dominant.features <- linf.dominant.features(Z)
table(dominant.features$label)

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

## ----load-data----------------------------------------------------------------
data(agp_gut)
str(agp_gut, max.level = 1)

dim(agp_gut$counts)   # samples x taxa

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

## ----filter-------------------------------------------------------------------
filt <- filter.asv(agp_gut$counts, min.lib = 1000, prev.prop = 0.05,
                   min.count = 2)
dim(filt$counts)
dim(filt$rel)

## ----normalize----------------------------------------------------------------
M <- normalize.linf(filt$counts)

# Verify: every row max is 1
stopifnot(all(abs(apply(M, 1, max) - 1) < 1e-10))

## ----dcst-depth1--------------------------------------------------------------
csts <- linf.csts(M, n0 = 30)

# dCST size distribution
sort(table(csts$lineage.label), decreasing = TRUE)

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

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