In georgeshirreff/multiNe: Multiple Ne estimators

Estimation of Ne from multiple independent genealogies

This tutorial implements an extension of the methods described in http://www.auai.org/uai2012/papers/310.pdf to infer effective population size trajectories from multiple indepenedent genealogies.

Example

We will first simulate 10 genealogies with 50 tips from a bottleneck demographic scenario.

#install.packages("devtools")
library("devtools")
install_github("georgeshirreff/multiNe")
library("multiNe")
library("INLA")

bottleneck_traj<-function(t){
  result=rep(0,length(t))
  result[t<=0.5]<-1
  result[t>0.5 & t<1]<-.1
  result[t>=1]<-1
  return(result)
}

my.out<-simulate.tree(n=50,N=10,simulator="thinning",Ne=bottleneck_traj,max=11)
my.out
par(mfrow=c(1,1))
plot(my.out$out[[1]],show.tip.label=FALSE)
ape::axisPhylo()

For the first genealogy, the following function implemented in phylodyn uses the INLA method to infer Ne at 100 points:

par(mfrow=c(1,2))
res_BNPR<-phylodyn::BNPR(data=my.out$out[[1]],lengthout=100,prec_alpha=0.01,prec_beta=0.01)
phylodyn::plot_BNPR(res_BNPR,traj=bottleneck_traj,ylim=c(5,0.05),main="from one genealogy")

res_BNPR = list()
res_BNPR[[1]] = phylodyn::BNPR(data = my.out$out[[1]], lengthout = 100, prec_alpha = 0.01, prec_beta = 0.01)
table_res<-data.frame(time=res_BNPR[[1]]$x,effpop=res_BNPR[[1]]$effpop,effpop025=res_BNPR[[1]]$effpop025,effpop975=res_BNPR[[1]]$effpop975)
for (i in 2:length(my.out$out))
{
  res_BNPR[[i]] = phylodyn::BNPR(data = my.out$out[[i]], lengthout = 100, prec_alpha = 0.01, prec_beta = 0.01)

  table_res<-rbind(table_res,data.frame(time=res_BNPR[[i]]$x,effpop=res_BNPR[[i]]$effpop,effpop025=res_BNPR[[i]]$effpop025,effpop975=res_BNPR[[i]]$effpop975))

}

conf.int.skyline(table_res,plot_type="linear",epsilon=0.01,xlim=c(0.7,0),log="y",main="Harmonic Mean")
points(seq(0,.7,.01),bottleneck_traj(seq(0,.7,.01)),lty=2,type="l")

The previous figure shows the harmonic mean estimator from the 10 independently estimated genealogies. In the following, we estimate $N_{e}$ using INLA, assuming the 10 genealogies are 10 independent realizations of the coalescent process.

res<-BNPR_Multiple(my.out$out, lengthout=100, ntrees=10)
par(mfrow=c(1,1))
phylodyn::plot_BNPR(res,traj=bottleneck_traj,main="INLA from multiple genealogies")

Let's compare it with Skyline plot

conf_int_obj<-Phy2Sky(my.out$out,output_type = "conf.int")
conf.int.skyline(conf_int_obj,xlim=c(0.7,0),main="Harmonic Mean")

georgeshirreff/multiNe documentation built on May 17, 2019, 1:15 a.m.