### this is FRRF data, too, but it looks at the parameters extracted for Iopt
### we then look at plassticity in these data, i.e. we make the slopes first 


library(data.table)
library(plyr)
library(tools)
library(zoo)
library(lattice)
library(ggplot2)
require(nlme)
require(minpack.lm)
library(reshape2)
#install.packages("devtools")
library(devtools)
#install_github('padpadpadpad/TeamPhytoplankton')
#install_github('padpadpadpad/nlsTools')
library(TeamPhytoplankton)
library(nlsLoop)
library(cowplot)
library(xlsx)
require(reshape2)
#require(gdata)
setwd("/Users/elisaschaum/Dropbox/Thally long fluc paper/Thally actual selection experiment/fastact by time and temp/t300 full curves")



f8<-read.csv('allfvfmdata_thally_t300_t000.csv')
colnames(f8)
str(f8)
pfcmelt<-melt(f8, id.vars=c("biorep","PAR","assayT","growthT","generation"), measure.vars=c( "Fo.F.." ,"Fm...", "Fv...","Fv.Fm.Fq..Fm..", "p" , "RSigma", "Sigma","Fa", "Tau","NPQ","NSV" , "C" ,"X.RCII.", "JPII", "JVPII","aLHII","X0","E","rP"), variable.name='variable')
pfcmelt[is.na(pfcmelt)] <- 0

pfcmeltnew <- subset(pfcmelt,value!=0)
str(pfcmeltnew)
summary(pfcmeltnew)

pfc2<-ddply(pfcmeltnew, .(PAR,biorep,variable,assayT,growthT,generation), function(x) x[nrow(x), ])
pfc2$PAR<-as.numeric(pfc2$PAR)
pfc2$variable<-as.character(pfc2$variable)
pfc2<-subset(pfc2,PAR<1600) #two measurements go up to 1900- made all the same length for now

# ##### plasticity in dark adapted fv/fm #####
dark <- subset(pfc2, PAR == 0 & variable == 'Fv.Fm.Fq..Fm..')
dark<-dark[,c(1,3:5,7)]
colnames(dark)<-c("biorep","assayT","growthT","generation","Fvfm")

str(dark)


#keep 22,32,FS FL and the Ancestor
dark$growthT<-as.character(dark$growthT)
dark$assayT<-as.numeric(dark$assayT)
dark$generation<-as.character(dark$generation)

dark$growthT[dark$generation == 't000'] <- 'Ancestor'
dark$assayT[dark$assayT == 20] <- 22
dark$assayT[dark$assayT == 25] <- 26
dark<-subset(dark, assayT!=26 & biorep!="0")

#plot


dark$growthT<-as.factor(dark$growthT)
dark$growthT <- factor(dark$growthT, levels = c('Ancestor', '22', '32', 'FS','FL'))

dark<-within(dark, id.plot<-as.character(as.factor(growthT):as.factor(biorep)))

darkplot<-qplot(as.factor(assayT),as.numeric(Fvfm), ylab="Fvfm dark", colour=as.factor(growthT), shape=biorep,data=dark) +
  theme_classic()+ theme(legend.position='top')+ 
  facet_wrap(~growthT)+
  geom_smooth(aes(group=id.plot, colour=growthT), method="lm", se=F)+
  scale_colour_manual(values=c("darkgrey","darkgreen","red","purple","black"))
darkplot

dark$Fvfm<-as.numeric(dark$Fvfm)

#now make the slopes 
slopefn<-function(d) { 
  m<-lm (Fvfm~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopedata<-ddply(dark,.(dark$growthT,dark$biorep),slopefn)
#slopedata[1,3]
slopedata
colnames(slopedata)[1]<-"growthT"
colnames(slopedata)[2]<-"biorep"

slopedata$growthT<-as.character(slopedata$growthT)
slopedata$biorep<-as.character(slopedata$biorep)

slopedata$growthT2<-ifelse(slopedata$growthT=="FL"& slopedata$biorep=="b1"|
                             slopedata$growthT=="FL"&slopedata$biorep=="b2"|
                             slopedata$growthT=="FL"&slopedata$biorep=="b5",
                           "FL2",
                           ifelse(slopedata$growthT=="FL"&slopedata$biorep=="b3"|
                                    slopedata$growthT=="FL"&slopedata$biorep=="b4"|
                                    slopedata$growthT=="FL"&slopedata$biorep=="b6",
                                  "FL1",as.character(slopedata$growthT)))

slopedata$growthT2<-as.factor(slopedata$growthT2)
slopedata$growthT2 <- factor(slopedata$growthT2, levels = c('Ancestor', '22', '32', 'FS','FL1','FL2'))

FVFMplasplot01<-qplot(as.factor(growthT2), slope, data=slopedata, geom="boxplot", fill=growthT2, ylab="Evolved plasticity (Slopes of reaction norms for FvFM)",xlab="Selection environment")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT2))+
  scale_fill_manual(values=c("darkgrey","darkgreen","darkred","purple","coral","royalblue"))+
  scale_colour_manual(values=c("lightgrey","lightgreen","red","magenta","gold","cornflowerblue"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
FVFMplasplot01


Mod0 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, method = 'ML', slopedata)
S1 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'REML', slopedata)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT2, adjust = "tukey") 

Mod0 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', slopedata)
S1 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', slopedata)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT, adjust = "tukey") 


##### now the plasticity or rP at I opt ####

rP_extract<-read.csv("rP_extracted_values.csv")
head(rP_extract)

rP_extract<-rP_extract[,c(2,6:8)]

str(rP_extract)


#keep 22,32,FS FL and the Ancestor
rP_extract$growthT<-as.character(rP_extract$growthT)
rP_extract$assayT<-as.numeric(rP_extract$assayT)


#plot
rP_extract$growthT<-as.factor(rP_extract$growthT)
rP_extract$growthT <- factor(rP_extract$growthT, levels = c('Ancestor', '22', '32', 'FS','FL'))

rP_extract<-within(rP_extract, id.plot<-as.character(as.factor(growthT):as.factor(biorep)))

rP_extract_plot<-qplot(as.factor(assayT),as.numeric(Pmax), ylab="rP", colour=as.factor(growthT), shape=biorep,data=rP_extract) +
  theme_classic()+ theme(legend.position='top')+ 
  facet_wrap(~growthT)+
  geom_smooth(aes(group=id.plot, colour=growthT), method="lm", se=F)+
  scale_colour_manual(values=c("darkgrey","darkgreen","red","purple","black"))
rP_extract_plot


#now make the slopes 
slopefn1<-function(d) { 
  m<-lm (Pmax~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopedata1<-ddply(rP_extract,.(rP_extract$growthT,rP_extract$biorep),slopefn1)
#slopedata[1,3]
slopedata1
colnames(slopedata1)[1]<-"growthT"
colnames(slopedata1)[2]<-"biorep"

slopedata1$growthT<-as.character(slopedata1$growthT)
slopedata1$biorep<-as.character(slopedata1$biorep)

slopedata1$growthT2<-ifelse(slopedata1$growthT=="FL"& slopedata1$biorep=="b1"|
                              slopedata1$growthT=="FL"&slopedata1$biorep=="b2"|
                              slopedata1$growthT=="FL"&slopedata1$biorep=="b5",
                           "FL2",
                           ifelse(slopedata1$growthT=="FL"&slopedata1$biorep=="b3"|
                                    slopedata1$growthT=="FL"&slopedata1$biorep=="b4"|
                                    slopedata1$growthT=="FL"&slopedata1$biorep=="b6",
                                  "FL1",as.character(slopedata1$growthT)))

slopedata1$growthT2<-as.factor(slopedata1$growthT2)
slopedata1$growthT2 <- factor(slopedata1$growthT2, levels = c('Ancestor', '22', '32', 'FS','FL1','FL2'))

slopedata1plas<-qplot(as.factor(growthT2), slope, data=slopedata1, geom="boxplot", fill=growthT2, ylab="Evolved plasticity (Slopes of reaction norms for FvFM)",xlab="Selection environment")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT2))+
  scale_fill_manual(values=c("darkgrey","darkgreen","darkred","purple","coral","royalblue"))+
  scale_colour_manual(values=c("lightgrey","lightgreen","red","magenta","gold","cornflowerblue"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
slopedata1plas

Mod0 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, method = 'ML', slopedata1)
S1 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata1)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata1)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'REML', slopedata1)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT2, adjust = "tukey") 

Mod0 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', slopedata1)
S1 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata1)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata1)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', slopedata1)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT, adjust = "tukey") 



###psi psiII plasticity (at Iopt) which is between 1200 and 1500 for all  ######

pfcnew<-subset(pfc2,PAR>1200 & PAR<1500)
#now make means 
pfcnew$value<-as.numeric(pfcnew$value)
pfcnewav<-ddply(pfcnew,c("growthT","assayT","generation","biorep","variable") , function(df) return(c(var.avg=mean(df$value),var.se=sd(df$value))))

Fv.FM<- subset(pfcnewav,var.avg>0&variable=="Fv.Fm.Fq..Fm.."& var.avg!="NA" &var.avg!="")
Fv.FM$growthT<-as.character(Fv.FM$growthT)
Fv.FM$generation<-as.character(Fv.FM$generation)
Fv.FM$growthT[Fv.FM$generation == "t000"] <- "Ancestor"
Fv.FM$assayT<-as.numeric(Fv.FM$assayT)
Fv.FM<-subset(Fv.FM, assayT!=26) 
Fv.FM<-subset(Fv.FM, growthT!="26") 


#plot
Fv.FM$growthT<-as.factor(Fv.FM$growthT)
Fv.FM$growthT <- factor(Fv.FM$growthT, levels = c('Ancestor', '22', '32', 'FS','FL'))

Fv.FM<-within(Fv.FM, id.plot<-as.character(as.factor(growthT):as.factor(biorep)))

Fv.FM_plot<-qplot(as.factor(assayT),as.numeric(var.avg), ylab=expression(psi~PSII), colour=as.factor(growthT), shape=biorep,data=Fv.FM) +
  theme_classic()+ theme(legend.position='top')+ 
  facet_wrap(~growthT)+
  geom_smooth(aes(group=id.plot, colour=growthT), method="lm", se=F)+
  scale_colour_manual(values=c("darkgrey","darkgreen","red","purple","black"))
Fv.FM_plot


#now make the slopes 
slopefn2<-function(d) { 
  m<-lm (var.avg~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopedata2<-ddply(Fv.FM,.(Fv.FM$growthT,Fv.FM$biorep),slopefn2)
#slopedata[1,3]
slopedata2
colnames(slopedata2)[1]<-"growthT"
colnames(slopedata2)[2]<-"biorep"

slopedata2$growthT<-as.character(slopedata2$growthT)
slopedata2$biorep<-as.character(slopedata2$biorep)

slopedata2$growthT2<-ifelse(slopedata2$growthT=="FL"& slopedata2$biorep=="b1"|
                              slopedata2$growthT=="FL"&slopedata2$biorep=="b2"|
                              slopedata2$growthT=="FL"&slopedata2$biorep=="b5",
                            "FL2",
                            ifelse(slopedata2$growthT=="FL"&slopedata2$biorep=="b3"|
                                     slopedata2$growthT=="FL"&slopedata2$biorep=="b4"|
                                     slopedata2$growthT=="FL"&slopedata2$biorep=="b6",
                                   "FL1",as.character(slopedata2$growthT)))

slopedata2$growthT2<-as.factor(slopedata2$growthT2)
slopedata2$growthT2 <- factor(slopedata2$growthT2, levels = c('Ancestor', '22', '32', 'FS','FL1','FL2'))

slopedata2plas<-qplot(as.factor(growthT2), slope, data=slopedata2, geom="boxplot", fill=growthT2, ylab="Evolved plasticity (Slopes of reaction norms for PsiPSII)",xlab="Selection environment")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT2))+
  scale_fill_manual(values=c("darkgrey","darkgreen","darkred","purple","coral","royalblue"))+
  scale_colour_manual(values=c("lightgrey","lightgreen","red","magenta","gold","cornflowerblue"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
slopedata2plas

slopedata <- slopedata2$slope[!is.na(slopedata2$slope)]
str(slopedata2)

slopedata2 <- na.omit(slopedata2)

Mod0 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, method = 'ML', slopedata2)
S1 <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata2)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'ML', slopedata2)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT2, random = ~ 1|biorep/growthT2, weights = varIdent( ~ 1|growthT2), method = 'REML', slopedata2)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT2, adjust = "tukey") 

Mod0 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, method = 'ML', slopedata2)
S1 <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata2)
S2 <- lme(fixed = slope ~ 1 , random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'ML', slopedata2)
anova(S1,S2) # growthT is siginficant effect
Sfinal <- lme(fixed = slope ~ 1 +  growthT, random = ~ 1|biorep/growthT, weights = varIdent( ~ 1|growthT), method = 'REML', slopedata2)
confint.mm <- data.frame(intervals(Sfinal, which='fixed')[1])
confint.mm
lsmeans(Sfinal, pairwise ~ growthT, adjust = "tukey") 

#####last one is for the thingies in closed state plasticity #####


C<- subset(pfc2,value>0&variable=="C"& value!="NA" &value!="")

C$growthT<-as.character(C$growthT)
C$generation<-as.character(C$generation)
C$growthT[C$generation == "t000"] <- "Ancestor"


C$assayT<-as.numeric(C$assayT)
C<-subset(C, assayT!=26) 
C<-subset(C, growthT!="26") 


#plot
C$growthT<-as.factor(C$growthT)
C$growthT <- factor(C$growthT, levels = c('Ancestor', '22', '32', 'FS','FL'))

C$value<-as.numeric(C$value)
C<-subset(C, value<=1)

Cavg<-ddply(C,c("growthT","assayT","generation","biorep","variable") , function(df) return(c(var.avg=mean(df$value),var.se=sd(df$value))))
Cavg<-within(Cavg, id.plot<-as.character(as.factor(growthT):as.factor(biorep)))


CPLOT<-qplot(as.factor(assayT),as.numeric(var.avg), ylab=expression(C), colour=as.factor(growthT), shape=biorep,data=Cavg) +
  theme_classic()+ theme(legend.position='top')+ 
  facet_wrap(~growthT)+
  geom_smooth(aes(group=id.plot, colour=growthT), method="lm", se=F)+
  scale_colour_manual(values=c("darkgrey","darkgreen","red","purple","black"))
CPLOT


#now make the slopes 
slopefn2<-function(d) { 
  m<-lm (var.avg~assayT,data=d) #can force through origin by setting offset to 0 
  sum.<-summary(m) 
  r2<-sum.$r.squared
  intercept<-m$coefficients[1]
  slope<-m$coefficients[2]
  output<-data.frame(slope,intercept,r2)
  output}
slopedata3<-ddply(Cavg,.(Cavg$growthT,Cavg$biorep),slopefn2)
#slopedata[1,3]
slopedata3
colnames(slopedata3)[1]<-"growthT"
colnames(slopedata3)[2]<-"biorep"

slopedata3$growthT<-as.character(slopedata3$growthT)
slopedata3$biorep<-as.character(slopedata3$biorep)

slopedata3$growthT2<-ifelse(slopedata3$growthT=="FL"& slopedata3$biorep=="b1"|
                              slopedata3$growthT=="FL"&slopedata3$biorep=="b2"|
                              slopedata3$growthT=="FL"&slopedata3$biorep=="b5",
                            "FL2",
                            ifelse(slopedata3$growthT=="FL"&slopedata3$biorep=="b3"|
                                     slopedata3$growthT=="FL"&slopedata3$biorep=="b4"|
                                     slopedata3$growthT=="FL"&slopedata3$biorep=="b6",
                                   "FL1",as.character(slopedata3$growthT)))

slopedata3$growthT2<-as.factor(slopedata3$growthT2)
slopedata3$growthT2 <- factor(slopedata3$growthT2, levels = c('Ancestor', '22', '32', 'FS','FL1','FL2'))

slopedata3plot<-qplot(as.factor(growthT2), slope, data=slopedata3, geom="boxplot", fill=growthT2, ylab="Evolved plasticity (Slopes of reaction norms for centres C in closed state)",xlab="Selection environment")+
  geom_point(size=2.5, position="jitter",aes(shape=biorep,colour=growthT2))+
  scale_fill_manual(values=c("darkgrey","darkgreen","darkred","purple","coral","royalblue"))+
  scale_colour_manual(values=c("lightgrey","lightgreen","red","magenta","gold","cornflowerblue"))+
  theme_classic(base_size = 14, base_family = 'Helvetica')
slopedata3plot
