In UMN-BarleyOatSilphium/GSSimTPUpdate: Code and data to accompany the publication [TITLE]

Introduction

Gather data from the simulation into a data frame for plotting

library(GSSimTPUpdate)
library(tidyverse)
library(stringr)

Load Data

Option 1: You have run the simulation yourself and want to use the files created

# Project directory, or the file path leading to the "GSSimTPUpdate" repository. You must
# edit this yourself
project_dir <- ""


# The plotdata directtory, or the directory in which the file ending with "plotdata.RData" is stored.
# You must edit this yourself
plotdata_dir <- file.path(project_dir, "inst/output/")

# Load the window data
filename <- file.path(output_dir, "simulation_results_MNxND_window_collective_plotdata.RData")
load(filename)

plot.list.window <- plot.list
plot.list.window <- plot.list.window[-12]

# Load the cumulative data
filename <- file.path(output_dir, "simulation_results_MNxND_cumulative_collective_plotdata.RData")
load(filename)

plot.list.cumulative <- plot.list

# Combine
plot.list <- mapply(plot.list.window, plot.list.cumulative, FUN = bind_rows)
# Remove missing
plot.list <- sapply(plot.list, na.omit)

plotting_data <- plot.list

Option 2: You want to use the data included in the package

# Load the data from the package
data("plotting_data")

Setup

Create other objects necessary for plotting

tp.change.factors <- as.factor(c(best = "Top", CDmean = "CDmean", nochange = "No Change", 
                       PEVmean = "PEVmean", random = "Random", tails = "Tails", worst = "Bottom") )


n.cycles <- plotting_data$df.acc$cycle %>% unique() %>% na.omit() %>% length()

Table 1 - Prediction Accuracy

df1 <- plotting_data$df.acc

# Determine value to display
df2 <- df1 %>% 
  mutate(disp = str_c(round(mean, 3), " (", round(mean - ci, 3), ", ", round(mean + ci, 3), ")"))

# Rearrange to make table 
df3 <- df2 %>% 
  select(exp_name, change, cycle, disp) %>% 
  spread(cycle, disp)

# Change names
df4 <- df3 %>%
  rename(Scenario = exp_name,
         Method = change) %>%
  mutate(Method = str_replace_all(Method, tp.change.factors))

Figure 3 - Prediction Accuracy

df1 <- plotting_data$df.acc

# Designate labels for the individual plots within facets
n.facets <- df1 %>% 
  select(exp_name, variable) %>%
  distinct() %>% 
  nrow()

gp <- df1 %>%
  ggplot(aes(x = cycle.offset, y = mean, col = change, shape = change)) +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), col = "black", width = 0.10) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  ylab(expression(Prediction~Accuracy~(italic(r)[MG]))) +
  xlab("Breeding Cycle") +
  scale_color_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_shape_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_x_continuous(breaks = seq(1, n.cycles + 1, 3)) +
  # Faceting
  facet_wrap("exp_name")

# Geom text data
gt.df <- ggplot_build(gp)$plot$data %>%
  group_by(variable, exp_name) %>%
  summarize(ymax = max(mean + ci)) %>%
  ungroup() %>%
  mutate(
    y = max(ymax) * 1.05,
    x = 1,
    label = LETTERS[seq_len(n.facets)])


# Modify fonts
gp1 <- gp + 
  theme(
    strip.text.x = element_text(face = "bold", size = 12),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.title = element_text(face = "bold", size = 12),
    legend.text = element_text(size = 10)) +
  geom_text(data = gt.df, aes(x = x, y = y, label = label, fontface = 2), 
            inherit.aes = FALSE, size = 5)

Figure 4 - Genotypic Value and Genetic Variance

plotting_data$df.genvar <- plotting_data$df.genvar %>% 
  mutate(variable = "Genetic Variance")
plotting_data$df.genval <- plotting_data$df.genval %>% 
  mutate(variable = "Genotypic Value")

df1 <- bind_rows(plotting_data$df.genvar, plotting_data$df.genval)

# Designate labels for the individual plots within facets
n.facets <- df1 %>% 
  select(exp_name, variable) %>%
  distinct() %>% 
  nrow()

gp <- df1 %>%
  ggplot(aes(x = cycle.offset, y = mean, col = change, shape = change)) +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), col = "black", width = 0.10) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  ylab("") +
  xlab("Breeding Cycle") +
  scale_color_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_shape_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_x_continuous(breaks = seq(1, n.cycles + 1, 3))

# Geom text data
gt.df <- ggplot_build(gp)$plot$data %>%
  group_by(variable, exp_name) %>%
  summarize(ymax = max(mean + ci)) %>%
  mutate(y = max(ymax) * 1.05,
         x = 1,
         label = LETTERS[seq_len(n.facets)])

gt.df <- ggplot_build(gp)$plot$data %>%
  group_by(variable, exp_name) %>%
  summarize(ymax = max(mean + ci)) %>%
  full_join(., summarize(., y = max(ymax) * 1.05)) %>%
  ungroup() %>%
  mutate(x = 1,
         label = LETTERS[seq_len(n.facets)])


# Modify fonts
gp1 <- gp + 
  theme(
    strip.text = element_text(face = "bold", size = 12),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.title = element_text(face = "bold", size = 12),
    legend.text = element_text(size = 10)) +
  facet_grid(variable ~ exp_name, scales = "free_y", switch = "y") +
  geom_text(data = gt.df, aes(x = x, y = y, label = label, fontface = 2), 
            inherit.aes = FALSE, size = 5)

Figure 5

Genomic relationship, persistence of phase, inbreeding, and number of QTL fixed for an allele

### Figure 5 will have the genomic relationship, persistence of phase, inbreeding,
### and number of QTL fixed for an allele

plotting_data$df.pers <- plotting_data$df.pers %>% 
  mutate(variable = "Persistence of LD Phase")
plotting_data$df.rel <- plotting_data$df.rel %>% 
  mutate(variable = "Average Relationship\n(Scaled to Base Population)")
plotting_data$df.inbred <- plotting_data$df.inbred %>% 
  mutate(variable = "Inbreeding\n(Scaled to Base Population)")
plotting_data$df.fixedqtl <- plotting_data$df.fixedqtl %>% 
  mutate(variable = "Number of Fixed QTL")


df1 <- bind_rows(plotting_data$df.pers, plotting_data$df.rel, plotting_data$df.inbred, plotting_data$df.fixedqtl)

# Designate labels for the individual plots within facets
n.facets <- df1 %>% 
  select(exp_name, variable) %>%
  distinct() %>% 
  nrow()


gp <- df1 %>%
  ggplot(aes(x = cycle.offset, y = mean, col = change, shape = change)) +
  geom_errorbar(aes(ymin = mean - ci, ymax = mean + ci), col = "black", width = 0.10) +
  geom_line(size = 1) +
  geom_point(size = 2) +
  ylab("") +
  xlab("Breeding Cycle") +
  scale_color_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_shape_discrete(name = "Update Method",
                       labels = as.character(tp.change.factors)) +
  scale_x_continuous(breaks = seq(1, n.cycles + 1, 3))

# Geom text data
gt.df <- ggplot_build(gp)$plot$data %>%
  group_by(variable, exp_name) %>%
  summarize(ymax = max(mean + ci)) %>%
  full_join(., summarize(., y = max(ymax) * 1.05)) %>%
  ungroup() %>%
  mutate(x = 1,
         label = LETTERS[seq_len(n.facets)])

gp1 <- gp + 
  theme(
    strip.text = element_text(face = "bold", size = 12),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10),
    legend.title = element_text(face = "bold", size = 12),
    legend.text = element_text(size = 10)) +
  facet_grid(variable ~ exp_name, scales = "free_y", switch = "y") +
  geom_text(data = gt.df, aes(x = x, y = y, label = label, fontface = 2), 
            inherit.aes = FALSE, size = 5)

UMN-BarleyOatSilphium/GSSimTPUpdate documentation built on May 9, 2019, 7:40 p.m.