###########################################################################################################
##calculate raw Ct technical variance and biological variance per gene for each pooling level
##calculate normalized Ct (to ref gene, aka dCt) technical variance and biological variance
#for each pooling level
##plot variance comparisons

#hld edited 12/26/17
###########################################################################################################

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)

qpcrPool.df <- read.table("2017-11-15_poolingQPCR_pooled.txt", header = TRUE)

#rm technical outliers
techVar.Pool <-  qpcrPool.df %>%
  group_by(GeneName, BioRep, Group) %>%
  mutate(techAve = mean(Cp), techSD = sd(Cp), techVar = abs(Cp - mean(Cp))) %>%
  filter(techVar <= 1.5*techSD, techVar <= 0.5) %>%
  ungroup()

#compute summary stats on tech reps
Ct.var <-  techVar.Pool %>%
  group_by(GeneName, BioRep, Group) %>%
  summarize(Count = n(), Cp_techAv = mean(Cp), Cp_techVar = var(Cp)) %>%
  group_by(GeneName, Group) %>%
  mutate(tech_avgVar = mean(Cp_techVar), bioRepAv = mean(Cp_techAv), bioRepVar = var(Cp_techAv),
         totalVar = (bioRepVar + tech_avgVar), percent_bioVar = 100*(bioRepVar/totalVar), 
         percent_techVar = 100*(tech_avgVar/totalVar)) %>%
  ungroup()

Ct.var.long <- Ct.var %>%
  group_by(GeneName, Group) %>%
  select(GeneName, Group, tech_avgVar, bioRepVar) %>%
  melt()

plt.Ct.var.Scx <- filter(Ct.var.long, GeneName == 'Scx') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Scx)),
       x = "Pooling Level", y = expression("C"[T] ~Variance ~(sigma^2))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12)) +
  scale_fill_manual(values = c("black", "grey"), name = "Source of Variance",
                    breaks = c("tech_avgVar", "bioRepVar"),
                      labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), name = "Source of Variance",
                      breaks = c("tech_avgVar", "bioRepVar"),
                      labels = c("technical", "biological")) +
  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.Ct.var.gapdh <- filter(Ct.var.long, GeneName == 'Gapdh') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Gapdh)),
       x = "Pooling Level", y = expression("C"[T] ~Variance ~(sigma^2))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12)) +
  scale_fill_manual(values = c("black", "grey"), name = "Source of Variance",
                    breaks = c("tech_avgVar", "bioRepVar"),
                    labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), name = "Source of Variance",
                      breaks = c("tech_avgVar", "bioRepVar"),
                      labels = c("technical", "biological")) +
  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.Ct.var <- ggarrange(plt.Ct.var.Scx, plt.Ct.var.gapdh, labels = c("A", "B"), 
                        font.label = list(size = 20),
                        common.legend = TRUE,
                        legend = "bottom")

png(file = "CtVariance_pooled_Scx-Gapdh.png", width = 2640, height = 3300,  res = 300)
plot(plt.Ct.var)
dev.off()

#same with stdev
Ct.SD.long <- Ct.var %>%
  group_by(GeneName, Group) %>%
  mutate(tech_avSD = sqrt(tech_avgVar), bioRepSD = sqrt(bioRepVar)) %>%
  select(GeneName, Group, tech_avSD, bioRepSD) %>%
  melt()

plt.Ct.SD.Scx <- filter(Ct.SD.long, GeneName == 'Scx') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Scx)),
       x = "Pooling Level", 
       y = expression("C"[T] ~"Standard Deviation" ~(sigma))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12)) +
  scale_fill_manual(values = c("black", "grey"), name = "Source of Variation",
                    breaks = c("tech_avSD", "bioRepSD"),
                    labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), name = "Source of Variation",
                      breaks = c("tech_avSD", "bioRepSD"),
                      labels = c("technical", "biological")) +
  expand_limits(y = c(0,11.5)) +
  scale_y_continuous(breaks = c(2,4,6,8,10,12)) +
  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.Ct.SD.gapdh <- filter(Ct.SD.long, GeneName == 'Gapdh') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Gapdh)),
       x = "Pooling Level", y = expression("C"[T] ~"Standard Deviation" ~(sigma))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12)) +
  scale_fill_manual(values = c("black", "grey"), name = "Source of Variation",
                    breaks = c("tech_avSD", "bioRepSD"),
                    labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), name = "Source of Variation",
                      breaks = c("tech_avSD", "bioRepSD"),
                      labels = c("technical", "biological")) +
  scale_y_continuous(breaks = c(2,4,6,8,10,12)) +
  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.Ct.SD <- ggarrange(plt.Ct.SD.Scx, plt.Ct.SD.gapdh, labels = c("A", "B"), 
                        font.label = list(size = 20),
                        common.legend = TRUE,
                        legend = "bottom")

png(file = "CtStdev_pooled_Scx-Gapdh.png", width = 2640, height = 3300,  res = 300)
plot(plt.Ct.SD)
dev.off()

#######################################################################################
#dCt
#######################################################################################
techAv.pool <-  techVar.Pool %>%
  group_by(GeneName, BioRep, Group, SampleType) %>%
  summarize(Count = n(), Cp_techAv = mean(Cp), Cp_techSD = sd(Cp), Cp_techVar = var(Cp))

#compute dCp
techAv.pool.gapdh <- Ct.var %>%
  filter(GeneName == "Gapdh")

dCp.gapdh.Pool <- Ct.var %>%
  inner_join(techAv.pool.gapdh, by = "BioRep") %>%
  rename(GeneName = GeneName.x, Group = Group.x, n.tech.goi = Count.x, Cp.techAv.goi = Cp_techAv.x,
         Cp_techVar.goi = Cp_techVar.x,
         RefGene = GeneName.y, n.tech.ref = Count.y, Cp.techAv.ref = Cp_techAv.y,
         Cp_techVar.ref = Cp_techVar.y) %>%
  filter(GeneName != "PPIA") %>%
  select(-bioRepAv.x, -bioRepVar.x, -bioRepAv.y, -bioRepVar.y, -percent_bioVar.x, -percent_bioVar.y,
         -totalVar.x, -totalVar.y, -percent_techVar.x, -percent_techVar.y) %>%
  mutate(dCp = (Cp.techAv.goi - Cp.techAv.ref)) %>%
  group_by(GeneName, Group) %>%
  summarize(av_dCp = mean(dCp), bioRepVar.goi = var(Cp.techAv.goi), bioRepVar.ref = var(Cp.techAv.ref),
            tech_avgVar.goi = mean(Cp_techVar.goi), tech_avgVar.ref = mean(Cp_techVar.ref),
            dCp.techVar.prop = sqrt((tech_avgVar.goi^2) + (tech_avgVar.ref^2)),
            var_dCp = var(dCp), dCp.bioVar.prop = sqrt((bioRepVar.ref^2) + (bioRepVar.goi^2))) %>%
  ungroup()

dCp.gapdh.long <- dCp.gapdh.Pool %>%
  group_by(GeneName, Group) %>%
  select(GeneName, Group, dCp.techVar.prop, dCp.bioVar.prop) %>%
  melt()

plt.dCt.var.Scx <- filter(dCp.gapdh.long, GeneName == 'Scx') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Scx)),
       x = "Pooling Level", y = expression(Delta ~"C"[T] ~Variance ~(sigma^2))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12), legend.position = c(0.1,0.9)) +
  scale_fill_manual(values = c("black", "grey"), guide = "legend", name = "Source of Variance",
                    breaks = c("dCp.techVar.prop", "dCp.bioVar.prop"),
                    labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), guide = "legend", name = "Source of Variance",
                      breaks = c("dCp.techVar.prop", "dCp.bioVar.prop"),
                      labels = c("technical", "biological")) +
  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"))

png(file = "deltaCtVariance_pooled_Scx-Gapdh.png", width = 3300, height = 2640,  res = 300)
plot(plt.dCt.var.Scx)
dev.off()

#same with stdev
dCp.gapdh.long.SD <- dCp.gapdh.Pool %>%
  group_by(GeneName, Group) %>%
  mutate(dCp.techSD.prop = sqrt(dCp.techVar.prop), dCp.bioSD.prop = sqrt(dCp.bioVar.prop)) %>%
  select(GeneName, Group, dCp.techSD.prop, dCp.bioSD.prop) %>%
  melt()

plt.dCt.SD.Scx <- filter(dCp.gapdh.long.SD, GeneName == 'Scx') %>%
  ggbarplot(x = 'Group', y = 'value',
            fill = 'variable', color = 'variable',
            label = FALSE) +
  labs(title = expression(italic(Scx)),
       x = "Pooling Level", y = expression(Delta ~"C"[T] ~"Standard Deviation" ~(sigma))) +
  theme(legend.background = element_rect(fill = NULL, color = NULL),
        legend.text = element_text(size = 12), legend.position = c(0.1,0.9)) +
  scale_fill_manual(values = c("black", "grey"), guide = "legend", name = "Source of Variation",
                    breaks = c("dCp.techSD.prop", "dCp.bioSD.prop"),
                    labels = c("technical", "biological")) +
  scale_colour_manual(values = c("black", "grey"), guide = "legend", name = "Source of Variation",
                      breaks = c("dCp.techSD.prop", "dCp.bioSD.prop"),
                      labels = c("technical", "biological")) +
  expand_limits(y = c(0,5)) +
  scale_y_continuous(breaks = c(1,2,3,4,5)) +
  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"))

png(file = "deltaCt_Stdev_pooled_Scx-Gapdh.png", width = 3300, height = 2640,  res = 300)
plot(plt.dCt.SD.Scx)
dev.off()