require(ggplot2) || install.packages(gglot2)
library(ggplot2)

#Working directory containing output from 2gate_classifer.pl script.
setwd("C:/analysis folder/analysis output/") #Please note that forward slashes are needed in the address
path <- "testfile.csv" #adjust to match chosen filename
data <- read.csv(file=path,header=TRUE)

#Thresholds for plot annotation (these are for graph annotation only, adjust to match thresholds used in Perl script data processing)
threshg1s <- c(1.5)
thresh780 <- c(0.004)

#Wells containing separate siRNA contextual data

NT <- data[data$Well %in% c("E5","F5","G5"),]
RB <- data[data$Well %in% c("E7","F7","G7"),]
cdk6 <- data[data$Well %in% c("B2","C2","D2"),]

#Commands to plot the data in Figures 6-8 are below. Commands for individual plots are separated empty lines.

#histogram plots to judge gating of S780 (Figure 6)

plotNT <- qplot(MeanRed, data = NT, geom = "histogram", fill = Label,  binwidth = 0.0003, xlab = "Nuclear P-S780 RB1 antibody intensity", xlim = c(0,0.02), main = "siRNA: Non-Targeting")
plotNT + scale_fill_manual(values=c("tomato","tomato","darkgrey","darkgrey"))+ guides(fill=FALSE)+ 
  geom_vline(aes(xintercept = thresh780),size = 1, data = NT)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

plotRB <- qplot(MeanRed, data = RB, geom = "histogram", fill = Label, binwidth = 0.0003, xlab = "Nuclear P-S780 RB1 antibody intensity", xlim = c(0,0.02), main = "siRNA: RB1")
plotRB + scale_fill_manual(values=c("tomato","tomato","darkgrey","darkgrey"))+ guides(fill=FALSE)+
  geom_vline(aes(xintercept = thresh780),size = 1, data = RB)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

plotCDK6 <- qplot(MeanRed, data = cdk6, geom = "histogram", fill = Label, binwidth = 0.0003, xlab = "Nuclear P-S780 RB1 antibody intensity", xlim = c(0,0.02), main = "siRNA: CDK6")
plotCDK6 + scale_fill_manual(values=c("tomato","tomato","darkgrey","darkgrey"))+ guides(fill=FALSE)+
  geom_vline(aes(xintercept = thresh780),size = 1, data = cdk6)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

#histogram plots to judge gating of GFP-CDK2 reporter (Figure 6)

plotNT <- qplot(GFP_Ratio, data = NT, geom = "histogram", fill = Label, binwidth = 0.125, xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", xlim = c(0,10), main = "siRNA: Non-Targeting")
plotNT + scale_fill_manual(values=c("darkgreen","darkgrey","darkgreen","darkgrey"))+ guides(fill=FALSE)+
  geom_vline(aes(xintercept = threshg1s),size = 1, data = NT)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

plotRB <- qplot(GFP_Ratio, data = RB, geom = "histogram", fill = Label, binwidth = 0.125, xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", xlim = c(0,10), main = "siRNA: RB1")
plotRB + scale_fill_manual(values=c("darkgreen","darkgrey","darkgreen","darkgrey"))+ guides(fill=FALSE)+
  geom_vline(aes(xintercept = threshg1s),size = 1, data = RB)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

plotcdk6 <- qplot(GFP_Ratio, data = cdk6, geom = "histogram", fill = Label, binwidth = 0.125, xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", xlim = c(0,10), main = "siRNA: CDK6")
plotcdk6 + scale_fill_manual(values=c("darkgreen","darkgrey","darkgreen","darkgrey"))+ guides(fill=FALSE)+
  geom_vline(aes(xintercept = threshg1s),size = 1, data = cdk6)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

#labelled xyplots used in Figures 7 and 8  (dependent on classification labels generated by 2gate_classifer.pl script)

NT_dim <- dim(NT)
ul <- 100*(sum(NT$Label == "P-S780+ nonG1")/NT_dim[1])
ll <- 100*(sum(NT$Label == "P-S780- nonG1")/NT_dim[1])
ur <- 100*(sum(NT$Label == "P-S780+ G1")/NT_dim[1])
lr <- 100*(sum(NT$Label == "P-S780- G1")/NT_dim[1])
NT_xy <- qplot(GFP_Ratio, MeanRed, data = NT, xlim = c(0,8), ylim = c(0,0.02), alpha = I(1/4), main = "siRNA: Non-targeting", xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", ylab = "Nuclear P-S780 RB1\nantibody intensity")
NT_xy + theme(panel.grid.minor = element_blank()) +
  geom_vline(aes(xintercept = threshg1s),size = 1.5, data = NT, colour = "darkgreen") +
  geom_hline(aes(yintercept = thresh780),size = 1.5, data = NT, colour = "tomato") +
  annotate("text", x = 0.25, y = 0.015,size = 8, label = signif(ul, digits = 3))+
  annotate("text", x = 0.25, y = 1e-04,size =8, label = signif(ll, digits = 3)) +
  annotate("text", x = 3, y = 0.015,size = 8, label = signif(ur, digits = 3)) +
  annotate("text", x = 3, y = 1e-04,size = 8, label = signif(lr, digits = 3))+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

cdk6_dim <- dim(cdk6)
ul <- 100*(sum(cdk6$Label == "P-S780+ nonG1")/cdk6_dim[1])
ll <- 100*(sum(cdk6$Label == "P-S780- nonG1")/cdk6_dim[1])
ur <- 100*(sum(cdk6$Label == "P-S780+ G1")/cdk6_dim[1])
lr <- 100*(sum(cdk6$Label == "P-S780- G1")/cdk6_dim[1])
cdk6_xy <- qplot(GFP_Ratio, MeanRed, data = cdk6, xlim = c(0,8), ylim = c(0,0.02), alpha = I(1/4), main = "siRNA: CDK6", xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", ylab = "Nuclear P-S780 RB1\nantibody intensity")
cdk6_xy + theme(panel.grid.minor = element_blank()) +
  geom_vline(aes(xintercept = threshg1s),size = 1.5, data = cdk6, colour = "darkgreen") +
  geom_hline(aes(yintercept = thresh780),size = 1.5, data = cdk6, colour = "tomato") +
  annotate("text", x = 0.25, y = 0.015,size = 8, label = signif(ul, digits = 4)) +
  annotate("text", x = 0.25, y = 1e-04,size =8, label = signif(ll, digits = 4)) +
  annotate("text", x = 3, y = 0.015,size = 8, label = signif(ur, digits = 4)) +
  annotate("text", x = 3, y = 1e-04,size = 8, label = signif(lr, digits = 4))+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

RB_dim <- dim(RB)
ul <- 100*(sum(RB$Label == "P-S780+ nonG1")/RB_dim[1])
ll <- 100*(sum(RB$Label == "P-S780- nonG1")/RB_dim[1])
ur <- 100*(sum(RB$Label == "P-S780+ G1")/RB_dim[1])
lr <- 100*(sum(RB$Label == "P-S780- G1")/RB_dim[1])
RB_xy <- qplot(GFP_Ratio, MeanRed, data = RB, xlim = c(0,8), ylim = c(0,0.02), alpha = I(1/4), main = "siRNA: RB1", xlab = "GFP-CDK2 Reporter (nuc/cyt) ratio", ylab = "Nuclear P-S780 RB1\nantibody intensity")
RB_xy + theme(panel.grid.minor = element_blank()) +
  geom_vline(aes(xintercept = threshg1s),size = 1.5, data = RB, colour = "darkgreen") +
  geom_hline(aes(yintercept = thresh780),size = 1.5, data = RB, colour = "tomato") +
  annotate("text", x = 0.25, y = 0.015,size = 8, label = signif(ul, digits = 3)) +
  annotate("text", x = 0.25, y = 1e-04,size =8, label = signif(ll, digits = 4)) +
  annotate("text", x = 3, y = 0.015,size = 8, label = signif(ur, digits = 3)) +
  annotate("text", x = 3, y = 1e-04,size = 8, label = signif(lr, digits = 4))+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))

#Histograms of integrated nuclear DNA intensity vs. single assay defined subclass

nonG1_cdk6 <- cdk6[cdk6$Label %in% c("P-S780+ nonG1","P-S780- nonG1"),]
plotall <- qplot(IntegBlue, data = nonG1_cdk6, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "non G1(siCDK6)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

G1_cdk6 <- cdk6[cdk6$Label %in% c("P-S780+ G1","P-S780- G1"),]
plotall <- qplot(IntegBlue, data = G1_cdk6, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "G1(siCDK6)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

nonP780_cdk6 <- cdk6[cdk6$Label %in% c("P-S780- G1","P-S780- nonG1"),]
plotall <- qplot(IntegBlue, data = nonP780_cdk6, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "non P-S780 RB1 (siCDK6)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

P780_cdk6 <- cdk6[cdk6$Label %in% c("P-S780+ G1","P-S780+ nonG1"),]
plotall <- qplot(IntegBlue, data = P780_cdk6, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "P-S780 RB1 (siCDK6)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

#Optional DNA profile analysis example for siNonTargeting cells (not used in manuscript but repeated here for comparison)

nonG1_NT <- NT[NT$Label %in% c("P-S780+ nonG1","P-S780- nonG1"),]
plotall <- qplot(IntegBlue, data = nonG1_NT, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main ="non G1(siNonTargeting)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

G1_NT <- NT[NT$Label %in% c("P-S780+ G1","P-S780- G1"),]
plotall <- qplot(IntegBlue, data = G1_NT, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "G1(siNonTargeting)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

nonP780_NT <- NT[NT$Label %in% c("P-S780- G1","P-S780- nonG1"),]
plotall <- qplot(IntegBlue, data = nonP780_NT, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "non P-S780 RB1 (siNT)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall

P780_NT <- NT[NT$Label %in% c("P-S780+ G1","P-S780+ nonG1"),]
plotall <- qplot(IntegBlue, data = P780_NT, geom = "histogram",fill = Label, binwidth = 0.125,xlim = c(0,15), xlab = "Integrated nuclear DNA intensity", main = "P-S780 RB1 (siNT)")+
  scale_fill_manual(values=c(rep("darkblue",2)))+ guides(fill=FALSE)+
  theme(axis.text.x = element_text(colour="grey20",size=20,angle=90,hjust=.5,vjust=.5,face="plain"),
        axis.text.y = element_text(colour="grey20",size=20,angle=0,hjust=1,vjust=0,face="plain"),  
        axis.title.x = element_text(size=20,angle=0,hjust=.5,vjust=0,face="bold"),
        axis.title.y = element_text(size=20,angle=90,hjust=.5,vjust=.5,face="bold"))+
  theme(plot.title = element_text(size = 20, lineheight=.8, face="bold"))+ 
  theme(strip.text.y = element_text(size = 20, angle = 90))
plotall