## ----knitr_options, include = FALSE-------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    message = FALSE
)

## ----setup--------------------------------------------------------------------
library(phylobar)   # phylogeny-aware bar plots
library(dplyr)      # data wrangling
library(ape)        # phylogenetic tree manipulation
library(stringr)    # string operations
library(zen4R)      # download from zenodo link

## ----data---------------------------------------------------------------------
download_zenodo("10.5281/zenodo.18791960", tempdir())
HFHSdata <- readRDS(str_c(tempdir(), "/HFHS-data.rds"))
normal <- HFHSdata$OTUdata_Normal
hfhs <- HFHSdata$OTUdata_HFHS
taxa <- HFHSdata$filtered_taxonomy

## ----tree_construction--------------------------------------------------------
tree <- taxa |>
    select(-X1, X1) |>
    mutate(
        across(
            everything(),
            ~if_else(str_ends(., "_"), NA, .)
        )
    ) |>
    taxonomy_to_tree()

checkValidPhylo(tree)

## ----ordering-----------------------------------------------------------------
# Extract day 7 samples and rename
normal_day7 <- normal[,, 4]
rownames(normal_day7) <- str_replace(rownames(normal_day7), "M_", "N")
hfhs_day7 <- hfhs[,, 4]
rownames(hfhs_day7) <- str_replace(rownames(hfhs_day7), "M_", "H")
# Combine the two datasets
all_day7 <- rbind(normal_day7, hfhs_day7)

# Order samples and species using hierarchical clustering
comp_norm <- normal_day7
comp_norm <- comp_norm / rowSums(comp_norm)
sample_order1 <- hclust(dist(comp_norm))$order
species_order <- hclust(dist(t(comp_norm)))$order

# Now order the HFHS samples
comp_hfhs <- hfhs_day7
comp_hfhs <- comp_hfhs / rowSums(comp_hfhs)
sample_order2 <- hclust(dist(comp_hfhs))$order

# Combine the sample orders
sample_order <- c(rownames(normal_day7)[sample_order1], rownames(hfhs_day7)[sample_order2])

# Reorder the combined data
x <- all_day7[sample_order, species_order]

## ----visualize----------------------------------------------------------------
comp <- x / rowSums(x)
colnames(comp) <- colnames(all_day7)[species_order]
phylobar(comp, tree, hclust_order = FALSE, sample_font_size = 9)

## ----fig1, echo=FALSE, fig.align='center', fig.cap="Phylobar plot of the mouse microbiome on day 7, comparing normal diet (N) and HFHS diet (H). The plot highlights the relative abundances of Bacteroidia class, which is known to decrease in response to an HFHS diet."----
knitr::include_graphics("Figures/fig1.png")

## ----fig2, echo=FALSE, fig.align='center', fig.cap="Phylobar plot of the mouse microbiome on day 7, comparing normal diet (N) and HFHS diet (H). The plot highlights the relative abundances of Firmicutes phylum, which is known to increase in response to an HFHS diet."----
knitr::include_graphics("Figures/fig2.png")

## ----melt---------------------------------------------------------------------
y <- hfhs[sample_order2, species_order, ] |>
    aperm(c(1, 3, 2)) |>
    matrix(nrow = nrow(hfhs)*dim(hfhs)[3] , ncol = ncol(hfhs))

rownames(y) <- apply(
    expand.grid(c(0, 1, 4, 7), str_replace(rownames(hfhs),"M_","")[sample_order2]),
    1, paste, collapse = "-"
)

## ----visualize_longitudinal---------------------------------------------------
comp <- y / rowSums(y)
colnames(comp) <- colnames(hfhs)[species_order]
phylobar(comp, tree, hclust_order = FALSE, sample_font_size = 9)

# p <- phylobar(
#     comp, tree, hclust_order = FALSE,
#     width = 1200, height = 800, sample_font_size = 8
# )
# htmlwidgets::saveWidget(
#     p, "interactive_scatter.html", selfcontained = TRUE
# )

## ----fig3, echo=FALSE, fig.align='center', fig.cap="Phylobar plot of the mouse microbiome for HFHS diet. The plot highlights the relative abundances of S24-7 family that belong to the Bacteroidetes phylum, which is known to decrease over time in response to an HFHS diet."----
knitr::include_graphics("Figures/fig3.png")

## ----longitudinal_reorder-----------------------------------------------------
time <- as.numeric(sub("-.*", "", rownames(comp)))
mouse <- factor(
    sub(".*-", "", rownames(comp)),
    levels = unique(sub(".*-", "", rownames(comp)))
)
comp_sort <- comp[order(time, mouse), ]
phylobar(comp_sort, tree, hclust_order = FALSE, sample_font_size = 9)

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