## VIRUS INOCULATION 

# Load packages
library(ggplot2)
library(plyr)
library(plotrix)
library(RColorBrewer)
library(gplots)

library(lme4)
library("lmerTest", lib.loc="~/Library/R/3.2/library")
library(plyr)

# Import data (cell densities 3 days after virus inoculation)
VirusRaw3d <- read.csv("~/Desktop/Figures & Analyses/VirusRaw3d.csv")
View(VirusRaw3d)

# Calculate means of flow cytometer repliactes (canto rep)
VirusData<-ddply(VirusRaw3d, .(Environment, Line, ResType, Treatment, Biorep, Techrep), summarize, CellsPerMl=mean(Cells))

# Calculate mean for each biorep
VirusPlot<-ddply(VirusData, .(Environment, Line, ResType, Treatment, Biorep), summarize, CellsPerMl=mean(CellsPerMl))

# Plot the figure
pd<-position_dodge(0.85)
(PlotVirus<- ggplot()+ geom_point(data=VirusPlot, aes(x=Treatment, y=MeanCells, shape=Line, col=ResType), size = 3, position=pd) + geom_errorbar(data = VirusPlot, aes(x=Treatment, y = MeanCells, ymin = MeanCells - std.error, ymax = MeanCells + std.error, col=ResType, shape = Line, panel.background="white"), size=0.3, width = 0.3, position=pd) +facet_grid(Environment~ResType) + geom_hline(aes(yintercept=100000), data=VirusPlot, linetype="dashed", size=0.5, col="black") + scale_shape_manual(values=c(15,16,17,18,19,25,0,1,2)) + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + xlab("OtV5 treatment") + ylab("Cells per ml") + theme(axis.title=element_text(size=15), axis.text=element_text(size=12), strip.text.x=element_text(size=13), strip.text.y=element_text(size=13), strip.background=element_rect(fill="white", colour="black", size=0.35)) +theme(axis.title.x=element_text(margin=margin(12))))

# Model
#Calculate means for canto replicates
VirusData<-ddply(VirusRaw3d, .(Environment, Line, ResType, Treatment, Biorep, Techrep, PopID), summarize, CellsPerMl=mean(Cells))

# Look at the data
plot(VirusData$Treatment, VirusData$CellsPerMl)
plot(VirusData$ResType, VirusData$CellsPerMl)
plot(VirusData$Environment, VirusData$CellsPerMl)

## Check whether the data are normally distributed
hist(VirusData$CellsPerMl, 100) #Not normally distributed
VirusLog<-ddply(VirusData, .(Environment, Line, ResType, Treatment, Biorep, Techrep, PopID, CellsPerMl), summarize, log10=log(CellsPerMl, base=10))
hist(VirusLog$log10, 100) #Still not normally distributed 
qqnorm(VirusLog$log10)
qqnorm(VirusLog$log10);qqline(VirusLog$log10, col = 2)
shapiro.test(VirusLog$log10)

# (it does not make a difference if the random effect included techrep)

#Mixed effects model
virus<-lmer(CellsPerMl~Environment*ResType*Treatment + (1|PopID), data=VirusData)
# Check the residuals
hist(resid(virus), 100) #Residuals are normally distributed
anova(virus)
summary(virus)

# Post hoc tests to examine where differences occurred between groups
lsmeans(virus, pairwise~Environment, adjust="Tukey")
lsmeans(virus, pairwise~ResType, adjust="Tukey")
lsmeans(virus, pairwise~Treatment, adjust="Tukey")


## Differences
# Difference in cell density as measured by not-inoculated mean minus inoculated mean
# Figure
Differences3days <- read.csv("~/Desktop/Figures & Analyses/Differences3Days.csv")
options(scipen = 999)
(PlotDifferences<- ggplot()+ geom_point(data=Differences3days, aes(x=Line, y=Difference, group=Line), size = 2) +facet_grid(Environment~ResType, scales="free_x") + geom_hline(aes(yintercept=0), data=Differences3days, linetype="dashed", size=0.5, col="black") + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + xlab("Line") + ylab("Differences in cell density") + theme(axis.title=element_text(size=15), axis.text=element_text(size=12), strip.text.x=element_text(size=13), strip.text.y=element_text(size=13), strip.background=element_rect(fill="white", colour="black", size=0.35)) +theme(axis.title.x=element_text(margin=margin(12))))

# Model
library(lme4)
library(plyr)
library("lmerTest", lib.loc="~/Library/R/3.2/library")
hist(Differences3days$Difference, 100)
difference<-lmer(Difference~Environment + ResType + (1|Biorep) +(1|Line), data=Differences3days)
anova(difference)


## Change in susceptibility of S at start and end of experiment
SusceptibleChange <- read.csv("~/Desktop/Figures & Analyses/SusceptibleChange.csv")
View(SusceptibleChange)
ChangeMean<-ddply(SusceptibleChange, .(Environment, Strain, Test), summarise, DifferenceMean=mean(Difference), std.error=std.error(Difference))
(SusceptibleChange<- ggplot()+ geom_point(data=ChangeMean, aes(x=Test, y=DifferenceMean), size = 2) +facet_grid(Environment~Strain) + geom_hline(aes(yintercept=0), data=ChangeMean, linetype="dashed", size=0.5, col="black")  + geom_errorbar(data = ChangeMean, aes(x=Test, y = DifferenceMean, ymin = DifferenceMean - std.error, ymax = DifferenceMean + std.error, panel.background="white"), size=0.3, width=0.15) +scale_x_discrete("OtV5 test", labels=c("First"="Start", "Last"="End")) + theme_bw() + theme(panel.grid.major=element_blank(), panel.grid.minor=element_blank()) + ylab("Change in cell density") + theme(axis.title=element_text(size=15), axis.text=element_text(size=12), strip.text.x=element_text(size=13), strip.text.y=element_text(size=13), strip.background=element_rect(fill="white", colour="black", size=0.35)) +theme(axis.title.x=element_text(margin=margin(12))))

# Model
SusceptibleChange <- read.csv("~/Desktop/Figures & Analyses/SusceptibleChange.csv")
change<-lmer(Difference~Environment + Test + (1|Biorep), data=SusceptibleChange)
anova(change)