diff --git a/lbm_phyloseq_test.R b/lbm_phyloseq_test.R index 918d16b..68ba345 100644 --- a/lbm_phyloseq_test.R +++ b/lbm_phyloseq_test.R @@ -3,23 +3,26 @@ library(ggplot2) library(sbm) library(biomformat) +source("utils.R") + # data("enterotype") +# the_data <- enterotype # data("mach") -data <- import_biom("data/chaillou/chaillou.biom") -net_table <- otu_table(data) +the_data <- import_biom("data/chaillou/chaillou.biom") +# the_data <- import_biom("data/mach/kinetic.biom") +# the_data <- import_biom("data/ravel/ravel.biom") +net_tables <- collapse_otu_at_taxo(the_data) +net_table <- as.matrix(net_tables[[3]]) 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"))) +res <- estimateBipartiteSBM(net_table, dimLabels = c("OTU", "Sample"), model = "poisson", estimOptions = list(plot = 0)) plot(res, type = "data") -plot(res_t, type = "data") - sdCol <- function(matrix) { sapply(seq_len(ncol(matrix)), function(col) sd(matrix[, col])) } @@ -37,17 +40,3 @@ ggplot(stats_df, aes(x = mean, y = var)) + 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)")