fEPI_metabolome.Rmd

---
title: 'Untargeted Analysis of the Serum Metabolome in Cats with Exocrine Pancreatic Insufficiency: A Pilot Study'
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, fig.align="center")
library(knitr)
library(dplyr)
library(kableExtra)
library(pheatmap)
library(fgsea)
library(factoextra)
library(EnhancedVolcano)
library(omu)
library(rstatix)
library(tidyr)
library(corrr)
library(RColorBrewer)
library(viridis)
library(visdat)
library(purrr)
library(ggpubr)
library(made4)
library(qvalue)
library(FactoMineR)
library(reshape2)
library(ggplus)

```

#**Load and Clean Data**

```{r}

#load scaled, imputed metabolite abundance data

metab <- read.csv("fEPI_metab_scaledimp.csv", header = T)

#relabel rows with arbitrary IDs for each metabolite 
#removes cumbersome chemical names

for (i in metab) {
  met_label <- paste("met", seq(1:length(i)), sep = "")
}

rownames(metab) <- met_label

#make an metabolite annotation table
annot <- metab[, 1:12]

#now isolate metabolite abundance data for analysis
metab <- metab[, 13:25]

#load sample metadata

meta <- read.csv("fEPI_metadata.csv", header = T)

str(meta)

meta$Group <- factor(meta$Group, levels = c("Healthy", "EPI"))

meta$Sex <- factor(meta$Sex, levels = c("MN", "FS"))

rownames(meta) <- meta[[1]]

#sanity check

names(metab) == meta[[1]]

rownames(meta) == names(metab)

sum(rownames(annot) == rownames(metab))

#list of microbial metabolites and co-metabolites from Metabolon Inc.

microbe <- read.csv("microbial_metabolites.csv", header = T)

#log-transform

metab_norm <- log10(metab)

```



#**Sunmmary Statistics for Sample Metadata**

```{r}

fTLI_stats <- meta %>%
  group_by(Group) %>%
  get_summary_stats(fTLI)

fTLI_stats

```

```{r}

hist(meta$fTLI)

shapiro.test(meta$fTLI)

```

```{r}

wilcox.test(fTLI~Group, data = meta)

ggplot(data = meta, aes(x = Group, y = fTLI, fill = Group)) + 
  geom_boxplot() + 
  geom_jitter() +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) +
  labs(fill = "Group", x = "Group", y = "Serum fTLI")

```

```{r}

Age <- meta %>%
  group_by(Group) %>%
  get_summary_stats(Age)

Age

```

```{r}
hist(meta$Age)

shapiro.test(meta$Age)
```

```{r}

wilcox.test(Age~Group, data = meta)

ggplot(data = meta, aes(x = Group, y = Age, fill = Group)) + 
  geom_boxplot() +
  geom_jitter() + 
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) +
  labs(fill = "Group", x = "Group", y = "Age")

```

# **Unsupervised Analysis of Metabolite Abundance Data**

## **Heatmap**

```{r}

annotation_colors <- list(Group = c(Healthy = "#00BFC4", EPI = "#F8766D"))

annotation_col <- as.data.frame(meta[, 2])

names(annotation_col)[1] <- "Group"

rownames(annotation_col) <- rownames(meta)

heatmap <- pheatmap(
  mat               = metab_norm,
  color             = inferno(10),
  border_color      = NA,
  show_colnames     = T,
  show_rownames     = FALSE,
  annotation_col    = annotation_col,
  annotation_colors = annotation_colors,
  cluster_cols = T,
  drop_levels       = TRUE,
  scale             = "row",
  fontsize          = 12,
  annotation_legend = F
)

heatmap

```


## **Principal Component Analysis**

```{r}

#make new df
dat_pca <- data.frame(t(metab_norm))
dat_pca$Group <- meta$Group

#perform PCA
pca <- PCA(dat_pca, quali.sup = 734, graph = F)

fviz_eig(pca)

#make a beautiful plot using ggplot2

dat_pca$pc1 <- pca$ind$coord[, 1] # indexing the first column

dat_pca$pc2 <- pca$ind$coord[, 2]  # indexing the second column

pca.vars <- pca$var$coord %>% data.frame

pca.vars$vars <- rownames(pca.vars)

pca.vars.m <- melt(pca.vars, id.vars = "vars")

pca_plot <- ggplot(data = dat_pca, aes(x = pc1, y = pc2, color = Group)) +
  scale_color_manual(values = c("#00BFC4", "#F8766D")) +
  geom_hline(yintercept = 0, lty = 2) +
  geom_vline(xintercept = 0, lty = 2) +
  guides(color = guide_legend(title = "Group")) +
  scale_shape_manual(values = c(15, 16, 16, 17, 18)) +
  geom_point(alpha = 0.8, size = 2) 


pca_plot <- pca_plot + stat_ellipse(geom="polygon", aes(fill = Group, color = Group), 
                      alpha = 0.2, 
                      show.legend = FALSE, 
                      level = 0.95) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) +
  xlab("PC 1 (29.34%)") + 
  ylab("PC 2 (17.85%)") +
  theme_minimal() +
  theme(panel.grid = element_blank(), 
        panel.border = element_rect(fill= "transparent")) +
  theme(legend.position="bottom") + 
  theme(legend.title=element_text(size=14), 
    legend.text=element_text(size=12))

pca_plot

```


# **Hypothesis Testing**

T-tests will be used on the scaled and normalized metabolite abundance values to test for differences between the healthy and fEPI groups. 

```{r, echo=FALSE}

metab_diff <- metab
metab_diff$Metabolite <- rownames(metab)

test <- omu_summary(count_data = metab_diff, 
                    metadata = meta, 
                    numerator = "EPI", 
                    denominator = "Healthy",
                    response_variable = "Metabolite",
                    Factor = "Group",
                    log_transform = T,
                    p_adjust = "fdr",
                    test_type = "welch")

test$FDR <- qvalue(test$pval)$qvalues

res <- test[, c(1:12, 26)]

res$Biochemical <- annot$BIOCHEMICAL[match(res$Metabolite, rownames(annot))]

res$SubPathway <- annot$SUB.PATHWAY[match(res$Metabolite, rownames(annot))]

res$SuperPathway <- annot$SUPER.PATHWAY[match(res$Metabolite, rownames(annot))]

```

Metabolites will be considered significant when:

* P<0.05
* Q<0.2 (local false discovery rate - FDR)
* |Log2Fold-Change| > 0.585 (equivalent to a |effect size| of 1.5)

```{r}

#threshold on p-value
sig <- res[res$pval < 0.05, ]

#threshold to select only results with an effect size of +/- 1.5 (log2FC +/- 0.585)
sig <- sig[sig$log2FoldChange < -0.585 | sig$log2FoldChange > 0.585, ]

#sig <- sig[, c(5, 7, 6, 2:4)]

sig <- sig[order(sig$SuperPathway, sig$log2FoldChange), ]

# sig$log2FoldChange <- round(sig$log2FoldChange, 3)
# sig$pval <- round(sig$pval, 3)
# sig$FDR <- round(sig$FDR, 3)

sig_inc <- sig[sig$log2FoldChange > 0, ]
rownames(sig_inc) <- c()

sig_dec <- sig[sig$log2FoldChange < 0, ]
rownames(sig_dec) <- c()

inc <- as.data.frame(table(sig_inc$SuperPathway))
names(inc)[1] <- "Super Pathway"
names(inc)[2] <- "Increased Metabolites"

sum(inc$`Increased Metabolites`)


dec <- as.data.frame(table(sig_dec$SuperPathway))
names(dec)[1] <- "Super Pathway"
names(dec)[2] <- "Decreased Metabolites"

sum(dec$`Decreased Metabolites`)

```

Plot significant results by super pathway

```{r}
inc2 <- inc
names(inc2)[2] <- "Significant_Changes"
inc2$colour <- "Increased"
inc2 <- inc2[!inc2$Significant_Changes == 0, ]

dec2 <- dec
names(dec2)[2] <- "Significant_Changes"
dec2$colour <- "Decreased"
dec2 <- dec2[!dec2$Significant_Changes == 0, ]

freq_tbl2 <- rbind(inc2, dec2)

freq_tbl2$colour <- factor(freq_tbl2$colour, 
                           levels = c("Increased", "Decreased"))

super_path.plot <- ggplot(data = freq_tbl2, aes(y = Significant_Changes, 
                                               x = as.factor(`Super Pathway`), 
                                               fill = colour)) + 
  geom_bar(stat = "identity", color = "black") + 
  labs(fill = "") +
  xlab("") +
  ylab("Number of Significant Metabolites") + 
  scale_fill_brewer(palette="Set1") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(legend.position = "right") + theme(legend.title.align=0.5) +
  theme(axis.text.x=element_text(angle=45,hjust=1))

super_path.plot
```

This will generate a pdf containing box plots of all metabolites found to vary significantly between cats with EPI and healthy controls
```{r, include=F}

sig_metab <- as.data.frame(t(metab))

names(sig_metab) <- annot$BIOCHEMICAL[match(names(sig_metab), rownames(annot))] 

sig_names <- as.character(sig$Biochemical)

sig_metab <- sig_metab[, sig_names]

sig_metab$Group <- meta$Group

cols <- c("#00BFC4", "#F8766D")

base_plot <-sig_metab %>%
  gather(Measure, Value, -Group) %>%
  ggplot(aes(x = factor(Group), y = Value, fill = Group)) +
  geom_boxplot() +
  geom_jitter()  +
  theme(plot.title = element_text(size = 6)) +
  scale_fill_manual(values = c("#00BFC4", "#F8766D")) +
  facet_wrap(~Measure, 
             scales = "free_y") +
  xlab("Group") +
  ylab("Metabolite Abundance (Median-Scaled))") 

pdf("S_Fig_sigBoxplots.pdf")
facet_multiple(base_plot, "Measure",
               nrow = ,
               ncol = 2,
               scales = "free_y")
dev.off()

```

# **Enrichment Analysis**

```{r}
##rank metabolites by the product of the log2FC and log(p-val)
#this generates a vector of metabolites ranked according a metric derived from significance testing
#the top end of the list are highly significant AND have high positive fold change
#the bottem end of the list are highly significant AND have negative fold change

diff_enrich <- res[res$pval < 0.05, ]

diff_enrich$rank <- -log(diff_enrich$pval)*diff_enrich$log2FoldChange

diff_enrich <- diff_enrich[order(diff_enrich$rank, decreasing = T), ]

#create ranking vector from product of log2FC*pval

rank <- diff_enrich[, 17]
names(rank) <- diff_enrich$Biochemical
# head(rank)
# tail(rank)

#Now generate lists of metabolites in each metabolic sub-pathway

pathways <- diff_enrich[, 14:15]

pathways <- pathways[order(pathways$SubPathway, pathways$Biochemical),]

#pathways$Biochemical <- as.character(pathways$Biochemical)
pathways$SubPathway <- as.character(pathways$SubPathway)

pathways <- split(pathways$Biochemical, as.character(pathways$SubPathway))

class(pathways)

fgseaRes <- fgsea(pathways = pathways, 
                  stats = rank,
                  minSize=2,
                  maxSize=500,
                  nperm=100000)

fgseaRes_out <- data.frame(fgseaRes)


```

Plot pathway enrichment
```{r}

topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=17), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))


topPathways <- unlist(topPathways)

fgseaRes_top20 <- fgseaRes[fgseaRes$pathway %in% topPathways, ]

fgseaRes_top20 <- fgseaRes_top20[order(fgseaRes_top20$NES), ]

fgseaRes_top20$pathway <- as.factor(fgseaRes_top20$pathway)

fgseaRes_top20$pathway <- factor(fgseaRes_top20$pathway, 
                                 levels = fgseaRes_top20$pathway[order(fgseaRes_top20$NES)])

gsea_plot <- ggplot(data = fgseaRes_top20, aes(y = fgseaRes_top20$NES, 
                                               x = fgseaRes_top20$pathway, 
                                               fill = ifelse(fgseaRes_top20$pval <0.05, 
                                                             "P<0.05", "P≥0.05"))) +
  geom_bar(stat = "identity", position = "dodge") + ylim(-3, 3) + coord_flip()

gsea_plot <- gsea_plot + 
  labs(fill = "") +
  xlab("") +
  ylab("Normalized Enrichment Score") + scale_fill_manual(values=c("orange", "grey")) + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(legend.position = "right") + theme(legend.title.align=0.5)
      
gsea_plot


```

Make Boxplots of leading edge metabolites in enriched pathways
```{r}
names_gsea <- fgseaRes[fgseaRes$padj < 0.2, ]

names_gsea <- sapply(names_gsea$leadingEdge, function(x) x[1])

metab_gsea <- as.data.frame(t(metab_norm))

annot_gsea <- annot[annot$BIOCHEMICAL %in% names_gsea, ]

metab_gsea <- metab_gsea[, names(metab_gsea) %in% rownames(annot_gsea)]

names(metab_gsea) <- annot_gsea$BIOCHEMICAL[match(rownames(annot_gsea), names(metab_gsea))]

metab_gsea$Group <- meta$Group


wrapper <- function(x, ...) {
  paste(strwrap(x, ...), collapse = " ")
}

boxplot_fun = function(data, x, y) {
  dat = data 
  ggplot(dat, aes(x = .data[[x]], y = .data[[y]], fill = Group) ) + 
    scale_fill_manual(values = c("#00BFC4", "#F8766D")) +
    geom_boxplot() + 
    geom_jitter() + 
    theme_bw() +
    labs(x = "", 
         y = "",
         fill = "Group") + 
  theme(legend.position="bottom") +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  theme(plot.title = element_text(hjust = 0.5)) + 
    ggtitle(wrapper(y, width = 8))  + 
    theme(plot.title = element_text(size = 8)) + 
    theme(legend.title=element_text(size=12), 
        legend.text=element_text(size=12))
}

names(metab_gsea) <- gsub("\\s*\\[[^\\)]+\\]","", as.character(names(metab_gsea)))
names(metab_gsea) <- gsub("\\*","", as.character(names(metab_gsea)))

length(unique(names(metab_gsea)))

length(names(metab_gsea))

metab_gsea$Group <- factor(metab_gsea$Group, levels = c("Healthy", "EPI"))

metab_plots <- map(names(metab_gsea[c(2, 3, 5, 6, 4, 7, 1)]), ~boxplot_fun(metab_gsea, "Group", .x))

metab_plots

msea.out <- data.frame(fgseaRes[order(fgseaRes$padj, fgseaRes$NES), ])

msea.out$leadingEdge <- sapply(msea.out$leadingEdge, function(x) x[1])

msea.out <- unnest(msea.out, leadingEdge)

```

# **Microbial Metabolites**

Metabolon Inc. maintains a library of compounds known to be microbial metabolites or microbe-host co-metabolites. This library was cross-referncef against the significant results to determine which metabolites that varied between groups were also microbial metabolites. 

Exclusive microbial metabolites:

```{r}

micro_metab <- microbe$Biochemical[which(microbe$Microbial == "microbe")]

comp <- comparelists(annot$BIOCHEMICAL, micro_metab)

comp2 <- comparelists(micro_metab, sig$Biochemical)
sig_microbe_metabs <- comp2$intersect

metab_microbe <- as.data.frame(t(metab_norm))

annot_microbe <- annot[annot$BIOCHEMICAL %in% sig_microbe_metabs, ]

metab_microbe <- metab_microbe[, names(metab_microbe) %in% rownames(annot_microbe)]

unique(annot_microbe$SUB.PATHWAY)

names(metab_microbe) <- annot_microbe$BIOCHEMICAL[match(names(metab_microbe), 
                                                     rownames(annot_microbe))]

metab_microbe$Group <- factor(meta$Group, levels = c("Healthy", "EPI"))

microbe_plots <- map(names(metab_microbe[1:8]), ~boxplot_fun(metab_microbe, "Group", .x))

microbe_plots

# annotate_figure(ggarrange(plotlist = microbe_plots, common.legend = T, legend = "bottom", labels = LETTERS[1:8]),
#                 left = text_grob("Normalized Metabolite Abundance", color = "black", rot = 90, size = 20)
#                )

```

# **Correlation Analysis**

## **Correlations Between Serum Metabolites and Serum fTLI**

```{r}

metab_cor <- data.frame(t(metab_norm))

metab_cor$fTLI <- meta$fTLI

cor.res <- metab_cor %>% cor_test(fTLI, method = "pearson")

cor.res$FDR <- qvalue(cor.res$p)$qvalues

cor.res$Biochemical <- annot$BIOCHEMICAL[match(cor.res$var2, rownames(annot))]

cor.res$SubPathway <- annot$SUB.PATHWAY[match(cor.res$var2, rownames(annot))]

cor.res$SuperPathway <- annot$SUPER.PATHWAY[match(cor.res$var2, rownames(annot))]

cor.res <- cor.res[order(cor.res$SubPathway, cor.res$cor, cor.res$p), ]


```


```{r}

cor_sig <- cor.res[cor.res$cor > 0.7 | cor.res$cor < -0.7, ]

cor_sig <-cor_sig[cor_sig$p < 0.05, ]

cor_sig$Biochemical <- annot$BIOCHEMICAL[match(cor_sig$var2, rownames(annot))]


```

```{r, warning=F, message=F}

names(metab_cor) <- annot$BIOCHEMICAL[match(names(metab_cor), rownames(annot))]

sig_names <- cor_sig$Biochemical

cor_df <- metab_cor[, sig_names]

oldnames <- names(cor_df)[1:79]

newnames <- as.character(cor_sig$Biochemical)

#length(oldnames) == length(newnames)

names(cor_df)[match(oldnames, names(cor_df))] <- newnames

cor_df$fTLI <- meta$fTLI

cor_plot <- cor_df[, 1:80] %>%
  gather(-fTLI, key = "var", value = "value") %>%
  ggplot(aes(x = value, y = fTLI)) + 
  facet_wrap(~var, scales = "free") +
  geom_point() +
  geom_smooth(method="lm") +
  xlab("Metabolite Anundance (Median-Scaled))") +
  ylab("Serum fTLI (µg/L)")

pdf("S_Fig_cor_plots.pdf")
facet_multiple(cor_plot, "var",
               nrow = 2, 
               ncol = 2, 
               scales = "free")
dev.off()
```