## Monday, March 5, 2018

### Some new functions for teaching phenotypic evolution

For the Evolution class I'm teaching this semester I've been working on adding some basic population genetic functions to my GitHub R package PopGen.

Today, the first one is a function to compute Hardy-Weinberg frequencies, but for an arbitrary number of alleles.

First, here is what the function looks like:

``````hardy.weinberg<-function(p=c(0.5,0.5),alleles=c("A","a"),as.matrix=FALSE){
if(sum(p)!=1) p<-p/sum(p)
nalleles<-max(length(p),length(alleles))
if(nalleles!=length(p)) p<-rep(1/nalleles,nalleles)
if(nalleles!=length(alleles)) {
if(length(p)<=26) alleles<-LETTERS[1:length(p)]
else alleles<-paste("A",1:length(p),sep="")
}
p<-setNames(p,alleles)
HW<-as.matrix(p)%*%t(as.matrix(p))
if(as.matrix) return(HW)
else {
hw<-vector()
k<-1
for(i in 1:length(alleles)){
for(j in i:length(alleles)){
hw[k]<-if(i==j) HW[i,j] else 2*HW[i,j]
names(hw)[k]<-paste(alleles[c(i,j)],collapse="")
k<-k+1
}
}
return(hw)
}
}
``````

Now let's use it. Obviously, it is straightforward in the biallelic case:

``````hardy.weinberg(p=c(0.7,0.3))
``````
``````##   AA   Aa   aa
## 0.49 0.42 0.09
``````
``````## or
barplot(hardy.weinberg(p=c(0.7,0.3)))
`````` but it also works nicely for multi-allelic situations in which Hardy-Weinberg frequencies are more laborious to compute:

``````hw<-hardy.weinberg(p=c(0.4,0.3,0.2,0.1))
hw
``````
``````##   AA   AB   AC   AD   BB   BC   BD   CC   CD   DD
## 0.16 0.24 0.16 0.08 0.09 0.12 0.06 0.04 0.04 0.01
``````
``````barplot(hw)
`````` The second function is designed merely to show that for multilocus genotypes in which allelic substitutions have additive effect, we quickly converge on a (approximately) normal phenotypic distribution for even a modest number of loci. Furthermore, that this is true regardless of our allele frequencies.

The way the function works is by computing all multi-locus genotypes and their HW frequencies, and then plotting the resultant phenotypic trait distribution.

Here's the function:

``````phenotype.freq<-function(nloci=10,p=NULL,effect=1){
if(is.null(p)) p<-rep(0.5,nloci)
genotypes<-t(apply(cbind(p,1-p),1,hardy.weinberg))
COMBN<-permutations(n=3,r=nloci,set=T,repeats.allowed=T)
PHEN<-rowSums(COMBN-2)*effect
FREQ<-rep(1,nrow(COMBN))
for(i in 1:nrow(COMBN)){
for(j in 1:nloci)   FREQ[i]<-FREQ[i]*genotypes[j,COMBN[i,j]]
}
phen<-unique(PHEN)
freq<-rep(0,length(phen))
for(i in 1:length(phen)) freq[i]<-sum(FREQ[which(PHEN==phen[i])])
plot(phen,freq,type="b",pch=21,bg="grey",
xlab="phenotypic trait value",ylab="relative frequency",
cex=1.5,xlim=range(phen)+c(-0.5,0.5)*effect)
for(i in 1:length(phen))
rect(phen[i]-0.4*effect,0,phen[i]+0.4*effect,
freq[i],border="grey",
col=make.transparent("blue",0.2))
}
``````

and here's what I mean using random allele frequencies at each locus (note that this function uses gtools & phytools):

``````library(gtools)
library(phytools)
## 1 locus
nloci<-1
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)
`````` ``````## 2 loci
nloci<-2
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)
`````` ``````## 4 loci
nloci<-4
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)
`````` ``````## 8 loci
nloci<-8
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)
`````` ``````## 12 loci
nloci<-12
phenotype.freq(nloci,p=runif(n=nloci),effect=1/nloci)
`````` Neat? I kind of think so.

Note that because the function finds all possible multilocus genotypes for `nloci` greater than 12 or 14 it is not going to work!