---
title: 'EEID Adult Newts Temperature Experiments: Microbiomes (16S & ITS)'
author: "Molly Bletz"
date: "2020"
output:
  html_document:
    smart_extension: no
    theme: sandstone
    highlight: tango
    toc: false
    toc_float: false
    df_print: paged
    code_folding: show
  pdf_document: default
editor_options:
  chunk_output_type: console
---

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

# {.tabset .tabset-fade}

>**Packages**

```{r 1, echo=TRUE}
library(tidyverse)
library(ggplot2)
library(ggthemes)
library(RColorBrewer)
library(lme4)
library(qiime2R)
library(ecodist)
library(vegan)
```


## 16S Bacterial Communities {.tabset .tabset-fade}

### Data Compiling {.tabset .tabset-fade}

> **Data Compiling**

**Metadata merging of antifungal function data**


```{r}
Metadata_AntiFungal_PCO_Alpha = read_tsv("16S_Metadata_AlphaDiv_AntiFungalPred_TimePre.txt")

```


### Alpha diversity {.tabset .tabset-fade}

*Alpha diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads*

>**Alpha Models**

**sOTU Richness models**

```{r}

# dist check
hist(Metadata_AntiFungal_PCO_Alpha$observed_otus)

# main model
mod_sOTU_1 <- glmer(observed_otus ~ Temperature + (1|CollectionLocation), family = poisson(link = "log"),data = Metadata_AntiFungal_PCO_Alpha)

car::Anova(mod_sOTU_1, type=c("III"), test.statistic=c("Chisq"))

# model fit
qqnorm(resid(mod_sOTU_1))
qqline(resid(mod_sOTU_1))

# post hoc test
emmeans::emmeans(mod_sOTU_1, list(pairwise ~ Temperature), adjust = "tukey")

```

**Phylogenetic diversity models**

```{r 9}

# dist check
hist(Metadata_AntiFungal_PCO_Alpha$faith_pd)

# main model
mod_PD_1 <-glmer(faith_pd ~ Temperature + (1|CollectionLocation), family = Gamma, data = Metadata_AntiFungal_PCO_Alpha)

car::Anova(mod_PD_1, type=c("III"), test.statistic=c("Chisq"))

qqnorm(resid(mod_PD_1))
qqline(resid(mod_PD_1))

# Post hoc test
emmeans::emmeans(mod_PD_1, list(pairwise ~ Temperature), adjust = "tukey")
```

>**Alpha diversity across temperatures plot**

```{r 7a}

trt_order <- c("6", "14", "22")

Metadata_AntiFungal_PCO_Alpha %>% 
  ggplot(aes(x = factor(Temperature, level = trt_order), y = observed_otus, color = factor(Temperature, level = trt_order))) +
  geom_violin()+
  geom_point(size = 3, alpha = 0.6) +
  scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  theme_classic() +
  xlab("Temperature (C)") +
  ylab("sOTU Richness") +
  ggtitle("Richness of adult newt bacterial communities across temperatures")

```

### Beta diversiy {.tabset .tabset-fade}


>**Beta diversity models**

*Beta diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads*

>**Weighted Unifrac models**


```{r}

#read in dist matrix
AdNewt_WUF<-read_qza("16S_weighted_unifrac_distance_matrix.qza")

AdNewt_WUF_mx <- as.matrix(AdNewt_WUF$data)

# filtering dist matrix and metadata to same samples
AdNewtT0_distList2<-colnames(AdNewt_WUF_mx)

Metadata_AntiFungal_PCO_Alpha_T0_BD2<- Metadata_AntiFungal_PCO_Alpha %>%
  filter(SampleID %in% AdNewtT0_distList2)

metaList<-as.vector(Metadata_AntiFungal_PCO_Alpha_T0_BD2$SampleID)
AdNewt_WUF_mx2 <- usedist::dist_subset(AdNewt_WUF_mx,metaList)
AdNewt_WUF_mx3 <-as.matrix(AdNewt_WUF_mx2)

# main model
adonis(AdNewt_WUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2)

```

**Pairwise models**
```{r}
# install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

# pairwise post hoc
WUFmodel <- pairwiseAdonis::pairwise.adonis(AdNewt_WUF_mx3, Metadata_AntiFungal_PCO_Alpha_T0_BD2$Temperature, p.adjust.m = "bonferroni")


```

>**Unweighted Unifrac models**

```{r}

# read in dist matrix

AdNewt_unWUF<-read_qza("16S_unweighted_unifrac_distance_matrix.qza")
AdNewt_unWUF_mx <- as.matrix(AdNewt_unWUF$data)

# filtering dist matrix and metadata to same samples

AdNewtT0_distList3<-colnames(AdNewt_unWUF_mx)
Metadata_AntiFungal_PCO_Alpha_T0_BD3<- Metadata_AntiFungal_PCO_Alpha %>% 
 filter(SampleID %in% AdNewtT0_distList3)

metaList3<-as.vector(Metadata_AntiFungal_PCO_Alpha_T0_BD3$SampleID)
AdNewt_unWUF_mx2 <- usedist::dist_subset(AdNewt_unWUF_mx,metaList3)
AdNewt_unWUF_mx3 <-as.matrix(AdNewt_unWUF_mx2)

# main model
adonis(AdNewt_unWUF_mx3 ~ Temperature, data = Metadata_AntiFungal_PCO_Alpha_T0_BD2)

```

**Pairwise models**

```{r}

# post hoc
pairwiseAdonis::pairwise.adonis(as.dist(AdNewt_unWUF_mx3), Metadata_AntiFungal_PCO_Alpha_T0_BD2$Temperature, p.adjust.m = "bonferroni")

```

>**PCoA Plots**

```{r}

trt_order <- c("6", "14", "22")

uWUF_16S_PCoA <- ape::pcoa(AdNewt_unWUF_mx3)

PCoA_Axes16S = as.data.frame(uWUF_16S_PCoA$vectors)
PCoA_Axes16S =tibble::rownames_to_column(PCoA_Axes16S, "SampleID")

PCoA_Axes16S %>%
  dplyr::select(SampleID, Axis.1, Axis.2) %>%
  left_join(Metadata_AntiFungal_PCO_Alpha) %>%
  filter(TimeWeek=="0")%>%
  ggplot(aes(x=Axis.1, y=Axis.2, color = as.factor(Temperature),size = Propor_TotalAntiFungal_BSALonly)) +
  geom_point(aes(color = as.factor(Temperature), alpha=0.5))+
   scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  stat_ellipse(type = "t",linetype = 2)+
   theme_classic() +
  xlab("PCo1 - 4.5%") +
  ylab("PCo2 - 2.8%")

```

### Taxa Plots

**Taxonomic Composition**

*Grouped sOTU tables with one column per temperature were calculated in QIIME2 with the feature-table group command using a individual sample OTU table rarified to 5000 reads and converted to a tab-delimited file using biom convert*

```{r}
TempTable2 = read_tsv("16S_TempGrouped_rarefied_table_r5000.txt")

Taxa = read_tsv("16S_Taxonomy.txt") %>% parse_taxonomy() %>%
  rownames_to_column("OTUID")

TempTable2Taxa = TempTable2 %>%
  left_join(.,Taxa, by = "OTUID")
```

```{r}

TempTable2Taxa_Order = TempTable2Taxa %>%
  group_by(Order) %>%
  summarize(AvgAbund6 = sum(t6),AvgAbund14 = sum(t14),AvgAbund22 = sum(t22), totalAbund = sum(t14+t6+t22))

## select rows that are greater than 300 in total abundance --> move to new table
TempTable2Taxa_Order_Top = TempTable2Taxa_Order %>%
  filter(totalAbund>=300) %>%
  select( -totalAbund)

TempTable2Taxa_Order_Top <- as.data.frame(TempTable2Taxa_Order_Top)

## select rows that are less than 300 and move to new table --> sum all to new table
TempTable2Taxa_Order_Low = TempTable2Taxa_Order %>%
  filter(totalAbund<300) %>%
  summarize(AvgAbund6 = sum(AvgAbund6),AvgAbund14 = sum(AvgAbund14),AvgAbund22 = sum(AvgAbund22)) %>%
  mutate(Order = "Other")
TempTable2Taxa_Order_Low <- as.data.frame(TempTable2Taxa_Order_Low)

## merge to tables and make it long format
TempTable2Taxa_Order_Summarize <- full_join(TempTable2Taxa_Order_Top,TempTable2Taxa_Order_Low)


TempTable2Taxa_Order_Summarize_long = TempTable2Taxa_Order_Summarize %>%
  gather(key = Temperature, value = Abundance, -Order)

TempTable2Taxa_Order_Summarize_long_prop = TempTable2Taxa_Order_Summarize_long %>% 
  group_by(Temperature) %>% 
   mutate(total = sum(Abundance), rel.freq = Abundance / total)

```

```{r}
trt_order2 <- c("AvgAbund6", "AvgAbund14", "AvgAbund22")

TempTable2Taxa_Order_Summarize_long_prop %>%
  ggplot(aes(x = factor(Temperature, level = trt_order2), y = rel.freq, fill = Order)) +
  geom_bar(stat = "identity")+
  scale_fill_stata() +
  ylab("Relative Abundance")+
  xlab("Temperature") +
  scale_x_discrete(labels=c("AvgAbund6" = "6", "AvgAbund22" = "22", "AvgAbund14" = "14"))

```


## ITS Fungal Communities {.tabset .tabset-fade}

### Alpha Diversity {.tabset .tabset-fade}


*Alpha diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads*

```{r}

ITSMeta_Alpha = read_tsv("ITS_Metadata_AlphaDiv_TimePre.txt")

```

> **Alpha Diversity models**

*Note: included Collection location as a random effect in alpha diversity models.*

```{r}

#dist check
hist(ITSMeta_Alpha$observed_otus)

#main model
mod_sOTU_ITS1 <- glmer(observed_otus ~ Temperature + (1|CollectionLocation), family = Gamma,data = ITSMeta_Alpha)

car::Anova(mod_sOTU_ITS1, type=c("III"), test.statistic=c("Chisq"))

qqnorm(resid(mod_sOTU_ITS1))
qqline(resid(mod_sOTU_ITS1))

# post hoc test
emmeans::emmeans(mod_sOTU_ITS1, list(pairwise ~ Temperature), adjust = "tukey")

```

> **Alpha Diversity plot**

```{r}

Temptrt_order <- c("6", "14", "22")

ITSMeta_Alpha %>% 
  filter(observed_otus<200) %>%
  ggplot(aes(x = factor(Temperature, level = Temptrt_order), y = observed_otus, color = factor(Temperature, level = Temptrt_order))) +
  geom_violin()+
  geom_point(size = 3 , alpha = 0.6) +
  scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  theme_classic() +
  xlab("Temperature (C)") +
  ylab("ITS sOTU Richness") +
  ggtitle("Richness of adult newt fungal communities across temperatures")

```

### Beta diversity {.tabset .tabset-fade}

**Beta diversity**

*Beta diversity metrics were calculated in QIIME2 with the diversity core-metrics-phylogenetic command with a rarefaction depth of 5000 reads*

```{r}

#read in dist matrix
AdNewt_JacITSmx<-read_qza("ITS_jaccard_distance_matrix.qza")

AdNewt_JacITSmx2 <- as.matrix(AdNewt_JacITSmx$data)

#filtering distance matrix and metadata

AdNewtT0_distList2<-colnames(AdNewt_JacITSmx2)
ITSMeta_ObsOTUS_AntiF_0b<- ITSMeta_Alpha %>%
  filter(SampleID %in% AdNewtT0_distList2)

metaList<-as.vector(ITSMeta_ObsOTUS_AntiF_0b$SampleID)

AdNewt_Jac_mx2c <- usedist::dist_subset(AdNewt_JacITSmx2,metaList)
AdNewt_JacITSmx3 <-as.matrix(AdNewt_JacITSmx2)

# main model
adonis(AdNewt_JacITSmx3 ~ Temperature, data = ITSMeta_Alpha)

# dispersion model
AdNewt_JacITSmx3_disp<-betadisper(as.dist(AdNewt_JacITSmx3), ITSMeta_Alpha$Temperature, type = c("median"), bias.adjust = FALSE, sqrt.dist = FALSE, add = FALSE)

anova(AdNewt_JacITSmx3_disp)

```

**Pairwise models**

```{r}

#devtools::install_github("pmartinezarbizu/pairwiseAdonis/pairwiseAdonis")

model <- pairwiseAdonis::pairwise.adonis(AdNewt_JacITSmx3, ITSMeta_Alpha$Temperature, p.adjust.m = "bonferroni")

```


**Beta diversity plot**
```{r 10a}

#read in dist matrix

PCoA_ITS <- ape::pcoa(AdNewt_JacITSmx3)

PCoA_AxesITS = as.data.frame(PCoA_ITS$vectors)
PCoA_AxesITS =tibble::rownames_to_column(PCoA_AxesITS, "SampleID")

PCoA_AxesITS %>%
  dplyr::select(SampleID, Axis.1, Axis.2) %>%
  left_join(ITSMeta_Alpha) %>%
  filter(TimeWeek=="0")%>%
  ggplot(aes(x=Axis.1, y=Axis.2, color = as.factor(Temperature))) +
  geom_point(aes(color = as.factor(Temperature), alpha=0.5))+
   scale_color_manual(values = c("#087CAF", "#E7B800", "#DB6E0B"),name = "Temperature (C)") +
  stat_ellipse(type = "t",linetype = 2)+
   theme_classic() +
  xlab("PCo1 - 4.5%") +
  ylab("PCo2 - 2.8%")

```

### Taxa Plots {.tabset .tabset-fade}

**ITS Taxonomic Composition**

*Grouped sOTU tables with one column per temperature were calculated in QIIME2 with the feature-table group command using a individual sample OTU table rarified to 1000 reads and converted to a tab-delimited file using biom convert*

```{r }
TempTable2ITS = read_tsv("ITS_TempGrouped_rarefied_table_r1000.txt")

TaxaITS = read_tsv("ITS_Taxonomy.txt") %>% parse_taxonomy() %>%
  rownames_to_column("OTUID")

TempTable2ITSTaxa = TempTable2ITS %>%
  left_join(.,TaxaITS, by = "OTUID")
```

```{r}

TempTable2ITSTaxa_Order = TempTable2ITSTaxa %>%
  group_by(Order) %>%
  summarize(AvgAbund6 = sum(t6),AvgAbund14 = sum(t14),AvgAbund22 = sum(t22), totalAbund = sum(t14+t6+t22))

## select rows that are greater than 100 in total abundance --> move to new table
TempTable2ITSTaxa_Order_Top = TempTable2ITSTaxa_Order %>%
  filter(totalAbund>=100) %>%
  select( -totalAbund)

TempTable2ITSTaxa_Order_Top <- as.data.frame(TempTable2ITSTaxa_Order_Top)

## select rows that are less than 100 and move to new table --> sum all to new table
TempTable2ITSTaxa_Order_Low = TempTable2ITSTaxa_Order %>%
  filter(totalAbund<100) %>%
  summarize(AvgAbund6 = sum(AvgAbund6),AvgAbund14 = sum(AvgAbund14),AvgAbund22 = sum(AvgAbund22)) %>%
  mutate(Order = "Other")
TempTable2ITSTaxa_Order_Low <- as.data.frame(TempTable2ITSTaxa_Order_Low)

## merge to tables and make it long format
TempTable2ITSTaxa_Order_Summarize <- full_join(TempTable2ITSTaxa_Order_Top,TempTable2ITSTaxa_Order_Low)


TempTable2ITSTaxa_Order_Summarize_long = TempTable2ITSTaxa_Order_Summarize %>%
  gather(key = Temperature, value = Abundance, -Order)

TempTable2ITSTaxa_Order_Summarize_long_prop = TempTable2ITSTaxa_Order_Summarize_long %>% 
  group_by(Temperature) %>% 
   mutate(total = sum(Abundance), rel.freq = Abundance / total) %>% 
  replace_na(list(Order = "Unknown Fungi"))

```

```{r}
trt_order2 <- c("AvgAbund6", "AvgAbund14", "AvgAbund22")

as.data.frame(TempTable2ITSTaxa_Order_Summarize_long_prop) %>%
  ggplot(aes(x = factor(Temperature, level = trt_order2), y = rel.freq, fill = Order)) +
  geom_bar(stat = "identity")+
  scale_fill_stata() +
  ylab("Relative Abundance")+
  xlab("Temperature") +
  scale_x_discrete(labels=c("AvgAbund6" = "6", "AvgAbund22" = "22", "AvgAbund14" = "14")) +
  theme_classic()

```