DAmixed.Rmd

---
title: "Differential Abundance for Repeated measures"
output: 
  html_document:
    number_sections: TRUE
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = F, warning = F, eval = T)
```

In a fair part of practical research designs including repeated sampling of the same individual occur. 
This is statistically a powerfull way to do research, as the effect of some pertubation (diet, medication, etc) is contrasted in the same individual. Usually, a measure at baseline is available, and then one or several after an intervention. Sometimes cross over designs is used, such that each individual is undergoing all treatments in the experiment. 

Modelling these data is in statistics refered to as mixed models. Mixed, due to the model having several terms modelling the random variation. For single univariate models, the _nlme_ or _lme4_ package has the functionallity. For multivariate analysis such as for the microbiome, there are dedicated methods which, as for the metagenomeSeq and DESeq2 methods, has parameter shrinkage to robustify the results. Here we use the [dream framework](https://www.bioconductor.org/packages/devel/bioc/vignettes/variancePartition/inst/doc/dream.html) (from Bioconductor). 

# Example data  {-}

The sausage intervention data is an example of repeated measures with baseline and after intervention data. 

```{r}
library(tidyverse)
library(phyloseq) 
library(variancePartition)
library(edgeR)
library(BiocParallel)

load('./data/Rats_inulin.RData')
phyX
metadata <- phyX %>% sample_data()
table(metadata$time,metadata$ID)
table(paste(metadata$Description,metadata$ID),metadata$time)
```
Here, the time==Start indicates pre-diet intervention, and hence here the rats should be similar regardless of diet. 

The model for the responses (the abundance of the microbiome) is:

$$ y_i = a(time_i) + b(time_i,diet_i) +\kappa(ID_i) + e_i$$
Here $a$ and $b$ describes the effect of time and time and diet in combination, while $\kappa$ reflects the difference between individual Rats, and $e$ reflects the residual uncertainty. 

## Setup the data  {-}

```{r}
phyX <- phyX %>% subset_samples(ID!='F1') # remove the F1 rat with no slut data. 

otutb <- as.data.frame(phyX@otu_table )
metadata <- phyX %>% sample_data() %>% data.frame() %>% 
  mutate(Diet = gsub('[[:punct:]]','',Description)) # we need to remove the + from the Description coloumn. 

GM = DGEList(otutb)
GM = calcNormFactors( GM)

# Specify parallel processing parameters
# this is used implicitly by dream() to run in parallel
param = SnowParam(4, "SOCK", progressbar=TRUE)
```

## Specify the model {-}

```{r}
# The variable to be tested must be a fixed effect
# form <- ~ time + time:Diet  + (1|ID)
form <- ~ 0 + time:Diet  + (1|ID)
```

## Run the model {-}


```{r}
# estimate weights using linear mixed model of dream
vobjDream = voomWithDreamWeights( GM, form, metadata, BPPARAM=param )

# Fit the dream model on each gene
# By default, uses the Satterthwaite approximation for the hypothesis test
fitmm = dream( vobjDream, form, metadata )
fitmm = eBayes(fitmm)
```

## Look at the results {-}

```{r}
topTable( fitmm) %>% data.frame() %>% rownames_to_column('asv')
```


This is overall effects. Obviously, diet changes the microbiome. 
However, what we are interested in is the contrasts. 

## Contrasts  {-}

This is refered to as contrasts. I.e. compare the effect from _Start_ to _Slut_ for Frankfurter with the same for chow; A change versus a change or also known as _Differences in Differences_

This is specified using the coloumn names from the output above.
Here the effect of Frankfurter is compare to the control diet (Chow). 

In the specification of the model above, we include 0 to remove the intercept, and : to avoid main effects in the interactionmodel. That produce an output from the the 2 x 3 design as just the flat means for each of the 6 groups (more or less). Establishing a difference from start to end of trial is hence a contrast of two terms, while a difference of differences is made up by four terms. 

```{r, fig.width=7}
L = makeContrastsDream( form, metadata, 
                        contrasts = c("timeSlut.DietFrankfurter - timeStart.DietFrankfurter - 
                                      timeSlut.DietChow + timeStart.DietChow", 
                                      
                                      "timeSlut.DietFrankfurterinulin - 
                                      timeStart.DietFrankfurterinulin - 
                                      timeSlut.DietChow + timeStart.DietChow"))

dimnames(L)[2]$Contrasts <-  c('C_FF','C_Inulin') # the labels are tooo long, so here shorter versions
plotContrasts(L)
```

## Fit the model

```{r}
fit = dream( vobjDream, form, metadata, L)
fit = eBayes(fit)
```

## Extract results and plot it

```{r}
res_Frankfurter <- topTable( fit, coef="C_FF", number=1000) # we just export all
res_Inulin <- topTable( fit, coef="C_Inulin", number=1000) # we just export all

TXtab <- phyX %>% tax_table() %>% data.frame() %>% # export TAxtable to glue onto results
  rownames_to_column('otu') 

res <- bind_rows(res_Frankfurter %>% mutate(contrast = 'Frankfurter vs Chow') %>% rownames_to_column('otu'), 
                 res_Inulin %>% mutate(contrast = 'FFinulin vs FF') %>% rownames_to_column('otu')) %>%
  left_join(TXtab, by = 'otu') # add Taxinfo


ggplot(data = res, aes(logFC, -log10(P.Value), label = otu, color = Rank2)) + 
  geom_point() + 
  geom_text(data = res[res$P.Value<0.000001,]) + 
  facet_wrap(~contrast)
```

... Lets plot one example OTU_6

```{r}

data.frame(metadata,otusel = as.numeric(otutb['OTU_6',])) %>% 
  ggplot(data =., aes(time,otusel + 1)) + geom_boxplot() + 
  geom_line(aes(group = ID)) + 
  scale_y_log10() + 
  facet_wrap(~Diet)

```


# Another Example

The study [Wastyk, H. C., Fragiadakis, G. K., Perelman, D., Dahan, D., Merrill, B. D., Feiqiao, B. Y., ... & Sonnenburg, J. L. (2021). Gut-microbiota-targeted diets modulate human immune status. Cell, 184(16), 4137-4153.](https://www.cell.com/cell/fulltext/S0092-8674(21)00754-6) has put all their data out on [github](https://github.com/SonnenburgLab/fiber-fermented-study)

try downloading the _phyloseq_obj_PilotStudy_log.rds_ from the data/16S folder and import it into R

```{r}
phyX <- readRDS('../MicrobiomeDataAnalysis/_site/data/phyloseq_obj_PilotStudy_log.rds')
phyX
```

... Take a look at the design

```{r}
table(paste(phyX@sam_data$Participant,phyX@sam_data$Group), phyX@sam_data$Timepoint)
```

It is a two armed longitudinal study with up to 9 timepoints. 


## Mixed model on Alpha-div

Before doing bug-per-bug analysis, the microbiome is condensed into one vector and used as response in an univariate mixed model. For the sake of the example, we use Faith phylogenetic diversity index as alpha diversity measure.

```{r}
library(picante)
FaithPD = picante::pd(samp = otu_table(phyX), tree = phy_tree(phyX), include.root = F)$PD
# combine with meta data
samdata <- data.frame(phyX@sam_data, FaithPD)

# plot it connecting the individuals with lines
ggplot(data = samdata, aes(Timepoint,FaithPD, group = Participant, color =Group)) + geom_point() + geom_line()


library(lme4)
library(lmerTest)
m <- lmer(data = samdata, FaithPD~factor(Timepoint)*Group + (1|Participant))
anova(m)
```

In line with the plot, there is not much of differences between groups and neither over time. 

```{r}
summary(m)
```

Looking at the random effect, we see that the variance related to individual is larger than the residual variation, meaning that there indeed is a lot of week to week variation conserved in the microbiome. 


Lets try DA on this dataset. Here, all code in one go. 

```{r}
phyXsel <- phyX %>% filter_taxa(function(x) sum(x > 0)>20, TRUE)

otutb <- as.data.frame(phyXsel@otu_table)
metadata <- phyXsel %>% sample_data() %>% data.frame() %>% 
  mutate(time = Timepoint %>% factor())

GM = otutb %>% t() %>% DGEList()
GM = calcNormFactors(GM)

# Specify parallel processing parameters
# this is used implicitly by dream() to run in parallel
param = SnowParam(4, "SOCK", progressbar=TRUE)


form <- ~ 0 + time:Group + (1|Participant)

# estimate weights using linear mixed model of dream
vobjDream = voomWithDreamWeights( GM, form, metadata, BPPARAM=param )

# Fit the dream model on each gene
# By default, uses the Satterthwaite approximation for the hypothesis test
fitmm = dream( vobjDream, form, metadata )
fitmm = eBayes(fitmm)
# topTable( fitmm) %>% data.frame() %>% rownames_to_column('asv')

L = makeContrastsDream( form, metadata, 
                        contrasts = c("time9.GroupFermented - time1.GroupFermented",
                                      "time9.GroupFiber - time1.GroupFiber",
                                      "time9.GroupFermented - time1.GroupFermented -
                                      time9.GroupFiber + time1.GroupFiber"))

dimnames(L)[2]$Contrasts <-  c('Eff_ferm','Eff_fiber','DDferm_fiber') # the labels are tooo long, so here shorter versions
# plotContrasts(L)

fit = dream( vobjDream, form, metadata, L)
fit = eBayes(fit)
```

Collect results and plot

```{r}
res_ferm <- topTable( fit, coef="Eff_ferm" ,number=1000) # we just export all
res_fiber <- topTable( fit, coef="Eff_fiber" ,number=1000) 
res_dd <- topTable( fit, coef="DDferm_fiber" ,number=1000) 

TXtab <- phyXsel %>% tax_table() %>% data.frame() %>% # export TAxtable to glue onto results
  rownames_to_column('otu') 

res <- bind_rows(res_ferm %>% mutate(contrast = 'Fermented') %>% rownames_to_column('otu'), 
                 res_fiber %>% mutate(contrast = 'Fiber') %>% rownames_to_column('otu'),
                 res_dd %>% mutate(contrast = 'Ferm - Fiber') %>% rownames_to_column('otu')) %>%
  left_join(TXtab, by = 'otu') # add Taxinfo


res %>% 
  mutate(Family2 = Family %>% factor() %>% fct_lump_n(5)) %>% 
  ggplot(data = ., aes(logFC, -log10(P.Value), label = otu, color = Family2)) + 
  geom_point() + 
  facet_wrap(~contrast)
```

```{r}
res %>% arrange(adj.P.Val) %>% 
  select(logFC,P.Value,adj.P.Val, B, contrast, Class:Species) %>% 
  head(10)
```

Plot the top flyers

```{r}
sel <- res %>% arrange(adj.P.Val) %>% head(3) # dig the strongest results out

# combine with meta data
xx <- otutb[,colnames(otutb) %in% sel$otu] %>% 
  data.frame() %>% 
  cbind(metadata)

# plot it
xx %>% 
  gather(otu,abu,-c(SampleID:time)) %>% 
  ggplot(data = ., aes(time,abu, fill = Group)) + 
  geom_boxplot() + scale_y_sqrt() + 
  facet_wrap(~otu)
```