## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                     fig.width = 6, fig.height = 4)
library(linf)
# Keep table rendering separate from the runnable analysis code.
knitr::opts_chunk$set(render = function(x, options) {
  rendered <- knitr::knit_print(x, options = options)
  if (inherits(x, "knitr_kable")) {
    knitr::asis_output(paste0(
      '<div class="linf-table" role="region" aria-label="Computed table; scroll horizontally if needed" tabindex="0">\n\n',
      rendered, '\n\n</div>\n'
    ))
  } else rendered
})

## ----inventory----------------------------------------------------------------
data(agp_gut, package = "linf")
data(valencia2k, package = "linf")
data(valencia13k_dcst_depth2_merged, package = "linf")
data(valencia13k_dcst_depth3_merged, package = "linf")
data(valencia_linf_hypercube_1k, package = "linf")

observed <- rbind(
  gut_counts = dim(agp_gut$counts),
  vaginal_rel = dim(valencia2k$rel),
  merged_depth2 = dim(valencia13k_dcst_depth2_merged$assignments),
  merged_depth3 = dim(valencia13k_dcst_depth3_merged$assignments),
  four_component_rel = dim(valencia_linf_hypercube_1k$rel4)
)
colnames(observed) <- c("Rows", "Columns")
stopifnot(identical(unname(observed),
                   matrix(c(766L, 2000L, 12881L, 12881L, 1000L,
                            314L, 178L, 6L, 7L, 4L), ncol = 2)))
knitr::kable(observed, caption = "Dimensions read from the bundled objects")

## ----provenance-location------------------------------------------------------
provenance <- system.file("DATA_PROVENANCE.md", package = "linf", mustWork = TRUE)
cat(readLines(provenance, warn = FALSE)[1:12], sep = "\n")

## ----alignment----------------------------------------------------------------
X <- agp_gut$counts
stopifnot(!anyDuplicated(agp_gut$meta$Run),
          !anyDuplicated(agp_gut$taxa))
sample.order <- match(rownames(X), agp_gut$meta$Run)
feature.order <- match(colnames(X), agp_gut$taxa)
stopifnot(!anyNA(sample.order), !anyNA(feature.order))
meta <- agp_gut$meta[sample.order, , drop = FALSE]
features <- data.frame(
  id = paste0("asv_", seq_len(ncol(X))),
  taxonomy = agp_gut$taxa[feature.order]
)
stopifnot(identical(rownames(X), meta$Run))

# After filtering, indices refer to the supplied X, in its original order.
filtered <- filter.asv(X, min.lib = 1000, prev.prop = 0.05,
                       min.count = 2, verbose = FALSE)
meta.filtered <- meta[filtered$kept.sample.idx, , drop = FALSE]
features.filtered <- features[filtered$kept.feature.idx, , drop = FALSE]
stopifnot(identical(rownames(filtered$counts), meta.filtered$Run),
          nrow(features.filtered) == ncol(filtered$counts))

## ----shared-source-join-------------------------------------------------------
merged <- valencia13k_dcst_depth3_merged
chart.meta <- valencia_linf_hypercube_1k$meta
stopifnot(!anyDuplicated(merged$assignments$source_row))
join <- match(chart.meta$source_row, merged$assignments$source_row)
stopifnot(!anyNA(join))
chart.depth3 <- merged$assignments$dcst_depth3[join]
stopifnot(identical(as.character(chart.meta$Val_CST),
                   as.character(merged$assignments$Val_CST[join])))
length(chart.depth3)
subset.join <- match(valencia2k$cst$source_row, merged$assignments$source_row)
stopifnot(!anyNA(subset.join),
          identical(as.character(valencia2k$cst$Val_CST),
                    as.character(merged$assignments$Val_CST[subset.join])))

## ----toy-edge-cases-----------------------------------------------------------
toy <- rbind(
  a1 = c(10, 3, 1), a2 = c(9, 2, 1),
  b1 = c(2, 10, 1), b2 = c(3, 9, 1),
  low_support = c(1, 2, 10), tie = c(4, 4, 0),
  zero = c(0, 0, 0), small = c(0.01, 0.005, 0)
)
colnames(toy) <- c("asv_1", "asv_2", "asv_3")
display <- c("Taxon A", "Taxon B", "Taxon C")
M <- normalize.linf(toy, tol = 0.02)
d <- linf.dominant.features(M, feature.labels = display)
fit <- linf.csts(M, feature.labels = display, n0 = 2)
edge.table <- data.frame(
  sample = rownames(M), maximum = apply(M, 1, max),
  dominant_ID = d$id, pure = fit$lineage.label.pure,
  absorb = fit$lineage.label.absorb
)
knitr::kable(edge.table, row.names = FALSE)
stopifnot(is.na(d$id[7]), d$id[6] == "asv_1", d$id[8] == "asv_1",
          identical(unname(M["small", ]), unname(toy["small", ])))

## ----random-tie---------------------------------------------------------------
set.seed(17)
random1 <- linf.dominant.features(M, tie.method = "random")$id
set.seed(17)
random2 <- linf.dominant.features(M, tie.method = "random")$id
stopifnot(identical(random1, random2))

## ----real-transfer------------------------------------------------------------
set.seed(20260914)
rows <- sample(seq_len(nrow(valencia2k$rel)), 160)
reference.rows <- rows[1:120]
query.rows <- rows[121:160]
stopifnot(length(intersect(reference.rows, query.rows)) == 0L)
reference <- valencia2k$rel[reference.rows, , drop = FALSE]
query <- valencia2k$rel[query.rows, , drop = FALSE]
reference.M <- normalize.linf(reference)
reference.fit <- linf.csts(reference.M, n0 = 5,
                           low.freq.policy = "absorb")
reference.fit <- refine.linf.csts(
  reference.M, reference.fit, n0 = 3, refinement.factor = 2,
  low.freq.policy = "absorb", verbose = FALSE
)

# Deliberately reverse query columns. Matching uses the unchanged feature IDs.
query <- query[, rev(seq_len(ncol(query))), drop = FALSE]
stopifnot(setequal(colnames(query), reference.fit$feature.ids))
transferred <- transfer.dcsts(query, reference.fit, match.by = "feature.ids")
knitr::kable(data.frame(
  depth = colnames(transferred$assignment),
  assigned = colSums(!is.na(transferred$assignment)),
  unassigned = colSums(is.na(transferred$assignment))
), row.names = FALSE)

## ----no-matching-feature------------------------------------------------------
unmatched <- matrix(5, nrow = 1, ncol = 1,
                    dimnames = list("unmatched_query", "not_in_reference"))
no.assignment <- transfer.dcsts(unmatched, reference.fit)$assignment
stopifnot(all(is.na(no.assignment)))
no.assignment

## ----static-chart, fig.cap="Two coordinates of the Li-reference chart. Colors describe dominance within the selected four taxa; overlapping points and the omitted third coordinate limit what this projection shows."----
four <- valencia_linf_hypercube_1k
embedding <- linf.hypercube.embedding(four$rel4, reference = "Li")
stopifnot(identical(dim(embedding), c(1000L, 3L)),
          all(is.finite(embedding)))
colors <- c(Li = "#0072B2", Lc = "#D55E00", Gv = "#009E73", Bv = "#CC79A7")
plot(embedding[, 1], embedding[, 2],
     col = colors[as.character(four$meta$dominant_component)],
     pch = 16, cex = 0.5, xlim = c(0, 1), ylim = c(0, 1),
     xlab = "Lc coordinate (reference Li)",
     ylab = "Gv coordinate (reference Li)")
legend("topright", legend = names(colors), col = colors, pch = 16,
       title = "Dominant component", bty = "n", cex = 0.8)