knitr::opts_chunk$set( echo = TRUE )
This is an R package with a number of functions to allow users to build up forward-in-time simulations to explore the dynamics of diverging genomes under various scenarios of gene flow, selection and genotype-phenotype maps.
It has been written by Claudio S. QuilodrĂ¡n, Eric C. Anderson, and Tim Coulson. The details of it appear in a manuscript submitted to Methods in Ecology and Evolution entitled, "The multiple population genetic and demographic routes to islands of genomic divergence." More information about the R package 'glads' can be found here.
To use it, you must:
r
dplyr,
magrittr,
pegas,
progress,
Rcpp (>= 0.12.7),
readr
{sh, eval=FALSE}
git clone https://github.com/eriqande/glads
R CMD INSTALL --no-multiarch --with-keep.source glads
If you prefer, after cloning it from GitHub, you can open the
RStudio project (glads/glads.Rproj
) within it and build and install
by using the "Install and Restart" button within RStudio.
Full documentation for the package can be found with:
help(package = "glads")
After that, to see examples of its use, you can run through the examples in
R-examples/Example1.R
and R-examples/Example2.R
within the glads
repository.
In fact, we will mirror R-examples/Example1.R
here:
#-------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------- # # Examples for submitted manuscript: # The many population genetic and demographic routes to islands of genomic divergence # # Example 1. Selection: recombination map and fitness function # Example 2. Neutral evolution: Mutation rate, average recombination, biallelic SNPs, # without fitness function. # Questions to: # Claudio S. QuilodrĂ¡n. Department of Zoology, University of Oxford. # Email: claudio.quilodran@zoo.ox.ac.uk; claudio.quilodran@unige.ch # Eric C. Anderson. Department of Ecology and Evolutionary Biology, University of California. # Email: eric.anderson@noaa.gov #-------------------------------------------------------------------------------------------- #-------------------------------------------------------------------------------------------- library(glads)
################################# ## Example 1: see Table 1 ## ## recombination map, ## ## with fitness function ## ################################# ######################## ### Generating data ### ######################## initial.population.size=400 #Initial population size n.loci=300 #Number of loci n.alleles.per.locus=20 #Number of alleles set.seed(2) start.1 <- start.2 <- initial.struct(initial.population.size,n.loci,n.alleles.per.locus) start <- list(start.1, start.1) #Initial populations ################### ### Parameters ### ################### loci.pos=1:n.loci #Loci position add.loci=50 #Number of additive loci fitness.pos <- 125:(125+(add.loci-1)) #Additive loci position sex.ratio <- 0.5 bvs = t(array( seq(0,1, length = n.alleles.per.locus) ,c(n.alleles.per.locus, add.loci))) #Breeding value of additive loci e.v=0.01 #Stochastic environmental variant d.v=1 #Stochastic demographic variant n.gens=100 #Number of generations recom.map= c(rep(0.5, 299)) #Recombination map Linked<-c(60:69, 150:159, 230:239) #Linked loci recom.map[Linked]<-0.0001 #Recombination map with linked loci #b0 Maximum generated offspring #b1 Phenotypic optima #b2 Variance of the fitness curve #b3 Density-dependent demographic variant ##Parameters of the fitness function param.w1 <- list(b0 = 6,b1 = 0.25, b2 = 0.5, b3 = 0.01, d.v = d.v, add.loci = add.loci) param.w2 <- list(b0 = 6,b1 = 0.75, b2 = 0.5, b3 = 0.005, d.v = d.v, add.loci = add.loci) ##Parameters of the phenotype function param.z1 <- list(sex.ratio, fitness.pos, bvs, add.loci, e.v) param.z2 <- param.z1
################### ### Simulations ### ################### example <- evolve( x = start, time = n.gens, type = "additive", recombination = "map", recom.rate = recom.map, param.z = list(param.z1, param.z2), param.w = list(param.w1, param.w2) )
################### ### Fst ### ################### ##Commputation of Fst values #We should first convert the output from class 'struct' to class 'loci' example.loci <- struct2pegas(example) #Estimating Fst with library 'pegas' FST <- pegas::Fst(example.loci)
#Figure plot(loci.pos, FST[,"Fst"], pch=16, bty="l", xlab="Loci position", ylab="Fst", ylim=c(0,1), type="l", xaxs="i", yaxs="i", las=1, lwd=2 )
################################# ## Example 2: Mutation rate, ## ## average recombination, ## ## biallelic SNPs ## ## without fitness function ## ################################# ######################## ### Generating data ### ######################## Nind = 200 #Number of individuals Nsnp = 1000 #Number of SNPs n.alleles = 2 #Random genetic identity set.seed(1) start.1 <- start.2 <- initial.struct(Nind, Nsnp, n.alleles) ################### ### Parameters ### ################### loci.pos<-sort(sample(1:12000000, 1000, replace=F) ) #Loci position chromo_mb<-12000000 #Chromosome size sex.ratio = 0.5 #Sex ratio mean.fitness = 3 #Offspring per breeding pair crossover <- 1.0/100000000.0 #Average recombination rate (cM/MB) mutation.rate=1.1*10^-8 #Mutation per site per generation n.gens=100 #Number of generations migration.rate=0.075 #Migration rate d.d=0.005 #Demographic effect to avoid exponential growth param.z0 <- list(sex.ratio = sex.ratio) param.w0 <- list(mean.fitness = mean.fitness, d.d = d.d) ################### ### Simulations ### ################### example2 <- evolve( x = list(start.1, start.2), time = n.gens, type = "dynamic", recombination = "average", recom.rate = crossover, migration.rate = migration.rate, mutation.rate = mutation.rate, loci.pos = loci.pos, chromo_mb = chromo_mb, param.z = list(param.z0, param.z0), param.w = list(param.w0, param.w0) ) ################### ### Fst ### ################### ##Commputation of Fst values #We should first convert the output from class 'struct' to class 'loci' example.loci2 <- struct2pegas(example2) #Estimating Fst with library 'pegas' FST <- pegas::Fst(example.loci2) #Figure plot(loci.pos, FST[,"Fst"], pch=16, bty="l", xlab="Loci position", ylab="Fst", ylim=c(0,1), type="l", xaxs="i", yaxs="i", las=1, lwd=2 )
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.