---
title: "NSF EEID Adult Newt Skin Secretion Experiments"
author: "Mitch Le Sage"
date: "2020"
output: html_document
---

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

## Load Packages & Data Prep

```{r prep, echo=TRUE}
library(ggplot2)
library(gridExtra)
library(bbmle)
library(emdbook)
library(dplyr)
library(tidyr)
library(tidyverse)
library(tibble)
library(knitr)
library(rmarkdown)
library(kableExtra)
library(stringr)
library(grid)
library(ggpubr)
library(dunn.test)
library(lawstat)

theme_mitch <- theme(panel.background=element_blank(), panel.grid.major=element_blank(), panel.grid.minor=element_blank(), panel.border=element_blank(), axis.line=element_line(color="black", size=1), text = element_text(colour = "black", size = 18),axis.text.x = element_text(color = "black", margin=unit(c(0,0.2,0.2,0.2), "cm")), axis.text.y = element_text(color = "black", margin=unit(c(0.2,0,0.2,0.2), "cm")), plot.title = element_text(hjust =))

AL <- read.csv(file.choose()) # choose "NOVISSData.csv" file that is provided
```

## Effects of Temperature on Total Skin Protein Fraction (Figure 6)

### Total Protein fraction recovery (Fig6a): µg of skin proteins recovered per gram body mass from adult NOVI at three temperatures

#### Data Prep
```{r 6a setup}
# TP Control Only Setup ####

TPCtlmaster <- filter(AL, Trtmt =="Control", Secretion=="Total Protein", InduceType!="Massage", InduceConc=="25.0nmol/g",  DPI=="000 dy")

# TP uggbw Control only####

TPTemp <- group_by(TPCtlmaster, Temp)
TPTempuggbw <- filter(TPTemp, !is.na(uggbw))
```
#### Data Summary
```{r 6a summary}
## summarize####

TPTempuggbwsum <- summarise(TPTempuggbw, count=n(), avg=mean(uggbw), standarddev=sd(uggbw), sem=sd(uggbw)/sqrt(n()))
F6aSum <- TPTempuggbwsum
F6aSum
```
#### Figure 6a
```{r 6a plot}
## plot means####

Fig6a <- ggplot(TPTempuggbwsum, aes(Temp, avg, fill=Temp)) + 
  geom_col() +   
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) +
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Protein Recovered(µg/gbw)") + 
  theme_mitch + 
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9))+
  scale_y_continuous(expand=c(0,0), limits=c(0,240))
Fig6a
```

#### Testing Assumptions & Transforming Data
```{r 6a assuptiontest}
## stats####

### test assumptions

levene.test(TPTempuggbw$uggbw, TPTempuggbw$Temp)
shapiro.test(TPTempuggbw$uggbw)

### transform data

TPTempuggbw$uggbw_log <- log10(TPTempuggbw$uggbw)

### test assumptions

levene.test(TPTempuggbw$uggbw_log, TPTempuggbw$Temp)
shapiro.test(TPTempuggbw$uggbw_log)
```
#### ANOVA
```{r 6a analyze}
### analyze data

TPTempuggbwanovalog <- aov(uggbw_log ~ Temp, data=TPTempuggbw)
F6aStats <- summary(TPTempuggbwanovalog)
F6aStats
```
#### PostHoc Test
```{r 6a posthoc}
F6aPostHoc <- TukeyHSD(TPTempuggbwanovalog)
F6aPostHoc
```

### Total Protein percent inhibition (Fig6b): percent of B. sal. growth inhibition by adult NOVI skin proteins (500µg/mL) at three temperatures measured using Cell Titre Glo assay

#### Data Prep
```{r 6b setup}
# TP % inhib Control only####

TPTempinhib <- filter(TPTemp, CellGloConc=="500", !is.na(CellGloPctInhib))
```
#### Data Summary
```{r 6b summary}
## summarize####

TPTempinhibsum <- summarise(TPTempinhib, count=n(), avg=mean(CellGloPctInhib), standarddev=sd(CellGloPctInhib), sem=sd(CellGloPctInhib)/sqrt(n()))
F6bSum <- TPTempinhibsum
F6bSum
```
#### Figure 6b
```{r 6b plot}
## plot means ####

Fig6b <- ggplot(TPTempinhibsum, aes(Temp, avg, fill=Temp)) + 
  geom_col() +  
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) + 
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Percent Inhibition") + 
  theme_mitch +
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9))+
  scale_y_continuous(expand=c(0,0), limits=c(0,100))
Fig6b
```

#### Testing Assumptions & Transforming Data
```{r 6b assumptiontest}
## stats####

### test assumptions

levene.test(TPTempinhib$CellGloPctInhib, TPTempinhib$Temp)
shapiro.test(TPTempinhib$CellGloPctInhib)
```
#### ANOVA
```{r 6b analyze}
### analyze data

TPTempinhibanova <- aov(CellGloPctInhib ~ Temp, data=TPTempinhib)
F6bStats <- summary(TPTempinhibanova)
F6bStats
```

### Total Protein Effectiveness (Fig6c): adult NOVI Total Protein effectiveness (µg produced * percent inhibition) at three temperatures

#### Data Prep
```{r 6c setup}
# TP Effectiveness Control only####

TPTempeffctl <- filter(TPTemp, CellGloConc=="500", !is.na(uggbw), !is.na(CellGloPctInhib))
TPTempeffctl$Eff <- TPTempeffctl$uggbw*TPTempeffctl$CellGloPctInhib
```
#### Data Summary
```{r 6c summary}
## summarize####

TPTempeffsumctl <- summarise(TPTempeffctl, count=n(), avg=mean(Eff), standarddev=sd(Eff), sem=sd(Eff)/sqrt(n()))
F6cSum <- TPTempeffsumctl
F6cSum
```
#### Figure 6c
```{r 6c plot}
### plot means####

Fig6c <- ggplot(TPTempeffsumctl, aes(Temp, avg, fill=Temp)) + 
  geom_col() +   
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) +
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Protein Effectiveness") + 
  theme_mitch + 
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9))+
  scale_y_continuous(expand=c(0,0), limits=c(0,12000))
Fig6c
```

#### Testing Assumptions & Transforming Data
```{r 6c assumptiontest}
## stats####

### test assumptions

levene.test(TPTempeffctl$Eff, TPTempeffctl$Temp)
shapiro.test(TPTempeffctl$Eff)

### transform data

TPTempeffctl$log_Eff <- log10(TPTempeffctl$Eff + 1)

### test assumptions

levene.test(TPTempeffctl$log_Eff, TPTempeffctl$Temp)
shapiro.test(TPTempeffctl$log_Eff)
```
#### Kruskal-Wallis Test
```{r 6c analyze}
### analyze data

F6cStats <- kruskal.test(Eff ~ Temp, data=TPTempeffctl)
F6cStats
```

## Effects of B.sal. exposure on Total Skin Protein and Hydrophobic Peptide production (Figure 7)

### Total Protein fraction recovery (Fig7a): µg of skin proteins recovered per gram body mass from adult NOVI exposed to different doses of B. sal. (all adult NOVI were kept at 16°C)

#### Data Prep
```{r 7a setup}
# TP 6C Only uggbw####
TPAllmaster <- filter(AL, Secretion=="Total Protein", InduceType!="Massage", InduceConc=="25.0nmol/g",  DPI=="000 dy")

# TP uggbw All Newts####

TPTempAll <- group_by(TPAllmaster, Temp)
TP6C <- filter(TPTempAll, Temp=="06°C", !is.na(uggbw))
TP6C <- group_by(TP6C, Trtmt)
```
#### Data Summary
```{r 7a summary}
## summarize####

TP6Csum <- summarize(TP6C, count=n(), avg=mean(uggbw), standarddev=sd(uggbw), sem=sd(uggbw)/sqrt(n()))
F7aSum <- TP6Csum
F7aSum
```
#### Figure 7a
```{r 7a plot}
## plot means####

Fig7a <- ggplot(TP6Csum, aes(Trtmt, avg, fill=Trtmt)) + 
  geom_col() + 
  theme(legend.position="none") +   
  scale_fill_manual(values=c("#B74641", "#730202", "#FB8383")) +
  xlab("Exposure Treatment(zsp)") +
  ylab("Recovered Proteins(µg/gbw)") + 
  theme_mitch + 
  geom_errorbar(aes(x=Trtmt, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9)) +
  scale_x_discrete(limits=c("Control", "5*10^3", "5*10^4"), labels=c("Control", "5,000", "50,000")) +
  scale_y_continuous(expand=c(0,0), limits=c(0,300))
Fig7a
```

#### Testing Assumptions & Transforming Data
```{r 7a assumptiontest}
## stats####

### test assumptions

levene.test(TP6C$uggbw, TP6C$Trtmt)
shapiro.test(TP6C$uggbw)

### transform data

TP6C$uggbw_log <- log10(TP6C$uggbw)

### test assumptions

levene.test(TP6C$uggbw_log, TP6C$Trtmt)
shapiro.test(TP6C$uggbw_log)
```
#### ANOVA
```{r 7a analyze}
### analyze data

TP6Canovalog <- aov(uggbw_log ~ Trtmt, data=TP6C)
F7aStats <- summary(TP6Canovalog)
F7aStats
```

### Hydrophobic Peptide fraction recovery (Fig7b): µg of hydrophobic peptides recovered per gram body mass from adult NOVI exposed to different doses of B. sal. (all adult NOVI were kept at 16°C)

#### Data Prep
```{r 7b setup}
# HP All Newts Setup ####

HPAllmaster <- filter(AL, Secretion=="Hydrophobic Peptide", InduceType!="Massage", InduceConc=="25.0nmol/g",  DPI=="000 dy")

# HP uggbw All Newts####

HPTempAll <- group_by(HPAllmaster, Temp)
HP6C <- filter(HPTempAll, Temp=="06°C", !is.na(uggbw))
HP6C <- group_by(HP6C, Trtmt)
```
#### Data Summary
```{r 7b summary}
## summarize####

HP6Csum <- summarize(HP6C, count=n(), avg=mean(uggbw), standarddev=sd(uggbw), sem=sd(uggbw)/sqrt(n()))
F7bSum <- HP6Csum
F7bSum
```
#### Figure 7b
```{r 7b plot}
## plot means####

Fig7b <- ggplot(HP6Csum, aes(Trtmt, avg, fill=Trtmt)) + 
  geom_col() +  
  theme(legend.position="none") +   
  scale_fill_manual(values=c("#B74641", "#730202", "#FB8383")) +
  xlab("Exposure Treatment(zsp)") +
  ylab("Recovered Peptides(µg/gbw)") + 
  theme_mitch + 
  geom_errorbar(aes(x=Trtmt, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9)) +
  scale_x_discrete(limits=c("Control", "5*10^3", "5*10^4"), labels=c("Control", "5,000", "50,000"))+
  scale_y_continuous(expand=c(0,0), limits=c(0,300))
Fig7b
```

#### Testing Assumptions & Transforming Data
```{r 7b assumptiontest}
## stats####

### test assumptions

levene.test(HP6C$uggbw, HP6C$Trtmt)
shapiro.test(HP6C$uggbw)

### transform data

HP6C$uggbw_log <- log10(HP6C$uggbw)

### test for normailty

levene.test(HP6C$uggbw_log, HP6C$Trtmt)
shapiro.test(HP6C$uggbw_log)
```
#### ANOVA
```{r 7b analyze}
### analyze data

HP6Canovalog <- aov(uggbw_log ~ Trtmt, data=HP6C)
F7bStats <- summary(HP6Canovalog)
F7bStats
```
#### PostHoc Test
```{r 7b posthoc}
F7bPostHoc <- TukeyHSD(HP6Canovalog)
F7bPostHoc
```

## Effects of Temperature on Hydrophobic Peptide Fraction (Figure S1)

### Hydrophobic Peptide fraction recovery (FigS1a): µg of hydrophobic peptides recovered per gram body mass from adult NOVI at three temperatures

#### Data Prep
```{r S1a setup}
# HP Control Only Setup #### 

HPCtlmaster <- filter(AL, Trtmt=="Control", Secretion=="Hydrophobic Peptide", InduceType!="Massage", InduceConc=="25.0nmol/g", DPI=="000 dy")

# HP uggbw Control only####

HPTemp <- group_by(HPCtlmaster, Temp)
HPTempuggbw <- filter(HPTemp, !is.na(uggbw))
```
#### Data Summary
```{r S1a summary}
## summarize####

HPTempuggbwsum <- summarise(HPTempuggbw, count=n(), avg=mean(uggbw), standarddev=sd(uggbw), sem=sd(uggbw)/sqrt(n()))
FS1aSum <- HPTempuggbwsum
FS1aSum
```
#### Figure S1a
```{r S1a plot}
## plot means####

FigS1a <- ggplot(HPTempuggbwsum, aes(Temp, avg, fill=Temp)) + 
  geom_col() +   
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) +
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Peptides Recovered(µg/gbw)") + 
  theme_mitch + 
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9)) +
  scale_y_continuous(expand=c(0,0), limits=c(0,240))
FigS1a
```

#### Testing Assumptions & Transforming Data
```{r S1a assumptiontest}
## stats####

### test assumptions

levene.test(HPTempuggbw$uggbw, HPTempuggbw$Temp)
shapiro.test(HPTempuggbw$uggbw)

### transform data

HPTempuggbw$uggbw_log <- log10(HPTempuggbw$uggbw)

### test assumptions

levene.test(HPTempuggbw$uggbw_log, HPTempuggbw$Temp)
shapiro.test(HPTempuggbw$uggbw_log)
```
#### ANOVA
```{r S1a analyze}
### analyze data

HPTempuggbwanovalog <- aov(uggbw_log ~ Temp, data=HPTempuggbw)
FS1aStats <- summary(HPTempuggbwanovalog)
FS1aStats
```

### Hydrophobic Peptide percent inhibition (FigS1b): percent of B. sal. growth inhibition by adult NOVI hydrophobic peptides (250µg/mL) at three temperatures measured using a 7-day Growth Inhibition assay

#### Data Prep
```{r S1b setup}
# HP % inhib @ 250ug/mL Control only####

HPTempinhib250 <- filter(HPTemp, !is.na(GIAPctInhib250))
```
#### Data Summary
```{r S1b summary}
## summarize####

HPTempinhibsum250 <- summarise(HPTempinhib250, count=n(), avg=mean(GIAPctInhib250), standarddev=sd(GIAPctInhib250), sem=sd(GIAPctInhib250)/sqrt(n()))
FS1bSum <- HPTempinhibsum250
FS1bSum
```
#### Figure S1b
```{r S1b plot}
## plot means####

FigS1b <- ggplot(HPTempinhibsum250, aes(Temp, avg, fill=Temp)) + 
  geom_col() +    
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) +
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Percent Inhibition") + 
  theme_mitch + 
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9)) +
  scale_y_continuous(expand=c(0,0), limits=c(0,100))
FigS1b
```

#### Testing Assumptions & Transforming Data
```{r S1b assumptiontest}
## stats####

### test for normality

levene.test(HPTempinhib250$GIAPctInhib250, HPTempinhib250$Temp)
shapiro.test(HPTempinhib250$GIAPctInhib250)
```
#### ANOVA
```{r S1b analyze}
### analyze

HPTempinhibanova250 <- aov(GIAPctInhib250 ~ Temp, data=HPTempinhib250)
FS1bStats <- summary(HPTempinhibanova250)
FS1bStats
```

### Hydrophobic Peptide Effectiveness (FigS1c): adult NOVI Hydrophobic Peptide effectiveness (µg produced * percent inhibition) at three temperatures

#### Data Prep
```{r S1c setup}
# HP Effectiveness Control only####

HPTempeffctrl <- filter(HPTemp, GIA250=="Y", !is.na(uggbw), !is.na(GIAPctInhib250))
HPTempeffctrl$Eff <- HPTempeffctrl$uggbw*HPTempeffctrl$GIAPctInhib250
```
#### Data Summary
```{r S1c summary}
## summarize####

HPTempeffsumctrl <- summarize(HPTempeffctrl, count=n(), avg=mean(Eff), standarddev=sd(Eff), sem=sd(Eff)/sqrt(n()))
FS1cSum <- HPTempeffsumctrl
FS1cSum
```
#### Figure S1c
```{r S1c plot}
## plot means####

FigS1c <- ggplot(HPTempeffsumctrl, aes(Temp, avg, fill=Temp)) +
  geom_col() +   
  scale_fill_manual(values=c("cornflowerblue", "goldenrod2", "firebrick")) +
  theme(legend.position="none") +
  xlab("Temperature(°C)") +
  ylab("Peptide Effectiveness") +
  theme_mitch +
  scale_x_discrete(limits=c("06°C", "14°C", "22°C"), labels=c("6", "14", "22")) +
  geom_errorbar(aes(x=Temp, ymin=avg-sem, ymax=avg+sem), width = 0.2, position=position_dodge(width=0.9)) +
  scale_y_continuous(expand=c(0,0), limits=c(0,12000))
FigS1c
```

#### Testing Assumptions & Transforming Data
```{r S1c assumptiontest}
## stats####

### test assumptions

levene.test(HPTempeffctrl$Eff, HPTempeffctrl$Temp)
shapiro.test(HPTempeffctrl$Eff)

### transform data

HPTempeffctrl$log_Eff <- log10(HPTempeffctrl$Eff + 1)

### test assumptions

levene.test(HPTempeffctrl$log_Eff, HPTempeffctrl$Temp)
shapiro.test(HPTempeffctrl$log_Eff)
```
#### ANOVA
```{r S1c analyze}
### analyze data

HPTempeffctrlanova <- aov(log_Eff ~ Temp, data=HPTempeffctrl)
FS1cStats <- summary(HPTempeffctrlanova)
FS1cStats
```