49 lines
1.3 KiB
R
49 lines
1.3 KiB
R
library(blockmodels)
|
|
library(phyloseq)
|
|
library(biomformat)
|
|
library(here)
|
|
source("utils.R")
|
|
|
|
options(parallelly.availableCores.custom = function() {
|
|
ncores <- min(parallelly::availableCores(), 64)
|
|
max(1L, ncores)
|
|
})
|
|
|
|
message(paste("Number of cores available:", parallelly::availableCores()))
|
|
args <- commandArgs(trailingOnly = TRUE)
|
|
if (length(args) == 0) {
|
|
author <- "mach"
|
|
} else {
|
|
author <- args[1]
|
|
}
|
|
|
|
message("For author: ", author)
|
|
the_data <- import_biom(here("data", author, paste0(author, ".biom")))
|
|
|
|
epoch <- as.integer(Sys.time())
|
|
level <- 2
|
|
message("At level: ", level, " starting at: ", epoch)
|
|
|
|
save_path <- here("results", "seq-effort", paste0(author, "_", level, "_", epoch, ".Rds"))
|
|
|
|
if (!dir.exists(save_path)) {
|
|
dir.create(save_path, recursive = TRUE)
|
|
}
|
|
|
|
per_taxa_networks <- collapse_otu_at_taxo(the_data)
|
|
|
|
covar_mat <- t(colSums(per_taxa_networks[[level]]) %*% t(rep(1, nrow(per_taxa_networks[[level]]))))
|
|
|
|
lcovar_mat <- log(covar_mat)
|
|
|
|
fitted_model <- BM_poisson_covariates("LBM",
|
|
adj = per_taxa_networks[[level]],
|
|
covariates = list(lcovar_mat),
|
|
autosave = save_path,
|
|
ncores = parallelly::availableCores()
|
|
)
|
|
fitted_model$estimate()
|
|
|
|
fitted_model$model_parameters[[which.max(fitted_model$ICL)]]
|
|
message("Saving the model results.")
|
|
saveRDS(fitted_model, save_path)
|