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

## ----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-cells----------------------------------------------------------------
cells <- linf.cells(Z)
table(cells$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$cell.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$cell.label), decreasing = TRUE)

## ----dcst-depth2--------------------------------------------------------------
csts2 <- refine.linf.csts(M, csts, n0 = 30)

# Show depth-2 dCSTs with >= 20 samples
tab2 <- sort(table(csts2$cell.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 = "dCST size distribution at depth 1. Bacteroides and Escherichia-Shigella dominate; rare dCSTs in the tail are collapsed into RARE_DOMINANT."----
tab1 <- sort(table(csts$cell.label), decreasing = TRUE)
par(mar = c(10, 4, 3, 1))
bp <- barplot(tab1,
        col = ifelse(names(tab1) == "RARE_DOMINANT", "#e74c3c", "#3498db"),
        las = 2, cex.names = 0.7,
        ylab = "Number of samples",
        main = "Gut dCST Size Distribution (depth 1)")
text(bp, tab1 + 5, labels = tab1, cex = 0.7, pos = 3)

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