diff --git a/utils.R b/utils.R
index ff4b4f8..3bc52e8 100644
--- a/utils.R
+++ b/utils.R
@@ -1,5 +1,7 @@
 library(knitr)
 library(kableExtra)
+library(dplyr)
+library(stringr)
 
 ellipse_table <- function(df, nr = 5, nc = 5, format = "latex") {
     nrows <- nrow(df)
@@ -56,3 +58,29 @@ ellipse_table <- function(df, nr = 5, nc = 5, format = "latex") {
     # Print with kable
     return(sub_df)
 }
+
+collapse_otu_at_taxo <- function(phylo_data, NArm = TRUE) {
+    tax_ranks <- colnames(tax_table(phylo_data))
+    tax_df <- tax_table(phylo_data)
+    otu_tables_by_rank <- lapply(tax_ranks, function(rank) {
+        ps_glom <- tax_glom(the_data, taxrank = rank, NArm = NArm)
+        otu <- as.data.frame(otu_table(ps_glom))
+
+        # Ensure taxa are rows
+        if (!taxa_are_rows(ps_glom)) {
+            otu <- t(otu)
+        }
+        otu_matrix <- as.matrix(otu)
+        id_in_taxo <- match(rownames(otu_matrix), rownames(tax_df))
+        taxo_names <- tax_df[id_in_taxo, seq(which(tax_ranks == rank))] %>%
+            as.data.frame() %>%
+            mutate_if(.predicate = is_character, .funs = partial(str_replace_all, pattern = "-", replacement = ".")) %>%
+            unite(col = "tax_full", sep = "-") %>%
+            unique()
+        rownames(otu_matrix) <- as.vector(taxo_names[, 1])
+
+        return(otu_matrix)
+    })
+    names(otu_tables_by_rank) <- tax_ranks
+    return(otu_tables_by_rank)
+}