###########################################################################################################
##data = 30 dpi inj vs uninj
#calculate ddCt for each biological replicate 
#plot technical mean +/- stdev for each biological replicate
#compute statistics on sample means - welch test ??

#hld edited 1/1/18
###########################################################################################################

library(dplyr)
library(tidyr)
library(reshape2)
library(ggpubr)
library(scales)
library(FSA)
library(ggrepel)

path <- "/Users/hdingwall/Documents/Research/gradSchool/GallowayLab/qPCR_data/LC480/poolingData"
setwd(path)

qpcrInj.df <- read.table("injuryDataForMethodPaper.txt", header = TRUE)

#rm technical outliers
techVar.inj <-  qpcrInj.df %>%
  group_by(GeneName, BioRep, SampleGroup, Group) %>%
  mutate(techAve = mean(Cp), techSD = sd(Cp), techVar = abs(Cp - mean(Cp))) %>%
  # filter(techVar <= 1.5*techSD, techVar <= 0.5) %>%
  ungroup()
  
techAv.inj <-  techVar.inj %>%
  group_by(GeneName, BioRep, SampleGroup, Group) %>%
  summarize(Count = n(), Cp_techAv = mean(Cp), Cp_techVar = var(Cp), Cp_techSD = sd(Cp)) %>%
  ungroup()

################################################################
#check bio rep var
bioVar.inj <- techAv.inj %>%
  group_by(GeneName, Group) %>%
  mutate(Cp_bioAv = mean(Cp_techAv), Cp_bioSD = sd(Cp_techAv),
         bioVar = abs(Cp_techAv - mean(Cp_techAv)))
################################################################

#compute dCt for each bio rep

techAv.inj.gapdh <- techAv.inj %>%
  filter(GeneName == "Gapdh")

dCp.inj.gapdh <- techAv.inj %>%
  inner_join(techAv.inj.gapdh, by = "BioRep") %>%
  rename(GOI = GeneName.x, SampleGroup = SampleGroup.x, Group = Group.x, n.tech.goi = Count.x,
         Cp.techAv.goi = Cp_techAv.x, techSD.goi = Cp_techSD.x, techVar.goi = Cp_techVar.x,
         RefGene = GeneName.y,n.tech.ref = Count.y, Cp.techAv.ref = Cp_techAv.y,
         techSD.ref = Cp_techSD.y, techVar.ref = Cp_techVar.y) %>%
  mutate(dCp = (Cp.techAv.goi - Cp.techAv.ref), dCp.genes.sd = sqrt((techSD.ref^2) + (techSD.goi^2)),
         dCp.genes.var = sqrt((techVar.ref^2)+(techVar.goi^2))) %>%
  select(GOI:Group, RefGene, dCp, dCp.genes.sd, dCp.genes.var) %>%
  ungroup()

#compute ddCt for each bio rep, normalize injured to internal uninjured control
dCp.gapdh.un <- dCp.inj.gapdh %>%
  group_by(GOI) %>%
  filter(Group == "un")

ddCp.gapdh.un <- dCp.inj.gapdh %>%
  group_by(GOI, SampleGroup) %>%
  inner_join(dCp.gapdh.un, by = c("GOI", "SampleGroup")) %>%
  filter(GOI != "Gapdh") %>%
  rename(Group = Group.x, BioRep = BioRep.x) %>%
  select(GOI:dCp.genes.var.x, Group.y:dCp.genes.var.y) %>%
  mutate(ddCp = (dCp.x - dCp.y)) %>%
  rename(control.group = Group.y, ddCp.prop.sd = dCp.genes.sd.x, ddCp.prop.var = dCp.genes.var.x) %>%
  mutate(ddCp.calc.sd = sqrt(ddCp.prop.var)) %>%
  select(GOI:Group, ddCp, ddCp.prop.var, ddCp.calc.sd) %>%
  ungroup()

#compute fold chance for each bio rep

ddCp.gapdh.un.FC <- ddCp.gapdh.un %>%
  mutate(CI.low = (ddCp + ddCp.calc.sd),
         CI.high = (ddCp - ddCp.calc.sd),
         relExp.ddCpFC = 2^(-ddCp),
         relExp.ddCpFC.low = 2^(-CI.low), relExp.ddCpFC.high = 2^(-CI.high))

ddCp.gapdh.un.FC$Group <- factor(ddCp.gapdh.un.FC$Group, levels = rev(levels(ddCp.gapdh.un.FC$Group)))

plt.Sox9.FC <- filter(ddCp.gapdh.un.FC, GOI == 'Sox9') %>%
  filter(SampleGroup != 'B') %>%
  ggplot(aes(x = SampleGroup, y = relExp.ddCpFC, fill = Group, order = Group)) +
  geom_bar(stat = 'identity', position = position_dodge()) +
  geom_errorbar(aes(ymin = relExp.ddCpFC.low, ymax = relExp.ddCpFC.high), width = 0.2, 
                position = position_dodge(0.9)) +
  scale_fill_grey(name = NULL, breaks = c("un", "inj"), labels = c("uninjured", "injured")) +
  scale_x_discrete(breaks = c("A", "C", "D"), labels = c("A", "B", "C")) +
  ggtitle(expression(italic("Sox9"))) +
  xlab("Biological Replicate") +
  ylab(expression("Fold change " (2^{~-Delta~Delta ~"C"[T]}))) +
  expand_limits(y = c(0,3.5)) +
  theme_pubr() +
  theme(plot.title = element_text(size = 18, hjust = 0.5, margin = margin(b = 15, unit = "pt"))) +
  theme(axis.title.x = element_text(size = 14, margin = margin(t = 10, unit = "pt"))) +
  theme(axis.title.y = element_text(size = 14)) +
  theme(plot.margin = unit(c(3,3,3,3), "pt"))
  
plt.Col1.FC <- filter(ddCp.gapdh.un.FC, GOI == 'ColI') %>%
  ggplot(aes(x = SampleGroup, y = relExp.ddCpFC, fill = Group, order = Group)) +
  geom_bar(stat = 'identity', position = position_dodge()) +
  geom_errorbar(aes(ymin = relExp.ddCpFC.low, ymax = relExp.ddCpFC.high), width = 0.2,
                position = position_dodge(0.9)) +
  scale_fill_grey(name = NULL, breaks = c("un", "inj"), labels = c("uninjured", "injured")) +
  scale_x_discrete(breaks = c("A", "B", "C"), labels = c("A", "B", "C")) +
  ggtitle(expression(italic("Col1a2"))) +
  xlab("Biological Replicate") +
  ylab(expression("Fold change " (2^{~-Delta~Delta ~"C"[T]}))) +
  expand_limits(y = c(0,6.5)) +
  theme_pubr() +
  theme(plot.title = element_text(size = 18, hjust = 0.5, margin = margin(b = 15, unit = "pt"))) +
  theme(axis.title.x = element_text(size = 14, margin = margin(t = 10, unit = "pt"))) +
  theme(axis.title.y = element_text(size = 14)) +
  theme(plot.margin = unit(c(3,3,3,3), "pt"))

plt.injury.FC <- ggarrange(plt.Sox9.FC, plt.Col1.FC, labels = c("A", "B"), 
                       font.label = list(size = 20),
                       common.legend = TRUE,
                       legend = "bottom")
png(file = "30dpi_injury_fig6", width = 2640, height = 3300,  res = 300)
plot(plt.injury.FC)
dev.off()