Add a qmd for stats panel

stats_panel.qmd

@@ -0,0 +1,98 @@
+---
+title: "Stats at Taxo"
+format: html
+---
+
+```{r}
+library(knitr)
+printer_tabset <- function(x, options, ...) {
+    if (is.null(names(x))) names(x) <- seq_along(x)
+
+    header <- ":::: {.panel-tabset}"
+    footer <- "::::"
+
+    res <- lapply(seq_along(x), function(i) {
+        knitr::knit_child(
+            text = c(
+                "##### `r names(x)[i]`",
+                "",
+                "```{r}",
+                "#| echo: false",
+                "x[[i]]",
+                "```"
+            ),
+            envir = environment(),
+            quiet = TRUE
+        )
+    })
+
+    out <- paste(c(header, res, footer), collapse = "\n\n")
+    knitr::asis_output(out)
+}
+registerS3method("knit_print", "list", printer_tabset)
+```
+
+```{r}
+source("utils.R")
+library(biomformat)
+library(phyloseq)
+library(ggplot2)
+library(tidyverse)
+library(R.utils)
+library(sbm)
+library(blockmodels)
+
+the_data <- import_biom("data/chaillou/chaillou.biom")
+
+taxa_names <- c("Kingdom", "Phylum", "Class", "Order", "Family", "Genus", "Species/OTU")
+
+taxo <- tax_table(the_data)
+per_taxa_network <- collapse_otu_at_taxo(the_data)
+names(per_taxa_network) <- taxa_names
+```
+
+```{r}
+#| freeze: true
+plot_mat <- function(matrix, rank) {
+    dat2 <-
+        matrix %>%
+        as_tibble(rownames = "Var1") %>%
+        pivot_longer(-Var1, names_to = "Var2", values_to = "value") %>%
+        mutate(Var1 = stringr::str_split(Var1, "-", simplify = TRUE)[, rank])
+    p <- ggplot(dat2, aes(Var1, Var2)) +
+        geom_tile(aes(fill = value)) +
+        # geom_text(aes(label = round(value, 1)), data = dat2 %>% filter(value > 0)) +
+        scale_fill_gradientn(
+            colours = c("white", "#FEECE3", "red"),
+            values = c(0, .Machine$double.eps, 1)
+        ) + # To set 0 in pure white and smallest non zero as off-white
+        theme_bw() +
+        labs(x = "Taxa", y = "Samples", title = paste0(length(unique(dat2$Var1)), " different taxa")) +
+        theme(
+            axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
+            axis.text.y = element_blank()
+        )
+
+    return(p)
+}
+otu_table_plot <- lapply(seq_along(per_taxa_network), function(idx) plot_mat(per_taxa_network[[idx]], rank = idx))
+names(otu_table_plot) <- taxa_names
+```
+
+# Count matrix plot
+```{r}
+otu_table_plot
+```
+
+# Standard deviation
+
+```{r}
+kables_sd <- lapply(per_taxa_network, function(matrix) {
+    sd_vec <- as.matrix(apply(matrix, 1, sd))
+    sd_vec <- sd_vec[order(rownames(sd_vec)), ]
+    tax_df_names <- do.call("rbind", rownames(matrix) %>% str_split("-")) %>% as.data.frame()
+    sd_df <- cbind(tax_df_names, sd_vec)
+    kable(sd_df, row.names = FALSE)
+})
+kables_sd
+```