## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
library(TaxSEA)

## ----matrix-example-----------------------------------------------------------
set.seed(42)
counts <- matrix(rpois(30 * 8, lambda = 10), nrow = 30, ncol = 8)
rownames(counts) <- paste0("Taxon_", seq_len(30))
colnames(counts) <- paste0("Sample_", seq_len(8))

# Plant a signal: raise one set's members in the last four samples
counts[1:6, 5:8] <- counts[1:6, 5:8] * 5

my_sets <- list(
  elevated_set = paste0("Taxon_", 1:6),
  control_set  = paste0("Taxon_", 20:27)
)

scores <- ssTaxSEA(counts, custom_db = my_sets, min_set_size = 3)
dim(scores)
round(scores, 2)

## ----matrix-signal------------------------------------------------------------
rowMeans(scores[, 5:8]) - rowMeans(scores[, 1:4])

## ----single-sample------------------------------------------------------------
ssTaxSEA(counts[, 1, drop = FALSE], custom_db = my_sets, min_set_size = 3)

## ----builtin, eval = FALSE----------------------------------------------------
# scores <- ssTaxSEA(counts)
# dim(scores)

## ----tse, message = FALSE, eval = requireNamespace("SummarizedExperiment", quietly = TRUE)----
library(SummarizedExperiment)

se <- SummarizedExperiment(assays = list(counts = counts))
scores_se <- ssTaxSEA(se, custom_db = my_sets, min_set_size = 3)

# Identical to the matrix result
all.equal(scores, scores_se)

## ----baselines----------------------------------------------------------------
rowMeans(scores)

## ----centering----------------------------------------------------------------
scaled <- t(scale(t(scores)))
round(scaled, 2)

## ----testing------------------------------------------------------------------
group <- rep(c("control", "case"), each = 4)
wilcox.test(scores["elevated_set", ] ~ group)$p.value

## ----pseudocount--------------------------------------------------------------
p1 <- ssTaxSEA(counts, custom_db = my_sets, min_set_size = 3, pseudocount = 1)
round(p1["elevated_set", ] - scores["elevated_set", ], 3)

## ----proportion---------------------------------------------------------------
pct <- sweep(counts, 2, colSums(counts), "/") * 100
from_pct <- ssTaxSEA(pct, custom_db = my_sets, min_set_size = 3,
                     scale = "proportion")
from_counts <- ssTaxSEA(counts, custom_db = my_sets, min_set_size = 3,
                        scale = "proportion")
all.equal(from_pct, from_counts)

## ----top_k--------------------------------------------------------------------
# Remove the rarest half of the detected taxa from one sample
depleted <- counts
rare <- order(depleted[, 1])[seq_len(nrow(depleted) / 2)]
depleted[rare, 1] <- 0

survivors <- list(survivors = rownames(counts)[-rare])
before <- ssTaxSEA(counts, custom_db = survivors, min_set_size = 3,
                   reference = "top_k", top_k = 10)
after <- ssTaxSEA(depleted, custom_db = survivors, min_set_size = 3,
                  reference = "top_k", top_k = 10)

# The survivors' score is unchanged with the top-k reference...
after[, 1] - before[, 1]
# ...but rises with the standard CLR
ssTaxSEA(depleted, custom_db = survivors, min_set_size = 3)[, 1] -
  ssTaxSEA(counts, custom_db = survivors, min_set_size = 3)[, 1]

## ----sessioninfo--------------------------------------------------------------
sessionInfo()