ANCOM-BC in BIOM

run_ancombc.R

library(biomformat)
library(ANCOMBC)
library(phyloseq)

tbl_file <- "data/shi_age_prediction/processed/processed_skin_tbl.biom"
tbl <- biomformat::read_biom(tbl_file)
tbl <- as.data.frame(as.matrix(biomformat::biom_data(tbl)))

md_file <- "data/shi_age_prediction/processed/processed_skin_md.tsv"
md <- read.table(md_file, sep="\t", row.names=1, header=T)
md$body_site <- relevel(as.factor(md$body_site), "head")
md$qiita_host_sex <- relevel(as.factor(md$qiita_host_sex), "male")

samples <- colnames(tbl)
md <- subset(md, rownames(md) %in% samples)
sample_order <- row.names(md)
tbl <- tbl[, sample_order]

taxa <- phyloseq::otu_table(tbl, taxa_are_rows=T)
meta <- phyloseq::sample_data(md)
physeq <- phyloseq::phyloseq(taxa, meta)

design_formula <- "body_site + qiita_host_sex"

ancombc.results <- ANCOMBC::ancombc(phyloseq=physeq, formula=design_formula,
                                    zero_cut=1.0)
results <- ancombc.results$res$beta

outfile <- "results/shi_age_prediction/ancombc_results_full.tsv"
write.table(results, file=outfile, sep="\t")