## ----install, eval=FALSE------------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#     install.packages("BiocManager")
# }
# BiocManager::install("strainspy")

## ----setup, message=FALSE, warning=FALSE--------------------------------------
library(strainspy)

meta_path <- system.file("extdata", "example_metadata.csv.gz", package = "strainspy")
sylph_path <- system.file("extdata", "example_sylph_profile.tsv.gz", package = "strainspy")
taxonomy_path <- system.file("extdata", "example_taxonomy.tsv.gz", package = "strainspy")

meta <- read.csv(meta_path)
meta$Case_status <- factor(meta$Case_status)

taxonomy <- read_taxonomy(taxonomy_path)
se <- read_sylph(sylph_path, meta_data = meta)

# Recommended default: keep strains observed in at least 10% of samples
min_nonzero <- ceiling(0.10 * ncol(se))
se <- filter_by_presence(se, min_nonzero = min_nonzero)

cat("Samples:", ncol(se), "\n")
cat("Retained strains:", nrow(se), "\n")
cat("min_nonzero used:", min_nonzero, "\n")

## ----fit, message=FALSE, warning=FALSE----------------------------------------
design <- as.formula("~ Case_status + Age_at_collection")

fit <- glmZiBFit(
  se,
  design,
  nthreads = 2L,
  return_vcov = FALSE
)

fit

## ----top-hits, message=FALSE, warning=FALSE-----------------------------------
hits <- top_hits(fit, coef = 2, alpha = 1)
head(hits[, c("Contig_name", "coefficient", "p_value", "p_adjust")])

## ----hadjust, message=FALSE, warning=FALSE------------------------------------
hier_hits <- hadjust(fit, coef = 2, taxonomy = taxonomy)
head(hier_hits[, c("Level", "Model", "Name", "p_adjust")])

## ----plot-data, message=FALSE, warning=FALSE----------------------------------
plot_volcano(fit, coef = 2)
plot_manhattan(fit, coef = 2, taxonomy = taxonomy)

## ----utilities, message=FALSE, warning=FALSE----------------------------------
rich <- estimate_sample_richness(se, taxonomy)
beta_dist <- estimate_beta_diversity(se, distance_only = TRUE)

pred_df <- prep_for_prediction(
  sy = se,
  outcome = "Case_status",
  contigs = rownames(se)[1:3],
  covariates = "Age_at_collection",
  use_genome_names = FALSE
)

length(rich$species_richness)
class(beta_dist)
dim(pred_df)

## ----other-inputs, message=FALSE, warning=FALSE-------------------------------
example_sourmash_path <- system.file("extdata", "example_sourmash.csv.gz", package = "strainspy")
example_metaphlan_path <- system.file("extdata", "metaphlan_merged.tsv.gz", package = "strainspy")

sm <- read_sourmash(example_sourmash_path, meta)
mp <- read_metaphlan(example_metaphlan_path, meta)

class(sm)
class(mp)

## ----session-info-------------------------------------------------------------
sessionInfo()