Existing data

library(arivale.data.interface)
library(mbtools)

classes <- c(public_client_id = "character", microbiome_id = "character")
cohort <- rbind(
    fread("no_weight_loss.csv", colClasses=classes),
    fread("successful_weight_loss.csv", colClasses=classes)
)[since_baseline == 0]
cohort[, "subset" := "controls"]
cohort[weight_change_relative < 0, "subset" := "weight loss"]
cohort[, subset := factor(subset)]
cohort[, "vendor_observation_id" := microbiome_id]

ps <- readRDS("/proj/arivale/microbiome/16S_processed/phyloseq.rds")
ps <- subset_samples(ps, vendor_observation_id %in% cohort[, vendor_observation_id])
ps

phyloseq-class experiment-level object
otu_table()   OTU Table:          [ 89022 taxa and 105 samples ]:
sample_data() Sample Data:        [ 105 samples by 24 sample variables ]:
tax_table()   Taxonomy Table:     [ 89022 taxa by 8 taxonomic ranks ]:
phy_tree()    Phylogenetic Tree:  [ 89022 tips and 82828 internal nodes ]:
taxa are columns

genus <- speedyseq::tax_glom(ps, "Genus")

sdata <- merge(
    as(sample_data(genus), "data.frame"),
    as.data.frame(cohort),
    by = c("vendor_observation_id", "sex")
)
rownames(sdata) <- sdata$id
sdata$age <- sdata$age.y
sdata$subset <- factor(sdata$subset)
sdata$sex <- factor(sdata$sex)
sample_data(genus) <- sdata

tests <- association(
    genus, 
    presence_threshold = 1,
    min_abundance = 10,
    in_samples = 0.5,
    confounders = c("age", "bmi", "sex"),
    variables = "subset",
    taxa_rank = "Genus",
    method = "deseq2",
    shrink=F
)
bmi <- association(
    genus, 
    presence_threshold = 1,
    min_abundance = 10,
    in_samples = 0.5,
    confounders = c("age", "sex"),
    variables = "bmi",
    taxa_rank = "Genus",
    method = "deseq2",
    shrink=F
)
tests <- rbind(tests, bmi)
tests <- merge.data.table(tests, as.data.frame(as(tax_table(genus), "matrix")), on="Genus")
fwrite(tests[order(padj)], "data/tests_16S_genus.csv")
tests

tests[, "t" := log2FoldChange / lfcSE]
wide <- dcast(tests, Genus ~ variable, value.var=c("t", "padj"), fill=0)
wide[, "sig" := "none"]
wide[padj_bmi < 0.05, "sig" := "BMI"]
wide[padj_subset < 0.05, "sig" := "weight loss"]
wide[padj_bmi < 0.05 & padj_subset < 0.05, "sig" := "both"]
wide[, "sig" := factor(sig, levels=c("none", "BMI", "weight loss", "both"))]
ggplot(wide, aes(x=t_bmi, y=t_subset, color=sig)) +
    geom_vline(xintercept = 0, color="gray30", lty="dashed") +
    geom_hline(yintercept = 0, color="gray30", lty="dashed") +
    geom_point() + stat_smooth(method="glm", aes(group = 1)) +
    labs(x = "t statistic BMI", y = "t statistic weight loss") +
    scale_color_manual(values=c(none = "black", BMI = "royalblue", `weight loss`="orange", both = "red"))

ggsave("figures/16S_genus_t.png", dpi=300, width=5, height=3)
cor.test(~ t_bmi + t_subset, data=wide, method="spearman")

    Spearman's rank correlation rho

data:  t_bmi and t_subset
S = 70764, p-value = 0.1301
alternative hypothesis: true rho is not equal to 0
sample estimates:
      rho 
0.1706048 

Overall explained variance

library(vegan)

m <- as(otu_table(genus), "matrix")
m[is.na(m)] <- min(as.numeric(m), na.rm=T)
perm <- adonis(m ~ age + sex + bmi + subset, data=as(sample_data(genus), "data.frame"), method="bray")
perm

Call:
adonis(formula = m ~ age + sex + bmi + subset, data = as(sample_data(genus),      "data.frame"), method = "bray") 

Permutation: free
Number of permutations: 999

Terms added sequentially (first to last)

           Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)    
age         1    0.3211 0.32105  2.0219 0.01876  0.026 *  
sex         1    0.1541 0.15408  0.9703 0.00900  0.429    
bmi         1    0.6281 0.62814  3.9558 0.03671  0.001 ***
subset      1    0.1283 0.12828  0.8078 0.00750  0.619    
Residuals 100   15.8789 0.15879         0.92802           
Total     104   17.1104                 1.00000           
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1