In ndukler/epiAllele: Tests for signals of variable rates of evolution between gene classes and lineages for probabalistic alleles

knitr::opts_chunk$set(echo = TRUE)

In order to simulate data we must first create a tree and set edge lengths to desired values

library(epiAllele)
library(ape)
library(cowplot)

set.seed(123)
### Create a tree
trTrue <- unroot(rtree(n = 5, br = runif))
trTrue=reorder.phylo(trTrue,"postorder")
trTrue$edge.length=sample(seq(0.1,1,by=0.1),size = length(trTrue$edge.length),replace = TRUE)
plot(trTrue)

Now we simulate from the true tree with a selected stationary distibution for allele frequency.

sData=simData(nSites = 100000,tr=trTrue,rate = 1,pi = c(0.3,0.7))

Now we create a copy of the true tree where all branch lengths are one for ease of comparing the fitted rates. We place the simulated data and the copy of the tree in an alleleData object and then a rateModel. We can check that the rate assignments are correct by plotting the rate model object.

trCopy=trTrue
trCopy$edge.length=rep(1,length(trCopy$edge.length))
ad=alleleData(data = sData$data,tree = trCopy)
eTab=getEdgeTable(ad)
eTab[,edgeGroup:=paste0("e",1:length(edgeID))]
fullMod=rateModel(data=ad,lineageTable=eTab)
plotTree(fullMod)

Now lets fit the model.

fittedTree=fit(fullMod,control = list(maxit=500))
plotFun = function() {plot(trTrue,show.tip.label = FALSE);tiplabels()}
g = plotTree(fittedTree$model,"index",TRUE)
plot_grid(plotFun,g)

We can also extract the scaled tree for another use if we want:

trScaled=exportScaledTree(fittedTree$model)
elCompare=matrix(0,ncol=2,nrow=nrow(trTrue$edge))
ind=1
for(i in trTrue$edge[,2]){
  elCompare[ind,1]=trTrue$edge.length[trTrue$edge[,2]==i]
  elCompare[ind,2]=trScaled[[1]]$edge.length[trScaled[[1]]$edge[,2]==i]
  ind=ind+1
}
r=cor(elCompare[,1],elCompare[,2])
plot(elCompare[,1],elCompare[,2],xlim = c(0,1),ylim = c(0,1))
abline(a=0,b=1)

ndukler/epiAllele documentation built on May 5, 2019, 4:50 p.m.