ResponsetoSelection

Validation for response to selection

We used the breeder’s equation to validate the implementation of selection on quantitative traits. This equation estimates the response to selection (R) from the heritability of the trait (h^2) and the selection differential (S).

                                         R=h^2 S

A quantitative trait with 10 QTLs having identical effect sizes (0.1) was simulated. All 10 QTLs were unlinked.
The response to selection for the simulations was estimated by subtracting the mean phenotypic value of the parental population from the mean phenotypic value of the offspring population, after one generation of truncating selection.

In total 100 independent simulations were performed for each scenario (different heritability and proportion of selected individuals). The simulated values were compared to the expected values from the breeder's equation given above.

In the above figure the x-axis corresponds to the proportion of selected individuals and y-axis to the response to selection (R). Three different plots were made for thee different values of heritability (0.3, 0.5, 0.8).

This demonstrates that selection on a quantitative trait is correctly implemented in MimicrEE2

Materials and Methods

All input files used to the validation can be found here.

In this section we will only describe how to produce the plot with heritability 0.5. If you want to produce the whole plot you just have to repeat the steps changing the heritability level used as input for MimicrEE2. The output is redirected from the terminal to a specified output file.

First, simulations with MimicrEE2 have to be performed for the different selection regimes.

:::bash
for i in {1..100}; do java -jar mim2.jar qt --haplotypes-g0 base_pops_45freq/basepop${i}.45freq.txt --recombination-rate rr_10loci.txt --effect-size ten_loci_chosen.txt --heritability 0.5 --replicate-runs 1 --snapshots 1 --selection-regime selection_regime_02.txt --threads 6 --output-sync response_to_selection.sync; done &> log_output_02selreg

for i in {1..100}; do java -jar mim2.jar qt --haplotypes-g0 base_pops_45freq/basepop${i}.45freq.txt --recombination-rate rr_10loci.txt --effect-size ten_loci_chosen.txt --heritability 0.5 --replicate-runs 1 --snapshots 1 --selection-regime selection_regime_04.txt --threads 6 --output-sync response_to_selection.sync; done &> log_output_04selreg

for i in {1..100}; do java -jar mim2.jar qt --haplotypes-g0 base_pops_45freq/basepop${i}.45freq.txt --recombination-rate rr_10loci.txt --effect-size ten_loci_chosen.txt --heritability 0.5 --replicate-runs 1 --snapshots 1 --selection-regime selection_regime_06.txt --threads 6 --output-sync response_to_selection.sync; done &> log_output_06selreg

for i in {1..100}; do java -jar mim2.jar qt --haplotypes-g0 base_pops_45freq/basepop${i}.45freq.txt --recombination-rate rr_10loci.txt --effect-size ten_loci_chosen.txt --heritability 0.5 --replicate-runs 1 --snapshots 1 --selection-regime selection_regime_08.txt --threads 6 --output-sync response_to_selection.sync; done &> log_output_08selreg

Next a Python script is provided which captures the summary statistics available in the MimicrEE2 output using regular expressions. The script further estimates R by subtracting the phenotypic mean of the offspring from the phenotypic mean of the base population. Use the --help flag to access further information about the script.

:::bash
python estimateR2.py --input log_output_02selreg --heritability 0.5 --selection-regime 0.2 > 02selreg_Rvalue

python estimateR2.py --input log_output_04selreg --heritability 0.5 --selection-regime 0.4 > 04selreg_Rvalue

python estimateR2.py --input log_output_06selreg --heritability 0.5 --selection-regime 0.6 > 06selreg_Rvalue

python estimateR2.py --input log_output_08selreg --heritability 0.5 --selection-regime 0.8 > 08selreg_Rvalue

Next we concatenate the output files in one file:

:::bash

cat 02selreg_Rvalue 04selreg_Rvalue 06selreg_Rvalue 08selreg_Rvalue > 05heritability_Rvalues

Finally, we load the output file in R to start with the plotting procedure.

:::R
setwd("/my/working/dir/")

her05_R<- read.table("05heritability_Rvalues")

colnames(her05_R)<- c("R", "Proportion_Selected", "heritability")

ggplot(her05_R, aes(x=Proportion_Selected, y=R, group=Proportion_Selected)) + geom_boxplot(width=0.01, fill="forestgreen")+ 
  theme_bw() +labs(x="Proportion Selected", y="R")

Next we describe how to estimate and plot the expected values for the response to selection. This is accomplished by the following R script which computes the expected response to selection (R) from the breeder's equation. In this example we will estimate R for the following proportions of selected individuals: 0.1, 0.2, 0.4, 0.6, 0.8 and 0.9.
The expected and observed values are plotted using ggplot2.

:::R

#Set your parameters

#heritability
h=0.5
#initial allele frequency
p=0.45
#effect size of the QTL
a=0.1
#proportion of selected individuals in the population
prop=c(0.1,0.2,0.4,0.6,0.8,0.9)




##estimate i using the hight of the ordinate of the truncation point in the normal ##distribution of the 
##phenotypic values

breeders<- function(h,p,a,prop){

get.i<-function(prop)
{
  z<-dnorm(qnorm(prop))
  # dnorm is height
  # quantile of the selected; sd = sigma_p
  i<-z/prop

  return(i)
}

i<- get.i(prop)


##estimate environmental variance
##for a given set of loci with heritability h, initial allele frequencies p
##and effect size a (vectors can be used)
get_ve<-function(h,p=c(),a=c()){
  q=1-p

  temp= 2*p*q*a^2
  vg=sum(temp)
  ve= (vg - vg*h) / h

  return(ve)
}

ve<-get_ve(h, rep(p, times=10), rep(a, times=10))


##estimate the standard phenotypic variance. Again p is the initial allele frequencies
##a the effect sizes and ve the environmental variance estimated above
get.sp<- function(p=c(),a=c(),ve){
  q=1-p

  temp= 2*p*q*a^2
  vg=sum(temp)
  sp= sqrt(ve+vg)

  return(sp)
}


sp<- get.sp(rep(p, times=10), rep(a, times=10), ve)



##estimate selection differential as a function of i(estimated above)
##and the standard phenotypic deviation
get.S<- function(i,sp){
  S=i*sp
  return(S)
}

S<- get.S(i,sp)


##estimate the response to selection as a function of the heritability and seelction ##differential
get.R<- function(h,S){
  R=h*S
  return(R)
}

R<- get.R(h, S)


print(R)
}

##How to run the function

R_theory=c()
for (i in 1:length(prop)){
  R_theory[i]<-as.numeric(breeders(h,p,a,prop[i]))
   }

R_theory<- data.frame(R=R_theory, Proportion_Selected=prop, heritability=0.5)


##Plot my results
ggplot(R_theory, aes(x=Proportion_Selected, y=R )) + geom_line() + 
  geom_boxplot(data=her05_R, aes(x=Proportion_Selected, y=R, group=Proportion_Selected), width=0.02, fill="forestgreen") + 
  facet_grid(.~heritability, scales="free_x", space = "free_x") + labs(x="Proportion Selected", y="R")+
    theme_bw()