library(multcomp)
library(ggplot2)
library(plyr)


setwd("C:/Users/awvater/Desktop")
qPCR_Tissue_161_343_Data <- read.csv("C:/Users/awvater/Desktop/Figure_4data.csv")
View(qPCR_Tissue_161_343_Data)

Mean_qPCR <- aggregate( LogDNA.ng~C_T_S + Day, qPCR_Tissue_161_343_Data, mean, na.rm=TRUE )
Mean_qPCR_DF <- data.frame(Mean_qPCR)
Mean_qPCR_DF$norm_counts = 10^Mean_qPCR_DF$LogDNA.ng
Mean_qPCR_DF$norm_counts[Mean_qPCR_DF$norm_counts == 1] <- 0
View(Mean_qPCR_DF)

#############161


Day161qPCR <- subset(qPCR_Tissue_161_343_Data, Day == 161)
LMqPCR <- lm(LogDNA.ng~ C_T_S, data = Day161qPCR)

hist(Day161qPCR$LogDNA.ng)

qqnorm(resid(LMqPCR))
shapiro.test (resid(LMqPCR))

res = resid(LMqPCR)

# get 1% and 99% for residual
q <- quantile(res, probs=c(.15, .85), na.rm = TRUE)

# perform windsorization 
for (i in 1:length(res)) {
  if (res[i] < q[1]) {
    Day161qPCR$LogDNA.ng[i] <- Day161qPCR$LogDNA.ng[i] - res[i] + q[1]
  } else if (res[i] > q[2]) {
    Day161qPCR$LogDNA.ng[i] <- Day161qPCR$LogDNA.ng[i] - res[i] + q[2]
  }
}

LMqPCR <- lm(LogDNA.ng~ C_T_S, data = Day161qPCR)
qqnorm(resid(LMqPCR))
shapiro.test (resid(LMqPCR))


#------------------343-----------------
Day343qPCR <- subset(qPCR_Tissue_161_343_Data, Day == 343)
LMqPCR <- lm(LogDNA.ng~ C_T_S, data = Day343qPCR)

hist(Day343qPCR$LogDNA.ng)

qqnorm(resid(LMqPCR))
shapiro.test (resid(LMqPCR))

ComparisonsqPCR<- glht (LMqPCR, mcp(C_T_S = "Tukey"))
summary(ComparisonsqPCR)


#---Done

# ------- construct data frame of the conditioned data -----
data <- rbind(Day161qPCR, Day343qPCR)

datasum <- ddply(data, c("Day", "C_T_S"), summarise,
                 N = length(LogDNA.ng), 
                 mean=mean(LogDNA.ng),
                 sd=sd(LogDNA.ng),
                 se=sd/sqrt(N) )

pd <- position_dodge(2) # move them .05 to the left and right

ggplot(datasum, aes(x=Day, y=mean, fill=C_T_S, group=C_T_S)) + 
  geom_bar(position=position_dodge(), stat="identity",
           colour="black",
           size=0.2) +
  geom_errorbar(aes(ymin=mean-se, ymax=mean+se), 
                colour="black", 
                width=20,
                position = position_dodge(175)) +
  xlab("Days") + 
  ylab("Log ( Ca. XC gene copy number / ng input DNA) )") +
  scale_fill_hue(name=" Experimental Treatment Groups",
                   labels = c("Ca. Xc - exposed, Ambient, Red",
                              "Ca. Xc - exposed, Ambient, White",
                              "Ca. Xc - exposed, Heated, Red",
                              "Ca. Xc - exposed, Heated, White",
                              "Naive, Heated, Red*",
                              "Naive, Heated, White*"))