Divergent responses of soil microbial community after amendment with thermally altered Pinus radiata needles

RScript written by Erinne Stirling 2018

R version 3.5.1 (2018-07-02) -- "Feather Spray"
Copyright (C) 2018 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin15.6.0 (64-bit)

library(readr)#for reading data
library(dplyr)#for manipulating data
library(ggplot2)#for plotting
library(RColorBrewer)#for palettes
library(agricolae)#for stastics
library(reshape2)#for melting
library(stringr)#for text extraction
library(nlme)

RESP_CUM<-read_csv("CSV/Cumulative_respiration.csv")

#----------
#Decay Curves
SCR_blank<-RESP_CUM%>%
  filter(TYPE=="BBB.C")

plot(SCR_blank$TIME,SCR_blank$CUM_RESP)

##----
rep_1<-SCR_blank%>%
  filter(REP=="1")

fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_1,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
rep_1$pred <- predict(fit_a)
plot(rep_1$TIME,rep_1$CUM_RESP)
lines(rep_1$TIME,rep_1$pred, col="blue",lty=2)

##----
rep_2<-SCR_blank%>%
  filter(REP=="2")

fit_b<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_2,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_b)
rep_2$pred <- predict(fit_b)
plot(rep_2$TIME,rep_2$CUM_RESP)
lines(rep_2$TIME,rep_2$pred, col="blue",lty=2)

##----
rep_3<-SCR_blank%>%
  filter(REP=="3")

fit_c<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_3,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_c)
rep_3$pred <- predict(fit_c)
plot(rep_3$TIME,rep_3$CUM_RESP)
lines(rep_3$TIME,rep_3$pred, col="blue",lty=2)

##----
rep_4<-SCR_blank%>%
  filter(REP=="4")

fit_d<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_4,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
d<-summary(fit_d)
rep_4$pred <- predict(fit_d)
plot(rep_4$TIME,rep_4$CUM_RESP)
lines(rep_4$TIME,rep_4$pred, col="blue",lty=2)

##----
rep_5<-SCR_blank%>%
  filter(REP=="5")

fit_e<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_5,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
e<-summary(fit_e)
rep_5$pred <- predict(fit_e)
plot(rep_5$TIME,rep_5$CUM_RESP)
lines(rep_5$TIME,rep_5$pred, col="blue",lty=2)

SCR_blank_ave<-data.frame(id=c(rep("Control",5)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1],d$parameters[1,1],e$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1],d$parameters[2,1],e$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1],d$parameters[3,1],e$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1],d$parameters[4,1],e$parameters[4,1]))

#----------
#Decay Curves
SCR_40C<-RESP_CUM%>%
  filter(TYPE=="040.C")

plot(SCR_40C$TIME,SCR_40C$CUM_RESP)

##----
rep_1<-SCR_40C%>%
  filter(REP=="1")

fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_1,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
rep_1$pred <- predict(fit_a)
plot(rep_1$TIME,rep_1$CUM_RESP)
lines(rep_1$TIME,rep_1$pred, col="blue",lty=2)

##----
rep_2<-SCR_40C%>%
  filter(REP=="2")

fit_b<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_2,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_b)
rep_2$pred <- predict(fit_b)
plot(rep_2$TIME,rep_2$CUM_RESP)
lines(rep_2$TIME,rep_2$pred, col="blue",lty=2)

##----
rep_3<-SCR_40C%>%
  filter(REP=="3")

fit_c<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_3,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_c)
rep_3$pred <- predict(fit_c)
plot(rep_3$TIME,rep_3$CUM_RESP)
lines(rep_3$TIME,rep_3$pred, col="blue",lty=2)

##----
rep_4<-SCR_40C%>%
  filter(REP=="4")

fit_d<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_4,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
d<-summary(fit_d)
rep_4$pred <- predict(fit_d)
plot(rep_4$TIME,rep_4$CUM_RESP)
lines(rep_4$TIME,rep_4$pred, col="blue",lty=2)

##----
rep_5<-SCR_40C%>%
  filter(REP=="5")

fit_e<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_5,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
e<-summary(fit_e)
rep_5$pred <- predict(fit_e)
plot(rep_5$TIME,rep_5$CUM_RESP)
lines(rep_5$TIME,rep_5$pred, col="blue",lty=2)

SCR_40C_ave<-data.frame(id=c(rep("40C",5)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1],d$parameters[1,1],e$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1],d$parameters[2,1],e$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1],d$parameters[3,1],e$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1],d$parameters[4,1],e$parameters[4,1]))

#----------
#Decay Curves
SCR_150C<-RESP_CUM%>%
  filter(TYPE=="150.C")

plot(SCR_150C$TIME,SCR_150C$CUM_RESP)

##----
rep_1<-SCR_150C%>%
  filter(REP=="1")

fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_1,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
rep_1$pred <- predict(fit_a)
plot(rep_1$TIME,rep_1$CUM_RESP)
lines(rep_1$TIME,rep_1$pred, col="blue",lty=2)

##----
rep_2<-SCR_150C%>%
  filter(REP=="2")

fit_b<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_2,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_b)
rep_2$pred <- predict(fit_b)
plot(rep_2$TIME,rep_2$CUM_RESP)
lines(rep_2$TIME,rep_2$pred, col="blue",lty=2)

##----
rep_3<-SCR_150C%>%
  filter(REP=="3")

fit_c<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_3,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_c)
rep_3$pred <- predict(fit_c)
plot(rep_3$TIME,rep_3$CUM_RESP)
lines(rep_3$TIME,rep_3$pred, col="blue",lty=2)

##----
rep_4<-SCR_150C%>%
  filter(REP=="4")

fit_d<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_4,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
d<-summary(fit_d)
rep_4$pred <- predict(fit_d)
plot(rep_4$TIME,rep_4$CUM_RESP)
lines(rep_4$TIME,rep_4$pred, col="blue",lty=2)

##----
rep_5<-SCR_150C%>%
  filter(REP=="5")

fit_e<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_5,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
e<-summary(fit_e)
rep_5$pred <- predict(fit_e)
plot(rep_5$TIME,rep_5$CUM_RESP)
lines(rep_5$TIME,rep_5$pred, col="blue",lty=2)

SCR_150C_ave<-data.frame(id=c(rep("150C",5)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1],d$parameters[1,1],e$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1],d$parameters[2,1],e$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1],d$parameters[3,1],e$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1],d$parameters[4,1],e$parameters[4,1]))

#----------
#Decay Curves
SCR_260C<-RESP_CUM%>%
  filter(TYPE=="260.C")

plot(SCR_260C$TIME,SCR_260C$CUM_RESP)

##----
rep_1<-SCR_260C%>%
  filter(REP=="1")

fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_1,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
rep_1$pred <- predict(fit_a)
plot(rep_1$TIME,rep_1$CUM_RESP)
lines(rep_1$TIME,rep_1$pred, col="blue",lty=2)

##----
rep_2<-SCR_260C%>%
  filter(REP=="2")

fit_b<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_2,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_b)
rep_2$pred <- predict(fit_b)
plot(rep_2$TIME,rep_2$CUM_RESP)
lines(rep_2$TIME,rep_2$pred, col="blue",lty=2)

##----
rep_3<-SCR_260C%>%
  filter(REP=="3")

fit_c<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_3,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_c)
rep_3$pred <- predict(fit_c)
plot(rep_3$TIME,rep_3$CUM_RESP)
lines(rep_3$TIME,rep_3$pred, col="blue",lty=2)

##----
rep_4<-SCR_260C%>%
  filter(REP=="4")

fit_d<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_4,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
d<-summary(fit_d)
rep_4$pred <- predict(fit_d)
plot(rep_4$TIME,rep_4$CUM_RESP)
lines(rep_4$TIME,rep_4$pred, col="blue",lty=2)

##----
rep_5<-SCR_260C%>%
  filter(REP=="5")

fit_e<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_5,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
e<-summary(fit_e)
rep_5$pred <- predict(fit_e)
plot(rep_5$TIME,rep_5$CUM_RESP)
lines(rep_5$TIME,rep_5$pred, col="blue",lty=2)

SCR_260C_ave<-data.frame(id=c(rep("260C",5)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1],d$parameters[1,1],e$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1],d$parameters[2,1],e$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1],d$parameters[3,1],e$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1],d$parameters[4,1],e$parameters[4,1]))

#----------
#Decay Curves
SCR_320C<-RESP_CUM%>%
  filter(TYPE=="320.C")

plot(SCR_320C$TIME,SCR_320C$CUM_RESP)

##----
rep_1<-SCR_320C%>%
  filter(REP=="1")

fit_a<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_1,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
a<-summary(fit_a)
rep_1$pred <- predict(fit_a)
plot(rep_1$TIME,rep_1$CUM_RESP)
lines(rep_1$TIME,rep_1$pred, col="blue",lty=2)

##----
rep_2<-SCR_320C%>%
  filter(REP=="2")

fit_b<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_2,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
b<-summary(fit_b)
rep_2$pred <- predict(fit_b)
plot(rep_2$TIME,rep_2$CUM_RESP)
lines(rep_2$TIME,rep_2$pred, col="blue",lty=2)

##----
rep_3<-SCR_320C%>%
  filter(REP=="3")

fit_c<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_3,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
c<-summary(fit_c)
rep_3$pred <- predict(fit_c)
plot(rep_3$TIME,rep_3$CUM_RESP)
lines(rep_3$TIME,rep_3$pred, col="blue",lty=2)

##----
rep_4<-SCR_320C%>%
  filter(REP=="4")

fit_d<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_4,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
d<-summary(fit_d)
rep_4$pred <- predict(fit_d)
plot(rep_4$TIME,rep_4$CUM_RESP)
lines(rep_4$TIME,rep_4$pred, col="blue",lty=2)

##----
rep_5<-SCR_320C%>%
  filter(REP=="5")

fit_e<-nls(CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME,data=rep_5,
           algorithm="port",
           start=c(C=10,k=1,a=1,m=0),lower=c(C=0,k=0,a=0,m=0), upper=c(C=Inf,k=Inf,a=Inf,m=Inf))
e<-summary(fit_e)
rep_5$pred <- predict(fit_e)
plot(rep_5$TIME,rep_5$CUM_RESP)
lines(rep_5$TIME,rep_5$pred, col="blue",lty=2)

SCR_320C_ave<-data.frame(id=c(rep("320C",5)),
                          C=c(a$parameters[1,1],b$parameters[1,1],c$parameters[1,1],d$parameters[1,1],e$parameters[1,1]),
                          k=c(a$parameters[2,1],b$parameters[2,1],c$parameters[2,1],d$parameters[2,1],e$parameters[2,1]),
                          a=c(a$parameters[3,1],b$parameters[3,1],c$parameters[3,1],d$parameters[3,1],e$parameters[3,1]),
                          m=c(a$parameters[4,1],b$parameters[4,1],c$parameters[4,1],d$parameters[4,1],e$parameters[4,1]))

#-------------
Curves<-rbind(SCR_blank_ave,
              SCR_40C_ave,SCR_150C_ave,SCR_260C_ave,SCR_320C_ave)

#----------STATS
stats<-Curves%>%
  group_by(id)

anova(aov(C~id, data=stats)) #***K
HSD.test(aov(C~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

anova(aov(k~id, data=stats)) #***r
HSD.test(aov(k~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

anova(aov(a~id, data=stats)) #***M0
HSD.test(aov(a~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

anova(aov(m~id, data=stats)) #***b
HSD.test(aov(m~id, data=stats), c("id"),
         console=TRUE,alpha=0.05)

####------------Plot

p_curves<-Curves%>%
  group_by(id)%>% 
  summarize_at(vars(C,k,a,m), funs(mean))

write.csv(p_curves,"CSV/p_curves.csv")

p_curves<-read.csv("CSV/p_curves.csv")%>%
  select(-X)

CUM_RESP~(a*C)/(a+(C-a)*exp(-k*TIME))+m*TIME


fun_control <- function(x) {
  (p_curves[1,4]*p_curves[1,2])/(p_curves[1,4]+(p_curves[1,2]-p_curves[1,4])*exp(-p_curves[1,3]*x))+p_curves[1,5]*x
}

fun_40C <- function(x) {
  (p_curves[2,4]*p_curves[2,2])/(p_curves[2,4]+(p_curves[2,2]-p_curves[2,4])*exp(-p_curves[2,3]*x))+p_curves[2,5]*x
}

fun_150C <- function(x) {
  (p_curves[3,4]*p_curves[3,2])/(p_curves[3,4]+(p_curves[3,2]-p_curves[3,4])*exp(-p_curves[3,3]*x))+p_curves[3,5]*x
}

fun_260C <- function(x) {
  (p_curves[4,4]*p_curves[4,2])/(p_curves[4,4]+(p_curves[4,2]-p_curves[4,4])*exp(-p_curves[4,3]*x))+p_curves[4,5]*x
}

fun_320C <- function(x) {
  (p_curves[5,4]*p_curves[5,2])/(p_curves[5,4]+(p_curves[5,2]-p_curves[5,4])*exp(-p_curves[5,3]*x))+p_curves[5,5]*x
}



Divstat_lines<-ggplot(data.frame(x=c(0,45)),aes(x=x))+
  stat_function(fun=fun_control, colour="grey")+
  stat_function(fun=fun_40C, colour="lightblue")+
  stat_function(fun=fun_150C, colour="orange")+
  stat_function(fun=fun_260C, colour="yellow")+
  stat_function(fun=fun_320C, colour="black")+
  scale_y_continuous(limits=c(0,420))+
  ggsave("Divstat_resp2a.svg", width=200, height=100, units="mm")

RESP_CUM%>%
  ggplot(aes(x=TIME,y=CUM_RESP))+
  geom_point(aes(colour=TYPE))+
  #scale_colour_manual(values=c("grey", "blue", "lightblue", "red", "orange"))+
  xlab("Time (days)") +
  ylab(expression(paste("Cumulative Respiration (",mu,"gCO2-C/gOC)"))) +
  scale_y_continuous(limits=c(0,420))+
  ggsave("Divstat_resp2b.svg", width=200, height=100, units="mm")


#-------
#testing

SCR_blank$pred<-fun_control(SCR_blank$TIME)
BLANK.lm = lm(CUM_RESP~pred, data=SCR_blank)
summary(BLANK.lm)$r.squared

SCR_40C$pred<-fun_40C(SCR_40C$TIME)
T40C.lm = lm(CUM_RESP~pred, data=SCR_40C)
summary(T40C.lm)$r.squared

SCR_150C$pred<-fun_150C(SCR_150C$TIME)
T150C.lm = lm(CUM_RESP~pred, data=SCR_150C)
summary(T150C.lm)$r.squared

SCR_260C$pred<-fun_260C(SCR_260C$TIME)
T260C.lm = lm(CUM_RESP~pred, data=SCR_260C)
summary(T260C.lm)$r.squared

SCR_320C$pred<-fun_320C(SCR_320C$TIME)
T320C.lm = lm(CUM_RESP~pred, data=SCR_320C)
summary(T320C.lm)$r.squared