Recombinastics_analysis.Rmd

---
title: "Recombinastics_analysis"
author: "Kenneth Matreyek and Nisha D. Kamath"
date: "initialized 6/17/2020 - last updated 9/22/2023"
output: github_document
---

```{r Set up the workspace}
rm(list = ls())
library(tidyverse)
library(ggrepel)
library(reshape)
library(abind)
library(gridExtra)
library(ggbeeswarm)
library(factoextra)
library(ggfortify)
library(patchwork)

theme_set(theme_bw())
theme_update(panel.grid.minor = element_blank())

R.Version()
```

## Look at the recombination method test data. Note: this is not part of the manuscript
```{r Testing recombination efficiency in two different plate formats, fig.height = 4, fig.width = 6}
methods_data <- read.csv(file = "data/Recombination_method_tests.csv", header = T)

methods_data_replicates <- ncol(methods_data) - 5
methods_data$mean <- rowMeans(methods_data[,c("Recombined_Rep1","Recombined_Rep2","Recombined_Rep3","Recombined_Rep4")], na.rm = T)
methods_data$sd <- sqrt((methods_data$Recombined_Rep1 - methods_data$mean)^2 + 
                        (methods_data$Recombined_Rep2 - methods_data$mean)^2 + 
                        (methods_data$Recombined_Rep3 - methods_data$mean)^2 +
                        (methods_data$Recombined_Rep4 - methods_data$mean)^2) 
methods_data$se <- methods_data$sd / sqrt(methods_data_replicates - 1)

methods_data$upper_conf <- methods_data$mean + methods_data$se * 1.96
methods_data$lower_conf <- methods_data$mean - methods_data$se * 1.96

methods_data2 <- methods_data %>% filter(!is.na(Recombined_Rep1) & Method != "none") %>% mutate(concat = paste("Method: ",Method,"\n","Bxb1: ",Bxb1,sep = ""))
methods_data2$WellSize <- factor(methods_data2$WellSize)
methods_data2$concat <- factor(methods_data2$concat, levels = methods_data2$concat[c(1,3,2,4)])

methods_data2[methods_data2$lower_conf < 0,"lower_conf"] <- 0

Recombination_method_plot <- ggplot() + 
  theme_classic() + theme(axis.text.x = element_text(angle = -90, vjust = 0.5)) +
  scale_y_continuous(expand = c(0,0), limits = c(0,65)) +
  ylab("% mCherry positive cells") +
  xlab(NULL) +
  geom_hline(yintercept = 0) + 
  geom_errorbar(data = methods_data2, aes(x = concat, ymin = lower_conf, ymax = upper_conf, color = WellSize), 
                width = 0.2, position=position_dodge(width=0.3), size = 0.3 ) +
  geom_point(data = methods_data2, aes(x= concat, y=mean, color = WellSize), 
             position=position_dodge(width=0.3), shape = 95, size = 2) +
  geom_jitter(data = methods_data2, aes(x= concat, y=Recombined_Rep1, color = WellSize), 
             size = 0.9, position=position_dodge(width=0.3), alpha = 0.4) +
  geom_jitter(data = methods_data2, aes(x= concat, y=Recombined_Rep2, color = WellSize), 
             size = 0.9, position=position_dodge(width=0.3), alpha = 0.4) +
  geom_jitter(data = methods_data2, aes(x= concat, y=Recombined_Rep3, color = WellSize), 
             size = 0.9, position=position_dodge(width=0.3), alpha = 0.4) +
  geom_jitter(data = methods_data2, aes(x= concat, y=Recombined_Rep4, color = WellSize), 
             size = 0.9, position=position_dodge(width=0.3), alpha = 0.4)
print(Recombination_method_plot)
```
## GT and GA Ortholognality
## Relevant to Figure 1
 This next chunk is for testing the orthogonality of the GT and GA Bxb1 sequences
```{r Initial testing of orthogonality of GT and GA recombination site pairs, fig.height = 2, fig.width = 3}
rep_1_frame <- read.csv(file = "data/Orthog_recomb/GT_vs_GA_200206.csv", header = T, stringsAsFactors = F) %>% arrange(attp, bxb1, mix) %>% filter(bxb1 != "none")
rep_2_frame <- read.csv(file = "data/Orthog_recomb/GT_vs_GA_200213.csv", header = T, stringsAsFactors = F) %>% arrange(attp, bxb1, mix) %>% filter(bxb1 != "none")
rep_3_frame <- read.csv(file = "data/Orthog_recomb/GT_vs_GA_200219.csv", header = T, stringsAsFactors = F) %>% arrange(attp, bxb1, mix) %>% filter(bxb1 != "none")

replicates <- 2
rep_1 <- log10(rep_1_frame[,5:8])
rep_2 <- log10(rep_2_frame[,5:8])
rep_3 <- log10(rep_3_frame[,5:8])
means <- (rep_1 + rep_2 + rep_3)/3

#https://stackoverflow.com/questions/32609926/performing-element-wise-standard-deviation-in-r-with-two-matrices
m <- abind(rep_1, rep_2, rep_3, along=3)
standard_devs <- data.frame(apply(m, 1:2, sd))
standard_errors <-  standard_devs / sqrt(replicates)

upper_conf <- means + standard_errors * 1.96
lower_conf <- means - standard_errors * 1.96

label_frame <- rep_1_frame[,1:4] %>% mutate(comboname = paste(bxb1, mix, sep = "\n"))

means2 <- cbind(label_frame,means)
upper_conf2 <- cbind(label_frame,upper_conf)
lower_conf2 <- cbind(label_frame,lower_conf)
#ga_data_filtered$attp <- factor(ga_data_filtered$attp, levels = c("gt","ga"))

upper_conf2_melted <- melt(upper_conf2[,c("comboname","attp","gfp","mcherry")], id = c("comboname","attp"))
lower_conf2_melted <- melt(lower_conf2[,c("comboname","attp","gfp","mcherry")], id = c("comboname","attp"))
means2_melted <- melt(means2[,c("comboname","attp","gfp","mcherry")], id = c("comboname","attp"))
means2_melted$geomean <- 10^means2_melted$value
means2_melted$geomean_rounded <- round(means2_melted$geomean,1)

green_red_colorscale <- c(gfp = "green", mcherry = "red")

GT_vs_GA_plot <- ggplot() + 
  theme_classic() + theme(axis.text.x = element_text(angle = -90, vjust = 0.5)) +
  ylab("% cells recombined") +
  xlab(NULL) +
  scale_y_log10() +
  scale_color_manual(values = green_red_colorscale) + 
  geom_errorbar(data = upper_conf2_melted, aes(x = comboname, color = variable, ymin = 10^lower_conf2_melted$value, ymax = 10^value),
                alpha = 0.5, width = 0.2, position = position_dodge(width = 0.5)) +
  geom_point(data = means2_melted, aes(x = comboname, y = geomean, color = variable), alpha = 0.5, position = position_dodge(width = 0.5)) +
  facet_wrap(~attp)
print(GT_vs_GA_plot)

ggsave(file = "Plots/GT_vs_GA_plot.pdf", GT_vs_GA_plot, height = 1.8, width = 3.8)
```
## Flanking recombinase sites to excise unwanted bacterial DNA sequences
## Relevant to Figure 2


```{r Flanking data - Flow Cyometry}
## Making a data frame to keep track of how many unexcised cells are observed with G718A
flanking_flow_df <- data.frame("flow" = NA, "frac_unexcised" = NA)
flank_red_pos_cutoff <- 3e3


## F84
f84_none <- read.csv(file = "data/flow/Flanking/F84/F84_None_unselected_Sample(1).csv.gz") %>% mutate("flow" = "F84", sample = "None")
f84_g718a <- read.csv(file = "data/flow/Flanking/F84/F84_G718A_hygro_Sample(4).csv.gz") %>% mutate("flow" = "F84", sample = "Flanked")
f84_g747a <- read.csv(file = "data/flow/Flanking/F84/F84_G747A_hygro_Sample(5).csv.gz") %>% mutate("flow" = "F84", sample = "Control")

flank_expt_cell_num <- 90000
flank_red_pos_cutoff <- 5e3

f84 <- rbind(f84_none[1:flank_expt_cell_num,], f84_g718a[1:flank_expt_cell_num,], f84_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f84, aes(x = YL2.A, y = BL1.A), alpha = 0.01) + facet_grid(rows = vars(sample))

f84$ratio <- f84$BL1.A / f84$YL2.A
f84_subset <- f84 %>% filter(YL2.A >= flank_red_pos_cutoff)

f84_control_95pct_interval <- c(quantile((f84_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f84_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f84_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f84_control_95pct_interval, linetype = 2)

f84_fraction_unexcised <- sum((f84_subset %>% filter(sample == "Flanked"))$ratio > f84_control_95pct_interval[1])/nrow((f84_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F84", "frac_unexcised" = f84_fraction_unexcised))

## F114
f114_g718a <- read.csv(file = "data/flow/Flanking/F114/F114_G718A_unselected_A2.csv.gz") %>% mutate("flow" = "F114", sample = "Flanked")
f114_g747a <- read.csv(file = "data/flow/Flanking/F114/F114_G747A_unselected_A3.csv.gz") %>% mutate("flow" = "F114", sample = "Control")

flank_expt_cell_num <- 45000
flank_red_pos_cutoff <- 3e3

f114 <- rbind(f114_g718a[1:flank_expt_cell_num,], f114_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f114, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f114$ratio <- f114$BL1.A / f114$YL2.A
f114_subset <- f114 %>% filter(YL2.A >= flank_red_pos_cutoff)

f114_control_95pct_interval <- c(quantile((f114_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f114_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f114_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f114_control_95pct_interval, linetype = 2)

f114_fraction_unexcised <- sum((f114_subset %>% filter(sample == "Flanked"))$ratio > f114_control_95pct_interval[1])/nrow((f114_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F114", "frac_unexcised" = f114_fraction_unexcised))


## F116
f116_g718a <- read.csv(file = "data/flow/Flanking/F116/F116_G718A_unselected_H11.csv.gz") %>% mutate("flow" = "F116", sample = "Flanked")
f116_g747a <- read.csv(file = "data/flow/Flanking/F116/F116_G747A_unselected_H12.csv.gz") %>% mutate("flow" = "F116", sample = "Control")
flank_expt_cell_num <- 25000
flank_red_pos_cutoff <- 5e3
f116 <- rbind(f116_g718a[1:flank_expt_cell_num,], f116_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f116, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f116$ratio <- f116$BL1.A / f116$YL2.A
f116_subset <- f116 %>% filter(YL2.A >= flank_red_pos_cutoff)

f116_control_95pct_interval <- c(quantile((f116_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f116_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f116_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f116_control_95pct_interval, linetype = 2)

f116_fraction_unexcised <- sum((f116_subset %>% filter(sample == "Flanked"))$ratio > f116_control_95pct_interval[1])/nrow((f116_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F116", "frac_unexcised" = f116_fraction_unexcised))


## F119
f119_none <- read.csv(file = "data/flow/Flanking/F119/F119_none_unselected_G11.csv.gz") %>% mutate("flow" = "F119", sample = "None")
f119_g718a <- read.csv(file = "data/flow/Flanking/F119/F119_G718A_unselected_H11.csv.gz") %>% mutate("flow" = "F119", sample = "Flanked")
f119_g747a <- read.csv(file = "data/flow/Flanking/F119/F119_G747A_unselected_H12.csv.gz") %>% mutate("flow" = "F119", sample = "Control")
flank_expt_cell_num <- 9000
flank_red_pos_cutoff <- 5e3
f119 <- rbind(f119_none[1:flank_expt_cell_num,], f119_g718a[1:flank_expt_cell_num,], f119_g747a[1:flank_expt_cell_num,]) %>% filter(!is.na(sample))

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f119, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f119$ratio <- f119$BL1.A / f119$YL2.A
f119_subset <- f119 %>% filter(YL2.A >= flank_red_pos_cutoff)

f119_control_95pct_interval <- c(quantile((f119_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f119_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f119_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f119_control_95pct_interval, linetype = 2)

f119_fraction_unexcised <- sum((f119_subset %>% filter(sample == "Flanked"))$ratio > f119_control_95pct_interval[1])/nrow((f119_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F119", "frac_unexcised" = f119_fraction_unexcised))

## F130
f130_none <- read.csv(file = "data/flow/Flanking/F130/F130_none_unselected_D10.csv.gz") %>% mutate("flow" = "F130", sample = "None")
f130_g718a <- read.csv(file = "data/flow/Flanking/F130/F130_G718A_unselected_D11.csv.gz") %>% mutate("flow" = "F130", sample = "Flanked")
f130_g747a <- read.csv(file = "data/flow/Flanking/F130/F130_G747A_unselected_D12.csv.gz") %>% mutate("flow" = "F130", sample = "Control")
flank_expt_cell_num <- 89000
flank_red_pos_cutoff <- 5e3
f130 <- rbind(f130_none[1:flank_expt_cell_num,], f130_g718a[1:flank_expt_cell_num,], f130_g747a[1:flank_expt_cell_num,]) %>% filter(!is.na(sample))

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f130, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f130$ratio <- f130$BL1.A / f130$YL2.A
f130_subset <- f130 %>% filter(YL2.A >= flank_red_pos_cutoff)

f130_control_95pct_interval <- c(quantile((f130_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f130_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f130_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f130_control_95pct_interval, linetype = 2)

f130_fraction_unexcised <- sum((f130_subset %>% filter(sample == "Flanked"))$ratio > f130_control_95pct_interval[1])/nrow((f130_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F130", "frac_unexcised" = f130_fraction_unexcised))

## F131
f131_none <- read.csv(file = "data/flow/Flanking/F131/F131_none_unselected_D4.csv.gz") %>% mutate("flow" = "F131", sample = "None")
f131_g718a <- read.csv(file = "data/flow/Flanking/F131/F131_G718A_unselected_D5.csv.gz") %>% mutate("flow" = "F131", sample = "Flanked")
f131_g747a <- read.csv(file = "data/flow/Flanking/F131/F131_G747A_unselected_D6.csv.gz") %>% mutate("flow" = "F131", sample = "Control")
flank_expt_cell_num <- 190000
flank_red_pos_cutoff <- 3e3
f131 <- rbind(f131_none[1:flank_expt_cell_num,], f131_g718a[1:flank_expt_cell_num,], f131_g747a[1:flank_expt_cell_num,]) %>% filter(!is.na(sample))

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f131, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f131$ratio <- f131$BL1.A / f131$YL2.A
f131_subset <- f131 %>% filter(YL2.A >= flank_red_pos_cutoff)

f131_control_95pct_interval <- c(quantile((f131_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f131_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f131_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f131_control_95pct_interval, linetype = 2)

f131_fraction_unexcised <- sum((f131_subset %>% filter(sample == "Flanked"))$ratio > f131_control_95pct_interval[1])/nrow((f131_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F131", "frac_unexcised" = f131_fraction_unexcised))

## F132
f132_g718a <- read.csv(file = "data/flow/Flanking/F132/F132_G718A_unselected_G12.csv.gz") %>% mutate("flow" = "F132", sample = "Flanked")
f132_g747a <- read.csv(file = "data/flow/Flanking/F132/F132_G747A_unselected_H1.csv.gz") %>% mutate("flow" = "F132", sample = "Control")
flank_expt_cell_num <- 180000
flank_red_pos_cutoff <- 3e3
f132 <- rbind(f132_g718a[1:flank_expt_cell_num,], f132_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f132, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f132$ratio <- f132$BL1.A / f132$YL2.A
f132_subset <- f132 %>% filter(YL2.A >= flank_red_pos_cutoff)

f132_control_95pct_interval <- c(quantile((f132_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f132_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f132_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f132_control_95pct_interval, linetype = 2)

f132_fraction_unexcised <- sum((f132_subset %>% filter(sample == "Flanked"))$ratio > f132_control_95pct_interval[1])/nrow((f132_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F132", "frac_unexcised" = f132_fraction_unexcised))

## f280
f280_g718a <- read.csv(file = "data/flow/Flanking/F280/F280_G718A_unselected.csv.gz") %>% mutate("flow" = "F280", sample = "Flanked")
f280_g747a <- read.csv(file = "data/flow/Flanking/F280/F280_G747A_unselected.csv.gz") %>% mutate("flow" = "F280", sample = "Control")
flank_expt_cell_num <- 180000
flank_red_pos_cutoff <- 3e3
f280 <- rbind(f280_g718a[1:flank_expt_cell_num,], f280_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f280, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f280$ratio <- f280$BL1.A / f280$YL2.A
f280_subset <- f280 %>% filter(YL2.A >= flank_red_pos_cutoff)

f280_control_95pct_interval <- c(quantile((f280_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f280_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f280_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f280_control_95pct_interval, linetype = 2)

f280_fraction_unexcised <- sum((f280_subset %>% filter(sample == "Flanked"))$ratio > f280_control_95pct_interval[1])/nrow((f280_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F280", "frac_unexcised" = f280_fraction_unexcised))


## f281
f281_g718a <- read.csv(file = "data/flow/Flanking/F281/F281_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F281", sample = "Flanked")
f281_g747a <- read.csv(file = "data/flow/Flanking/F281/F281_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F281", sample = "Control")
flank_expt_cell_num <- 310000
flank_red_pos_cutoff <- 3e3
f281 <- rbind(f281_g718a[1:flank_expt_cell_num,], f281_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f281, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f281$ratio <- f281$BL1.A / f281$YL2.A
f281_subset <- f281 %>% filter(YL2.A >= flank_red_pos_cutoff)

f281_control_95pct_interval <- c(quantile((f281_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f281_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f281_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f281_control_95pct_interval, linetype = 2)

f281_fraction_unexcised <- sum((f281_subset %>% filter(sample == "Flanked"))$ratio > f281_control_95pct_interval[1])/nrow((f281_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F281", "frac_unexcised" = f281_fraction_unexcised))


## f281b
f281b_g718a <- read.csv(file = "data/flow/Flanking/F281/F281b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F281b", sample = "Flanked")
flank_expt_cell_num <- 230000
flank_red_pos_cutoff <- 3e3
f281b <- rbind(f281b_g718a[1:flank_expt_cell_num,], f281_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f281b, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f281b$ratio <- f281b$BL1.A / f281b$YL2.A
f281b_subset <- f281b %>% filter(YL2.A >= flank_red_pos_cutoff)

f281b_control_95pct_interval <- c(quantile((f281b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f281b_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f281b_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f281b_control_95pct_interval, linetype = 2)

f281b_fraction_unexcised <- sum((f281b_subset %>% filter(sample == "Flanked"))$ratio > f281b_control_95pct_interval[1])/nrow((f281b_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F281b", "frac_unexcised" = f281b_fraction_unexcised))


## f282
f282_g718a <- read.csv(file = "data/flow/Flanking/F282/F282_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F282", sample = "Flanked")
f282_g747a <- read.csv(file = "data/flow/Flanking/F282/F282_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F282", sample = "Control")
flank_expt_cell_num <- 280000
flank_red_pos_cutoff <- 3e3
f282 <- rbind(f282_g718a[1:flank_expt_cell_num,], f282_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f282, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f282$ratio <- f282$BL1.A / f282$YL2.A
f282_subset <- f282 %>% filter(YL2.A >= flank_red_pos_cutoff)

f282_control_95pct_interval <- c(quantile((f282_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f282_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f282_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample), scales="free_y") + geom_vline(xintercept  = f282_control_95pct_interval, linetype = 2)

f282_fraction_unexcised <- sum((f282_subset %>% filter(sample == "Flanked"))$ratio > f282_control_95pct_interval[1])/nrow((f282_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F282", "frac_unexcised" = f282_fraction_unexcised))


## f282b
f282b_g718a <- read.csv(file = "data/flow/Flanking/F282/F282b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F282b", sample = "Flanked")
flank_expt_cell_num <- 930000
flank_red_pos_cutoff <- 3e3
f282b <- rbind(f282b_g718a[1:flank_expt_cell_num,], f282_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f282b, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f282b$ratio <- f282b$BL1.A / f282b$YL2.A
f282b_subset <- f282b %>% filter(YL2.A >= flank_red_pos_cutoff)

f282b_control_95pct_interval <- c(quantile((f282b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f282b_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f282b_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample), scales="free_y") + geom_vline(xintercept  = f282b_control_95pct_interval, linetype = 2)

f282b_fraction_unexcised <- sum((f282b_subset %>% filter(sample == "Flanked"))$ratio > f282b_control_95pct_interval[1])/nrow((f282b_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F282b", "frac_unexcised" = f282b_fraction_unexcised))


## f283
f283_g718a <- read.csv(file = "data/flow/Flanking/F283/F283_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F283", sample = "Flanked")
f283_g747a <- read.csv(file = "data/flow/Flanking/F283/F283_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F283", sample = "Control")
flank_expt_cell_num <- 305000
flank_red_pos_cutoff <- 3e3
f283 <- rbind(f283_g718a[1:flank_expt_cell_num,], f283_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) +geom_point(data = f283, aes(x = YL2.A, y = BL1.A), alpha = 0.1) +facet_grid(rows = vars(sample))

f283$ratio <- f283$BL1.A / f283$YL2.A
f283_subset <- f283 %>% filter(YL2.A >= flank_red_pos_cutoff)

f283_control_95pct_interval <- c(quantile((f283_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f283_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f283_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample), scales="free_y") + geom_vline(xintercept  = f283_control_95pct_interval, linetype = 2)

f283_fraction_unexcised <- sum((f283_subset %>% filter(sample == "Flanked"))$ratio > f283_control_95pct_interval[1])/nrow((f283_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F283", "frac_unexcised" = f283_fraction_unexcised))


## f283b
f283b_g718a <- read.csv(file = "data/flow/Flanking/F283/F283b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F283b", sample = "Flanked")
flank_expt_cell_num <- 80000
flank_red_pos_cutoff <- 3e3
f283b <- rbind(f283b_g718a[1:flank_expt_cell_num,], f283_g747a[1:flank_expt_cell_num,])

#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f283b, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))

f283b$ratio <- f283b$BL1.A / f283b$YL2.A
f283b_subset <- f283b %>% filter(YL2.A >= flank_red_pos_cutoff)

f283b_control_95pct_interval <- c(quantile((f283b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f283b_subset %>% filter(sample == "Control"))$ratio,0.975))

#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f283b_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample), scales="free_y") + geom_vline(xintercept  = f283b_control_95pct_interval, linetype = 2)

f283b_fraction_unexcised <- sum((f283b_subset %>% filter(sample == "Flanked"))$ratio > f283b_control_95pct_interval[1])/nrow((f283b_subset %>% filter(sample == "Flanked")))
flanking_flow_df <- rbind(flanking_flow_df, data.frame("flow" = "F283b", "frac_unexcised" = f283b_fraction_unexcised))
```

```{r Some more flanking data - Flow Cyometry}
## Some summary graphs
## An example set of samples
f131_example <- rbind(f131_none[1:20000,], f131_g718a[1:20000,], f131_g747a[1:20000,]) %>% filter(!is.na(sample))
f131_example$sample <- factor(f131_example$sample, levels = c("None", "Flanked", "Control"))
F131_example_scatterplot <- ggplot() + theme(panel.grid.minor = element_blank()) + 
  scale_x_log10(limits = c(1e1,1e5), expand = c(0,0), breaks = c(1e2,1e4)) + 
  scale_y_log10(limits = c(1e1,1e4), expand = c(0,0), breaks = c(1e2,1e3)) + 
  geom_point(data = f131_example, aes(x = YL2.A, y = BL1.A), alpha = 0.1, size = 0.5) +
  facet_grid(cols = vars(sample))
ggsave(file = "plots/Flanking_example_scatterplot.pdf", F131_example_scatterplot, height = 1.6, width = 4.5)
F131_example_scatterplot

## To respond to reviewer comments, compare the fraction of red cells in the control and flanked transfections. 
ggplot() + scale_x_log10() +
  geom_histogram(data = f131_example, aes(x = YL2.A)) + facet_grid(rows = vars(sample))

paste("F131 flanked recomb percentage:",round(nrow(subset(f131_g718a, YL2.A > 1000)) / nrow(f131_g718a) * 100,1))
paste("F131 control recomb percentage:",round(nrow(subset(f131_g747a, YL2.A > 1000)) / nrow(f131_g747a) * 100,1))

f131_example <- rbind(f131_none[1:194000,], f131_g718a[1:194000,], f131_g747a[1:194000,]) %>% filter(!is.na(sample))
f131_example$sample <- factor(f131_example$sample, levels = c("None", "Flanked", "Control"))
f131_example$ratio <- (f131_example$BL1.A + 250) / (f131_example$YL2.A + 250)
f131_example_subset <- f131_example %>% filter(YL2.A >= flank_red_pos_cutoff & sample != "None")

f131_example_control_95pct_interval <- c(quantile((f131_example_subset %>% filter(sample == "Control"))$ratio,0.05))

F131_ratio_histogram <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), panel.grid.major.y = element_blank()) + 
  labs(x = "Green / Red ratio", y = "Number\nof cells") +
  scale_x_continuous(limits = c(0, 0.6), breaks = c(0, 0.2, 0.4, 0.6), expand = c(0,0)) + scale_y_continuous(breaks = c(0,100)) +
  geom_histogram(data = f131_example_subset, aes(x = ratio), binwidth = 0.01) + facet_grid(rows = vars(sample)) + geom_vline(xintercept  = f131_example_control_95pct_interval, linetype = 2)
ggsave(file = "Plots/F131_ratio_histogram.pdf", F131_ratio_histogram, height = 1.4, width = 1.8)
F131_ratio_histogram

## To respond to reviewer comments, compare the fraction of pared ratio cells to unexcised
paste("F131 flanked and unexcised:",round(nrow(subset(f131_example_subset, sample == "Flanked" & YL2.A > 1000 & ratio > 0.15)) / nrow(subset(f131_example_subset, sample == "Flanked" & YL2.A > 1000)) * 100,1))
paste("F131 control and unexcised:",round(nrow(subset(f131_example_subset, sample == "Control" & YL2.A > 1000 & ratio > 0.15)) / nrow(subset(f131_example_subset, sample == "Control" & YL2.A > 1000)) * 100,1))

```

```{r Flanking data averages}
## Looking at fraction excised over time
flanking_flow_days <- read.csv(file = "data/Flanking_expt_timepoints.csv", header = T, stringsAsFactors = F)
flanking_flow_df2 <- merge(flanking_flow_df, flanking_flow_days, by = "flow")

flanking_flow_df2$replicate <- as.factor(flanking_flow_df2$replicate)

Flanking_timeplot <- ggplot() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_y_log10(limits = c(0.001,1)) + scale_x_log10() +
  labs(x = "Days after\ntransfection", y = "Fraction of\ncells unexcised") +
  geom_point(data = flanking_flow_df2 %>% filter(instrument == "Attune"), aes(x = days, y = frac_unexcised, color = replicate), color = "red")
ggsave(file = "plots/Flanking_timeplot.pdf", Flanking_timeplot, height = 1.6, width = 2)
Flanking_timeplot


## Looking at the reproducibility of recombination
flank_indep_recombs <- flanking_flow_df2 %>% filter(frac_unexcised != 0) %>% group_by(replicate) %>% summarize(frac_unexcised = 10^mean(log10(frac_unexcised)))
flank_flow_geomean <- 10^mean(log10(flank_indep_recombs$frac_unexcised))
flank_flow_upper_conf <- 10^(mean(log10(flank_indep_recombs$frac_unexcised)) + (sd(log10(flank_indep_recombs$frac_unexcised))/sqrt(nrow(flank_indep_recombs)) * 1.96))
flank_flow_lower_conf <- 10^(mean(log10(flank_indep_recombs$frac_unexcised)) - (sd(log10(flank_indep_recombs$frac_unexcised))/sqrt(nrow(flank_indep_recombs)) * 1.96))

Repeated_excisions_plot <- ggplot() + theme(panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), axis.text.x = element_blank(), axis.ticks.x = element_blank()) + 
  labs(x = NULL, y = "Fraction of cells unexcised") + 
  scale_y_log10(limits = c(0.001,1)) + #scale_x_continuous(limits = c(-0.1,0.1)) +
  geom_point(data = flank_indep_recombs, aes(x = 0, y = frac_unexcised), alpha = 0.5, color = "Red") +
  geom_point(aes(x = 0, y = flank_flow_geomean), shape = 95, size = 10, color = "Black") +
  geom_errorbar(aes(x = 0, ymax = flank_flow_upper_conf, ymin = flank_flow_lower_conf), width = 0.01, alpha = 0.4)
ggsave(file = "plots/Repeated_excisions_plot.pdf", Repeated_excisions_plot, height = 1.1, width = 1)
Repeated_excisions_plot
```

```{r Effect of adding AP1903 to flanking recombinations}
## f280
f280_g718a <- read.csv(file = "data/flow/Flanking/F280/F280_G718A_unselected.csv.gz") %>% mutate("flow" = "F280", sample = "Flanked", treatment = "none")
f280_g747a <- read.csv(file = "data/flow/Flanking/F280/F280_G747A_unselected.csv.gz") %>% mutate("flow" = "F280", sample = "Control", treatment = "none")
f280_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F280/F280_G718A_AP1903.csv.gz") %>% mutate("flow" = "F280", sample = "Flanked", treatment = "AP1903")
f280_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F280/F280_G747A_AP1903.csv.gz") %>% mutate("flow" = "F280", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 374000
flank_red_pos_cutoff <- 3e3
f280 <- rbind(f280_g718a[1:flank_expt_cell_num,], f280_g747a[1:flank_expt_cell_num,], f280_g718a_ap1903[1:flank_expt_cell_num,], f280_g747a_ap1903[1:flank_expt_cell_num,])
#ggplot() + theme_bw() + scale_x_log10(limits = c(1e1,1e5), expand = c(0,0)) + scale_y_log10(limits = c(1e1,1e5), expand = c(0,0)) + geom_vline(xintercept = flank_red_pos_cutoff, linetype = 2) + geom_point(data = f280, aes(x = YL2.A, y = BL1.A), alpha = 0.1) + facet_grid(rows = vars(sample))
f280$ratio <- f280$BL1.A / f280$YL2.A
f280_subset <- f280 %>% filter(YL2.A >= flank_red_pos_cutoff)
f280_control_95pct_interval <- c(quantile((f280_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f280_subset %>% filter(sample == "Control"))$ratio,0.975))
f280_subset$label <- paste0(f280_subset$sample,"_",f280_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f280_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f280_control_95pct_interval, linetype = 2)
f280_subset_summary1 <- f280_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f280_control_95pct_interval))
f280_subset_summary2 <- f280_subset %>% group_by(label) %>% count()
f280_subset_summary <- merge(f280_subset_summary1, f280_subset_summary2, by = "label")
f280_subset_summary$frac_unexcised <- f280_subset_summary$unexcised / f280_subset_summary$n
f280_subset_summary$rep <- "F280"

## F281
f281_g718a <- read.csv(file = "data/flow/Flanking/F281/F281_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F281", sample = "Flanked", treatment = "none")
f281_g747a <- read.csv(file = "data/flow/Flanking/F281/F281_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F281", sample = "Control", treatment = "none")
f281_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F281/F281_G718A_AP1903_B1.csv.gz") %>% mutate("flow" = "F281", sample = "Flanked", treatment = "AP1903")
f281_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F281/F281_G747A_AP1903_B2.csv.gz") %>% mutate("flow" = "F281", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 287000
flank_red_pos_cutoff <- 3e3
f281 <- rbind(f281_g718a[1:flank_expt_cell_num,], f281_g747a[1:flank_expt_cell_num,], f281_g718a_ap1903[1:flank_expt_cell_num,], f281_g747a_ap1903[1:flank_expt_cell_num,])
f281$ratio <- f281$BL1.A / f281$YL2.A
f281_subset <- f281 %>% filter(YL2.A >= flank_red_pos_cutoff)
f281_control_95pct_interval <- c(quantile((f281_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f281_subset %>% filter(sample == "Control"))$ratio,0.975))
f281_subset$label <- paste0(f281_subset$sample,"_",f281_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f281_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f281_control_95pct_interval, linetype = 2)
f281_subset_summary1 <- f281_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f281_control_95pct_interval))
f281_subset_summary2 <- f281_subset %>% group_by(label) %>% count()
f281_subset_summary <- merge(f281_subset_summary1, f281_subset_summary2, by = "label")
f281_subset_summary$frac_unexcised <- f281_subset_summary$unexcised / f281_subset_summary$n
f281_subset_summary$rep <- "F281"

## F281b
f281b_g718a <- read.csv(file = "data/flow/Flanking/F281/F281b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F281b", sample = "Flanked", treatment = "none")
f281b_g747a <- read.csv(file = "data/flow/Flanking/F281/F281b_G747A_unselected_A4.csv.gz") %>% mutate("flow" = "F281b", sample = "Control", treatment = "none")
f281b_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F281/F281b_G718A_AP1903_B3.csv.gz") %>% mutate("flow" = "F281b", sample = "Flanked", treatment = "AP1903")
f281b_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F281/F281b_G747A_AP1903_B4.csv.gz") %>% mutate("flow" = "F281b", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 157000
flank_red_pos_cutoff <- 3e3
f281b <- rbind(f281b_g718a[1:flank_expt_cell_num,], f281b_g747a[1:flank_expt_cell_num,], f281b_g718a_ap1903[1:flank_expt_cell_num,], f281b_g747a_ap1903[1:flank_expt_cell_num,])
f281b$ratio <- f281b$BL1.A / f281b$YL2.A
f281b_subset <- f281b %>% filter(YL2.A >= flank_red_pos_cutoff)
f281b_control_95pct_interval <- c(quantile((f281b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f281b_subset %>% filter(sample == "Control"))$ratio,0.975))
f281b_subset$label <- paste0(f281b_subset$sample,"_",f281b_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f281b_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f281b_control_95pct_interval, linetype = 2)
f281b_subset_summary1 <- f281b_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f281b_control_95pct_interval))
f281b_subset_summary2 <- f281b_subset %>% group_by(label) %>% count()
f281b_subset_summary <- merge(f281b_subset_summary1, f281b_subset_summary2, by = "label")
f281b_subset_summary$frac_unexcised <- f281b_subset_summary$unexcised / f281b_subset_summary$n
f281b_subset_summary$rep <- "F281b"


## F282
f282_g718a <- read.csv(file = "data/flow/Flanking/F282/F282_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F282", sample = "Flanked", treatment = "none")
f282_g747a <- read.csv(file = "data/flow/Flanking/F282/F282_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F282", sample = "Control", treatment = "none")
f282_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F282/F282_G718A_AP1903_B1.csv.gz") %>% mutate("flow" = "F282", sample = "Flanked", treatment = "AP1903")
f282_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F282/F282_G747A_AP1903_B2.csv.gz") %>% mutate("flow" = "F282", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 190000
flank_red_pos_cutoff <- 3e3
f282 <- rbind(f282_g718a[1:flank_expt_cell_num,], f282_g747a[1:flank_expt_cell_num,], f282_g718a_ap1903[1:flank_expt_cell_num,], f282_g747a_ap1903[1:flank_expt_cell_num,])
f282$ratio <- f282$BL1.A / f282$YL2.A
f282_subset <- f282 %>% filter(YL2.A >= flank_red_pos_cutoff)
f282_control_95pct_interval <- c(quantile((f282_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f282_subset %>% filter(sample == "Control"))$ratio,0.975))
f282_subset$label <- paste0(f282_subset$sample,"_",f282_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f282_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f282_control_95pct_interval, linetype = 2)
f282_subset_summary1 <- f282_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f282_control_95pct_interval))
f282_subset_summary2 <- f282_subset %>% group_by(label) %>% count()
f282_subset_summary <- merge(f282_subset_summary1, f282_subset_summary2, by = "label")
f282_subset_summary$frac_unexcised <- f282_subset_summary$unexcised / f282_subset_summary$n
f282_subset_summary$rep <- "F282"

## F282b
f282b_g718a <- read.csv(file = "data/flow/Flanking/F282/F282b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F282b", sample = "Flanked", treatment = "none")
f282b_g747a <- read.csv(file = "data/flow/Flanking/F282/F282b_G747A_unselected_A4.csv.gz") %>% mutate("flow" = "F282b", sample = "Control", treatment = "none")
f282b_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F282/F282b_G718A_AP1903_B3.csv.gz") %>% mutate("flow" = "F282b", sample = "Flanked", treatment = "AP1903")
f282b_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F282/F282b_G747A_AP1903_B4.csv.gz") %>% mutate("flow" = "F282b", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 93000
flank_red_pos_cutoff <- 3e3
f282b <- rbind(f282b_g718a[1:flank_expt_cell_num,], f282b_g747a[1:flank_expt_cell_num,], f282b_g718a_ap1903[1:flank_expt_cell_num,], f282b_g747a_ap1903[1:flank_expt_cell_num,])
f282b$ratio <- f282b$BL1.A / f282b$YL2.A
f282b_subset <- f282b %>% filter(YL2.A >= flank_red_pos_cutoff)
f282b_control_95pct_interval <- c(quantile((f282b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f282b_subset %>% filter(sample == "Control"))$ratio,0.975))
f282b_subset$label <- paste0(f282b_subset$sample,"_",f282b_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f282b_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f282b_control_95pct_interval, linetype = 2)
f282b_subset_summary1 <- f282b_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f282b_control_95pct_interval))
f282b_subset_summary2 <- f282b_subset %>% group_by(label) %>% count()
f282b_subset_summary <- merge(f282b_subset_summary1, f282b_subset_summary2, by = "label")
f282b_subset_summary$frac_unexcised <- f282b_subset_summary$unexcised / f282b_subset_summary$n
f282b_subset_summary$rep <- "F282b"


## F283
f283_g718a <- read.csv(file = "data/flow/Flanking/F283/F283_G718A_unselected_A1.csv.gz") %>% mutate("flow" = "F283", sample = "Flanked", treatment = "none")
f283_g747a <- read.csv(file = "data/flow/Flanking/F283/F283_G747A_unselected_A2.csv.gz") %>% mutate("flow" = "F283", sample = "Control", treatment = "none")
f283_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F283/F283_G718A_AP1903_B1.csv.gz") %>% mutate("flow" = "F283", sample = "Flanked", treatment = "AP1903")
f283_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F283/F283_G747A_AP1903_B2.csv.gz") %>% mutate("flow" = "F283", sample = "Control", treatment = "AP1903")
f283_g718a_hygro <- read.csv(file = "data/flow/Flanking/F283/F283_G718A_Hygro_C1.csv.gz") %>% mutate("flow" = "F283b", sample = "Flanked", treatment = "Hygro")
f283_g747a_hygro <- read.csv(file = "data/flow/Flanking/F283/F283_G747A_Hygro_C2.csv.gz") %>% mutate("flow" = "F283b", sample = "Control", treatment = "Hygro")
flank_expt_cell_num <- 12000
flank_red_pos_cutoff <- 3e3
f283 <- rbind(f283_g718a[1:flank_expt_cell_num,], f283_g747a[1:flank_expt_cell_num,], f283_g718a_ap1903[1:flank_expt_cell_num,], f283_g747a_ap1903[1:flank_expt_cell_num,], f283_g718a_hygro[1:flank_expt_cell_num,], f283_g747a_hygro[1:flank_expt_cell_num,])
f283$ratio <- f283$BL1.A / f283$YL2.A
f283_subset <- f283 %>% filter(YL2.A >= flank_red_pos_cutoff)
f283_control_95pct_interval <- c(quantile((f283_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f283_subset %>% filter(sample == "Control"))$ratio,0.975))
f283_subset$label <- paste0(f283_subset$sample,"_",f283_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f283_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f283_control_95pct_interval, linetype = 2)
f283_subset_summary1 <- f283_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f283_control_95pct_interval))
f283_subset_summary2 <- f283_subset %>% group_by(label) %>% count()
f283_subset_summary <- merge(f283_subset_summary1, f283_subset_summary2, by = "label")
f283_subset_summary$frac_unexcised <- f283_subset_summary$unexcised / f283_subset_summary$n
f283_subset_summary$rep <- "F283"

## F283b
f283b_g718a <- read.csv(file = "data/flow/Flanking/F283/F283b_G718A_unselected_A3.csv.gz") %>% mutate("flow" = "F283b", sample = "Flanked", treatment = "none")
f283b_g747a <- read.csv(file = "data/flow/Flanking/F283/F283b_G747A_unselected_A4.csv.gz") %>% mutate("flow" = "F283b", sample = "Control", treatment = "none")
f283b_g718a_ap1903 <- read.csv(file = "data/flow/Flanking/F283/F283b_G718A_AP1903_B3.csv.gz") %>% mutate("flow" = "F283b", sample = "Flanked", treatment = "AP1903")
f283b_g747a_ap1903 <- read.csv(file = "data/flow/Flanking/F283/F283b_G747A_AP1903_B4.csv.gz") %>% mutate("flow" = "F283b", sample = "Control", treatment = "AP1903")
flank_expt_cell_num <- 83000
flank_red_pos_cutoff <- 3e3
f283b <- rbind(f283b_g718a[1:flank_expt_cell_num,], f283b_g747a[1:flank_expt_cell_num,], f283b_g718a_ap1903[1:flank_expt_cell_num,], f283b_g747a_ap1903[1:flank_expt_cell_num,])
f283b$ratio <- f283b$BL1.A / f283b$YL2.A
f283b_subset <- f283b %>% filter(YL2.A >= flank_red_pos_cutoff)
f283b_control_95pct_interval <- c(quantile((f283b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f283b_subset %>% filter(sample == "Control"))$ratio,0.975))
f283b_subset$label <- paste0(f283b_subset$sample,"_",f283b_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f283b_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f283b_control_95pct_interval, linetype = 2)
f283b_subset_summary1 <- f283b_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f283b_control_95pct_interval))
f283b_subset_summary2 <- f283b_subset %>% group_by(label) %>% count()
f283b_subset_summary <- merge(f283b_subset_summary1, f283b_subset_summary2, by = "label")
f283b_subset_summary$frac_unexcised <- f283b_subset_summary$unexcised / f283b_subset_summary$n
f283b_subset_summary$rep <- "F283b"

## F284
f284_g718a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F284/F284_G718A_AP1903_Hygro_A1.csv.gz") %>% mutate("flow" = "F284", sample = "Flanked", treatment = "AP1903_Hygro")
f284_g747a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F284/F284_G747A_AP1903_Hygro_A2.csv.gz") %>% mutate("flow" = "F284", sample = "Control", treatment = "AP1903_Hygro")
f284_g718a_hygro <- read.csv(file = "data/flow/Flanking/F284/F284_G718A_Hygro_B1.csv.gz") %>% mutate("flow" = "F284", sample = "Flanked", treatment = "Hygro")
f284_g747a_hygro <- read.csv(file = "data/flow/Flanking/F284/F284_G747A_Hygro_B2.csv.gz") %>% mutate("flow" = "F284", sample = "Control", treatment = "Hygro")
flank_expt_cell_num <- 51000
flank_red_pos_cutoff <- 3e3
f284 <- rbind(f284_g718a_ap1903_hygro[1:flank_expt_cell_num,], f284_g747a_ap1903_hygro[1:flank_expt_cell_num,], f284_g718a_hygro[1:flank_expt_cell_num,], f284_g747a_hygro[1:flank_expt_cell_num,])
f284$ratio <- f284$BL1.A / f284$YL2.A
f284_subset <- f284 %>% filter(YL2.A >= flank_red_pos_cutoff)
f284_control_95pct_interval <- c(quantile((f284_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f284_subset %>% filter(sample == "Control"))$ratio,0.975))
f284_subset$label <- paste0(f284_subset$sample,"_",f284_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f284_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f284_control_95pct_interval, linetype = 2)
f284_subset_summary1 <- f284_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f284_control_95pct_interval))
f284_subset_summary2 <- f284_subset %>% group_by(label) %>% count()
f284_subset_summary <- merge(f284_subset_summary1, f284_subset_summary2, by = "label")
f284_subset_summary$frac_unexcised <- f284_subset_summary$unexcised / f284_subset_summary$n
f284_subset_summary$rep <- "F284"

## F284b
f284b_g718a_hygro <- read.csv(file = "data/flow/Flanking/F284/F284b_G718A_Hygro_B3.csv.gz") %>% mutate("flow" = "F284b", sample = "Flanked", treatment = "Hygro")
f284b_g747a_hygro <- read.csv(file = "data/flow/Flanking/F284/F284b_G747A_Hygro_B4.csv.gz") %>% mutate("flow" = "F284b", sample = "Control", treatment = "Hygro")
f284b_g718a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F284/F284b_G718A_AP1903_Hygro_A3.csv.gz") %>% mutate("flow" = "F284b", sample = "Flanked", treatment = "AP1903_Hygro")
f284b_g747a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F284/F284b_G747A_AP1903_Hygro_A4.csv.gz") %>% mutate("flow" = "F284b", sample = "Control", treatment = "AP1903_Hygro")
flank_expt_cell_num <- 35000
flank_red_pos_cutoff <- 3e3
f284b <- rbind(f284b_g718a_hygro[1:flank_expt_cell_num,], f284b_g747a_hygro[1:flank_expt_cell_num,], f284b_g718a_ap1903_hygro[1:flank_expt_cell_num,], f284b_g747a_ap1903_hygro[1:flank_expt_cell_num,])
f284b$ratio <- f284b$BL1.A / f284b$YL2.A
f284b_subset <- f284b %>% filter(YL2.A >= flank_red_pos_cutoff)
f284b_control_95pct_interval <- c(quantile((f284b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f284b_subset %>% filter(sample == "Control"))$ratio,0.975))
f284b_subset$label <- paste0(f284b_subset$sample,"_",f284b_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f284b_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f284b_control_95pct_interval, linetype = 2)
f284b_subset_summary1 <- f284b_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f284b_control_95pct_interval))
f284b_subset_summary2 <- f284b_subset %>% group_by(label) %>% count()
f284b_subset_summary <- merge(f284b_subset_summary1, f284b_subset_summary2, by = "label")
f284b_subset_summary$frac_unexcised <- f284b_subset_summary$unexcised / f284b_subset_summary$n
f284b_subset_summary$rep <- "F284b"


## F285
f285_g718a_hygro <- read.csv(file = "data/flow/Flanking/F285/F285_G718A_Hygro_A3.csv.gz") %>% mutate("flow" = "F285", sample = "Flanked", treatment = "Hygro")
f285_g747a_hygro <- read.csv(file = "data/flow/Flanking/F285/F285_G747A_Hygro_A4.csv.gz") %>% mutate("flow" = "F285", sample = "Control", treatment = "Hygro")
flank_expt_cell_num <- 235000
flank_red_pos_cutoff <- 3e3
f285 <- rbind(f285_g718a_hygro[1:flank_expt_cell_num,], f285_g747a_hygro[1:flank_expt_cell_num,])
f285$ratio <- f285$BL1.A / f285$YL2.A
f285_subset <- f285 %>% filter(YL2.A >= flank_red_pos_cutoff)
f285_control_95pct_interval <- c(quantile((f285_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f285_subset %>% filter(sample == "Control"))$ratio,0.975))
f285_subset$label <- paste0(f285_subset$sample,"_",f285_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f285_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f285_control_95pct_interval, linetype = 2)
f285_subset_summary1 <- f285_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f285_control_95pct_interval))
f285_subset_summary2 <- f285_subset %>% group_by(label) %>% count()
f285_subset_summary <- merge(f285_subset_summary1, f285_subset_summary2, by = "label")
f285_subset_summary$frac_unexcised <- f285_subset_summary$unexcised / f285_subset_summary$n
f285_subset_summary$rep <- "F285"

## F286
f286_g718a_hygro <- read.csv(file = "data/flow/Flanking/F286/F286_G718A_Hygro_B1.csv.gz") %>% mutate("flow" = "F286", sample = "Flanked", treatment = "Hygro")
f286_g747a_hygro <- read.csv(file = "data/flow/Flanking/F286/F286_G747A_Hygro_B2.csv.gz") %>% mutate("flow" = "F286", sample = "Control", treatment = "Hygro")
flank_expt_cell_num <- 88000
flank_red_pos_cutoff <- 3e3
f286 <- rbind(f286_g718a_hygro[1:flank_expt_cell_num,], f286_g747a_hygro[1:flank_expt_cell_num,])
f286$ratio <- f286$BL1.A / f286$YL2.A
f286_subset <- f286 %>% filter(YL2.A >= flank_red_pos_cutoff)
f286_control_95pct_interval <- c(quantile((f286_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f286_subset %>% filter(sample == "Control"))$ratio,0.975))
f286_subset$label <- paste0(f286_subset$sample,"_",f286_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f286_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f286_control_95pct_interval, linetype = 2)
f286_subset_summary1 <- f286_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f286_control_95pct_interval))
f286_subset_summary2 <- f286_subset %>% group_by(label) %>% count()
f286_subset_summary <- merge(f286_subset_summary1, f286_subset_summary2, by = "label")
f286_subset_summary$frac_unexcised <- f286_subset_summary$unexcised / f286_subset_summary$n
f286_subset_summary$rep <- "F286"

## F286b
f286b_g718a_hygro <- read.csv(file = "data/flow/Flanking/F286/F286b_G718A_Hygro_B3.csv.gz") %>% mutate("flow" = "F286b", sample = "Flanked", treatment = "Hygro")
f286b_g747a_hygro <- read.csv(file = "data/flow/Flanking/F286/F286b_G747A_Hygro_B4.csv.gz") %>% mutate("flow" = "F286b", sample = "Control", treatment = "Hygro")
f286b_g718a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F286/F286b_G718A_AP1903_Hygro_A1.csv.gz") %>% mutate("flow" = "F286b", sample = "Flanked", treatment = "AP1903_Hygro")
f286b_g747a_ap1903_hygro <- read.csv(file = "data/flow/Flanking/F286/F286b_G747A_AP1903_Hygro_A2.csv.gz") %>% mutate("flow" = "F286b", sample = "Control", treatment = "AP1903_Hygro")
flank_expt_cell_num <- 51000
flank_red_pos_cutoff <- 3e3
f286b <- rbind(f286b_g718a_hygro[1:flank_expt_cell_num,], f286b_g747a_hygro[1:flank_expt_cell_num,], f286b_g718a_ap1903_hygro[1:flank_expt_cell_num,], f286b_g747a_ap1903_hygro[1:flank_expt_cell_num,])
f286b$ratio <- f286b$BL1.A / f286b$YL2.A
f286b_subset <- f286b %>% filter(YL2.A >= flank_red_pos_cutoff)
f286b_control_95pct_interval <- c(quantile((f286b_subset %>% filter(sample == "Control"))$ratio,0.05),quantile((f286b_subset %>% filter(sample == "Control"))$ratio,0.975))
f286b_subset$label <- paste0(f286b_subset$sample,"_",f286b_subset$treatment)
#ggplot() + theme_bw() + scale_x_continuous(limits = c(-0.05, 0.5)) + geom_histogram(data = f286b_subset, aes(x = ratio), binwidth = 0.01) + facet_wrap(~label, scales = "free_y") + geom_vline(xintercept  = f286b_control_95pct_interval, linetype = 2)
f286b_subset_summary1 <- f286b_subset %>% group_by(label) %>% summarize(unexcised = sum(ratio > f286b_control_95pct_interval))
f286b_subset_summary2 <- f286b_subset %>% group_by(label) %>% count()
f286b_subset_summary <- merge(f286b_subset_summary1, f286b_subset_summary2, by = "label")
f286b_subset_summary$frac_unexcised <- f286b_subset_summary$unexcised / f286b_subset_summary$n
f286b_subset_summary$rep <- "F286b"


## Now combining everything
flank_ap1903_summary <- rbind(f280_subset_summary, f281_subset_summary, f281b_subset_summary, f282_subset_summary, f282b_subset_summary, f283_subset_summary, f283b_subset_summary, f284_subset_summary, f284b_subset_summary, f285_subset_summary, f286_subset_summary, f286b_subset_summary)

flank_ap1903_summary2 <- flank_ap1903_summary %>% group_by(label) %>% summarize(geomean = 10^mean(log10(frac_unexcised)))

flank_ap1903_summary$label <- factor(flank_ap1903_summary$label, levels = c("Control_none", "Control_AP1903", "Control_Hygro", "Control_AP1903_Hygro", "Flanked_none", "Flanked_AP1903", "Flanked_Hygro", "Flanked_AP1903_Hygro"))
flank_ap1903_summary2$label <- factor(flank_ap1903_summary2$label, levels = c("Control_none", "Control_AP1903", "Control_Hygro", "Control_AP1903_Hygro", "Flanked_none", "Flanked_AP1903", "Flanked_Hygro", "Flanked_AP1903_Hygro"))

Flanking_AP1903_plot <- ggplot() + theme(panel.grid.major.x = element_blank(), axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1)) + 
  scale_y_log10(breaks = c(0.01,0.1,1)) + 
  labs(x = NULL, y = "Fraction of cells\nunexcised") + 
  geom_beeswarm(data = flank_ap1903_summary, aes(x = label, y = frac_unexcised), color = "red", alpha = 0.5, size = 1) +
  geom_point(data = flank_ap1903_summary2, aes(x = label, y = geomean), shape = 95, size = 8, color = "Black", alpha = 0.4) +
  NULL
ggsave(file = "plots/Flanking_AP1903_plot.pdf", Flanking_AP1903_plot, height = 2.25, width = 2.25)
Flanking_AP1903_plot

```

```{r}
## To respond to reviewer comments, compare the fraction of red cells in the control and flanked transfections. 
ggplot() + scale_x_log10() +
  geom_histogram(data = f131_example, aes(x = YL2.A)) + facet_grid(rows = vars(sample))

flanked_unselected_recomb_rates <- data.frame(rbind(
c("F130",nrow(subset(f130_g718a, YL2.A > 1000)) / nrow(f130_g718a) * 100,nrow(subset(f130_g747a, YL2.A > 1000)) / nrow(f130_g747a) * 100),
c("F131",nrow(subset(f131_g718a, YL2.A > 1000)) / nrow(f131_g718a) * 100,nrow(subset(f131_g747a, YL2.A > 1000)) / nrow(f131_g747a) * 100),
c("F132",nrow(subset(f132_g718a, YL2.A > 1000)) / nrow(f132_g718a) * 100,nrow(subset(f132_g747a, YL2.A > 1000)) / nrow(f132_g747a) * 100),
c("F280",nrow(subset(f280_g718a, YL2.A > 1000)) / nrow(f280_g718a) * 100,nrow(subset(f280_g747a, YL2.A > 1000)) / nrow(f280_g747a) * 100),
c("F282",nrow(subset(f282_g718a, YL2.A > 1000)) / nrow(f282_g718a) * 100,nrow(subset(f282_g747a, YL2.A > 1000)) / nrow(f282_g747a) * 100),
c("F283",nrow(subset(f283_g718a, YL2.A > 1000)) / nrow(f283_g718a) * 100,nrow(subset(f283_g747a, YL2.A > 1000)) / nrow(f283_g747a) * 100)))

colnames(flanked_unselected_recomb_rates) <- c("sample","flanked_rate","control_rate")

flanked_unselected_recomb_rates$flanked_rate <- as.numeric(flanked_unselected_recomb_rates$flanked_rate)
flanked_unselected_recomb_rates$control_rate <- as.numeric(flanked_unselected_recomb_rates$control_rate)

flanked_unselected_recomb_rates$fold_diff <- 0
flanked_unselected_recomb_rates$fold_diff <- flanked_unselected_recomb_rates$control_rate / flanked_unselected_recomb_rates$flanked_rate

mean(c(flanked_unselected_recomb_rates$flanked_rate,flanked_unselected_recomb_rates$control_rate))

##

control_red_not_green <- data.frame(rbind(
c("F130",nrow(subset(f130_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f130_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f130_subset, sample == "Control" & YL2.A > 1000))),
c("F131",nrow(subset(f131_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f131_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f131_subset, sample == "Control" & YL2.A > 1000))),
c("F132",nrow(subset(f132_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f132_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f132_subset, sample == "Control" & YL2.A > 1000))),
c("F280",nrow(subset(f280_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f280_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f280_subset, sample == "Control" & YL2.A > 1000))),
c("F282",nrow(subset(f282_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f282_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f282_subset, sample == "Control" & YL2.A > 1000))),
c("F283",nrow(subset(f283_subset, sample == "Control" & YL2.A > 1000 & ratio < quantile((f283_subset %>% filter(sample == "Flanked"))$ratio,0.95)))/nrow(subset(f283_subset, sample == "Control" & YL2.A > 1000)))
))
colnames(control_red_not_green) <- c("sample","control_redpos_grnneg")
control_red_not_green$control_redpos_grnneg <- as.numeric(control_red_not_green$control_redpos_grnneg)

10^mean(log10(control_red_not_green$control_redpos_grnneg))
```



```{r Flanking sequencing data}
flanked_seq <- read.csv(file = "data/Flanked_sequencing.csv", header = T, stringsAsFactors = F)

flanked_seq_geomean <- flanked_seq %>% group_by(sample, result) %>% summarize(geomean_fraction = 10^(mean(log10(fraction))))

Flanking_sequencing <- ggplot() + theme(panel.grid.major.x = element_blank(), legend.position = "top") + 
  labs(x = NULL, y = "Fraction of reads") +
  geom_point(data = flanked_seq, aes(x = sample, y = fraction, color = result), position = position_dodge(width = 0.5), alpha = 0.4) +
  geom_point(data = flanked_seq_geomean, aes(x = sample, y = geomean_fraction, color = result), position = position_dodge(width = 0.5), shape = 95, size = 10)
ggsave(file = "plots/Flanking_sequencing.pdf", Flanking_sequencing, height = 2, width = 2.2)
Flanking_sequencing

## Reviewer 2 asked how many reads underlied this figure. To answer this question, I created the table below.

flanked_read_counts <- flanked_seq %>% group_by(replicate, plasmid) %>% summarize(reads_sum = sum(reads))
print(flanked_read_counts$reads_sum)
```


## Creating and initially testing the double landing pad
## This is relevant to Figure 3

```{r Initial validation data for G542Ac3 and G783Ac2 cells}
c2_orig_none_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone2_none.csv.gz")
c2_orig_none <- c2_orig_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_orig_none) <- c("fsc","ssc","blu","grn","red","nir")
c2_orig_none$clone <- "c2"; c2_orig_none$treatment <- "none"
c2_orig_recomb_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone2_recomb.csv.gz")
c2_orig_recomb <- c2_orig_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_orig_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c2_orig_recomb$clone <- "c2"; c2_orig_recomb$treatment <- "recomb"
c2_orig_ap1903_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone2_ap1903.csv.gz")
c2_orig_ap1903 <- c2_orig_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_orig_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c2_orig_ap1903$clone <- "c2"; c2_orig_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c2_orig_none[1:cell_number,], c2_orig_recomb[1:cell_number,], c2_orig_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c2_orig_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c2_orig_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c2_orig_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c2_orig_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c2_orig_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c2_orig_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c2_orig_gr_plot

c2_orig_arranged_plot <- grid.arrange(c2_orig_bn_plot, c2_orig_bg_plot, c2_orig_br_plot, c2_orig_ng_plot, c2_orig_nr_plot, c2_orig_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201002/c2_orig_arranged_plot.png", c2_orig_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c2_orig_recomb$grn >= 3e4 & c2_orig_recomb$red >= 3e4) / nrow(c2_orig_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c2_orig_ap1903$grn >= 3e4 & c2_orig_ap1903$red >= 3e4) / nrow(c2_orig_ap1903) * 100,2))
```

```{r Initial validation data for G542Ac3 and G783Ac6 cells}
c6_orig_none_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone6_none.csv.gz")
c6_orig_none <- c6_orig_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_orig_none) <- c("fsc","ssc","blu","grn","red","nir")
c6_orig_none$clone <- "c6"; c6_orig_none$treatment <- "none"
c6_orig_recomb_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone6_recomb.csv.gz")
c6_orig_recomb <- c6_orig_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_orig_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c6_orig_recomb$clone <- "c6"; c6_orig_recomb$treatment <- "recomb"
c6_orig_ap1903_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone6_ap1903.csv.gz")
c6_orig_ap1903 <- c6_orig_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_orig_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c6_orig_ap1903$clone <- "c6"; c6_orig_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c6_orig_none[1:cell_number,], c6_orig_recomb[1:cell_number,], c6_orig_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c6_orig_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c6_orig_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c6_orig_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c6_orig_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c6_orig_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c6_orig_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c6_orig_gr_plot

c6_orig_arranged_plot <- grid.arrange(c6_orig_bn_plot, c6_orig_bg_plot, c6_orig_br_plot, c6_orig_ng_plot, c6_orig_nr_plot, c6_orig_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201002/c6_orig_arranged_plot.png", c6_orig_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c6_orig_recomb$grn >= 3e4 & c6_orig_recomb$red >= 3e4) / nrow(c6_orig_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c6_orig_ap1903$grn >= 3e4 & c6_orig_ap1903$red >= 3e4) / nrow(c6_orig_ap1903) * 100,2))
```

```{r Initial validation data for G542Ac3 and G783Ac11 cells}
c11_orig_none_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone11_none.csv.gz")
c11_orig_none <- c11_orig_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_orig_none) <- c("fsc","ssc","blu","grn","red","nir")
c11_orig_none$clone <- "c11"; c11_orig_none$treatment <- "none"
c11_orig_recomb_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone11_recomb.csv.gz")
c11_orig_recomb <- c11_orig_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_orig_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c11_orig_recomb$clone <- "c11"; c11_orig_recomb$treatment <- "recomb"
c11_orig_ap1903_raw <- read.csv(file = "data/flow/201002_F69_G783A_Clone_Comparisons/Clone11_ap1903.csv.gz")
c11_orig_ap1903 <- c11_orig_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_orig_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c11_orig_ap1903$clone <- "c11"; c11_orig_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c11_orig_none[1:cell_number,], c11_orig_recomb[1:cell_number,], c11_orig_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c11_orig_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c11_orig_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c11_orig_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c11_orig_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c11_orig_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c11_orig_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c11_orig_gr_plot

c11_orig_arranged_plot <- grid.arrange(c11_orig_bn_plot, c11_orig_bg_plot, c11_orig_br_plot, c11_orig_ng_plot, c11_orig_nr_plot, c11_orig_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201002/c11_orig_arranged_plot.png", c11_orig_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c11_orig_recomb$grn >= 3e4 & c11_orig_recomb$red >= 3e4) / nrow(c11_orig_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c11_orig_ap1903$grn >= 3e4 & c11_orig_ap1903$red >= 3e4) / nrow(c11_orig_ap1903) * 100,2))
```

Below is now repeat data for testing these three clones to confirm which is best

```{r Initial validation data for G542Ac3 and G783Ac2 cells repeated}
c2_none_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone2_none.csv.gz")
c2_none <- c2_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_none) <- c("fsc","ssc","blu","grn","red","nir")
c2_none$clone <- "c2"; c2_none$treatment <- "none"
c2_recomb_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone2_recomb.csv.gz")
c2_recomb <- c2_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c2_recomb$clone <- "c2"; c2_recomb$treatment <- "recomb"
c2_ap1903_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone2_ap1903.csv.gz")
c2_ap1903 <- c2_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c2_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c2_ap1903$clone <- "c2"; c2_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c2_none[1:cell_number,], c2_recomb[1:cell_number,], c2_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c2_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c2_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c2_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c2_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c2_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c2_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c2_gr_plot

c2_arranged_plot <- grid.arrange(c2_bn_plot, c2_bg_plot, c2_br_plot, c2_ng_plot, c2_nr_plot, c2_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201016/c2_arranged_plot.png", c2_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c2_recomb$grn >= 3e4 & c2_recomb$red >= 3e4) / nrow(c2_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c2_ap1903$grn >= 3e4 & c2_ap1903$red >= 3e4) / nrow(c2_ap1903) * 100,2))
```

```{r Initial validation data for G542Ac3 and G783Ac6 cells repeated}
c6_none_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone6_none.csv.gz")
c6_none <- c6_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_none) <- c("fsc","ssc","blu","grn","red","nir")
c6_none$clone <- "c6"; c6_none$treatment <- "none"
c6_recomb_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone6_recomb.csv.gz")
c6_recomb <- c6_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c6_recomb$clone <- "c6"; c6_recomb$treatment <- "recomb"
c6_ap1903_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone6_ap1903.csv.gz")
c6_ap1903 <- c6_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c6_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c6_ap1903$clone <- "c6"; c6_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c6_none[1:cell_number,], c6_recomb[1:cell_number,], c6_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c6_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c6_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c6_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c6_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c6_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c6_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c6_gr_plot

c6_arranged_plot <- grid.arrange(c6_bn_plot, c6_bg_plot, c6_br_plot, c6_ng_plot, c6_nr_plot, c6_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201016/c6_arranged_plot.png", c6_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c6_recomb$grn >= 3e4 & c6_recomb$red >= 3e4) / nrow(c6_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c6_ap1903$grn >= 3e4 & c6_ap1903$red >= 3e4) / nrow(c6_ap1903) * 100,2))
```

```{r Initial validation data for G542Ac3 and G783Ac11 cells repeated}
c11_none_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone11_none.csv.gz")
c11_none <- c11_none_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_none) <- c("fsc","ssc","blu","grn","red","nir")
c11_none$clone <- "c11"; c11_none$treatment <- "none"
c11_recomb_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone11_recomb.csv.gz")
c11_recomb <- c11_recomb_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_recomb) <- c("fsc","ssc","blu","grn","red","nir")
c11_recomb$clone <- "c11"; c11_recomb$treatment <- "recomb"
c11_ap1903_raw <- read.csv(file = "data/flow/201016_F74_G783A_Clone_Comparisons/Clone11_ap1903.csv.gz")
c11_ap1903 <- c11_ap1903_raw[,c("FSC.A","SSC.A","VL1.A","BL1.A","YL2.A","RL1.A")]; colnames(c11_ap1903) <- c("fsc","ssc","blu","grn","red","nir")
c11_ap1903$clone <- "c11"; c11_ap1903$treatment <- "ap1903"


cell_number <- 2000
combined_data <- rbind(c11_none[1:cell_number,], c11_recomb[1:cell_number,], c11_ap1903[1:cell_number,])
combined_data$treatment <- factor(combined_data$treatment, levels = c("recomb","none","ap1903"))

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

custom_color_scale <- c("none" = "magenta", "recomb" = "black", "ap1903" = "cyan")

c11_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = nir, color = treatment), alpha = plot_alpha)

c11_bg_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = grn, color = treatment), alpha = plot_alpha)

c11_br_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = blu, y = red, color = treatment), alpha = plot_alpha)

c11_ng_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = grn, color = treatment), alpha = plot_alpha)

c11_nr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = nir, y = red, color = treatment), alpha = plot_alpha)

c11_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data, aes(x = grn, y = red, color = treatment), alpha = plot_alpha) +
  geom_rect(mapping = aes(xmin = 3e4, xmax = 1e6, ymin = 3e4, ymax = 1e6), alpha = 0, color = "black", linetype = 2)
c11_gr_plot

c11_arranged_plot <- grid.arrange(c11_bn_plot, c11_bg_plot, c11_br_plot, c11_ng_plot, c11_nr_plot, c11_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201016/c11_arranged_plot.png", c11_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c11_recomb$grn >= 3e4 & c11_recomb$red >= 3e4) / nrow(c11_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c11_ap1903$grn >= 3e4 & c11_ap1903$red >= 3e4) / nrow(c11_ap1903) * 100,2))
```

## We're moving forward with clone 11, so show that data here 
```{r G542Ac3 and G783Ac11 cell plots for manuscript}

plot_alpha <- 0.2
axis_limits <- c(10,1e6)

c11_none_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "none"), aes(x = blu, y = nir), alpha = plot_alpha)
c11_none_bn_plot

c11_recomb_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "recomb") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "recomb"), aes(x = blu, y = nir), alpha = plot_alpha)
c11_recomb_bn_plot

c11_ap1903_bn_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "ap1903") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "ap1903"), aes(x = blu, y = nir), alpha = plot_alpha)
c11_ap1903_bn_plot

c11_none_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "none"), aes(x = grn, y = red), alpha = plot_alpha)
c11_none_gr_plot

c11_recomb_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "recomb") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "recomb"), aes(x = grn, y = red), alpha = plot_alpha)
c11_recomb_gr_plot

c11_ap1903_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "ap1903") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = combined_data %>% filter(treatment == "ap1903"), aes(x = grn, y = red), alpha = plot_alpha)
c11_ap1903_gr_plot


c11_arranged_plot <- grid.arrange(c11_bn_plot, c11_bg_plot, c11_br_plot, c11_ng_plot, c11_nr_plot, c11_gr_plot, ncol=6, nrow=1)

ggsave(file = "plots/201016/c11_arranged_plot.png", c11_arranged_plot, height = 3, width = 18)

paste("Percent double positive before selection:", round(sum(c11_recomb$grn >= 3e4 & c11_recomb$red >= 3e4) / nrow(c11_recomb) * 100,2))
paste("Percent double positive after selection:", round(sum(c11_ap1903$grn >= 3e4 & c11_ap1903$red >= 3e4) / nrow(c11_ap1903) * 100,2))


## Based on Reviewer #1s comment, lets figure out what the color profile of the green+/red- cells in the ap1903 treated sample are

c11_ap1903_red_partial <- combined_data %>% filter(treatment == "ap1903" & grn > 3e4 & red < 3e4)
c11_ap1903_red_partial_melt <- melt(c11_ap1903_red_partial[,c("blu","grn","red","nir")])

ggplot() + scale_x_log10() + 
  geom_histogram(data = c11_ap1903_red_partial_melt, aes(x = value)) +
  facet_grid(rows = vars(variable))

paste("Fraction of selected cells that were GFP-high and mCherry-low",nrow(combined_data %>% filter(treatment == "ap1903" & grn > 3e4 & red < 1e4))/nrow(combined_data %>% filter(treatment == "ap1903")))

## Looks like it's cells where the second, GA landing pad is not fully transcribing. There are definitely untransfected double landing pad cells that are Blue-high but NIR-low. I'll try seeing what percentage of the starting cells are that population.

c11_melt <- melt(subset(combined_data, treatment == "none")[,c("blu","grn","red","nir")])

ggplot() + scale_x_log10() + 
  geom_histogram(data = c11_melt, aes(x = value)) +
  facet_grid(rows = vars(variable))

c11_nir_partial <- combined_data %>% filter(treatment == "none" & blu > 5e2 & nir < 2e3)
c11_nir_partial_melt <- melt(c11_nir_partial[,c("blu","grn","red","nir")])

ggplot() + scale_x_log10() + 
  geom_histogram(data = c11_nir_partial_melt, aes(x = value)) +
  facet_grid(rows = vars(variable))

paste("Fraction of starting cells that were BFP-high and NIR-low",nrow(combined_data %>% filter(treatment == "none" & blu > 5e2 & nir < 2e3))/nrow(combined_data %>% filter(treatment == "none")))

## After all that, what percentage of AP1903 treated cells were Green-high and Red-low?

c11_ap1903_green_red <- combined_data %>% filter(treatment == "ap1903" & grn > 3e4 & red > 3e4)

paste("Fraction of selected cells that were GFP-high and mCherry-high",nrow(c11_ap1903_green_red)/nrow(combined_data %>% filter(treatment == "ap1903")))


```

Plotting the double landing pad PCA data for all relevant colors

```{r Principal component analysis}
# http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/118-principal-component-analysis-in-r-prcomp-vs-princomp/

combined_data <- rbind(c11_none[1:2000,], c11_recomb[1:2000,], c11_ap1903[1:2000,])

ggplot() + theme_bw() + scale_x_log10() + scale_y_log10() + geom_point(data = combined_data, aes(x = red, y = nir))

combined_data$zir <- combined_data$nir - combined_data$red * 0.01
ggplot() + theme_bw() + scale_x_log10() + scale_y_log10() + geom_point(data = combined_data, aes(x = red, y = zir))

combined_for_pca <- combined_data[,c("blu","grn","red","zir")]

clone2.pca <- prcomp(combined_for_pca, center = TRUE,scale. = TRUE)
fviz_eig(clone2.pca)

autoplot(clone2.pca, alpha = 0.2) +
  scale_x_continuous(limits = c(-0.3,0.3)) + scale_y_continuous(limits = c(-0.3,0.3))

fviz_pca_var(clone2.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )

clone2.pca$rotation

pca2 <- data.frame("pc1" = clone2.pca$x[,1], "pc2" = clone2.pca$x[,2])

pca2$treatment <- combined_data[,"treatment"]

pca2$treatment <- factor(pca2$treatment, levels = c("none","recomb","ap1903"))
custom_color_scale <- c("none" = "magenta", "ap1903" = "cyan", "recomb" = "green")

ggplot() + theme_bw() + theme(panel.grid.minor = element_blank()) + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_continuous(limits = c(-5,5)) +
  scale_y_continuous(limits = c(-3,2)) +
  geom_point(data = pca2, aes(x = pc1, y = pc2, color = treatment), alpha = 0.05)

ggplot() + theme_bw() + theme(panel.grid.minor = element_blank()) + 
  scale_color_manual(values = custom_color_scale) +
  scale_x_continuous(limits = c(-5,5)) +
  scale_y_continuous(limits = c(-3,2)) +
  geom_point(data = pca2, aes(x = pc1, y = pc2), alpha = 0.05) +
  facet_grid(rows = vars(treatment))

pca_details <- cbind(combined_for_pca,pca2)

pca_details_melted <- melt(pca_details[,c("treatment","blu","grn","red","zir")], id. = "treatment")

ggplot() + theme_bw() +
  scale_x_log10() + 
  geom_histogram(data = pca_details_melted, aes(x = value)) +
  facet_grid(cols = vars(variable))

blu_cutoff <- 1e3
grn_cutoff <- 1e4
red_cutoff <- 1e4
zir_cutoff <- 1e3

ggplot() + theme_bw() + scale_x_log10() + geom_histogram(data = pca_details, aes(x = blu)) + geom_vline(xintercept = blu_cutoff)
ggplot() + theme_bw() + scale_x_log10() + geom_histogram(data = pca_details, aes(x = grn)) + geom_vline(xintercept = grn_cutoff)
ggplot() + theme_bw() + scale_x_log10() + geom_histogram(data = pca_details, aes(x = red)) + geom_vline(xintercept = red_cutoff)
ggplot() + theme_bw() + scale_x_log10() + geom_histogram(data = pca_details, aes(x = zir)) + geom_vline(xintercept = zir_cutoff)

pca_details$blu2 <- 0
pca_details$grn2 <- 0
pca_details$red2 <- 0
pca_details$zir2 <- 0

for(x in 1:nrow(pca_details)){
  if(pca_details$blu[x] > blu_cutoff){pca_details$blu2[x] <- 1}
  if(pca_details$grn[x] > grn_cutoff){pca_details$grn2[x] <- 1}
  if(pca_details$red[x] > red_cutoff){pca_details$red2[x] <- 1}
  if(pca_details$zir[x] > zir_cutoff){pca_details$zir2[x] <- 1}
}

pca_details$bgrn <- paste(pca_details$blu2,pca_details$grn2,pca_details$red2,pca_details$zir2,sep="")

bgrn_table <- data.frame(table(pca_details$bgrn)) %>% arrange(desc(Freq))

pca_details$label <- "Unannotated"

for(x in 1:nrow(pca_details)){
  if(pca_details$bgrn[x] == "1001"){pca_details$label[x] <- "Unrecombined"}
  if(pca_details$bgrn[x] == "0110"){pca_details$label[x] <- "Doubly_recombined"}
  if(pca_details$bgrn[x] == "1010"){pca_details$label[x] <- "GA_only"}
  if(pca_details$bgrn[x] == "0101"){pca_details$label[x] <- "GT_only"}
}

custom_color_scale2 <- c("Unannotated" = "black", "Unrecombined" = "blue", "Doubly_recombined" = "purple", "GA_only" = "red", "GT_only" = "green")

Clone11_PCA_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank()) + 
  scale_color_manual(values = custom_color_scale2) +
  scale_x_continuous(limits = c(-5,5)) +
  scale_y_continuous(limits = c(-3,2)) +
  geom_point(data = pca_details, aes(x = pc1, y = pc2, color = label), alpha = 0.1, size = 0.5) +
  facet_grid(rows = vars(treatment))
Clone11_PCA_plot
ggsave(file = "plots/Clone11_PCA_plot.pdf", Clone11_PCA_plot, height = 3.4, width = 4.5)

```

```{r Individual flow plots}
plot_alpha = 0.02
point_size = 0.25

c11_none_bz_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "none"), aes(x = blu, y = zir), alpha = plot_alpha, size = point_size)
c11_none_bz_plot

c11_recomb_bz_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "recomb"), aes(x = blu, y = zir), alpha = plot_alpha, size = point_size)
c11_recomb_bz_plot

c11_ap1903_bz_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "ap1903"), aes(x = blu, y = zir), alpha = plot_alpha, size = point_size)
c11_ap1903_bz_plot

c11_none_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "none") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "none"), aes(x = grn, y = red), alpha = plot_alpha, size = point_size)
c11_none_gr_plot

c11_recomb_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "recomb") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "recomb"), aes(x = grn, y = red), alpha = plot_alpha, size = point_size)
c11_recomb_gr_plot

c11_ap1903_gr_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank(), legend.position = "ap1903") + 
  scale_x_log10(limits = axis_limits) + scale_y_log10(limits = axis_limits) +
  geom_point(data = pca_details %>% filter(treatment == "ap1903"), aes(x = grn, y = red), alpha = plot_alpha, size = point_size)
c11_ap1903_gr_plot

c11_arranged_plot <- grid.arrange(c11_none_bz_plot, c11_recomb_bz_plot, c11_ap1903_bz_plot, c11_none_gr_plot, c11_recomb_gr_plot, c11_ap1903_gr_plot, ncol=3, nrow=2)

ggsave(file = "plots/c11_arranged_plot.pdf", c11_arranged_plot, height = 2.2, width = 4.1)

#paste("Percent double positive before selection:", round(sum(c11_recomb$grn >= 3e4 & c11_recomb$red >= 3e4) / nrow(c11_recomb) * 100,2))
#paste("Percent double positive after selection:", round(sum(c11_ap1903$grn >= 3e4 & c11_ap1903$red >= 3e4) / nrow(c11_ap1903) * 100,2))

pca_vectorplot <- fviz_pca_var(clone2.pca,
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             )
pca_vectorplot

ggsave(file = "plots/pca_vectorplot.pdf", pca_vectorplot, height = 2, width = 2)
```
## Since reviewer 2 asked, I am now also collecting summary statistics on the double landing pad flow cytometry

```{r}
c2_orig_none <- c2_orig_none %>% mutate(replicate = 1)
c2_orig_recomb <- c2_orig_recomb %>% mutate(replicate = 1)
c2_orig_ap1903 <- c2_orig_ap1903 %>% mutate(replicate = 1)
c6_orig_none <- c6_orig_none %>% mutate(replicate = 1)
c6_orig_recomb <- c6_orig_recomb %>% mutate(replicate = 1)
c6_orig_ap1903 <- c6_orig_ap1903 %>% mutate(replicate = 1)
c11_orig_none <- c11_orig_none %>% mutate(replicate = 1)
c11_orig_recomb <- c11_orig_recomb %>% mutate(replicate = 1)
c11_orig_ap1903 <- c11_orig_ap1903 %>% mutate(replicate = 1)
c2_none <- c2_none %>% mutate(replicate = 2)
c2_recomb <- c2_recomb %>% mutate(replicate = 2)
c2_ap1903 <- c2_ap1903 %>% mutate(replicate = 2)
c6_none <- c6_none %>% mutate(replicate = 2)
c6_recomb <- c6_recomb %>% mutate(replicate = 2)
c6_ap1903 <- c6_ap1903 %>% mutate(replicate = 2)
c11_none <- c11_none %>% mutate(replicate = 2)
c11_recomb <- c11_recomb %>% mutate(replicate = 2)
c11_ap1903 <- c11_ap1903 %>% mutate(replicate = 2)

dlp_analysis <- rbind(c2_orig_none,c2_orig_recomb,c2_orig_ap1903,c6_orig_none,c6_orig_recomb,c6_orig_ap1903,c11_orig_none,c11_orig_recomb,c11_orig_ap1903,c2_none,c2_recomb,c2_ap1903,c6_none,c6_recomb,c6_ap1903,c11_none,c11_recomb,c11_ap1903)

ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = blu))
ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = grn))
ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = red))
ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = nir))
```


```{r}
blu_stat_density_plot <- ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = blu))
blu_stat_density_values <- blu_stat_density_plot %>% 
  ggplot_build() %>% 
  `[[`("data") %>% 
  `[[`(1) %>% select(x,y)

ggplot() + geom_point(data = blu_stat_density_values, aes(x = x, y = y))
blu_stat_density_miny <- as.numeric(blu_stat_density_values %>% filter(x > 2.75 & x < 3.25) %>% summarize(miny = min(y)))
blu_stat_density_minx <- blu_stat_density_values[blu_stat_density_values$y == blu_stat_density_miny,"x"]
ggplot() + geom_point(data = blu_stat_density_values, aes(x = x, y = y)) + geom_vline(xintercept = blu_stat_density_minx)

dlp_analysis_all <- dlp_analysis %>% group_by(clone, treatment, replicate) %>% tally()
dlp_analysis_blupos <- dlp_analysis %>% filter(blu >= 10^blu_stat_density_minx) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_summary <- left_join(dlp_analysis_all,dlp_analysis_blupos,c("clone", "treatment", "replicate"))
colnames(dlp_analysis_summary)[4] <- "all"
colnames(dlp_analysis_summary)[5] <- "blu_pos"
dlp_analysis_summary$blupos_pct <- dlp_analysis_summary$blu_pos/dlp_analysis_summary$all
dlp_analysis_summary$treatment <- factor(dlp_analysis_summary$treatment, levels = c("none","recomb","ap1903"))
DLP_summary_blupos_plot <- ggplot() + labs(x = NULL, y = "Percent of cells", title = "Blue+") + scale_y_continuous(limits = c(0,1)) +
  geom_quasirandom(data = dlp_analysis_summary, aes(x = treatment, y = blupos_pct, color = clone), alpha = 0.5) 
ggsave(file = "plots/DLP_summary_blupos_plot.pdf", DLP_summary_blupos_plot, height = 1.25, width = 2.75)

## We first looked at blue positive cells. Now lets next look at NIR
nir_stat_density_plot <- ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = nir))
nir_stat_density_values <- nir_stat_density_plot %>% 
  ggplot_build() %>% 
  `[[`("data") %>% 
  `[[`(1) %>% select(x,y)

ggplot() + geom_point(data = nir_stat_density_values, aes(x = x, y = y))
nir_stat_density_miny <- as.numeric(nir_stat_density_values %>% filter(x > 3.25 & x < 4) %>% summarize(miny = min(y)))
nir_stat_density_minx <- nir_stat_density_values[nir_stat_density_values$y == nir_stat_density_miny,"x"]
ggplot() + geom_point(data = nir_stat_density_values, aes(x = x, y = y)) + geom_vline(xintercept = nir_stat_density_minx)

dlp_analysis_nirpos <- dlp_analysis %>% filter(nir >= 10^nir_stat_density_minx) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_summary <- left_join(dlp_analysis_summary,dlp_analysis_nirpos,c("clone", "treatment", "replicate"))
colnames(dlp_analysis_summary)[7] <- "nir_pos"
dlp_analysis_summary$nirpos_pct <- dlp_analysis_summary$nir_pos/dlp_analysis_summary$all
DLP_summary_nirpos_plot <- ggplot() + labs(x = NULL, y = "Percent of cells", title = "NIR+") + scale_y_continuous(limits = c(0,1)) +
  geom_quasirandom(data = dlp_analysis_summary, aes(x = treatment, y = nirpos_pct, color = clone), alpha = 0.5)
ggsave(file = "plots/DLP_summary_nirpos_plot.pdf", DLP_summary_nirpos_plot, height = 1.25, width = 2.75)

## Next look at green cells
grn_stat_density_plot <- ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = grn))
grn_stat_density_values <- grn_stat_density_plot %>% 
  ggplot_build() %>% 
  `[[`("data") %>% 
  `[[`(1) %>% select(x,y)

ggplot() + geom_point(data = grn_stat_density_values, aes(x = x, y = y))
grn_stat_density_miny <- as.numeric(grn_stat_density_values %>% filter(x > 2 & x < 4) %>% summarize(miny = min(y)))
grn_stat_density_minx <- grn_stat_density_values[grn_stat_density_values$y == grn_stat_density_miny,"x"]
ggplot() + geom_point(data = grn_stat_density_values, aes(x = x, y = y)) + geom_vline(xintercept = grn_stat_density_minx)

## Next look at green cells
red_stat_density_plot <- ggplot() + scale_x_log10() + geom_density(data = dlp_analysis, aes(x = red))
red_stat_density_values <- red_stat_density_plot %>% 
  ggplot_build() %>% 
  `[[`("data") %>% 
  `[[`(1) %>% select(x,y)

ggplot() + geom_point(data = red_stat_density_values, aes(x = x, y = y))
red_stat_density_miny <- as.numeric(red_stat_density_values %>% filter(x > 2 & x < 4) %>% summarize(miny = min(y)))
red_stat_density_minx <- red_stat_density_values[red_stat_density_values$y == red_stat_density_miny,"x"]
ggplot() + geom_point(data = red_stat_density_values, aes(x = x, y = y)) + geom_vline(xintercept = red_stat_density_minx)

## Now look at green only, red only, and green and red cells
dlp_analysis_grnpos_redneg <- dlp_analysis %>% filter(grn >= 10^grn_stat_density_minx & red <= 10^red_stat_density_minx) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_summary <- left_join(dlp_analysis_summary,dlp_analysis_grnpos_redneg,c("clone", "treatment", "replicate"))
colnames(dlp_analysis_summary)[9] <- "grnpos_redneg"
dlp_analysis_summary$grnpos_redneg_pct <- dlp_analysis_summary$grnpos_redneg/dlp_analysis_summary$all
DLP_summary_grnpos_redneg_plot <- ggplot() + labs(x = NULL, y = "Percent of cells", title = "Green+ Red-") + scale_y_continuous(limits = c(0,1)) +
  geom_quasirandom(data = dlp_analysis_summary, aes(x = treatment, y = grnpos_redneg_pct, color = clone), alpha = 0.5) 
ggsave(file = "plots/DLP_summary_grnpos_redneg_plot.pdf", DLP_summary_grnpos_redneg_plot, height = 1.25, width = 2.75)

dlp_analysis_grnpos_redpos <- dlp_analysis %>% filter(grn >= 10^grn_stat_density_minx & red >= 10^red_stat_density_minx) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_summary <- left_join(dlp_analysis_summary,dlp_analysis_grnpos_redpos,c("clone", "treatment", "replicate"))
colnames(dlp_analysis_summary)[11] <- "grnpos_redpos"
dlp_analysis_summary$grnpos_redpos_pct <- dlp_analysis_summary$grnpos_redpos/dlp_analysis_summary$all
DLP_summary_grnpos_redpos_plot <- ggplot() + labs(x = NULL, y = "Percent of cells", title = "Green+ Red+") + scale_y_continuous(limits = c(0,1)) +
  geom_quasirandom(data = dlp_analysis_summary, aes(x = treatment, y = grnpos_redpos_pct, color = clone), alpha = 0.5) 
ggsave(file = "plots/DLP_summary_grnpos_redpos_plot.pdf", DLP_summary_grnpos_redpos_plot, height = 1.25, width = 2.75)

dlp_analysis_grnneg_redpos <- dlp_analysis %>% filter(grn <= 10^grn_stat_density_minx & red >= 10^red_stat_density_minx) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_summary <- left_join(dlp_analysis_summary,dlp_analysis_grnneg_redpos,c("clone", "treatment", "replicate"))
colnames(dlp_analysis_summary)[13] <- "grnneg_redpos"
dlp_analysis_summary$grnneg_redpos_pct <- dlp_analysis_summary$grnneg_redpos/dlp_analysis_summary$all
DLP_summary_grnneg_redpos_plot <- ggplot() + labs(x = NULL, y = "Percent of cells", title = "Green- Red+") + scale_y_continuous(limits = c(0,1)) +
  geom_quasirandom(data = dlp_analysis_summary, aes(x = treatment, y = grnneg_redpos_pct, color = clone), alpha = 0.5) 
ggsave(file = "plots/DLP_summary_grnneg_redpos_plot.pdf", DLP_summary_grnneg_redpos_plot, height = 1.25, width = 2.75)

## Print out some statistics for writing into the paper
recombined_dlp_only <- dlp_analysis_summary %>% filter(treatment == "recomb")

paste("Mean of singly recombined cells by green-red analysis:", round(mean(c((recombined_dlp_only$grnpos_redneg_pct),(recombined_dlp_only$grnneg_redpos_pct))),3))

paste("Mean of doubly recombined cells by green-red analysis:", round(mean(c((recombined_dlp_only$grnpos_redpos_pct))),3))

selected_dlp_only <- dlp_analysis_summary %>% filter(treatment == "ap1903" & clone != "c6")

paste("Mean of doubly recombined cells after selection by green-red analysis:", round(mean(c((selected_dlp_only$grnpos_redpos_pct))),3))
```

```{r PCA on all datapoints from all replicates from all clones}

dlp_analysis_for_pca <- dlp_analysis[,c("blu","grn","red","nir")]
dlp_analysis_for_pca2 <- cbind(dlp_analysis_for_pca, scale(dlp_analysis_for_pca, clone2.pca$center, clone2.pca$scale) %*% clone2.pca$rotation)
dlp_analysis_for_pca3 <- cbind(dlp_analysis_for_pca2,dlp_analysis[,c("clone", "treatment", "replicate")])

dlp_analysis_for_pca3_c6 <- dlp_analysis_for_pca3 %>% mutate(PC2B = PC2 - 2, label = paste0(clone,"_",replicate), pc1_2b_ratio = PC1/PC2B)
dlp_analysis_for_pca3_c6$treatment <- factor(dlp_analysis_for_pca3_c6$treatment, levels = c("none","recomb","ap1903"))

custom_color_scale3 <- c("c2_1" = "green", "c2_2" = "green", "c6_1" = "blue", "c6_2" = "blue", "c11_1" = "red", "c11_2" = "red")

All_data_PCA_histogram <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank()) + 
  scale_fill_manual(values = custom_color_scale3) +
  scale_x_continuous(limits = c(-1.1,1.5)) +
  geom_density(data = dlp_analysis_for_pca3_c6, aes(x = pc1_2b_ratio, fill = label), alpha = 0.2, size = 0.3) +
  geom_vline(xintercept = -0.2, linetype = 2, alpha = 0.4, color = "red") +
  geom_vline(xintercept = 0.5, linetype = 2, alpha = 0.4, color = "purple") +
  facet_grid(rows = vars(treatment), scales = "free_y")
All_data_PCA_histogram
ggsave(file = "plots/All_data_PCA_histogram.pdf", All_data_PCA_histogram, height = 2.4, width = 3)


dlp_analysis_for_pca3_c6_summary_all <- dlp_analysis_for_pca3_c6 %>% group_by(clone, treatment, replicate) %>% tally()
dlp_analysis_for_pca3_c6_summary_none <- dlp_analysis_for_pca3_c6 %>% filter(pc1_2b_ratio < -0.2) %>% group_by(clone, treatment, replicate) %>% tally()
dlp_analysis_for_pca3_c6_summary_single <- dlp_analysis_for_pca3_c6 %>% filter(pc1_2b_ratio > -0.2 & pc1_2b_ratio < 0.5) %>% group_by(clone, treatment, replicate) %>% tally()
dlp_analysis_for_pca3_c6_summary_double <- dlp_analysis_for_pca3_c6 %>% filter(pc1_2b_ratio > 0.5) %>% group_by(clone, treatment, replicate) %>% tally()

dlp_analysis_for_pca3_c6_summary_combined <- left_join(dlp_analysis_for_pca3_c6_summary_all,dlp_analysis_for_pca3_c6_summary_none, c("clone","treatment","replicate"))
dlp_analysis_for_pca3_c6_summary_combined2 <- left_join(dlp_analysis_for_pca3_c6_summary_combined, dlp_analysis_for_pca3_c6_summary_single, c("clone","treatment","replicate")) 
dlp_analysis_for_pca3_c6_summary_combined3 <- left_join(dlp_analysis_for_pca3_c6_summary_combined2, dlp_analysis_for_pca3_c6_summary_double, c("clone","treatment","replicate")) 

colnames(dlp_analysis_for_pca3_c6_summary_combined3)[4:7] <- c("all","none","single","double")
dlp_analysis_for_pca3_c6_summary_combined3 <- dlp_analysis_for_pca3_c6_summary_combined3 %>% mutate(none_ratio = none / all, single_ratio = single / all, double_ratio = double / all, label = paste0(clone,"_",replicate))

library(reshape2)
dlp_analysis_for_pca3_c6_summary_combined3_melt <- melt(dlp_analysis_for_pca3_c6_summary_combined3[,c("treatment","label","none_ratio","single_ratio","double_ratio")], "id" = c("treatment","label"))

PCA_recombination_plot <- ggplot() + theme_bw() + theme(panel.grid.minor.x = element_blank(), panel.grid.major.x = element_blank()) + 
  scale_color_manual(values = custom_color_scale3) +
  scale_y_continuous(breaks = c(0,0.5,1)) +
  geom_quasirandom(data = dlp_analysis_for_pca3_c6_summary_combined3_melt, aes(x = variable, y = value, color = label), alpha = 0.25) +
  facet_grid(rows = vars(treatment))
PCA_recombination_plot
ggsave(file = "plots/PCA_recombination_plot.pdf", PCA_recombination_plot, height = 2.45, width = 2.6)


## Print out some statistics for writing into the paper
recombined_dlp_pca_only <- dlp_analysis_for_pca3_c6_summary_combined3_melt %>% filter(treatment == "recomb")

paste("Mean of singly recombined cells by pca-based analysis:", round(mean(subset(recombined_dlp_pca_only,variable == "single_ratio")$value),3))

paste("Mean of doubly recombined cells by pca-based analysis:", round(mean(subset(recombined_dlp_pca_only,variable == "double_ratio")$value),3))

selected_dlp_pca_only <- dlp_analysis_for_pca3_c6_summary_combined3_melt %>% filter(treatment == "ap1903" & label != "c6_1" & label != "c6_2")

paste("Mean of doubly recombined cells after selection by pca-based analysis:", round(mean(subset(selected_dlp_pca_only,variable == "double_ratio")$value),3))
```


## Double landing pad cells for biosensors
## GCaMP and STIM1 analysis relevant to figure 4

```{r First import the GCaMP data}
rep1 = read.csv("Data/jgcamp_data/rep1_signal.csv")
rep1 = rep1[-1]
rep2 = read.csv("Data/jgcamp_data/rep2_signal.csv")
rep2 = rep2[-1]

## Shifting times so addition of drug is time zero
rep1$ND.T <- rep1$ND.T - 30
rep2$ND.T <- rep2$ND.T - 30

#separate based on stimulus and variant
rep1_mch_cch = rep1 %>% filter(Name == "mCherry_cch")
rep1_005_cch = rep1 %>% filter(Name == "pNK005_cch")
rep1_013_cch = rep1 %>% filter(Name == "pNK013_cch")

rep1_mch_iono = rep1 %>% filter(Name == "mCherry_ionomycin")
rep1_005_iono = rep1 %>% filter(Name == "pNK005_ionomycin")
rep1_013_iono = rep1 %>% filter(Name == "pNK013_ionomycin")

rep2_mch_cch = rep2 %>% filter(Name == "mCherry_cch")
rep2_005_cch = rep2 %>% filter(Name == "pNK005_cch")
rep2_013_cch = rep2 %>% filter(Name == "pNK013_cch")

rep2_mch_iono = rep2 %>% filter(Name == "mCherry_ionomycin")
rep2_005_iono = rep2 %>% filter(Name == "pNK005_ionomycin")
rep2_013_iono = rep2 %>% filter(Name == "pNK013_ionomycin")

```

Figure 4B - Representative plot of MFI vs time for fluorescence microscopy measurements - add carbachol and use neg control
```{r Representative GcaMP plot with negative control}
rep_4b = ggplot() +
  scale_x_continuous(limits = c(-30,150), breaks = c(-30,0,30,60,90,120,150,180)) + 
  scale_y_continuous(limits = c(0,2500))+
  geom_point(data = rep1_mch_cch, aes(x = ND.T, y = signal, color = "100uM Carbachol"), size =3, alpha = 0.2, color = "black")+
  geom_point(data = rep1_mch_iono, aes(x = ND.T, y = signal, color = "2uM Ionomycin"), size =3, alpha = 0.2, color = "orange")+
  labs(x = "Time (seconds)", y = "Green MFI", color = "Condition") +
  theme(legend.position = "top", panel.grid.major = element_blank())
rep_4b
ggsave("Plots/rep_4b.pdf", rep_4b, width = 3.5, height = 1.75, useDingbats = F)
```

Supplementary Figure 3 - un-normalized fluorescence measurement curves for variants in carbachol and ionomycin
```{r Raw GCaMP curves for supplement}
#fluorescence measurements for cch of each replicate
cch_rep1 = bind_rows(rep1_mch_cch, rep1_005_cch, rep1_013_cch) 
cch_rep1 = cch_rep1 %>% mutate(rep = '1')
cch_rep2 = bind_rows(rep2_mch_cch, rep2_005_cch, rep2_013_cch)
cch_rep2 = cch_rep2 %>% mutate(rep = '2')
cch_reps = bind_rows(cch_rep1, cch_rep2)
cch_reps[cch_reps == "mCherry_cch"] = "Negative control"
cch_reps[cch_reps == "pNK005_cch"] = "WT STIM1"
cch_reps[cch_reps == "pNK013_cch"] = "R429C STIM1"

cch_reps$Name <- factor(cch_reps$Name, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

b = ggplot(data = cch_reps, aes(x = ND.T, y = signal, color = Name))+ theme(panel.grid.major = element_blank(), legend.position = "top") +
  scale_y_continuous(limits = c(0,2500))+ 
  scale_x_continuous(limits= c(-30,150), breaks = c(0,60,120))+
  labs(title = "100uM Carbachol", x = "Time (s)", y = "Green MFI", color = "Variant")+ 
  geom_point(size = 1, alpha = 0.2) +
  facet_wrap(~rep)
b
ggsave("Plots/cch_reps.pdf", b, width = 7, height = 4)


#fluorescence measurements for ionomycin of each replicate
iono_rep1 = bind_rows(rep1_mch_iono, rep1_005_iono, rep1_013_iono) 
iono_rep1 = iono_rep1 %>% mutate(rep = '1')
iono_rep2 = bind_rows(rep2_mch_iono, rep2_005_iono, rep2_013_iono)
iono_rep2 = iono_rep2 %>% mutate(rep = '2')
iono_reps = bind_rows(iono_rep1, iono_rep2)
iono_reps[iono_reps == "mCherry_ionomycin"] = "Negative control"
iono_reps[iono_reps == "pNK005_ionomycin"] = "WT STIM1"
iono_reps[iono_reps == "pNK013_ionomycin"] = "R429C STIM1"

iono_reps$Name <- factor(iono_reps$Name, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

c = ggplot(data = iono_reps, aes(x = ND.T, y = signal, color = Name))+ theme(panel.grid.major = element_blank(), legend.position = "top") +
  scale_y_continuous(limits = c(0,6000))+ scale_x_continuous(limits= c(-30,150), breaks = c(0,60,120))+ 
  labs(title = "2uM Ionomycin", x = "Time (s)", y = "Green MFI")+ 
  geom_point(size = 1, alpha = 0.2) +
  facet_wrap(~rep)
c
ggsave("Plots/iono_reps.pdf", c, width = 7, height = 4)

raw_gcamp_plots <- b | c

raw_gcamp_plots

ggsave(file = "Plots/Raw_gcamp_plots.pdf", raw_gcamp_plots, height = 3, width = 5)
```

Figure 4C - Differences in initial vs peak fluorescence for each variant for carbachol and ionomycin
```{r GCaMP initial and peak comparisons}
#100uM carbachol 
#replicate 1
comp_rep1_cch = data.frame(matrix(ncol = 3, nrow = 3))
colnames(comp_rep1_cch) = c("variant", "Start value", "Peak value")
comp_rep1_cch[1,1] = "Negative control"
comp_rep1_cch[2,1] = "WT STIM1"
comp_rep1_cch[3,1] = "R429C STIM1"
comp_rep1_cch[1,2] = mean(rep1_mch_cch$signal[1:30]) #min(rep1_mch_cch$signal)
comp_rep1_cch[2,2] = mean(rep1_005_cch$signal[1:30]) #min(rep1_005_cch$signal)
comp_rep1_cch[3,2] = mean(rep1_013_cch$signal[1:30]) #min(rep1_013_cch$signal)
comp_rep1_cch[1,3] = quantile(rep1_mch_cch$signal, 0.95) #max(rep1_mch_cch$signal)
comp_rep1_cch[2,3] = quantile(rep1_005_cch$signal, 0.95) #max(rep1_005_cch$signal)
comp_rep1_cch[3,3] = quantile(rep1_013_cch$signal, 0.95) #max(rep1_013_cch$signal)
comp_rep1_cchlg = comp_rep1_cch %>% melt() %>% mutate(replicate = '1')

#replicate 2
comp_rep2_cch = data.frame(matrix(ncol = 3, nrow = 3))
colnames(comp_rep2_cch) = c("variant", "Start value", "Peak value")
comp_rep2_cch[1,1] = "Negative control"
comp_rep2_cch[2,1] = "WT STIM1"
comp_rep2_cch[3,1] = "R429C STIM1"
comp_rep2_cch[1,2] = mean(rep2_mch_cch$signal[1:30]) #min(rep2_mch_cch$signal)
comp_rep2_cch[2,2] = mean(rep2_005_cch$signal[1:30]) #min(rep2_005_cch$signal)
comp_rep2_cch[3,2] = mean(rep2_013_cch$signal[1:30]) #min(rep2_013_cch$signal)
comp_rep2_cch[1,3] = quantile(rep2_mch_cch$signal, 0.95) #max(rep2_mch_cch$signal)
comp_rep2_cch[2,3] = quantile(rep2_005_cch$signal, 0.95) #max(rep2_005_cch$signal)
comp_rep2_cch[3,3] = quantile(rep2_013_cch$signal, 0.95) #max(rep2_013_cch$signal)
comp_rep2_cchlg = comp_rep2_cch %>% melt() %>% mutate(replicate = '2')

comp_cch = bind_rows(comp_rep1_cchlg, comp_rep2_cchlg)

cch_4c = ggplot(data = comp_cch, aes(x =variable, y= value, color = variant))+ 
  geom_point(position=position_dodge(width = 0.25)) + 
  stat_summary(geom = "point", shape = 8,  position = position_dodge(width = 0.25))+
  scale_x_discrete(labels = c("Initial", "Peak"))+ 
  labs(x = "Time point", y = "Green MFI", color = "Variant")+ 
  theme(text=element_text(size=20), legend.position = "none")
cch_4c
ggsave("Plots/cch_4c.pdf", cch_4c, width = 4, height = 3)

#ionomycin
#replicate 1
comp_rep1_iono = data.frame(matrix(ncol = 3, nrow = 3))
colnames(comp_rep1_iono) = c("variant", "Start value", "Peak value")
comp_rep1_iono[1,1] = "Negative control"
comp_rep1_iono[2,1] = "WT STIM1"
comp_rep1_iono[3,1] = "R429C STIM1"
comp_rep1_iono[1,2] = mean(rep1_mch_iono$signal[1:30]) #min(rep1_mch_iono$signal)
comp_rep1_iono[2,2] = mean(rep1_005_iono$signal[1:30]) #min(rep1_005_iono$signal)
comp_rep1_iono[3,2] = mean(rep1_013_iono$signal[1:30]) #min(rep1_013_iono$signal)
comp_rep1_iono[1,3] =  #max(rep1_mch_iono$signal)
comp_rep1_iono[2,3] = quantile(rep1_005_iono$signal, 0.95) #max(rep1_005_iono$signal)
comp_rep1_iono[3,3] = quantile(rep1_013_iono$signal, 0.95) #max(rep1_013_iono$signal)
comp_rep1_ionolg = comp_rep1_iono %>% melt() %>% mutate(replicate = '1')

#replicate 2
comp_rep2_iono = data.frame(matrix(ncol = 3, nrow = 3))
colnames(comp_rep2_iono) = c("variant", "Start value", "Peak value")
comp_rep2_iono[1,1] = "Negative control"
comp_rep2_iono[2,1] = "WT STIM1"
comp_rep2_iono[3,1] = "R429C STIM1"
comp_rep2_iono[1,2] = mean(rep2_mch_iono$signal[1:30]) #min(rep2_mch_iono$signal)
comp_rep2_iono[2,2] = mean(rep2_005_iono$signal[1:30]) #min(rep2_005_iono$signal)
comp_rep2_iono[3,2] = mean(rep2_013_iono$signal[1:30]) #min(rep2_013_iono$signal)
comp_rep2_iono[1,3] = quantile(rep2_mch_iono$signal, 0.95) #max(rep2_mch_iono$signal)
comp_rep2_iono[2,3] = quantile(rep2_005_iono$signal, 0.95) #max(rep2_005_iono$signal)
comp_rep2_iono[3,3] = quantile(rep2_013_iono$signal, 0.95) #max(rep2_013_iono$signal)
comp_rep2_ionolg = comp_rep2_iono %>% melt() %>% mutate(replicate = '2')

comp_iono = bind_rows(comp_rep1_ionolg, comp_rep2_ionolg)

iono_4c = ggplot(data = comp_iono, aes(x =variable, y= value, color = variant))+ scale_y_continuous(limits = c(0,6000))+
 geom_point(position=position_dodge(width = 0.25)) + 
  stat_summary(geom = "point", shape = 8,  position = position_dodge(width = 0.25))+
  scale_x_discrete(labels = c("Initial", "Peak"))+ 
  labs(x = "Time point", y = "Green MFI", color = "Variant")+ 
  theme(text=element_text(size=20), legend.position = "none")
iono_4c
ggsave("Plots/iono_4c.pdf", iono_4c, width = 4, height = 3)


## Another version of the plots, from the way I do it (KAM)

comp_cch$variant <- factor(comp_cch$variant, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

comp_iono$variant <- factor(comp_iono$variant, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

cch4c_version2 <- ggplot() + theme(legend.position = "top", panel.grid.major.x = element_blank(), axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) + 
  labs(title = "Carbachol", x = NULL, y = "Mean fluorescent intensity") +
  scale_y_continuous(limits = c(-200,6000)) +
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_point(data = comp_cch, aes(x = variable, y = value, color = variant), position = position_dodge(width = 0.6)) +
  geom_point(data = comp_cch %>% group_by(variant, variable) %>% summarize(value = mean(value)), aes(x = variable, y = value, color = variant), position = position_dodge(width = 0.6), shape = 95, size = 5)
cch4c_version2

iono4c_version2 <- ggplot() + theme(legend.position = "top", panel.grid.major.x = element_blank(), axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) + 
  labs(title = "Ionomycin", x = NULL, y = "Mean fluorescent intensity") +
  scale_y_continuous(limits = c(-200,6000)) +
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_point(data = comp_iono, aes(x = variable, y = value, color = variant), position = position_dodge(width = 0.6)) +
  geom_point(data = comp_iono %>% group_by(variant, variable) %>% summarize(value = mean(value)), aes(x = variable, y = value, color = variant), position = position_dodge(width = 0.6), shape = 95, size = 5)
#iono4c_version2

peak_figure <- cch4c_version2 | iono4c_version2
peak_figure
ggsave(file = "Plots/Peak_value.pdf", peak_figure, height = 2.5, width = 3)
```

Figure 4D - normalized curves to peak values 
```{r Normalized GCaMP data for plot}
#carbachol
#rep 1
rep1_cch_mch_decay = rep1_mch_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_cch_005_decay = rep1_005_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_cch_013_decay = rep1_013_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_cch_decay_normpeak = bind_rows(rep1_cch_mch_decay, rep1_cch_005_decay, rep1_cch_013_decay)

#rep 2
rep2_cch_mch_decay = rep2_mch_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_cch_005_decay = rep2_005_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_cch_013_decay = rep2_013_cch %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_cch_decay_normpeak = bind_rows(rep2_cch_mch_decay, rep2_cch_005_decay, rep2_cch_013_decay)

rep_normpeak_decay = bind_rows(rep1_cch_decay_normpeak, rep2_cch_decay_normpeak)
rep_normpeak_decay[rep_normpeak_decay == "mCherry_cch"] = "Negative control"
rep_normpeak_decay[rep_normpeak_decay == "pNK005_cch"] = "WT STIM1"
rep_normpeak_decay[rep_normpeak_decay == "pNK013_cch"] = "R429C STIM1"

cch_normpeak_4d = ggplot(data = rep_normpeak_decay, aes(x = ND.T, y = norm_peak, color = Name)) +
  theme_bw() + 
  scale_y_continuous(breaks = c(0,0.5,1)) + 
  scale_x_continuous(breaks = c(-30,0,30,60,90,120,150))+ 
  labs(x = "Time (s)", y = "Green MFI", color = "Variant")+ theme(text = element_text(size = 20)) + 
  geom_line(size = 1) +
  facet_wrap(~rep)
cch_normpeak_4d

## Kenny mods

rep_normpeak_decay2 <- rep_normpeak_decay %>% group_by(Name, ND.T) %>% summarize(mean_normpeak_decay = mean(norm_peak))

rep_normpeak_decay2$Name <- factor(rep_normpeak_decay2$Name, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

cch_normpeak_average_decay <- ggplot() + theme(legend.position = "top") +
  labs(x = "Time (seconds)", y = "Normalized fluorescence", title = "Carbachol") +
  scale_x_continuous(breaks = c(-30,0,30,60,90,120,150))+ 
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_point(data = rep_normpeak_decay2, aes(x = ND.T, y = mean_normpeak_decay, color = Name), alpha = 0.2)
cch_normpeak_average_decay


#ionomycin
#rep 1
rep1_iono_mch_decay = rep1_mch_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_iono_005_decay = rep1_005_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_iono_013_decay = rep1_013_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep1_iono_decay_normpeak = bind_rows(rep1_iono_mch_decay, rep1_iono_005_decay, rep1_iono_013_decay)

#rep 2
rep2_iono_mch_decay = rep2_mch_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_iono_005_decay = rep2_005_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_iono_013_decay = rep2_013_iono %>% mutate(norm_peak = signal/quantile(signal, 0.95))
rep2_iono_decay_normpeak = bind_rows(rep2_iono_mch_decay, rep2_iono_005_decay, rep2_iono_013_decay)

rep_normpeak_decay_iono = bind_rows(rep1_iono_decay_normpeak, rep2_iono_decay_normpeak)
rep_normpeak_decay_iono[rep_normpeak_decay_iono$Name == "mCherry_ionomycin","Name"] = "Negative control"
rep_normpeak_decay_iono[rep_normpeak_decay_iono$Name == "pNK005_ionomycin","Name"] = "WT STIM1"
rep_normpeak_decay_iono[rep_normpeak_decay_iono$Name == "pNK013_ionomycin","Name"] = "R429C STIM1"

iono_normpeak_4d = ggplot(data = rep_normpeak_decay_iono, aes(x = ND.T, y = norm_peak, color = Name))+theme_bw()+ scale_y_continuous(breaks = c(0,0.5,1))+ 
  scale_x_continuous(breaks = c(-30,0,30,60,90,120,150))+ 
  labs(x = "Time (s)", y = "Green MFI", color = "Variant")+ theme(text = element_text(size = 20)) + 
  geom_line(size = 1) +
  facet_wrap(~rep)
iono_normpeak_4d

## Kenny mods
rep_normpeak_decay_iono2 <- rep_normpeak_decay_iono %>% group_by(Name, ND.T) %>% summarize(mean_normpeak_decay = mean(norm_peak))
rep_normpeak_decay_iono2$Name <- factor(rep_normpeak_decay_iono2$Name, levels = c("Negative control", "WT STIM1", "R429C STIM1"))

iono_normpeak_average_decay <- ggplot() + theme(legend.position = "top") +
  labs(x = "Time (seconds)", y = "Normalized fluorescence", title = "Ionomycin") +
  scale_x_continuous(breaks = c(-30,0,30,60,90,120,150))+ 
  geom_hline(yintercept = 0, alpha = 0.5) + 
  geom_point(data = rep_normpeak_decay_iono2, aes(x = ND.T, y = mean_normpeak_decay, color = Name), alpha = 0.2)
iono_normpeak_average_decay

# ggsave("cch_normpeak_4d.pdf", cch_normpeak_4d, width = 7, height = 4)
# ggsave("iono_normpeak_4d.pdf", iono_normpeak_4d, width = 7, height = 4)

grid.arrange(cch_normpeak_4d, iono_normpeak_4d)

Decay_curves <- cch_normpeak_average_decay | iono_normpeak_average_decay
Decay_curves
ggsave(file = "Plots/Decay_curves.pdf", Decay_curves, height = 2.5, width = 5)
```

## Double landing pad with protein interaction pairs
## Relevant to Figure 5 Vif <-> Apobec experiments

```{r Import the vif orientations data}

vif_apobec <- read.csv(file = "data/flow/Vif_Apobec/Orientations.csv")
vif_apobec_fusions <- vif_apobec %>% filter(gene != "None")
vif_apobec_egfp <- vif_apobec %>% filter(gene == "None")


Apobec_fusions_plot <- ggplot() + theme_bw() + theme(panel.grid.minor = element_blank()) +
  labs(x = "Geometric mean of MFI (Apobec only)", y = "Geometric mean of MFI\n(HIV-1 vif coexpressed)") +
  scale_x_log10() + scale_y_log10() +
  geom_abline(slope = 1, alpha = 0.2, linetype = 2, alpha = 0.5, size = 1.5) +
  geom_point(data = vif_apobec_fusions, aes(x = stop, y = wt, color = gene, shape = orientation), size = 3, alpha = 0.7) +
  geom_point(data = vif_apobec_egfp, aes(x = stop, y = wt), size = 3, shape = 15, alpha = 0.7)
Apobec_fusions_plot
ggsave(file = "plots/Apobec_fusions_plot.pdf", Apobec_fusions_plot, height = 2.2, width = 4)

```

```{r Vif Apobec heatmap}
va_rep1 <- read.csv(file = "Data/Vif_apobec/Vif_apobec_rep1.csv", header = T)
for(x in 1:nrow(va_rep1)){
  for(y in 2:ncol(va_rep1)){
    va_rep1[x,y] <- va_rep1[x,y] / va_rep1[6,y]
  }
}
va_rep1_melt <- melt(va_rep1) %>% filter(variable != "EGFP.A3F")
colnames(va_rep1_melt) <- c("vif","apobec","normalized")
va_rep1_melt$rep <- 1

va_rep2 <- read.csv(file = "Data/Vif_apobec/Vif_apobec_rep2.csv", header = T)
for(x in 1:nrow(va_rep2)){
  for(y in 2:ncol(va_rep2)){
    va_rep2[x,y] <- va_rep2[x,y] / va_rep2[8,y]
  }
}
va_rep2_melt <- melt(va_rep2)
colnames(va_rep2_melt) <- c("vif","apobec","normalized")
va_rep2_melt$rep <- 2

va_rep3 <- read.csv(file = "Data/Vif_apobec/Vif_apobec_rep3.csv", header = T)
for(x in 1:nrow(va_rep3)){
  for(y in 2:ncol(va_rep3)){
    va_rep3[x,y] <- va_rep3[x,y] / va_rep3[8,y]
  }
}
va_rep3_melt <- melt(va_rep3)
colnames(va_rep3_melt) <- c("vif","apobec","normalized")
va_rep3_melt$rep <- 2

va_rep_combined <- rbind(va_rep1_melt, va_rep2_melt, va_rep3_melt)

va_rep_combined_summary <- va_rep_combined %>% group_by(vif, apobec) %>% summarize(mean = mean(normalized), sd = sd(normalized)) %>% mutate(mean_log2 = log2(mean), sd_log2 = log2(sd))

va_rep_combined_summary$apobec <- factor(va_rep_combined_summary$apobec, levels = c("EGFP","A3G.EGFP","EGFP.A3D","EGFP.A3H2"))
va_rep_combined_summary$vif <- factor(va_rep_combined_summary$vif, levels = rev(c("WT","STOP","W89R","A149Y","W11A","H43E","F39A","G84K")))

vif_apobec_heatmap1 <- ggplot() + theme(axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "EGFP or Apobec fusion", y = "HIV-1 Vif variant") +
  #scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "brown4", limits = c(0.25,1.3), midpoint = 1) + 
  scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "green4", limits = c(log2(0.25),log2(1.3)), midpoint = 1) + 
  geom_tile(data = va_rep1_melt, aes(x = apobec, y = vif, fill = log2(normalized)))
vif_apobec_heatmap1
ggsave(file = "Plots/Vif_apobec_heatmap1.pdf", vif_apobec_heatmap1, height = 2, width = 3.5)

vif_apobec_heatmap2 <- ggplot() + theme(axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "EGFP or Apobec fusion", y = "HIV-1 Vif variant") +
  #scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "brown4", limits = c(0.25,1.3), midpoint = 1) + 
  scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "green4", limits = c(log2(0.18),log2(1.32)), midpoint = 1) + 
  geom_tile(data = va_rep2_melt, aes(x = apobec, y = vif, fill = log2(normalized)))
vif_apobec_heatmap2
ggsave(file = "Plots/Vif_apobec_heatmap2.pdf", vif_apobec_heatmap2, height = 2, width = 3.5)

vif_apobec_heatmap3 <- ggplot() + theme(axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "EGFP or Apobec fusion", y = "HIV-1 Vif variant") +
  #scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "brown4", limits = c(0.25,1.3), midpoint = 1) + 
  scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "green4", limits = c(log2(0.2),log2(1.62)), midpoint = 1) + 
  geom_tile(data = va_rep3_melt, aes(x = apobec, y = vif, fill = log2(normalized)))
vif_apobec_heatmap3
ggsave(file = "Plots/Vif_apobec_heatmap3.pdf", vif_apobec_heatmap3, height = 2, width = 3.5)

vif_apobec_heatmap_combined <- ggplot() + theme(axis.text.x.bottom = element_text(angle = 45, hjust = 1, vjust = 1)) +
  labs(x = "EGFP or Apobec fusion", y = "HIV-1 Vif variant") +
  #scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "brown4", limits = c(0.25,1.3), midpoint = 1) + 
  scale_fill_gradient2(low = "white", mid = "darkgreen", high = "red", na.value = "green4", limits = c(-2.15,0.4), midpoint = 1) + 
  geom_tile(data = va_rep_combined_summary, aes(x = apobec, y = vif, fill = mean_log2))
vif_apobec_heatmap_combined
ggsave(file = "Plots/vif_apobec_heatmap_combined.pdf", vif_apobec_heatmap_combined, height = 2, width = 3)
```

```{r Individual Vif-Apobec data for supplementary figure}

va_rep_combined$vif <- factor(va_rep_combined$vif, levels = c("WT","STOP","W89R","A149Y","W11A","H43E","F39A","G84K"))
va_rep_combined_summary$vif <- factor(va_rep_combined_summary$vif, levels = c("WT","STOP","W89R","A149Y","W11A","H43E","F39A","G84K"))

EGFP_plot <- ggplot() + theme(panel.grid.major.x = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  scale_y_log10(limits = c(0.18,1.62), breaks = c(0.25,0.5,1)) + labs(title = "EGFP", x = NULL) +
  geom_hline(yintercept = 1, alpha = 0.4) +
  geom_point(data = va_rep_combined %>% filter(apobec == "EGFP"), aes(x = vif, y = normalized), alpha = 0.4, color = "red") +
  geom_point(data = va_rep_combined_summary %>% filter(apobec == "EGFP"), aes(x = vif, y = mean), shape = 95, size = 6)
EGFP_plot

A3G.EGFP_plot <- ggplot() + theme(panel.grid.major.x = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  scale_y_log10(limits = c(0.18,1.62), breaks = c(0.25,0.5,1)) + labs(title = "A3G-EGFP", x = NULL) +
  geom_hline(yintercept = 1, alpha = 0.4) +
  geom_point(data = va_rep_combined %>% filter(apobec == "A3G.EGFP"), aes(x = vif, y = normalized), alpha = 0.4, color = "red") +
  geom_point(data = va_rep_combined_summary %>% filter(apobec == "A3G.EGFP"), aes(x = vif, y = mean), shape = 95, size = 6)
A3G.EGFP_plot

EGFP.A3D_plot <- ggplot() + theme(panel.grid.major.x = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  scale_y_log10(limits = c(0.18,1.62), breaks = c(0.25,0.5,1)) + labs(title = "EGFP-A3D", x = NULL) +
  geom_hline(yintercept = 1, alpha = 0.4) +
  geom_point(data = va_rep_combined %>% filter(apobec == "EGFP.A3D"), aes(x = vif, y = normalized), alpha = 0.4, color = "red") +
  geom_point(data = va_rep_combined_summary %>% filter(apobec == "EGFP.A3D"), aes(x = vif, y = mean), shape = 95, size = 6)
EGFP.A3D_plot

EGFP.A3H2_plot <- ggplot() + theme(panel.grid.major.x = element_blank(), plot.title = element_text(hjust = 0.5), axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) + 
  scale_y_log10(limits = c(0.18,1.62), breaks = c(0.25,0.5,1)) + labs(title = "EGFP-A3H2", x = NULL) +
  geom_hline(yintercept = 1, alpha = 0.4) +
  geom_point(data = va_rep_combined %>% filter(apobec == "EGFP.A3H2"), aes(x = vif, y = normalized), alpha = 0.4, color = "red") +
  geom_point(data = va_rep_combined_summary %>% filter(apobec == "EGFP.A3H2"), aes(x = vif, y = mean), shape = 95, size = 6)
EGFP.A3H2_plot

Apobec_vif_supp_fig <- EGFP_plot | A3G.EGFP_plot | EGFP.A3D_plot | EGFP.A3H2_plot
Apobec_vif_supp_fig
ggsave(file = "Plots/Apobec_vif_supp_fig.pdf", Apobec_vif_supp_fig, height = 2, width = 8)
```

## Testing additional GN dinucleotide pairs in combination
## Relevant to Figure 6

```{r Importing GN sequencing data}

attB_plasmid_mixture <- read.delim(file = "data/GN_recombinations/attBmix_outfile.tsv", header = T, stringsAsFactors = F)
attB_plasmid_mixture_table <- data.frame(table(attB_plasmid_mixture$attB))
attB_plasmid_mixture_table$Freq2 <- attB_plasmid_mixture_table$Freq / sum(attB_plasmid_mixture_table$Freq)
colnames(attB_plasmid_mixture_table) <- c("attB","start_count","start_freq")

## Rep 1
rep1_ga <- read.delim(file = "data/GN_recombinations/Rep1_GA_outfile.tsv", header = T, stringsAsFactors = F)
rep1_ga2 <- data.frame(table(rep1_ga$attB)) %>% mutate(attP = "GA")
rep1_gt <- read.delim(file = "data/GN_recombinations/Rep1_GT_outfile.tsv", header = T, stringsAsFactors = F)
rep1_gt2 <- data.frame(table(rep1_gt$attB)) %>% mutate(attP = "GT")
rep1_gc <- read.delim(file = "data/GN_recombinations/Rep1_GC_outfile.tsv", header = T, stringsAsFactors = F)
rep1_gc2 <- data.frame(table(rep1_gc$attB)) %>% mutate(attP = "GC")
rep1_gg <- read.delim(file = "data/GN_recombinations/Rep1_GG_outfile.tsv", header = T, stringsAsFactors = F)
rep1_gg2 <- data.frame(table(rep1_gg$attB)) %>% mutate(attP = "GG")

rep1_comb <- rbind(rep1_ga2, rep1_gt2, rep1_gc2, rep1_gg2)
colnames(rep1_comb) <- c("attB","Freq","attP")

gn_rep1_results_frame <- data.frame("attB" = NA, "attP" = NA, "final_freq" = NA, "start_freq" = NA)
for( x in c("GT","GA","GC","GG")){
  temp_subset <- subset(rep1_comb, attP == x)
  temp_subset$final_freq <- temp_subset$Freq / sum(temp_subset$Freq)
  temp_subset2 <- merge(temp_subset[,c("attP","attB","final_freq")], attB_plasmid_mixture_table[,c("attB","start_freq")])
  gn_rep1_results_frame <- rbind(gn_rep1_results_frame, temp_subset2)
}
gn_rep1_results_frame <- gn_rep1_results_frame %>% filter(!is.na(attB)) %>% mutate(rep = 1)

## Rep 2
gn_rep2 <- read.delim(file = "data/GN_recombinations/Rep2_outfile.tsv", header = T, stringsAsFactors = F)
gn_rep2_table <- data.frame(with(gn_rep2,table(attP,attB)))

gn_rep2_results_frame <- data.frame("attB" = NA, "attP" = NA, "final_freq" = NA, "start_freq" = NA)
for( x in c("GT","GA","GC","GG")){
  temp_subset <- subset(gn_rep2_table, attP == x)
  temp_subset$final_freq <- temp_subset$Freq / sum(temp_subset$Freq)
  temp_subset2 <- merge(temp_subset[,c("attP","attB","final_freq")], attB_plasmid_mixture_table[,c("attB","start_freq")])
  gn_rep2_results_frame <- rbind(gn_rep2_results_frame, temp_subset2)
}
gn_rep2_results_frame <- gn_rep2_results_frame %>% filter(!is.na(attB)) %>% mutate(rep = 2)

## Rep 3
gn_rep3 <- read.delim(file = "data/GN_recombinations/Rep3_outfile.tsv", header = T, stringsAsFactors = F)
gn_rep3_table <- data.frame(with(gn_rep3,table(attP,attB)))

gn_rep3_results_frame <- data.frame("attB" = NA, "attP" = NA, "final_freq" = NA, "start_freq" = NA)
for( x in c("GT","GA","GC","GG")){
  temp_subset <- subset(gn_rep3_table, attP == x)
  temp_subset$final_freq <- temp_subset$Freq / sum(temp_subset$Freq)
  temp_subset2 <- merge(temp_subset[,c("attP","attB","final_freq")], attB_plasmid_mixture_table[,c("attB","start_freq")])
  gn_rep3_results_frame <- rbind(gn_rep3_results_frame, temp_subset2)
}
gn_rep3_results_frame <- gn_rep3_results_frame %>% filter(!is.na(attB)) %>% mutate(rep = 3)

## Combine the results of all three replicates

gn_combined <- rbind(gn_rep1_results_frame, gn_rep2_results_frame, gn_rep3_results_frame)

gn_combined$enrichment <- gn_combined$final_freq/gn_combined$start_freq

write.csv(file = "data/GN_recombination_results.csv", gn_combined, row.names = FALSE)

gn_combined_summary <- gn_combined %>% group_by(attP, attB) %>% summarize(mean = mean(final_freq), sd = sd(final_freq))

gn_combined_summary2 <- rbind(gn_combined_summary, data.frame("attP" = "Input", "attB" = attB_plasmid_mixture_table$attB, "mean" = attB_plasmid_mixture_table$start_freq, "sd" = NA))


```


```{r Recombination mixture experiment plot}
gn_combined$attP <- factor(gn_combined$attP, levels = c("GT","GA","GC","GG"))
gn_combined_summary2$attP <- factor(gn_combined_summary2$attP, levels = c("Input", "GT","GA","GC","GG"))

Orthogonal_plasmid_recombinations <- ggplot() + scale_y_log10() + theme(panel.grid.major.x = element_blank()) +
  labs(x = "AttP plasmid (individually tested)", y = "Fraction of\nrecombinants", color = "AttB\nplasmid\n(mixed)") + 
  geom_point(data = gn_combined_summary2, aes(x = attP, y = mean, color = attB), position = position_dodge(width = 0.6), shape = 95, size = 4) +
  geom_point(data = gn_combined, aes(x = attP, y = final_freq, color = attB), position = position_dodge(width = 0.6), alpha = 0.4) +
  #geom_errorbar(data = gn_combined_summary2, aes(x = attP, ymin = mean - sd/sqrt(2)*1.96, ymax = mean + sd/sqrt(3)*1.96, color = attB), position = position_dodge(width = 0.6), width = 0.5, alpha = 0.5) +
  geom_point(data = gn_combined_summary2, aes(x = attP, y = mean, color = attB), position = position_dodge(width = 0.6), shape = 95, size = 4) +
  NULL
ggsave(file = "plots/Orthogonal_plasmid_recombinations.pdf", Orthogonal_plasmid_recombinations, height = 2, width = 4.2)
Orthogonal_plasmid_recombinations

Orthogonal_plasmid_recombinations_linear <- ggplot() + theme(panel.grid.major.x = element_blank()) +
  labs(x = "AttP plasmid (individually tested)", y = "Fraction of\nrecombinants", color = "AttB\nplasmid\n(mixed)") + 
  geom_point(data = gn_combined_summary2, aes(x = attP, y = mean, color = attB), position = position_dodge(width = 0.6), shape = 95, size = 4) +
  geom_point(data = gn_combined, aes(x = attP, y = final_freq, color = attB), position = position_dodge(width = 0.6), alpha = 0.4) +
  #geom_errorbar(data = gn_combined_summary2, aes(x = attP, ymin = mean - sd/sqrt(2)*1.96, ymax = mean + sd/sqrt(3)*1.96, color = attB), position = position_dodge(width = 0.6), width = 0.5, alpha = 0.5) +
  geom_point(data = gn_combined_summary2, aes(x = attP, y = mean, color = attB), position = position_dodge(width = 0.6), shape = 95, size = 4) +
  NULL
ggsave(file = "plots/Orthogonal_plasmid_recombinations_linear.pdf", Orthogonal_plasmid_recombinations_linear, height = 2, width = 4.2)
Orthogonal_plasmid_recombinations_linear
```

```{r Modeling predicted pairwise specificities based on the GN data}
recomb_mixture_test <- gn_combined_summary[1:16,c("attP","attB")]
colnames(recomb_mixture_test) <- c("lp1","lp2")

simulation_test_size = 100000
recomb_mixture_test$on_target <- NA

for(x in 1:nrow(recomb_mixture_test)){
  recomb_mix_lp1 <- recomb_mixture_test$lp1[x]
  recomb_mix_lp2 <- recomb_mixture_test$lp2[x]
  if(recomb_mix_lp1 == recomb_mix_lp2){
    lp1_sampling <- sample(x = c(1,0), prob = c(0.5,0.5) , size = simulation_test_size, replace = T)
    lp2_sampling <- sample(x = c(1,0), prob = c(0.5,0.5) , size = simulation_test_size, replace = T)
    results_frame <- data.frame("lp1_on_target" = lp1_sampling == 1, "lp2_on_target" = lp2_sampling == 1)
    results_frame$sum <- results_frame$lp1_on_target + results_frame$lp2_on_target
    recomb_mixture_test$on_target[x] <- sum(results_frame$sum == 1)
  }
  if(recomb_mix_lp1 != recomb_mix_lp2){
    lp1_probs <- gn_combined_summary %>% filter(attP == recomb_mix_lp1 & attB %in% c(recomb_mix_lp1, recomb_mix_lp2))
    lp1_probs$norm <- lp1_probs$mean/sum(lp1_probs$mean)
    lp2_probs <- gn_combined_summary %>% filter(attP == recomb_mix_lp2 & attB %in% c(recomb_mix_lp1, recomb_mix_lp2))
    lp2_probs$norm <- lp2_probs$mean/sum(lp2_probs$mean)
    lp1_sampling <- sample(x = lp1_probs$attB, prob = lp1_probs$norm , size = simulation_test_size, replace = T)
    lp2_sampling <- sample(x = lp2_probs$attB, prob = lp2_probs$norm , size = simulation_test_size, replace = T)
    results_frame <- data.frame("lp1_on_target" = lp1_sampling == recomb_mix_lp1, "lp2_on_target" = lp2_sampling == recomb_mix_lp2)
    results_frame$sum <- results_frame$lp1_on_target + results_frame$lp2_on_target
    recomb_mixture_test$on_target[x] <- sum(results_frame$sum == 2) + sum(results_frame$sum == 0)
  }
}
recomb_mixture_test$frac_on_target <- recomb_mixture_test$on_target / simulation_test_size

Double_LP_heatmap <- ggplot() + labs(x = "First landing pad", y = "Second landing pad", fill = "Fraction\nwith unique\nrecombinant\nplasmid pair") + 
  scale_fill_gradient2(low = "white", high = "red", midpoint = 0.8) + 
  geom_tile(data = recomb_mixture_test, aes(x = lp1, y = lp2, fill = frac_on_target)) +
  geom_tile(data = recomb_mixture_test %>% filter(lp1 == "GC" | lp2 == "GC"), aes(x = lp1, y = lp2, fill = frac_on_target), color = "blue", alpha = 0, size = 0.5) +
  geom_text(data = recomb_mixture_test, aes(x = lp1, y = lp2, label = round(frac_on_target,2)), size = 3)
ggsave(file = "plots/Double_LP_heatmap.pdf", Double_LP_heatmap, height = 2, width = 3.3)
Double_LP_heatmap
```