#loading data into r
#install.packages("gdata",dependencies = TRUE)
library(gdata)
frogdat<-read.xls("Sustainable_Tissues_RawData.xlsx",sheet=1,method="csv",header=TRUE)


#Whole data analysis
boxplot(frogdat$Mass_mg, xlab="Mass in mg")

boxplot(frogdat$DNA_ng.uL, xlab="DNA concentration")

databreaks <- c(.5:20.5)
hist(frogdat$Mass_mg, 
     xlab="Mass in mg",
     breaks=databreaks,
     col="steelblue1",
     freq=T,
     main="Samples per mass (mg)"
     )

hist(frogdat$DNA_ng.uL, xlab="DNA concentration",breaks=50)



#Exp1 correlation

hist(frogdat$Mass_mg[frogdat$Exp==1],xlab="Mass",breaks=10)
hist(frogdat$DNA_ng.uL[frogdat$Exp==1],xlab="DNA concentration",breaks=10)

cormethod <- c("pearson", "kendall", "spearman")
frog1cor <- matrix(nrow=3, ncol=3)

for(i in seq_along(cormethod)){
  frog1cor[1,i]<- cormethod[i]
  frog1cor[2,i] <- cor(frogdat$Mass_mg[frogdat$Exp==1],frogdat$DNA_ng.uL[frogdat$Exp==1], method=cormethod[i])
  frog1cor[3,i] <- cor.test(frogdat$Mass_mg[frogdat$Exp==1],frogdat$DNA_ng.uL[frogdat$Exp==1], method=cormethod[i])$p.value

}
plot(x=frogdat$Mass_mg[frogdat$Exp==1], y=frogdat$DNA_ng.uL[frogdat$Exp==1], xlab="Mass (mg)", ylab="[DNA] (ug)")

cor.test(frogdat$Mass_mg[frogdat$Exp==1],frogdat$DNA_ng.uL[frogdat$Exp==1])

#Exp2 scatterplot
plot(frogdat$Mass_mg[frogdat$Exp==2],frogdat$DNA_ng.uL[frogdat$Exp==2],
     pch=19,
     col="steelblue1",
     xlab="Mass of Tissue in mg",
     ylab="[DNA] in ng/uL",
     main="Effect of Tissue Mass on DNA Concentration",
     xaxt="n"
     )

axis(1,at=c(1,2,4,8,10,12,14,16,20),labels=T)
axis(2,at=seq(from=0,to=250,by=25))


#Exp2 scatterplot w/ model
  #install.packages("drc")
  #install.packages("ggplot2")
  library(drc)
  library(ggplot2)
frogmodel <- drm(frogdat$DNA_ng.uL[frogdat$Exp==2]~frogdat$Mass_mg[frogdat$Exp==2], fct=MM.2())
plot.new()
plot(frogmodel,
     pch=19,
     col="steelblue1",
     xlab="Mass of Tissue in mg",
     ylab="[DNA] in ng/uL",
     main="Effect of Tissue Mass on DNA Concentration",
     lab=c(1,6,7)
)
axis(1,at=c(0,2,4,6,8,12,14,16,20),labels=T)

frogcurve <- function(x)36.523*log(x)+48.021

plot(x=frogdat$Mass_mg[frogdat$Exp==2],y=frogdat$DNA_ng.uL[frogdat$Exp==2], xlab="Mass in mg", ylab="[DNA] in ng/uL", main="DNA Concentrations in Experiment 2")
plot.function(frogcurve, from=0, to=20, col="red",add=T)





#Exp2 sep boxplots
colorful <- colorRampPalette(colors=c("steelblue1","darkolivegreen2"))
boxcol <- c(colorful(9))

frogbox<- boxplot(frogdat$DNA_ng.uL[frogdat$Class_No=="1"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="2"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="4"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="8"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="10"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="12"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="14"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="16"&frogdat$Exp==2],
                  frogdat$DNA_ng.uL[frogdat$Class_No=="20"&frogdat$Exp==2],
                  xlab="Mass in mg",
                  ylab="[DNA] in ng/uL",
                  col=boxcol,
                  main="Effect of Tissue Mass on DNA Concentration"
                  )

axis(1,at=c(0:9),labels=c(0,1,2,4,8,10,12,14,16,20))
axis(2,at=seq(from=0,to=400,by=25))

boxlab <- vector()
for(i in 1:9){
  x <- frogbox$stats[3,i]
  boxlab[i] <- x
}
text(y=(frogbox$stats[5,]-6),x=seq(from=.85,to=8.85,by=1),labels=boxlab)

##
colorful <- colorRampPalette(colors=c("red","blue"))
boxcol <- c(colorful(9))

frogbox<- boxplot(frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="1"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="2"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="4"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="8"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="10"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="12"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="14"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="16"&frogdat$Exp==2],
                  frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="20"&frogdat$Exp==2],
                  xlab="Mass in mg",
                  ylab="DNA.Mass_ng.mg",
                  col=boxcol,
                  main="Effect of Tissue Mass on DNA Concentration per mg"
)

axis(1,at=c(0:9),labels=c(0,1,2,4,8,10,12,14,16,20))
axis(2,at=seq(from=0,to=400,by=25))

boxlab <- vector()
for(i in 1:9){
  x <- frogbox$stats[3,i]
  boxlab[i] <- x
}
text(y=(frogbox$stats[5,]-6),x=seq(from=.85,to=8.85,by=1),labels=boxlab)


box8 <- boxplot( frogdat$DNA.Mass_ng.mg[frogdat$Class_No=="8"&frogdat$Exp==2])

#Exp2 ANOVA

exp2anova <- aov(formula=frogdat$DNA_ng.uL[frogdat$Exp=="2"]~frogdat$Class_No[frogdat$Exp=="2"], data=frogdat[frogdat$Exp=="2",])
summary(exp2anova)
exp2tukey <- TukeyHSD(exp2anova)

colnames(exp2tukey$`frogdat$Class_No[frogdat$Exp == "2"]`)

exp2tpval <- data.frame(exp2tukey$`frogdat$Class_No[frogdat$Exp == "2"]`[,"p adj"], row.names = rownames(exp2tukey$`frogdat$Class_No[frogdat$Exp == "2"]`))

m1 <- do.call(rbind, strsplit(rownames(exp2tukey$`frogdat$Class_No[frogdat$Exp == "2"]`),"-"))
exp2tpval['first'] <- m1[,1]
exp2tpval['second'] <- m1[,2]


massnames <- c(1,2,4,8,10,12,14,16,20)

masslist <- list(massnames,massnames)

exp2vars <- matrix(nrow=9,ncol=9, dimnames = masslist)


exp2sdtest <- bartlett.test(formula=frogdat$DNA_ng.uL[frogdat$Exp=="2"]~frogdat$Class_No[frogdat$Exp=="2"], data=frogdat[frogdat$Exp=="2",])

for(i in seq_along(massnames)){
  exp2vars[i,1] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="1"])$p.value
  exp2vars[i,2] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="2"])$p.value
  exp2vars[i,3] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="4"])$p.value
  exp2vars[i,4] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="8"])$p.value
  exp2vars[i,5] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="10"])$p.value
  exp2vars[i,6] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="12"])$p.value
  exp2vars[i,7] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="14"])$p.value
  exp2vars[i,8] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="16"])$p.value
  exp2vars[i,9] <- var.test(frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No==massnames[i]],frogdat$DNA_ng.uL[frogdat$Exp=="2"&frogdat$Class_No=="20"])$p.value
}





#Exp4 correlation

hist(frogdat$Mass_mg[frogdat$Exp==4],xlab="Mass",breaks=10)
hist(frogdat$DNA_ng.uL[frogdat$Exp==4],xlab="DNA concentration",breaks=10)

cormethod <- c("pearson", "kendall", "spearman")
frog4cor <- matrix(nrow=2, ncol=3)

for(i in seq_along(cormethod)){
  frog4cor[1,i]<- cormethod[i]
  frog4cor[2,i] <- cor(frogdat$Mass_mg[frogdat$Exp==4],frogdat$DNA_ng.uL[frogdat$Exp==4], method=cormethod[i])
}

plot(x=frogdat$Year[frogdat$Exp==4], y=frogdat$DNA_ng.uL[frogdat$Exp==4], xlab="Collection Date", ylab= "[DNA] in ng/uL", main="DNA Concentration of Archival Specimens")

cor.test(frogdat$Mass_mg[frogdat$Exp==4],frogdat$DNA_ng.uL[frogdat$Exp==4],method="pearson")


#Analyzing data
summary(frogdat)

  #install.packages("psych")
  library(psych)

describe(frogdat)

prelim<-read.csv("PrelimCor.csv",header=TRUE)
cor(prelim$Mass..mg.,prelim$X.DNA.)

cor(frogdat$Mass_mg,frogdat$DNA_ng.uL)

frogaov <- aov(DNA_ng.uL~as.character(Class_No), frogdat)
summary(frogaov)

cor(frogdat$Year[frogdat$Year!=2016],frogdat$DNA_ng.uL[frogdat$Year!=2016])

classes <- c(1,2,4,8,10,12,14,16,20)
ranking <- matrix(nrow=9,ncol=3,dimnames=list(classes,list("Mean [DNA] ng/uL", "Standard Deviation","Ranking Number")))
for(i in seq_along(classes)){
  ranking[i,1] <-mean(frogdat$DNA_ng.uL[frogdat$Exp==2&frogdat$Class_No==classes[i]])
  ranking[i,2] <- sd(frogdat$DNA_ng.uL[frogdat$Exp==2&frogdat$Class_No==classes[i]])
  ranking[i,3] <- ranking[i,1]-ranking[i,2]
}

exp3rank <- mean(frogdat$DNA_ng.uL[frogdat$Exp==3])-sd(frogdat$DNA_ng.uL[frogdat$Exp==3])


#summary table
sumnames <- list("Target Mass mg","Replications","Mean [DNA] ng/uL","St.Dev [DNA]","Standard Error")

frogsum <- matrix(nrow=9, ncol=5,dimnames=list(list(),sumnames))

for(i in seq_along(classes)){
  frogsum[i,1] <- classes[i]
  frogsum[i,2] <- sum(frogdat$Class_No==classes[i],na.rm=T)
  frogsum[i,3] <- mean(frogdat$DNA_ng.uL[frogdat$Class_No==classes[i]],na.rm=T)
  frogsum[i,4] <- sd(frogdat$DNA_ng.uL[frogdat$Class_No==classes[i]],na.rm=T)
  frogsum[i,5] <- frogsum[i,4]/sqrt(frogsum[i,2])
}





#Exp 3 ANOVA

exp3anova <- aov(formula=frogdat$DNA_ng.uL[frogdat$Exp=="3"]~frogdat$Field_No[frogdat$Exp=="3"], data=frogdat[frogdat$Exp=="3",])
summary(exp3anova)
TukeyHSD(exp3anova)

exp3tukey <- TukeyHSD(exp3anova)

colnames(exp3tukey$`frogdat$Field_No[frogdat$Exp == "3"]`)

exp3tpval <- data.frame(exp3tukey$`frogdat$Field_No[frogdat$Exp == "3"]`[,"p adj"], row.names = rownames(exp3tukey$`frogdat$Field_No[frogdat$Exp == "3"]`))

m1 <- do.call(rbind, strsplit(rownames(exp3tukey$`frogdat$Field_No[frogdat$Exp == "3"]`),"-"))
exp3tpval['first'] <- m1[,1]
exp3tpval['second'] <- m1[,2]

plot(TukeyHSD(exp3anova))



###########################################################################################################################
###########################################################################################################################
###########################################################################################################################

#Final figures

#Figure 1:Tissue mass in Experiment 1
plot(x=frogdat$Mass_mg[frogdat$Exp==1], 
     y=frogdat$DNA_ng.uL[frogdat$Exp==1]*100, 
     xlab="Mass (mg)", 
     ylab="Total DNA yield (ng)",
     main="Total DNA yield vs tissue mass in Experiment 1",
     pch=16,
     cex=1.5,
     cex.axis=1.5,
     cex.lab=2,
     cex.main=2.5,
     mgp=c(2,0.5,0),
     tcl=.5)

abline(lm(frogdat$DNA_ng.uL[frogdat$Exp==1]*100~frogdat$Mass_mg[frogdat$Exp==1]),col="red",lwd=2)

#Figure 2:Non-linear relationship between tissue mass and concentration, best fit 
         #Non-linear relatioship between tissue mass and total DNA yield
plot(x=frogdat$Mass_mg[frogdat$Exp==2], 
     y=frogdat$DNA_ng.uL[frogdat$Exp==2]*100, 
     xlab="Mass (mg)", 
     ylab="DNA yield (ng)",
     main="DNA yield vs tissue mass in Experiment 2",
     pch=16,
     cex=1.5,
     cex.axis=1.5,
     cex.lab=2,
     cex.main=2.5,
     mgp=c(2,0.5,0),
     tcl=.5,
     lab=c(1,6,7))

frogcurve <- function(x)3317.2*log(x)+5030.3
plot.function(frogcurve, from=0, to=20, col="red",add=T,lwd=2)

frogcurve2 <- function(x)-1656*log(x)+5436.3
plot.function(frogcurve2, from=0, to=20, col="blue",add=T,lwd=2)

axis(1,at=c(0,1,2,4,8,10,12,14,16),labels=T, cex.axis=1.5, tcl=.5, mgp=c(2,0.5,0))

legend(x=1,
       y=25000,
       legend=c("Total DNA (ng)","DNA per unit mass (ng/mg)"),
       col=c("red","blue"), 
       lty=1:1,
       lwd=2,
       box.lty=0,
       text.width=.5)


#Figure 3:Tissue age vs DNA yield
plot(x=frogdat$Year[frogdat$Exp==4], 
     y=frogdat$DNA_ng.uL[frogdat$Exp==4]*100, 
     xlab="Year collected", 
     ylab="DNA yield (ng)",
     main="DNA yield vs tissue age in Experiment 4",
     pch=16,
     cex=1.5,
     cex.axis=1.5,
     cex.lab=2,
     cex.main=2.5,
     mgp=c(2,0.5,0),
     tcl=.5,
     lab=c(1,6,7))

axis(1,at=c(1984,1986,1988,1990,1992,1994,1996,1998,2000,2002),labels=T, cex.axis=1.5, tcl=.5, mgp=c(2,0.5,0))

abline(lm(frogdat$DNA_ng.uL[frogdat$Exp==4]*100~frogdat$Year[frogdat$Exp==4]),col="red",lwd=2)