1 Setup

library(here)
source(here::here("code", "scripts", "fecundity", "source.R"))

2 Inbreeding

Data saved here (copied from Felix’s table in the paper): data/fecundity/20210212_inbreeding.csv

2.1 Read in data

raw_dat = here::here("data", "fecundity", "20210212_inbreeding.csv")
inbreed_df = readr::read_csv(raw_dat)
## Rows: 12 Columns: 8
## ── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): MIKK, Inbreeding (Generation), N of strains lost
## dbl (5): not productive, infertile male, no males, no females, all cause death
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
total_crosses = 253

inbreed_df_tidy = inbreed_df %>%
  # Remove final row, which makes `N of strains lost` character
  dplyr::slice(-c(11,12)) %>% 
  # Change column names
  dplyr::select(GENERATION = `Inbreeding (Generation)`,
                N_LOST_TOTAL = `N of strains lost`,
                everything()) %>% 
  # Make N_LOST integer
  dplyr::mutate(N_LOST_TOTAL = as.integer(N_LOST_TOTAL)) %>% 
  # Remove MIKK column
  dplyr::select(-MIKK) %>% 
  # Get cumulative lost and N_SURVIVING
  dplyr::mutate(N_LOST_CUM = cumsum(N_LOST_TOTAL),
                N_SURVIVING = total_crosses - N_LOST_CUM) %>% 
  # Replace 
  dplyr::mutate(dplyr::across(c(`not productive`,
                                `infertile male`,
                                `no males`,
                                `no females`,
                                `all cause death`),
                              ~tidyr::replace_na(.x, 0))) %>% 
  # Get cumulative lost for each category
  dplyr::mutate(dplyr::across(c(`not productive`,
                                `infertile male`,
                                `no males`,
                                `no females`,
                                `all cause death`),
                              ~cumsum(.x))) %>% 
  # Amend GENERATION column
  dplyr::mutate(GENERATION = stringr::str_extract(GENERATION, "(?:F\\d{1,2})")) %>% 
  # Add first row
  rbind(., c("F1", rep(0,7), total_crosses)) %>% 
  # Add rows for F12-14
  rbind(., dplyr::filter(., GENERATION == "F11") %>%
            dplyr::slice(rep(1:n(), each = 3)) %>%
            dplyr::mutate(GENERATION = c("F12", "F13", "F14"))) %>% 
  # Factorise GENERATION
  dplyr::mutate(GENERATION = factor(GENERATION, levels = paste("F", seq(1:14), sep = ""))) %>% 
  # Order
  .[order(.$GENERATION), ] %>% 
  # Gather 
  tidyr::pivot_longer(cols = c(`not productive`,
                               `infertile male`,
                               `no males`, 
                               `no females`,
                               `all cause death`,
                               N_SURVIVING),
                      names_to = "STATUS",
                      values_to = "N") %>% 
  # Convert `N` to integer
  dplyr::mutate(N = as.integer(N))

knitr::kable(head(inbreed_df_tidy))
GENERATION N_LOST_TOTAL N_LOST_CUM STATUS N
F1 0 0 not productive 0
F1 0 0 infertile male 0
F1 0 0 no males 0
F1 0 0 no females 0
F1 0 0 all cause death 0
F1 0 0 N_SURVIVING 253

2.2 Plot

recode_vec = c("Infertile male",
               "No males",
               "No females", 
               "Not productive",
               "All-cause death", 
               "Reproducing")
names(recode_vec) = c("infertile male", 
                      "no males", 
                      "no females", 
                      "not productive", 
                      "all cause death", 
                      "N_SURVIVING")

test = inbreed_df_tidy %>% 
  # geom_area() seems to need `x` to be numeric
  dplyr::mutate(GENERATION = as.integer(gsub("F", "", GENERATION))) %>%
  # recode STATUS values and order
  dplyr::mutate(STATUS = dplyr::recode(STATUS, !!!recode_vec),
                STATUS = factor(STATUS, levels = recode_vec))# %>% 
out2 = ggplot(test) +
    geom_area(aes(x=GENERATION,
                  y=N,
                  fill = STATUS)) +
#    geom_text(aes(x=GENERATION, y=N, label = STATUS)) +
    scale_fill_viridis_d(option = "magma") +
    scale_x_continuous(breaks = unique(test$GENERATION),
                       labels = paste("F", unique(test$GENERATION), sep = "")) +
    theme_cowplot(font_size = 10) +
    xlab("Generation") +
    ylab("Number of MIKK panel lines") +
    labs(fill = "Status") +
    theme(panel.grid = element_blank())
#    theme(legend.position=c(0.2, 0.3),
#          legend.box.background = element_blank())

out2

2.2.1 New palette

reviewer_pal = c("#FC4E07", "#FFBF00", "#0BC166", "#00AFBB", "#360568", "#D84797")

out3 = ggplot(test) +
    geom_area(aes(x=GENERATION,
                  y=N,
                  fill = STATUS)) +
    scale_fill_manual(values = reviewer_pal) +
    scale_x_continuous(breaks = unique(test$GENERATION),
                       labels = paste("F", unique(test$GENERATION), sep = "")) +
    theme_cowplot(font_size = 16) +
    xlab("Generation") +
    ylab("Number of MIKK panel lines") +
    labs(fill = "Status") +
    theme(panel.grid = element_blank())

out3

3 Fecundity

3.1 Read in data

in_file = here::here("data", "fecundity", "20210205_semiquantitative.xlsx")
df_semi = readxl::read_xlsx(in_file, range = "A1:C81")
out_file = here::here("data", "fecundity", "20210208_semiquant.csv")

# Get strain levels
strain_levels = unique(df_semi$Pair)

# Create recode vector
date_recode = c("Feb 2019", "Jul 2020")
names(date_recode) = c("2/19", "7/20")
recode_vec_1 = c(0, 1, 2, 3, 4, 5)
names(recode_vec_1) = c(0, "o", "x", "x/", "xx", "xxx")
recode_vec_2 = c("Not producing",
                 "<3 eggs; not every day",
                 "<5 eggs; not every day",
                 "1-3 eggs per day",
                 "1-5 eggs per day",
                 "5-10 eggs per day")
names(recode_vec_2) = c(0, 1, 2, 3, 4, 5)
recode_vec_3 = gsub("; ", ";\n", recode_vec_2)
recode_vec_3 = gsub("eggs per", "eggs\nper", recode_vec_3)

# Tidy

semi_out = df_semi %>% 
  # pivot fecundity
  tidyr::pivot_longer(cols = contains("fecundity"), 
                      names_to = "DATE",
                      names_prefix = "fecundity ",
                      values_to = "FECUNDITY") %>% 
  # recode fecundity measures
  dplyr::mutate(DATE = dplyr::recode(DATE, !!!date_recode),
                FECUNDITY = dplyr::recode(FECUNDITY, "xx/" = "xx"),
                FECUNDITY = dplyr::na_if(FECUNDITY, "do not prod. Yet"),
                FECUNDITY = ifelse(is.na(FECUNDITY), 0, FECUNDITY),
                FECUNDITY = dplyr::recode(FECUNDITY, !!!factor(recode_vec_1)),
                KEY = dplyr::recode(FECUNDITY, !!!recode_vec_2)) %>%
  # rename STRAIN
  dplyr::rename(STRAIN = Pair) %>% 
  # factorise
  dplyr::mutate(STRAIN = factor(STRAIN, levels = strain_levels)) %>% 
  # write to file
  readr::write_csv(out_file, na = "")

knitr::kable(head(semi_out))
STRAIN DATE FECUNDITY KEY
4-1 Jul 2020 2 <5 eggs; not every day
4-1 Feb 2019 2 <5 eggs; not every day
4-2 Jul 2020 2 <5 eggs; not every day
4-2 Feb 2019 4 1-5 eggs per day
5-1 Jul 2020 2 <5 eggs; not every day
5-1 Feb 2019 4 1-5 eggs per day

3.2 Plot

# Process data
final_hor_df = semi_out %>% 
  dplyr::mutate(KEY = gsub("; ", ";\n", KEY),
                KEY = gsub("eggs per", "eggs\nper", KEY),
                KEY = factor(KEY, levels = recode_vec_3)) %>% 
  dplyr::filter(DATE == "Jul 2020") %>% 
  dplyr::mutate(STRAIN = factor(STRAIN, levels=unique(STRAIN[order(KEY,decreasing = T)]))) %>% 
  # remove "Not producing"
  dplyr::filter(KEY != "Not producing")
  # order by category

# Generate histogram
fec_count = final_hor_df %>% 
  ggplot() +
    geom_histogram(aes(KEY, fill = KEY),
                   stat = "count") +
    scale_fill_viridis_d() +
    theme_cowplot(font_size = 16) +
    xlab("Fecundity") +
    ylab("Number of MIKK panel lines") +
    guides(fill = F) +
    theme(axis.text = element_text(size = 14))

fec_count

4 Add eye measurement data

4.1 Read in data

eye_dat_file = here::here("data/fecundity/eye_rel_measurements_line_ids.txt")

eye_df = readr::read_delim(eye_dat_file,
                           delim = "\t",
                           col_types = c("cdc"))

eye_df %>% 
  DT::datatable()

4.2 Explore data

# How many of each line and sex?
eye_df %>% 
  dplyr::count(line, gender)
## # A tibble: 152 × 3
##    line  gender     n
##    <chr> <chr>  <int>
##  1 10-1  Female     2
##  2 10-1  Male       3
##  3 104-1 Female     2
##  4 104-1 Male       1
##  5 106-1 Female     2
##  6 106-1 Male       1
##  7 106-2 Female     2
##  8 106-2 Male       1
##  9 11-1  Female     2
## 10 11-1  Male       1
## # … with 142 more rows
# Exclude males and get means
eye_df_filt = eye_df %>% 
  #dplyr::filter(gender == "Female") %>% 
  dplyr::group_by(line, gender) %>% 
  dplyr::summarise(MEAN_EYE = mean(eye_rel)) %>% 
  dplyr::arrange(MEAN_EYE)
## `summarise()` has grouped output by 'line'. You can override using the `.groups` argument.
eye_df_filt %>% 
  ggplot() +
    geom_histogram(aes(MEAN_EYE), bins = 40) +
    facet_wrap(vars(gender)) +
    theme_cowplot() +
    guides(fill = "none", colour = "none")

gender_pal_1 = c("#0D1B2A", "#1B263B")
gender_pal_2 = c("#FA9702", "#922294")

eye_dens = eye_df_filt %>% 
  ggplot() +
    geom_density(aes(MEAN_EYE, fill = gender, colour = gender)) +
    scale_fill_manual(values = gender_pal_2) +
    scale_colour_manual(values = darker(gender_pal_2, 70)) +
    facet_wrap(vars(gender)) +
    theme_cowplot(font_size = 16) +
    guides(fill = "none", colour = "none") +
    xlab("Mean relative eye size") +
    ylab("Density")

eye_dens

4.3 Plot

# order 
eye_df_filt_f = eye_df_filt %>% 
  dplyr::filter(gender == "Female") %>% 
  dplyr::arrange(MEAN_EYE) 

line_order = eye_df_filt_f$line

eye_df_filt_f %>% 
  dplyr::mutate(line = factor(line, levels = line_order)) %>% 
  ggplot() +
    geom_col(aes(line, MEAN_EYE, fill = line)) +
    scale_fill_viridis_d(option = "magma") +
    theme_cowplot() +
    guides(fill = "none") +
    theme(axis.text.x = element_text(size = 5, angle = 45, hjust = 1))

4.4 Compile

final = ggdraw() +
  draw_plot(out3, x = 0, y = 0.45, width = .55, height = 0.55) +
  draw_plot(fec_count, x = .55, y = 0.45, width = .45, height = 0.55) +
  draw_plot(eye_dens, x = 0, y = 0, width = 1, height = 0.45) +
  draw_plot_label(label = c("A", "B", "C"), size = 24,
                  x = c(0, 0.55, 0), y = c(1, 1, 0.45))

final

# SVG
ggsave(here::here("docs/plots/fecundity/20220209_final_figure.pdf"),
       final,
       device = "pdf",
       units = "in",
       height = 12,
       width = 16)

# PNG
ggsave(here::here("docs/plots/fecundity/20220209_final_figure.png"),
       final,
       device = "png",
       dpi = 400,
       units = "in",
       height = 12,
       width = 16)