metabolome.rmd

---
title: "metabolome"
output: html_notebook
---

```{r}
library(mbtools)
library(phyloseq)
library(SummarizedExperiment)
theme_set(theme_minimal())
```

# meteome analysis

Let's start by loading some previously preprocessed metabolome data.

```{r}
mets <- readRDS("/proj/gibbons/toadnet/data/metabolome.rds")
mets
```

This includes data for all of Arivale so let's filter the sample we are interested in:

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

manifest <- cohort
met_samples <- colData(mets) %>% as.data.table()
matched <- met_samples[manifest, on = c("public_client_id", "days_in_program"), nomatch=0]
matched[, "subset" := "controls"]
matched[weight_change_relative < 0, subset:="weight loss"]
sdata <- as.data.frame(matched)
sdata$subset <- factor(sdata$subset)
sdata$sex <- factor(sdata$sex)
rownames(sdata) <- sdata$sample_id
sdata$high_bmi <- factor(as.numeric(sdata$bmi > 30))
mets <- mets[, mets$sample_id %in% matched$sample_id]
taxa <- matrix(rownames(mets), ncol=1)
colnames(taxa) <- "metabolite"
rownames(taxa) <- rownames(mets)
mat <- log2(t(assay(mets)))
mat[is.na(mat)] <- -9
ps <- phyloseq(
    otu_table(mat, taxa_are_rows = FALSE),
    sample_data(sdata),
    tax_table(taxa)
)
```

```{r}
tests <- association(
    ps, 
    presence_threshold = -5,
    min_abundance = -3,
    in_samples = 0.5,
    confounders = c("age", "bmi", "sex"),
    variables = "subset",
    taxa_rank = NA,
    method = "lm"
)
bmi <- association(
    ps, 
    presence_threshold = -5,
    min_abundance = -3,
    in_samples = 0.5,
    confounders = c("age", "sex"),
    variables = "bmi",
    taxa_rank = NA,
    method = "lm"
)
tests <- rbind(tests, bmi)
rowData(mets)$metabolite <- rownames(mets)
tests[, "metabolite" := variant]
tests <- merge.data.table(tests, rowData(mets), on="metabolite")
fwrite(tests[order(padj)], "data/tests_metabolites.csv")
tests
```

Volcano plots

```{r}

ggplot(tests, aes(log2FoldChange, y=-log10(pvalue), 
                  shape=padj<0.1, col=variable, size=padj<0.1)) +
    geom_point() + scale_size_discrete(range=c(2, 3)) +
    labs(shape="FDR < 0.1", size="FDR < 0.1")
ggsave("figures/metabolome_volcano.png", dpi=300, width=4, height=3)
```

Corr

```{r, fig.width=5, fig.height=3}
wide <- dcast(tests, metabolite ~ 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/metabolome_t.png", dpi=300, width=5, height=3)
cor.test(~ t_bmi + t_subset, data=wide)
```

## Overall explained variance

```{r}
library(vegan)

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