---
title: "R Notebook"
output: html_notebook
---


```{r}
# This clears all stored variables from memory
rm(list = ls())


# Library imports
# If any of the libraries are missing, run the following command to install them: install.packages('')
library(ggplot2)
library(ggh4x, include.only = 'facet_nested')
library(Hmisc, include.only = 'binconf')
library(ggbeeswarm, include.only = 'geom_quasirandom')
library(quantreg)


# This function creates a PDF file of a ggplot and saves it to a specified directory
# we create to version of each, one with the legend overlaid, and one without
PDF_me <- function(Fig, name, directory, x, y) {
  # Create a PDF file with the specified file name, width, and height
  cairo_pdf(
    filename  = paste(paste(directory, '/', sep = ''), name, ".pdf", sep = ''),
    width = x,
    height = y
  )
  
  # Create the plot with R code
  # We position the legend over the graph so it doesn't affect the size
  # This needs to be manually moved in software like inkscape
  print(Fig + theme(legend.position = c(0.5, 0.5)))
  
  # Run dev.off() to create the file!
  dev.off()
  
  
  # Create a PDF file with the specified file name, width, and height
  # In this one we removed the legend
  cairo_pdf(
    filename  = paste(
      paste(directory, '/', sep = ''),
      paste(name, '_no_legend'),
      ".pdf",
      sep = ''
    ),
    width = x,
    height = y
  )
  
  # Create the plot with R code
  print(Fig + theme(legend.position = "none"))
  
  # Run dev.off() to create the file!
  dev.off()
  
  
  return()
}


# This adds some standard aesthetics to ensure consistency between graphs
S_A <- function(ggplot_Element) {
  # Aesthetics common to all graphs
  ggplot_Element <- ggplot_Element +
    scale_color_brewer(palette = "Accent") +
    theme_linedraw() +
    scale_shape_manual(values = c(0, 1)) +
    scale_y_continuous(
      labels = function(x)
        paste0(x * 100, ""),
      limits = c(0, 1),
      breaks = c(0, 0.25, 0.5, 0.75, 1)
    ) +
    scale_size_continuous(
      limits = c(1, 300),
      breaks = c(1, 75, 150, 300),
      range = c(.000000001, 3),
      name = expression('Number of I'[2] * '\n larvae scored')
    ) +
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(), text = element_text(size = 12))
  
  return(ggplot_Element)
}

# A conversion list so match symbols to sex
Sex_Symbol <- list('\u2642', '\u2640') # mars and venus symbols
names(Sex_Symbol) <- c('Male', 'Female')

# Setting the working directory, change for your computer
# This is for the folder where the csv file can be found
wd <- paste(getwd(), "/", sep = "")
```



```{r}


# Here we ingest the data and calculate outcomes per parent and per cross (pool)

# Read in the specified csv file
raw_data = read.csv(file = "Cross_Data.csv")

# Later we will produce a graph that shows the somatic phenotype rate for
# combinations of gRNA+/- and Cas9+/- inheritance
# "Cas9_gRNA_Pheno_Reps" will hold the number of progeny in that group, as opposed to 'Reps' which holds the overall number of progeny from that parent/pool
raw_data$Cas9_gRNA_Pheno_Reps <- NA

# Convert the Cross and Parent columns to factors
raw_data$Cross <- as.factor(raw_data$Cross)
raw_data$Parent <- as.factor(raw_data$Parent)

# Initialize an empty data frame to store summary statistics
Sum_Stats <- raw_data[0,]

# Here we loop over all the individual crosses in the database
# i is every unique cross
for (i in unique(raw_data$Cross)) {
  # We add 'pool' as a unique parent, using it to sum the outcomes of all parent together
  for (ii in unique(c(as.character(raw_data[raw_data$Cross == i,]$Parent), 'Pool'))) {
    # If the parent is 'Pool' we sum over all parents
    if (ii == 'Pool') {
      current_cross_and_parent <- raw_data[raw_data$Cross == i,]
      
      current_cross_and_parent$Parent <- 'Pool'
      
    } else {
      # If the parent is not 'Pool', we create an entry for the individual parent
      current_cross_and_parent <-
        raw_data[raw_data$Cross == i & raw_data$Parent == ii,]
    }
    
    # Here we check for errors in the database (values that should be the same that are not)
    # We check that withing the cross anything not an F2 etor/parent is the same
    if (length(unique(current_cross_and_parent[, names(current_cross_and_parent)[!names(current_cross_and_parent) %in%
                                                                                 c('F2.Sex',
                                                                                   'F2.Som_Pheno',
                                                                                   'F2.gRNA',
                                                                                   'F2.Nuclease',
                                                                                   'Parent',
                                                                                   'Reps')]])[, 1]) > 1) {
      print(unique(current_cross_and_parent[, names(current_cross_and_parent)[!names(current_cross_and_parent) %in%
                                                                                c('F2.Sex',
                                                                                  'F2.Som_Pheno',
                                                                                  'F2.gRNA',
                                                                                  'F2.Nuclease',
                                                                                  'Parent',
                                                                                  'Reps')]]))
      
      stop('Factors within a single cross that should be the same are not (most anything not F2)')
    }
    
    # Initialize a data frame for summary statistics for the current cross and parent
    # This will be the summary, so we only copy a single row
    Sum_Stats_Current <- current_cross_and_parent[1,]
    
    # Overall A inheritance
    # Calculate the number of positive reps (F2s with gRNA == 'Yes')
    Sum_Stats_Current$Pos_Reps <-
      sum(current_cross_and_parent[current_cross_and_parent$F2.gRNA == 'Yes',]$Reps)
    
    
    # Overall F2s in this cross
    Sum_Stats_Current$Reps <- sum(current_cross_and_parent$Reps)
    
    #Overall A inheritance, ignoring individual crosses
    Sum_Stats_Current$Pos_Frac <-
      Sum_Stats_Current$Pos_Reps / Sum_Stats_Current$Reps
    
    # Here we list the status of the other factors
    Sum_Stats_Current$F2.Som_Pheno <- 'Either'
    Sum_Stats_Current$Outcome <- 'gRNA_gene'
    Sum_Stats_Current$F2.Nuclease <- 'Either'
    Sum_Stats_Current$F2.gRNA <- 'Yes'
    
    # If there are any reps, add the current summary statistics to the main data frame
    if (!is.na(sum(Sum_Stats_Current$Reps))) {
      Sum_Stats <- rbind(Sum_Stats, Sum_Stats_Current)
    }
    
    # We do the same, but for overall non-WT somatic phenotype
    
    # Initialize a data frame for summary statistics for the current cross and parent
    # This will be the summary, so we only copy a single row
    Sum_Stats_Current <- current_cross_and_parent[1,]
    
    # Overall no-WT phenotype
    Sum_Stats_Current$Pos_Reps <-
      sum(current_cross_and_parent[current_cross_and_parent$F2.Som_Pheno != 'WT',]$Reps)
    
    
    # Overall F2s in this cross
    Sum_Stats_Current$Reps <- sum(current_cross_and_parent$Reps)
    
    #Overall phenotype rate, ignoring individual crosses
    Sum_Stats_Current$Pos_Frac <-
      Sum_Stats_Current$Pos_Reps / Sum_Stats_Current$Reps
    
    
    # Here we list the status of the other factors
    Sum_Stats_Current$F2.Som_Pheno <- '!non WT/Not Scored'
    Sum_Stats_Current$Outcome <- 'nonWT_Pheno'
    Sum_Stats_Current$F2.Nuclease <- 'Either'
    Sum_Stats_Current$F2.gRNA <- 'Either'
    
    
    # If there are any reps, add the current summary statistics to the main data frame
    if (!is.na(sum(Sum_Stats_Current$Reps))) {
      Sum_Stats <- rbind(Sum_Stats, Sum_Stats_Current)
    }
    
    # Here we create summary statistics for the
    # different combinations of F2s gRNA/Cas9 status
    # We will create a separate row for each combination of factors
    
    # gRNA inheritance status
    D_status_options <- c('Yes', 'Yes', 'No', 'No')
    # Nuclease inheritance status
    N_status_options <- c('Yes', 'No', 'Yes', 'No')
    
    # Here we loop over the three phenotype options
    Pheno_status_options <- c('WT', 'Mosaic/WE', 'Not Scored')
    
    # Here we loop over the 3 combinations of phenotypes
    # And the 4 combinations of gRNA+/- and Cas9+/-
    for (p in 1:length(Pheno_status_options)) {
      for (f in 1:length(N_status_options)) {
        # Initialize a data frame for summary statistics for the current cross and parent
        # This will be the summary, so we only copy a single row
        Sum_Stats_Current <- current_cross_and_parent[1,]
        
        # We record the outcome as:
        Sum_Stats_Current$Outcome <- 'Cas9_gRNA_Pheno'
        
        # The phenotype is set to the current one we are looping over
        Sum_Stats_Current$F2.Som_Pheno <- Pheno_status_options[p]
        
        # The gRNA status is set to the current one we are looping over
        Sum_Stats_Current$F2.gRNA <- D_status_options[f]
        
        # The Cas9 status is set to the current one we are looping over
        Sum_Stats_Current$F2.Nuclease <- N_status_options[f]
        
        
        # Here we record the overall number of reps in this cross and parent
        Sum_Stats_Current$Reps <-
          sum(current_cross_and_parent$Reps)
        
        # Here we record the number of reps in this cross and parent and gRNA/Cas9 status
        Sum_Stats_Current$Cas9_gRNA_Pheno_Reps <-
          sum(current_cross_and_parent[(
            current_cross_and_parent$F2.gRNA ==  D_status_options[f] &
              current_cross_and_parent$F2.Nuclease ==  N_status_options[f]
          ),]$Reps)
        
        
        # Here we record the number of reps in this cross and parent and gRNA/Cas9 status and phenotype status
        Sum_Stats_Current$Pos_Reps <-
          sum(current_cross_and_parent[(
            current_cross_and_parent$F2.Som_Pheno == Pheno_status_options[p] &
              current_cross_and_parent$F2.gRNA ==  D_status_options[f] &
              current_cross_and_parent$F2.Nuclease ==  N_status_options[f]
          ),]$Reps)
        
        # If there are any reps, add the current summary statistics to the main data frame
        if (!is.na(sum(Sum_Stats_Current$Reps)) &
            sum(Sum_Stats_Current$Reps) > 0) {
          Sum_Stats_Current$Pos_Frac <-
            Sum_Stats_Current$Pos_Reps / Sum_Stats_Current$Reps
          
          Sum_Stats <- rbind(Sum_Stats, Sum_Stats_Current)
        }
      }
    }
  }
}

# We add columns for the p-value and p-value symbol
Sum_Stats$Ftestp_value <- NA
Sum_Stats$Ftestsymbol <- NA
Sum_Stats$Lower_binconf <- NA
Sum_Stats$Upper_binconf <- NA


# Loop over the rows of the data frame with Outcome == 'gRNA_gene'
# Fisher’s Exact tests to the control

for (i in 1:length(Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene', ][, 1])) {
  f_test <- fisher.test(matrix(c(
    c(Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene', ][i, ]$Pos_Reps,
      Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene', ][i, ]$Reps),
    c(sum(Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene' &
                      Sum_Stats$Nuclease.Promoter == 'Ctrl', ]$Pos_Reps),
      sum(Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene' &
                      Sum_Stats$Nuclease.Promoter == 'Ctrl', ]$Reps))
  ), nrow = 2, ncol = 2))
  
  Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Ftestp_value <-
    f_test$p.value
  Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Ftestsymbol <-
    ifelse(f_test$p.value <= 0.001,
           "***",
           ifelse(
             f_test$p.value <= 0.01,
             "**",
             ifelse(f_test$p.value <= 0.05, "*", 'ns')
           ))
  
  # Here we calculate the binomial confidence intervals
  Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Lower_binconf <-
    binconf(
      Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Pos_Reps,
      Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Lower']
  
  Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Upper_binconf <-
    binconf(
      Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Pos_Reps,
      Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',][i,]$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Upper']
  
}



#  gRNA inheritance control value is used to set the y-intercept of the dotted line
Ctrl <-
  Sum_Stats[Sum_Stats$Nuclease.Promoter == 'Ctrl' &
            Sum_Stats$Outcome == 'gRNA_gene',]$Pos_Reps / Sum_Stats[Sum_Stats$Nuclease.Promoter == 'Ctrl' &
                                                                    Sum_Stats$Outcome == 'gRNA_gene',]$Reps

# Here we output the summary statistics table
Sum_Stats

```


```{r}


# Create a new data frame with only the phenotype data from Sum_Stats
DF_Rows_UI_Phe <-
  Sum_Stats[Sum_Stats$Outcome == 'Cas9_gRNA_Pheno',]

# Here we reorder the phenotype factors so it look better in the graphs
DF_Rows_UI_Phe$F2.Som_Pheno <-
  ordered(DF_Rows_UI_Phe$F2.Som_Pheno,
          levels = c("Not Scored", "WT", "Mosaic/WE"))

# Here we reorder the rows of the phenotype data so it look better in the graphs
DF_Rows_UI_Phe <-
  DF_Rows_UI_Phe[with(DF_Rows_UI_Phe, rev(order(Cross, F2.Nuclease, F2.gRNA, F2.Som_Pheno))), ]

# We output the table in the markdown document
DF_Rows_UI_Phe


Sum_Stats <- Sum_Stats[Sum_Stats$Outcome == 'gRNA_gene',]

# Write the data frame to a csv file
write.csv(Sum_Stats, 'gRNA_Inheritance.csv')


# Write the data frame to a csv file
write.csv(DF_Rows_UI_Phe, 'Phenotypes.csv')

```


```{r}
# Here we create a figure for the male P0 crosses
# Fig 1 C and D

# We copy the main data table and then subset it
DF_Rows_UI <- Sum_Stats

# We select out the relevant crosses
DF_Rows_UI <-
  DF_Rows_UI[DF_Rows_UI$F0.Male == 'Nuc' | DF_Rows_UI$F0.Male == 'N/A',] # We select the crosses where the P0 parent carried only Cas9, and the control cross
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Outcome == 'gRNA_gene',]
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Parent == 'Pool',]
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Pos_Reps > 0,]


# We copy the main phenotype data table and then subset it
DF_Rows_UI_Phe_n <- DF_Rows_UI_Phe

# We select out the relevant crosses
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$F0.Male == 'Nuc'| DF_Rows_UI_Phe_n$F0.Male == 'N/A',]
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$Parent == 'Pool',]
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$Pos_Reps > 0,]


# Here we define the core data for the Fig 1C

ggplot_fig <- ggplot(data = DF_Rows_UI,
               aes(
                 x = as.character(Sex_Symbol[Drive.F1.sex]),
                 # X position is P1 sex
                 y = Pos_Reps / Reps,
                 # Y position is gRNA inheritance rate
                 stroke = 2,
                 color = Drive.F1.sex # circle color is P1 sex
               )) +
  geom_hline(yintercept = Ctrl,
             # We draw a line for the control value
             linetype = 2,
             color = 'grey') +
  geom_errorbar(aes(
    ymin = binconf(
      Pos_Reps,
      Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Lower'],
    ymax = binconf(
      Pos_Reps,
      Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Upper']
  ),
  width = .2,
  color = "black")

# Here we add a whole bunch of aesthetics
ggplot_fig <- ggplot_fig +
  geom_point(stat = "identity",
             alpha = 1,
             stroke = 1.25) +
  geom_text(
    label = DF_Rows_UI$Ftestsymbol,
    # We add the significance symbol
    size = 3.5,
    nudge_x = 0,
    nudge_y = binconf(
      DF_Rows_UI$Pos_Reps,
      DF_Rows_UI$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 3] - DF_Rows_UI$Pos_Frac + 0.035,
    check_overlap = F,
    show.legend  = F
  ) +
  geom_text(
    colour = 'black',
    size = 2.5,
    label = DF_Rows_UI$Reps,
    # We add the number of progeny scored to the figure. Some values get cut-off and need to be added manually
    nudge_x = 0,
    nudge_y = binconf(
      DF_Rows_UI$Pos_Reps,
      DF_Rows_UI$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 3] - DF_Rows_UI$Pos_Frac - 0.08, # change - to + to replicate graphs
    check_overlap = F,
    show.legend  = F
  ) +
  facet_nested(
    cols = vars(Nuclease.Promoter,
                InsertionVariant),
    # We organised the data by nuclease promoter and insertion variant
    drop = TRUE,
    scales = "free",
    space = "free"
  ) +
  xlab(expression('Trans-het P'[1] * ' sex')) +
  ylab(expression('kmo' ^ 'sgRNAs' * ' inheritance rate in P' [2] * ' larvae (%)'))

ggplot_fig$labels$colour <- expression('Trans-het P' [1] * ' sex')

name <- 'Fig1C'
x <- 9.75 # length of pdf file
y <- 3.25  # height of pdf file

Figure <- S_A(ggplot_fig) +  theme(
  panel.spacing = unit(0, "lines"),
  strip.background = element_rect(color = "black", fill =
                                    "white"),
                                  strip.text = element_text(color = "black"),
  axis.ticks.x = element_blank()
)

# We output the figure in the markdown document
Figure

# We create a pdf with and without the legend
PDF_me(Figure,
       name,
       wd,
       x,
       y)

# Here we define the core data for the Fig1D
Fig1p <- ggplot(
  data = DF_Rows_UI_Phe_n,
  aes(
    x = as.character(Sex_Symbol[Drive.F1.sex]),
    # X position is P1 sex
    y = Pos_Reps / Reps,
    # Y height is the fraction of P2 progeny with the particular phenotype and Cas9/gRNA status
    fill = interaction(F2.Nuclease, F2.gRNA),
    # Bars are colored according to Cas9/gRNA status
    linetype = F2.Som_Pheno,
    # Line type is set according to the phenotype
    alpha = F2.Som_Pheno,
    # alpha is set according to the phenotype
  )
) +
  geom_bar(
    position = "fill",
    stat = "identity",
    size = 0.5,
    color = 'black'
  ) +
  facet_nested(
    cols = vars(Nuclease.Promoter, # We organise the data by nuclease promoter and insertion variant
                InsertionVariant),
    drop = TRUE,
    scales = "free",
    space = "free"
  ) +
  scale_fill_manual(values = c("#808080ff",
                                          "#ff5555ff",
                                          '#2a7fffff',
                                          "#9955ffff")) + # Here we set the gRNA/Cas9 colors
                                            scale_alpha_manual(values = c(0.9, 1, 0.4), # Here we set somatic phenotype alpha level
                                                               guide = F) +
  xlab(expression('Trans-het P'[1] * ' sex')) +
  ylab('mosaic/white-eyed larvae (%)')


Fig1p <-
  S_A(Fig1p) + scale_linetype_manual(values = c("dashed", "solid", 'solid')) # Here we set somatic phenotype border type
Fig1p$labels$fill <- expression('(+/-)Cas9;(+/-)gRNA')
Fig1p$labels$colour <- expression('Scored or not')

name <- 'Fig1D'

Figure <- Fig1p +  theme(
  panel.spacing = unit(0, "lines"),
  strip.background = element_rect(color = "black", fill =
                                    "white"),
                                  strip.text = element_text(color = "black"),
  axis.ticks.x = element_blank()
)

# Fig 1D does not have the control data
# We make the figure slightly smaller to keep them aligned
x <- x - 0.4

# We output the figure in the markdown document
Figure

# We create a pdf with and without the legend
PDF_me(Figure,
       name,
       wd,
       x,
       y)

```





```{r}


# Here we create a figure for the female M0 crosses
# Fig 2 B and C


# We copy the main data table and then subset it
DF_Rows_UI <- Sum_Stats

# We select out the relevant crosses
DF_Rows_UI <-
  DF_Rows_UI[DF_Rows_UI$F0.Female == 'Nuc;gRNA',] # M0 parent is double het
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Outcome == 'gRNA_gene',]
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Parent == 'Pool',]
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Pos_Reps > 0,]

# We copy the main phenotype data table and then subset it
DF_Rows_UI_Phe_n <- DF_Rows_UI_Phe

# We select out the relevant crosses
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$F0.Female == 'Nuc;gRNA',]
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$Parent == 'Pool',]
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$Pos_Reps > 0,]


# Here we define the core data for Fig 2B
ggplot_fig <-
  ggplot(data = DF_Rows_UI,
         aes(
           x = as.character(Sex_Symbol[Drive.F1.sex]),
           # X position is M1 sex
           y = Pos_Reps / Reps,
           # Y position is gRNA inheritance rate
           stroke = 2,
           color = Drive.F1.sex # circle color is M1 sex
         )) +
  geom_hline(yintercept = Ctrl,
             # We draw a line for the control value
             linetype = 2,
             color = 'grey') +
  geom_errorbar(aes(
    ymin = binconf(
      Pos_Reps,
      Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Lower'],
    ymax = binconf(
      Pos_Reps,
      Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 'Upper']
  ),
  width = .2,
  color = "black")

# Here we add a whole bunch of aesthetics
ggplot_fig <- ggplot_fig +
  geom_point(stat = "identity",
             alpha = 1,
             stroke = 1.25) +
  geom_text(
    label = DF_Rows_UI$Ftestsymbol,
    # We add the significance symbol
    size = 3.5,
    nudge_x = 0,
    nudge_y = binconf(
      DF_Rows_UI$Pos_Reps,
      DF_Rows_UI$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 3] - DF_Rows_UI$Pos_Frac + 0.035,
    check_overlap = F,
    show.legend  = F
  ) +
  geom_text(
    colour = 'black',
    size = 2.5,
    label = DF_Rows_UI$Reps,
    # We add the number of progeny scored to the figure. Some values get cut-off and need to be added manually
    nudge_x = 0,
    nudge_y = binconf(
      DF_Rows_UI$Pos_Reps,
      DF_Rows_UI$Reps,
      alpha = 0.05,
      method = c("wilson"),
      include.x = FALSE,
      include.n = FALSE,
      return.df = FALSE
    )[, 3] - DF_Rows_UI$Pos_Frac - 0.08, # change - to + to replicate graphs
    check_overlap = F,
    show.legend  = F
  ) +
  facet_nested(
    cols = vars(Nuclease.Promoter, # We organise the data by nuclease promoter and insertion variant 
                InsertionVariant),
    drop = TRUE,
    scales = "free",
    space = "free"
  ) +
  xlab(expression('Trans-het M'[1] * ' sex')) +
  ylab(expression('kmo' ^ 'sgRNAs' * ' inheritance rate in M' [2] * ' larvae (%)'))

ggplot_fig$labels$colour <- expression('Trans-het M' [1] * ' sex')

name <- 'Fig2B'
x <- 3 # length of pdf file
y <- 3.5  # height of pdf file

Figure <- S_A(ggplot_fig) + theme(
    panel.spacing = unit(0, "lines"),
    strip.background = element_rect(color = "black",
                                    fill = "white"),
    strip.text = element_text(color = "black"),
    axis.ticks.x = element_blank()
  )

# We output the figure in the markdown document
Figure

# We create a pdf with and without the legend
PDF_me(
  Figure,
  name,
  wd,
  x,
  y
)


# Here we define the core data for Fig2C
Fig1p <- ggplot(
  data = DF_Rows_UI_Phe_n,
  aes(
    x = as.character(Sex_Symbol[Drive.F1.sex]), # X position is M1 sex
    y = Pos_Reps/Reps, # Y height is the fraction of M2 progeny with the particular phenotype and Cas9/gRNA status
    fill = interaction(F2.Nuclease, F2.gRNA), # Bars are colored according to Cas9/gRNA status
    linetype = F2.Som_Pheno, # Line type is set according to the phenotype
    alpha = F2.Som_Pheno # Alpha is set according to the phenotype
  )
) +
  geom_bar(position = "fill",
           stat = "identity",
           size = 0.5,
           color = 'black') +
  facet_nested(
    cols = vars(Nuclease.Promoter, # We organise the data by nuclease promoter and insertion variant
                InsertionVariant),
    drop = TRUE,
    scales = "free",
    space = "free"
  ) +
  scale_fill_manual(values = c("#808080ff",
                               "#ff5555ff",
                               '#2a7fffff',
                               "#9955ffff")) + # Here we set the gRNA/Cas9 colors
  scale_alpha_manual(values = c(0.9,
                                1,
                                0.4), # Here we set somatic phenotype alpha level
                     guide = F) +
  xlab(expression('Trans-het M'[1] * ' sex')) +
  ylab('mosaic/white-eyed larvae (%)') 
  # ylab(expression('kmo' ^ 's' * 'M' [2] * 'mosaic/white-eyed larvae (%)')) # used to make the B and C fig align better vertically


Fig1p <-
  S_A(Fig1p) + scale_linetype_manual(values = c("dashed", "solid", 'solid')) # Here we set somatic phenotype border type
Fig1p$labels$fill <- expression('(+/-)Cas9;(+/-)gRNA')
Fig1p$labels$colour <- expression('Scored or not')

name <- 'Fig2C'

Figure <- Fig1p + theme(
    panel.spacing = unit(0,
                         "lines"),
    strip.background = element_rect(color = "black",
                                    fill = "white"),
    strip.text = element_text(color = "black"),
    axis.ticks.x = element_blank()
  )


# We output the figure in the markdown document
Figure

# We create a pdf with and without the legend
PDF_me(
  Figure,
  name,
  wd,
  x,
  y
)

```



```{r}

# Here we create a figure for the individual crosses
# Fig3B and Fig3C


# We copy the main data table and then subset it
DF_Rows_UI <- Sum_Stats

# We select out the relevant crosses
DF_Rows_UI <- DF_Rows_UI[DF_Rows_UI$Outcome == 'gRNA_gene',]
DF_Rows_UI <-
  DF_Rows_UI[DF_Rows_UI$Parent != 'Pool',] #We select the non-pool data

# We copy the main phenotype data table and then subset it
DF_Rows_UI_Phe_n <- DF_Rows_UI_Phe

# We select out the relevant crosses
DF_Rows_UI_Phe_n <-
  DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$Outcome == 'Cas9_gRNA_Pheno' &
                     DF_Rows_UI_Phe_n$Parent != 'Pool',] #We select the non-pool data

# Here we define the core data for Fig 3B
fig <-
  ggplot(data = DF_Rows_UI) +
  geom_boxplot(
    aes(
      x = as.character(Sex_Symbol[Drive.F1.sex]),
      # X position is I1 sex
      y = Pos_Reps / Reps,
      # Y position is gRNA inheritance rate
      weight = Reps,
      # We weight the boxplot calculation by the number of reps of each observation
    ),
    outlier.shape = NA,
    width = 0.5
  ) +
  geom_quasirandom(
    aes(
      x = as.character(Sex_Symbol[Drive.F1.sex]),
      # X position is I1 sex
      y = Pos_Reps / Reps,
      # Y position is gRNA inheritance rate
      size = Reps,
      # The circle size is the number of I2 in that individual cross
      color = as.character(Sex_Symbol[Drive.F1.sex]) # circle color is I1 sex
    ),
    fill = NA,
    varwidth = TRUE,
    stroke = 1,
    shape = 1,
    width = 0.25
  ) +
  geom_hline(yintercept = Ctrl,
             # We draw a line for the control value
             linetype = 2,
             color = 'grey') +
  xlab(expression('Trans-het I'[1] * ' sex')) +
  ylab(expression('kmo' ^ 'sgRNAs' * ' inheritance rate in I' [2] * ' larvae (%)'))

fig$labels$size <- expression('Number of I'[2] * ' larvae scored')
fig$labels$colour <- expression('Trans-het I'[1] * ' sex')
#fig$labels$shape <- expression('Genotype/phenotype scored in I'[2] * ' larvae')

fig <-
  fig + geom_text(
    # We add the number of progeny scored to the figure.
    aes(x = as.character(Sex_Symbol['Male']),
        y = 0.70),
    colour = 'black',
    size = 2.5,
    label = sum(DF_Rows_UI[DF_Rows_UI$Outcome == 'gRNA_gene' &
                             DF_Rows_UI$Drive.F1.sex ==  'Male', ]$Reps),
    nudge_x = 0,
    nudge_y = 0,
    check_overlap = F,
    show.legend  = F
  ) +
  geom_text(
    # We add the number of progeny scored to the figure.
    aes(x = as.character(Sex_Symbol['Female']),
        y = 0.70),
    colour = 'black',
    size = 2.5,
    label = sum(DF_Rows_UI[DF_Rows_UI$Outcome == 'gRNA_gene' &
                             DF_Rows_UI$Drive.F1.sex ==  'Female', ]$Reps),
    nudge_x = 0,
    nudge_y = 0,
    check_overlap = F,
    show.legend  = F
  )


name <- 'Fig3B'
x <- 3.5 # length of pdf file
y <- 3.5 # height of pdf file


figure <- S_A(fig) + theme(
  panel.spacing = unit(3, "lines"),
  strip.background = element_rect(color = "black", fill =
                                    "white"),
                                  strip.text = element_text(color = "black"),
  axis.ticks.x = element_blank()
)

# We output the figure in the markdown document
figure

# We create a pdf with and without the legend
PDF_me(figure,
       name,
       wd,
       x,
       y)

# Below er create the phenotype data figure for the individual crosses
# The data are organised by 3 factors on the X-axis
# Here we simplify the labeling so it fits on the figure

# gRNA Yes/No is set to g/-
DF_Rows_UI_Phe_n$F2.gRNA[DF_Rows_UI_Phe_n$F2.gRNA == 'Yes'] <- 'g'
DF_Rows_UI_Phe_n$F2.gRNA[DF_Rows_UI_Phe_n$F2.gRNA == 'No'] <- '-'

# Cas9 Yes/No is set to N/-
DF_Rows_UI_Phe_n$F2.Nuclease[DF_Rows_UI_Phe_n$F2.Nuclease == 'Yes'] <-
  'N'
DF_Rows_UI_Phe_n$F2.Nuclease[DF_Rows_UI_Phe_n$F2.Nuclease == 'No'] <-
  '-'

# Here we reorder the gRNA factors
DF_Rows_UI_Phe_n$F2.gRNA <-
  factor(DF_Rows_UI_Phe_n$F2.gRNA, levels = c("g", "-"))

# Here we reorder the Cas9 factors
DF_Rows_UI_Phe_n$F2.Nuclease <-
  factor(DF_Rows_UI_Phe_n$F2.Nuclease, levels = c("N", "-"))

name <- 'Fig3C'

# Here we define the core data for Fig3C
fig <-
  ggplot(data = DF_Rows_UI_Phe_n[DF_Rows_UI_Phe_n$F2.Som_Pheno == 'Mosaic/WE',]) +
  geom_boxplot(
    aes(
      x = interaction(F2.Nuclease, interaction(F2.gRNA, as.character(Sex_Symbol[Drive.F1.sex]))),
      # X position is I2 sex, gRNA, Cas9 status
      y = Pos_Reps / Cas9_gRNA_Pheno_Reps,
      # Y height is the fraction of I2 progeny with the particular phenotype and Cas9/gRNA status
      weight = Cas9_gRNA_Pheno_Reps # We weight the boxplot calculation by the number of reps of each observation
    ),
    outlier.shape = NA
  ) +
  geom_quasirandom(
    aes(
      x = interaction(F2.Nuclease, interaction(F2.gRNA, as.character(Sex_Symbol[Drive.F1.sex]))),
      y = Pos_Reps / Cas9_gRNA_Pheno_Reps,
      size = Cas9_gRNA_Pheno_Reps,
      color = Drive.F1.sex
    ),
    fill = NA,
    varwidth = TRUE,
    stroke = 0.7,
    shape = 2,
    width = 0.25
  ) +
  xlab(expression('Trans-het I'[1] * ' sex, sgRNA inheritance, and Cas9 inheritance')) +
  #ylab(expression('mosaic/white-eyed I' [2] * ' larvae (%)'))
  ylab(expression('kmo' ^ 's' * 'M' [2] * 'mosaic/white-eyed larvae (%)')) # used to make the B and C figs align better vertically


figure <- S_A(fig) + scale_shape_manual(values = c(2)) +
  scale_shape_manual(values =
                       c(1, 2)) +
  theme(
    panel.spacing = unit(0, "lines"),
    strip.background = element_rect(color =
                                      "black", fill = "white"),
                                    strip.text = element_text(color =
                                                                "black"),
                                                              axis.ticks.x = element_blank(),
    legend.position = c(0.5, 0.5)
  )

# We output the figure in the markdown document
figure

# We create a pdf with and without the legend
PDF_me(figure,
       name,
       wd,
       x,
       y)

```
