polarolouis/human-microbiome-compendium
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

lbm on phyloseq

Browse source
This commit is contained in:
Louis 2026-01-28 15:55:12 +01:00
parent 05e9bca613
commit 2fc1cf4e65
1 changed files with 9 additions and 20 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

29
lbm_phyloseq_test.R
View file

@ -3,23 +3,26 @@ library(ggplot2)
library(sbm) library(sbm)
library(biomformat) library(biomformat)
source("utils.R")
# data("enterotype") # data("enterotype")
# the_data <- enterotype
# data("mach") # data("mach")
data <- import_biom("data/chaillou/chaillou.biom") the_data <- import_biom("data/chaillou/chaillou.biom")
net_table <- otu_table(data) # 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) dim(net_table)
nb_inter <- sum(net_table > 0) nb_inter <- sum(net_table > 0)
all_possible_inter <- nrow(net_table) * ncol(net_table) all_possible_inter <- nrow(net_table) * ncol(net_table)
(density <- nb_inter / all_possible_inter) (density <- nb_inter / all_possible_inter)
res <- estimateBipartiteSBM(net_table, dimLabels = c("OTU", "Sample"), model = "poisson") res <- estimateBipartiteSBM(net_table, dimLabels = c("OTU", "Sample"), model = "poisson", estimOptions = list(plot = 0))
res_t <- estimateBipartiteSBM(t(net_table), dimLabels = rev(c("OTU", "Sample")))
plot(res, type = "data") plot(res, type = "data")
plot(res_t, type = "data")
sdCol <- function(matrix) { sdCol <- function(matrix) {
sapply(seq_len(ncol(matrix)), function(col) sd(matrix[, col])) 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_y_log10() +
scale_x_log10() + scale_x_log10() +
ggtitle("OTU for log(mean) = log(var)") 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)")