# qPCR analysis
# copy number of (some) FACS samples versus fluorescence

# 20180802 qPCR#4 (H5, with CFP and rbsB primers, 5 orders of M dilutions)
# biorad qpcr machine BFB_1st floor

#H5 standart in parallel 
# $Content = Std (Standart), NTC (samples)
# $Samples: # 1-6; 8-9 == D8; #7 == H5 ctr.; H5==Standart  --> 3 replicates

#use Pfaffl method, efficiency slightly below 90% 

##  CFP intensity of macroscope picture (for_qPCR-636247363150899708-CFP) -> 4 replicates


######## read in data 
setwd("/Users/itomanek/Documents/promoter_evolution/experiments/qPCR/qPCR_20180802")
results_table= read.csv2("BollbackGroup_2018-08-02 16-49-05_qPCR -  Quantification Cq Results_0.csv", header=TRUE, sep= ",",stringsAsFactors=FALSE)
results_table$Cq=as.numeric(results_table$Cq)

#change name: Starting.Quantity..SQ.  to nanogram (ng)
names(results_table)[11] <-c("ng")
results_table$ng <- as.numeric(results_table$ng) 
#delete empty well with no content (also - ctr; to remove NA)
results_table=results_table[-78,]
results_table=results_table[-68,]
results_table[56,]$ng = 0.001  #manually correct conc. of std. (see labbook)

library("ggplot2")

#plot data # overview (x-axis not too meaningful)
ggplot(data=results_table)+
  geom_point(aes(x=(as.numeric(row.names(results_table))), y=Cq, col=Sample, shape=Target),size=2)

#plot data #with ng information
ggplot(data=(results_table))+      
  geom_point(aes(x=log10(ng), y=Cq, col=Target, shape=Target),size=2)+ #the "outlier" of E-2 CFP STD is just mislabelled (see well content; compare labbook)
  geom_smooth(data=subset(results_table, Target == "CFP"& Sample == "H5"),aes(x=log10(ng), y=Cq),color = "dark green", method = "lm",size=0.4, se = FALSE)+
  geom_smooth(data=subset(results_table, Target == "rbsB"& Sample == "H5"),aes(x=log10(ng), y=Cq),color = "dark blue", method = "lm",size=0.4, se = FALSE)+
  ggtitle("Standart Curve") +
  theme_bw()+  #no grey background
  theme(panel.grid.minor = element_blank(),panel.grid.major=element_blank(),text=element_text(size=20))+ 
  labs(colour = "Locus") #change legend title
  #both primer efficiency (slopes look very similar)
   # position="jitter"


##### primer efficiency :  ###################################
## amplified molecules should double every cycle:
##10-fold ng should give 3.3 times smaller Cq value (=100% efficiency)

#E =-1+10^(-1/slope).

#EFFICIENCY H5 - CFP primers
H5_CFP=subset(results_table, Sample == "H5" & Target == "CFP")
plot((H5_CFP$ng),H5_CFP$Cq) #nonlinear relationship
plot(log10(H5_CFP$ng),H5_CFP$Cq) 
H5_CFP_lm=lm((H5_CFP$Cq)~log10(H5_CFP$ng))
abline(H5_CFP_lm)
slope=H5_CFP_lm$coefficients[2]
E =-1+10^(-1/slope)  #perfect would be 2 (double the template after 1 cycle)
E*100 #Efficiency; 88.5%
E =10^(-1/slope)  #efficiency for Pfaffl method!
E

#EFFICIENCY H5 - rbsB primers
H5_rbsB=subset(results_table, Sample == "H5" & Target == "rbsB")
plot((H5_rbsB$ng),H5_rbsB$Cq) #nonlinear relationship
plot(log10(H5_rbsB$ng),H5_rbsB$Cq) 
H5_rbsB_lm=lm((H5_rbsB$Cq)~log10(H5_rbsB$ng))
abline(H5_rbsB_lm)
slope=H5_rbsB_lm$coefficients[2]
E =-1+10^(-1/slope)  #perfect would be 2 (double the template after 1 cycle)
E*100 #Efficiency; 89.5% 
E =10^(-1/slope)  #efficiency for Pfaffl method!
E =10^(-1/slope)  #efficiency for Pfaffl method!
E
## > efficiency of primers varies by 1%. Almost 90%. Use Pfaffl method to correct for efficiency <90%.
## this reproduces the result of qPCR_20180801_analysis.R !

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

### Pfaffl 2001 method (take into account efficiency)
E =10^(-1/slope)  #efficiency for Pfaffl method!
E_target=1.885 #CFP efficiency
E_ref=1.895 #rbsB efficiency

samplelist=c(1:9)
ratios=c()

for (i in 1:length(samplelist)) { 
CFP_cq_D8 = mean(subset(results_table, Sample == paste(i) & Target == "CFP")$Cq)
CFP_cq_H5 = mean(subset(results_table, Sample =="7" & Target == "CFP")$Cq)  #this is H5 = ctr.
rbsB_cq_H5 = mean(subset(results_table, Sample =="7" & Target == "rbsB")$Cq)
rbsB_cq_D8 = mean(subset(results_table, Sample == paste(i) & Target == "rbsB")$Cq)
ratio=E_target^(CFP_cq_H5 - CFP_cq_D8)/E_ref^(rbsB_cq_H5-rbsB_cq_D8)

ratios[i]=ratio
}
ratios #fold change relative to H5 ctr ("7")
barplot(ratios)


#### delta delta Ct method #######################  (just to check)

#d_cq_D8= CFP_cq - rbs_cq
#d_d_cq= d_cq_D8 - d_cq_H5
# relative normalized expression = 2^(-d_d_cq)
#D8
CFP_cq = mean(subset(results_table, Sample =="6" & Target == "CFP")$Cq)
rbs_cq=mean(subset(results_table, Sample == "6" & Target == "rbsB")$Cq)
d_cq_D8=CFP_cq - rbs_cq
#
#H5
CFP_cq = mean(subset(results_table, Sample =="7" & Target == "CFP")$Cq)
rbs_cq=mean(subset(results_table, Sample == "7" & Target == "rbsB")$Cq)
d_cq_H5=CFP_cq - rbs_cq
# # 
d_d_cq=d_cq_D8 - d_cq_H5
d_d_cq  
2^(-d_d_cq) # 6.695 fold change :-) 

#################################################
#what is the fluorescence-fold change ? RFU pin macroscope pix (fiji)
CFP=read.csv2("Results_CFP.csv",  header=TRUE, sep= ";" ,stringsAsFactors=FALSE)
YFP=read.csv2("Results_YFP.csv",  header=TRUE, sep= ";" ,stringsAsFactors=FALSE)

head(CFP)
CFP$Sample=as.character(CFP$Sample)
names(CFP)
names(CFP)= c("colony","Label","area","meanCFP","sdCFP","medianCFP","Sample")
names(YFP)= c("colony","Label_y","area_y","meanYFP","sdYFP","medianYFP")
unique(CFP$colony)

## pair CFP with YFP
library("dplyr")
results <- inner_join(CFP, YFP, by="colony", copy=TRUE)


## pair CFP with results_table

results_2 <- inner_join(results_table, results, by="Sample", copy=TRUE)
results<-results_2 #rename to work with below code

results$meanCFP=as.numeric(results$meanCFP)

CFPmeans=tapply(as.numeric(results$meanCFP),as.factor(results$Sample),mean)
YFPmeans=tapply(as.numeric(results$meanYFP),as.factor(results$Sample),mean)

CFPsds=tapply(as.numeric(results$meanCFP),as.factor(results$Sample),sd) ##  standard devs.for each Sample (four replicates)
YFPsds=tapply(as.numeric(results$meanYFP),as.factor(results$Sample),sd) ##  standard devs.for each Sample (four replicates)



plot(CFPmeans,ratios, ylab="relative copy #  -- as determined by qPCR", xlab="CFP fluorescene of colony [RFU]", main="CFP fluorescence as proxy for copy nr.")
abline(lm(ratios~CFPmeans))
summary((lm(ratios~CFPmeans)))

plot(CFPmeans,YFPmeans, ylab="relative copy #  -- as determined by qPCR", xlab="CFP fluorescene of colony [RFU]", main="CFP fluorescence as proxy for copy nr.")
abline(lm(ratios~YFPmeans))  #promoter mutant H5 is an outlier (as expected)

#same for errorbars! (calculate ratio for )
### repeat Pffaffl method for individual wells
samplelist=c(1:9)
ratios_all=matrix(0,9,3)

for (i in 1:length(samplelist)) { 
  CFP_cq_D8 = (subset(results_table, Sample == paste(i) & Target == "CFP")$Cq)
  CFP_cq_H5 = mean(subset(results_table, Sample =="7" & Target == "CFP")$Cq)  #this is H5 = ctr.
  rbsB_cq_H5 = mean(subset(results_table, Sample =="7" & Target == "rbsB")$Cq)
  rbsB_cq_D8 = (subset(results_table, Sample == paste(i) & Target == "rbsB")$Cq)
  ratio=E_target^(CFP_cq_H5 - CFP_cq_D8)/E_ref^(rbsB_cq_H5-rbsB_cq_D8)
  
  ratios_all[i,]=ratio
}
ratios_all #fold change relative to H5 ctr ("7")

library(matrixStats)
ratiosds=rowSds(ratios_all)

plotdf=data.frame(samplelist,ratios,CFPmeans,ratiosds,CFPsds,YFPmeans,YFPsds)

#plot in AMGET paper
ggplot(data=plotdf)+
  geom_point(aes(x=CFPmeans,y=ratios), size=2.5)+
  geom_smooth(data=plotdf,aes(x=CFPmeans, y=ratios), method = "lm",size=0.4, se = FALSE)+
  ggtitle("CFP as a proxy for copy number") +
  theme_bw()+ 
  theme(panel.grid.minor = element_blank(),panel.grid.major=element_blank(),text=element_text(size=20))+
  geom_errorbarh(aes(x=CFPmeans,y=ratios,xmin=CFPmeans-CFPsds,xmax=CFPmeans+CFPsds),size=0.5,width=0.5)+
  geom_errorbar(aes(x=CFPmeans,y=ratios,ymin=ratios-ratiosds,ymax=ratios+ratiosds),size=0.5,width=0.5)+
  labs(x="CFP fluorescence [a.u.]",y="relative copy number (from qPCR)") 

#hindrance paper
#CFP (switch axes to above plot)
a=ggplot(data=plotdf)+
   geom_point(aes(x=ratios,y=CFPmeans), size=2.5)+
  geom_point(data=subset(plotdf,samplelist=="7"),aes(x=ratios,y=CFPmeans), size=2.5,color="green")+
  geom_smooth(data=plotdf,aes(x=ratios, y=CFPmeans), method = "lm",size=1, se = FALSE,color="blue")+
  ggtitle("CFP as a proxy for copy number") +
  theme_bw()+ 
  theme(panel.grid.minor = element_blank(),panel.grid.major=element_blank(),text=element_text(size=20))+
  geom_errorbarh(aes(x=ratios,y=CFPmeans,xmin=ratios-ratiosds,xmax=ratios+ratiosds),size=0.5,width=0.5)+
  geom_errorbar(aes(x=ratios,y=CFPmeans,ymin=CFPmeans-CFPsds,ymax=CFPmeans+CFPsds),size=0.5,width=0.5)+
  labs(x="relative copy number (from qPCR)",y="CFP fluorescence [a.u.]") 


#YFP
#YFP (switch axes to above plot)
b=ggplot(data=plotdf)+
  geom_point(aes(x=ratios,y=YFPmeans), size=2.5)+
  geom_point(data=subset(plotdf,samplelist=="7"),aes(x=ratios,y=YFPmeans), size=2.5,color="green")+
  #exclude promoter mut.from stats, plot individually 
  geom_smooth(data=subset(plotdf,samplelist!="7"),aes(x=ratios, y=YFPmeans), method = "lm",size=1, se = FALSE,color="orange")+
  ggtitle("YFP as a proxy for copy number") +
  theme_bw()+ 
  theme(panel.grid.minor = element_blank(),panel.grid.major=element_blank(),text=element_text(size=20))+
  geom_errorbarh(aes(x=ratios,y=YFPmeans,xmin=ratios-ratiosds,xmax=ratios+ratiosds),size=0.5,width=0.5)+
  geom_errorbar(aes(x=ratios,y=YFPmeans,ymin=YFPmeans-YFPsds,ymax=YFPmeans+YFPsds),size=0.5,width=0.5)+
  labs(x="relative copy number (from qPCR)",y="YFP fluorescence [a.u.]") 

#YFP and CFP   
c=ggplot(data=plotdf)+
  geom_point(aes(x=CFPmeans,y=YFPmeans), size=2.5)+
  geom_point(data=subset(plotdf,samplelist=="7"),aes(x=CFPmeans,y=YFPmeans), size=2.5,color="green")+
  #exclude promoter mut.from stats, plot individually 
  geom_smooth(data=subset(plotdf,samplelist!="7"),aes(x=CFPmeans, y=YFPmeans), method = "lm",size=1, se = FALSE,color="dark green")+
  ggtitle("colony phenotype") +
  theme_bw()+ 
  theme(panel.grid.minor = element_blank(),panel.grid.major=element_blank(),text=element_text(size=20))+
  geom_errorbarh(aes(x=CFPmeans,y=YFPmeans,xmin=CFPmeans-CFPsds,xmax=CFPmeans+CFPsds),size=0.5,width=0.5)+
  geom_errorbar(aes(x=CFPmeans,y=YFPmeans,ymin=YFPmeans-YFPsds,ymax=YFPmeans+YFPsds),size=0.5,width=0.5)+
  labs(x="CFP fluorescence [a.u.]",y="YFP fluorescence [a.u.]") 
a
b
c
multiplot(a,b,c,cols=1)
library("Rmisc")
