---
title: "MD3F Application to Kostic Microbiome"
author: "Alexander V. Alekseyenko"
date: "`r date()`"
output: html_document
---
1. Load packages & Data.
```{r}
library(plyr)
library(phyloseq)
library(dplyr)
library(glmnet)
library(ggplot2)
library(gridExtra)
```

```{r}
load("kostic_withsampda.Rdata")
```

2. Attach sample data to species level kostic phyloseq object. 
```{r}
kostic_2_full = merge_phyloseq(kostic_species, sample_data(kostsamptab))
#Report number of sample fore each subject included. 
table(sample_data(kostic_2_full)$Subject_ID)
```
3. For convenience in analysis, create index for samples.
```{r}
kostic_2_full = subset_samples(kostic_2_full)
sample_data(kostic_2_full)$index = 1:nsamples(kostic_2_full)
```

4. Examine sample sums for each sample in study (across all subjects). 
```{r}
plot(sample_sums(kostic_2_full))
```

5. Filter the kostic samples to sums greater than 0.8. 
```{r}
kostic_2_full = subset_samples(kostic_2_full, sample_sums(kostic_2_full)>0.8)
```

6. MD3F Application to Kostic samples. 
```{r}
subjects = unique(sample_data(kostic_2_full)$Subject_ID)
plots = list()
fits = list()
results = data.frame()

for(s in subjects){
  subj.data = subset_samples(kostic_2_full, 
                   Subject_ID == s)
  ns = nsamples(subj.data)
  if(ns<3) #Each subject must have a minimum of 3 longitudinal samples collected. 
    next
#Create Distance matrix. 
  EuclidDistCLRmatrix =
  as.matrix(
  distance(transform_sample_counts(subj.data,
                                  function(x) {y=log(x/sum(x)+10^-4); y-mean(y)}),
          method="euclidean"))
  
  to = order(sample_data(subj.data)$Age_at_Collection, decreasing = F)
  
  dT = sample_data(subj.data)$Age_at_Collection[to[-1]] - sample_data(subj.data)$Age_at_Collection[to[1]]
  dX = EuclidDistCLRmatrix[to[1], to[-1]]

  # Limit to first year of life
  keep = dT<365
  dT = dT[keep]
  dX = dX[keep]
  
  resName = s
  fits[[resName]] = lm(dX^2 ~ dT^2)

  
    plots[[resName]] = 
      ggplot(data.frame(dTSquared=dT^2, dXSquared=dX^2), aes(x=dTSquared, y=dXSquared)) + 
      geom_point(show.legend = FALSE) +
      geom_smooth(method="lm")+
      ylab("Normalized distance squared") +
      xlab("Time interval squared") +
      ggtitle(paste(s, " (", sample_data(subj.data)$Case_Control[1],")", sep=""))
    
    results = rbind(results,
                    data.frame(
                      resname = resName,
                      subject = s,
                      caseStatus = sample_data(subj.data)$Case_Control[1],
                      speed = coef(fits[[resName]])[2],
                      speed.low = confint(fits[[resName]])[2,1],
                      speed.high = confint(fits[[resName]])[2,2],
                      p.value = summary(fits[[resName]])$coefficients[2,4]
                    ))

}

with(results, tapply(speed, caseStatus, mean))
```










