R Notebook

library(tidyverse)
library(ggjoy)
library(viridis)
library(cowplot)

Munge!

Simulated data from \(N_e=10,000\) with bottleneck and growth following Beissinger et al. 2016. Munged data will be neutral \(\pi\) and \(\xi_1\) every 100 generations in 25 4Kb windows centered around “gene” which has neutral and deleterious mutations. Deleterious mutations are drawn from a gamma with mean \(s=-0.0005 (N_es\approx5)\).

read_sim<-function(pidata,xidata){
  #read data
  pi<-as.tibble(read.csv(pidata,header=F))
  pi<-mutate(pi,stat=rep("pi",length(pi[,1])))
  xi<-as.tibble(read.csv(xidata,header=F)) 
  xi<-mutate(xi,stat=rep("xi",length(xi[,1])))

  #munge into single tibble of means
  data<-rbind(xi,pi) %>%
    `colnames<-`(c("gen", 1:25, "sim", "stat")) %>%
    gather(window,value,2:26) %>% 
    mutate(gen=gen-100000) %>%
    group_by(gen,stat,window) %>%
    summarize(mean_val=mean(value)) %>%
    mutate(window=(as.numeric(window)-13)*4000*3E-5,mean_val=if_else(window==13,mean_val,mean_val/4))
  
  return(data)
}

simplot<-function(somedata,mystat,label){
  ggplot(filter(somedata,gen %% 100 <1,stat==mystat), aes(x=window, y=mean_val, group = gen,color=gen)) + 
  scale_color_viridis() +   
  background_grid(major = "xy", minor = "none") +
  geom_line() + 
  ylab(label)  + xlab("cM") +
  scale_x_continuous(breaks = seq(-1.5,1.5,0.5))
}

Visualize

\(s=5E-4\)

Maize

maize<-read_sim("~/Desktop/sim_neutralPi_s5e-4.csv","~/Desktop/sim_singleton_s5e-4.csv")

## Warning: package 'bindrcpp' was built under R version 3.3.2

\(\pi\)

simplot(maize,"pi",expression(bar(pi)))

\(\xi_1\)

simplot(maize,"xi",expression(bar(xi[i])))

Europeans

euros<-read_sim("~/Desktop/sim_neutralPi_European.csv","~/Desktop/sim_singleton_European.csv")

\(\pi\)

simplot(euros,"pi",expression(bar(pi)))

\(\xi_1\)

simplot(euros,"xi",expression(bar(xi[i])))

\(s=5E-3\)

Maize

maize<-read_sim("~/Desktop/sim_neutralPi_s5e-3.csv","~/Desktop/sim_singleton_s5e-3.csv")

\(\pi\)

simplot(maize,"pi",expression(bar(pi)))

\(\xi_1\)

simplot(maize,"xi",expression(bar(xi[i])))

\(s=5E-3\)

Maize

maize<-read_sim("~/Desktop/sim_neutralPi_s5e-5.csv","~/Desktop/sim_singleton_s5e-5.csv")

\(\pi\)

simplot(maize,"pi",expression(bar(pi)))

\(\xi_1\)

simplot(maize,"xi",expression(bar(xi[i])))