diff --git a/lbm_phyloseq_test.R b/lbm_phyloseq_test.R
new file mode 100644
index 0000000..918d16b
--- /dev/null
+++ b/lbm_phyloseq_test.R
@@ -0,0 +1,53 @@
+library(phyloseq)
+library(ggplot2)
+library(sbm)
+library(biomformat)
+
+# data("enterotype")
+# data("mach")
+data <- import_biom("data/chaillou/chaillou.biom")
+net_table <- otu_table(data)
+dim(net_table)
+
+nb_inter <- sum(net_table > 0)
+all_possible_inter <- nrow(net_table) * ncol(net_table)
+(density <- nb_inter / all_possible_inter)
+
+res <- estimateBipartiteSBM(net_table, dimLabels = c("OTU", "Sample"), model = "poisson")
+res_t <- estimateBipartiteSBM(t(net_table), dimLabels = rev(c("OTU", "Sample")))
+
+plot(res, type = "data")
+
+plot(res_t, type = "data")
+
+sdCol <- function(matrix) {
+    sapply(seq_len(ncol(matrix)), function(col) sd(matrix[, col]))
+}
+
+t_net_table <- t(net_table)
+std_otu <- sdCol(t_net_table)
+mean_otu <- rowMeans(net_table)
+
+stats_df <- data.frame(otu = rownames(net_table), mean = mean_otu, std = std_otu, var = std_otu^2)
+
+ggplot(stats_df, aes(x = mean, y = var)) +
+    geom_line() +
+    geom_point() +
+    geom_abline(slope = 1, color = "red") +
+    scale_y_log10() +
+    scale_x_log10() +
+    ggtitle("OTU for log(mean) = log(var)")
+
+ggplot(stats_df, aes(x = mean, y = var / mean)) +
+    geom_line() +
+    geom_point() +
+    geom_hline(yintercept = 1, color = "red") +
+    ggtitle("OTU for f(mean) = var/mean = 1")
+
+ggplot(stats_df, aes(x = 1 / mean, y = var / mean^2)) +
+    geom_line() +
+    geom_point() +
+    geom_abline(slope = 1, color = "red") +
+    scale_y_log10() +
+    scale_x_log10() +
+    ggtitle("OTU for log(1/mean) = log(var/mean^2)")