The second quantitative trait code post. The set up is the same as in the last QT post. I simulate parents with L genotypes drawn binomially from the population, and then make kids from these parents.

This one makes the parental mid-point vs child phenotype regression plot. You get to specify the number of loci (L) and environmental variance (environ.var), and genetic.var=1. The slope of the regression of child_pheno ~ (mum_pheno+dad_pheno)/2 is h^2, i.e. genetic.var/(environ.var+genetic.var). What I find really quite neat is how few loci you need for the h^2 regression to work.

Ideas for this year: Use height loci? Add in dominance and other relationships, or get the students to do this. They struggle to understand the dominance variance (and to be honest I always do too).

This slideshow requires JavaScript.

##Paste this in first
par.off.corr<-function(L=20, environ.var,print.slope=FALSE){
##Quantitative genetics sims
allele.freq<-0.5   ###each locus is assumed to have the same allele frequencies. This is just to simplify the coding, in reality these results work when each locus has its own frequency (and the coding wouldn't be too much harder). 
 
Num_inds=1000 

##MAKE A MUM
## For each mother, at each locus we draw an allele (either 0 or 1) from the population allele frequency. 
##We do this twice for each mother two represent the two haplotypes in the mother 
mum.hap.1<-replicate(Num_inds, rbinom(L,1,allele.freq) )
mum.hap.2<-replicate(Num_inds, rbinom(L,1,allele.freq) )
##type mum.hap.1[,1] to see the 1st mothers 1st haplotype

##Each mothers genotype at each locus is either 0,1,2
mum.geno<-mum.hap.1+mum.hap.2

additive.genetic<-colSums(mum.geno)
genetic.sd<-sd(additive.genetic)
mean.genetic<-mean(additive.genetic)

additive.genetic<-additive.genetic / sd(additive.genetic)
mum.pheno<- additive.genetic + rnorm(Num_inds,sd=sqrt(environ.var))
mum.pheno<-mum.pheno-mean(mum.pheno)



##MAKE A DAD (same code as make a mum, only said in a deeper voice)
dad.hap.1<-replicate(Num_inds, rbinom(L,1,allele.freq) )
dad.hap.2<-replicate(Num_inds, rbinom(L,1,allele.freq) )
dad.geno<-dad.hap.1+dad.hap.2


additive.genetic<-colSums(dad.geno)
additive.genetic<-additive.genetic / sd(additive.genetic)
dad.pheno<- additive.genetic + rnorm(Num_inds,sd=sqrt(environ.var))
dad.pheno<-dad.pheno-mean(dad.pheno)

### Make a child
child.geno<-dad.hap.1+mum.hap.1 ##1 haplotype from mum 1 haplotype from dad

additive.genetic<-colSums(child.geno)
additive.genetic<-additive.genetic / sd(additive.genetic)
child.pheno<- additive.genetic + rnorm(Num_inds,sd=sqrt(environ.var))
child.pheno<-child.pheno-mean(child.pheno)



##Calculate midpoints, linear model and plots

parental.midpoint<-(mum.pheno+dad.pheno)/2 ##avg. parents

lm(child.pheno~parental.midpoint) ##linear model between child and mid point
							##the slope of this line is the narrow sense heritability.

# plot parental midpoint against offsprings phenotype.
layout(1) ###done in case this is run after the code with 3 plots
plot(parental.midpoint,child.pheno,xlab="Parental midpoint",ylab="Child's phenotype",cex=1.5,main=paste("L =",L,"environ.var=",environ.var,"genetic.var=1"))
## plot the regression in red
abline(h=0,col="grey",lwd=2)
abline(v=0,col="grey",lwd=2)
 abline(lm(child.pheno~parental.midpoint),col="red",lwd=2)

if(print.slope) text(x=-2,y=3,label=paste("slope= ",format(lm(child.pheno~parental.midpoint)$coeff[2],digit=3)),col="red",lwd=4)

abline(0,1,col="blue",lwd=3)
 
 }
 
for(environ.var in c(0.00001,0.5,1,2)){
png(file=paste("parent_offspring_regression_varg_1_vare_",format(environ.var,dig=1),".png",sep=""))
 par.off.corr(L=200, environ.var=environ.var,print.slope=TRUE)
dev.off()
 }
 

This work is licensed under a Creative Commons Attribution 3.0 Unported License.