Published
Edited
Feb 2, 2021
Importers
coalesce =function(height,nodes){
const parentNode = {height:height,children:nodes,id:uuid.v4()}
return parentNode;
}
class PopulationDemographic{
constructor({N0}){
this.N0=N0;
}
getT(tau){
throw Error("don't call base class methods")
}
getIntegral(t0,t1){
throw Error("don't call base class methods")
}
getDemographic(t){
throw Error("don't call base class methods")
}
}
class ExponetialPopulation extends PopulationDemographic{
constructor(options){
super(options);
this.g= options.g;
}
getT(tau){
if(this.g===0){
return tau
}
return (1/this.g)*(Math.log(1+this.g*tau))
}
getTau(t){
return (1/this.g) * (Math.exp(this.g*t)-1)
}
getIntegral(t0,t1){
if(this.g===0){
return (t1-t0)*(this.N0)
}else{
return 1/(this.N0*this.g)*(Math.exp(this.g*t1)-Math.exp(this.g*t0))
}
}
getDemographic(t){
return this.N0*Math.exp(-this.g*t)
}
}
class ConstantPopulation extends PopulationDemographic{
constructor(options){
super(options);
}
getT(tau){
return tau;
}
getTau(t){
return t;
}
getIntegral(t0,t1){
return (t1-t0)*(1/this.N0)
}
getDemographic(t){
return this.N0;
}
}
simulation=populationDemographicModel=>{
let currentHeight=0;
const helper = function(population){
//Filter population to just nodes that are in play and can coalesce
const activePopulation = population.filter(n=>n.height<=currentHeight);
// these are tips that have yet to be sampled.
const toBeSampled = population.filter(n=>n.height>currentHeight)
// this is the next tip to be sampled (if it exists)
const nextSample = toBeSampled.length>0?toBeSampled[d3_array.minIndex(toBeSampled,d=>d.height)]:null;
// If the population only has one node return the node. it is the root. We are done.
if(population.length===1){
const tree= new Figtree.Tree(population[0],{lengthsKnown:false,heightsKnown:true});
//convert to new times
tree.nodeList.forEach(n=>n.height=populationDemographicModel.getT(n.height));
return tree;
}
//If there is only 1 active population then we can't coalesce to move up to the next sample time
if(activePopulation.length===1 && nextSample){
currentHeight=nextSample.height;
return helper(population)
}
// if we get this far then there are nodes that can coalsce
const k = activePopulation.length;
const Ne = populationDemographicModel.N0;
// draw the wait time for the next coalescent event. the event will happen at currentHeight+tau
const tau = sampleExponetial(k*(k-1)/(2*Ne));
// If there is a node to sample still and this coalecent event happend after that node would have been sampled reject the event, move up to the sample time and try again
if(nextSample && currentHeight+tau>nextSample.height){
//didn't happen in time
currentHeight=nextSample.height;
return helper(population);
}else{
// there is a coalsenct event
// pick nodes to coalesce
const luckyDucks = sampleLineages(activePopulation);
// get time of coalescent
currentHeight+=tau;
// make parent node
const parent = coalesce(currentHeight,luckyDucks);
// replace the children nodes with their parent in the population and go again
const updatedPopulation = population.filter(n=>!luckyDucks.includes(n)).concat(parent);
return helper(updatedPopulation);
}
}
return helper;
}
function simulateCoalescentTree({demographic,maxTipHeight,numberOfTips}){
return compose(simulation(demographic),makeTips(demographic)(maxTipHeight))(numberOfTips);
}
function sampleLineages(lineages){
const samples = [lineages[Math.round(Math.random()*(lineages.length-1))]];
const remainingLineages = lineages.filter(l=>!samples.includes(l))
return samples.concat(remainingLineages[Math.round(Math.random()*(remainingLineages.length-1))])
}
sampleExponetial=(lambda)=>{
return -1*Math.log(Math.random())/lambda
}
makeTips=populationDemographic=>(maxSampleHeight)=>numberofTips=>{
function makeTipsHelper(tips){
const heightDraw = Math.random()*maxSampleHeight;
// get Tau because the simulation will be run in the imaginary time tau and then converted to real
// time. This why the maxSampleHieght can be given in real time.
const newTips = tips.concat({height:populationDemographic.getTau(heightDraw),name:uuid.v4()})
if(newTips.length<numberofTips){
return makeTipsHelper(newTips)
}else{
return newTips;
}
}
return makeTipsHelper([{height:0,name:uuid.v4()}])
}
simtree = simulateCoalescentTree({demographic:new ExponetialPopulation({N0:100,g:2}),maxTipHeight:1,numberOfTips:10})
figure = {
const treeSVG = document.getElementById('figtree');
const layout = new Figtree.RectangularLayout(simtree);
const figTree = new Figtree.FigTree(treeSVG, layout, { top: 10, bottom: 70, left: 10, right: 150},
{ hoverBorder: 5, backgroundBorder: 2,
baubles: [ new Figtree.CircleBauble({r:20,vertexFilter:v=>v.degree===1})],
styles:{"nodes":{"fill":(d)=>colorRange[2]},
"branches":{"stroke":(d)=>colorRange[9],
"stroke-width":d=>4}}},
)
figTree.draw()
return figTree;
}
getIntervalList=(tree)=>{
let k=tree.nodeList.filter(n=>n.height===0).length
return tree.nodeList.sort((a,b)=>a.height-b.height).map((n,i,nodes)=> {if(n!==tree.rootNode){
const interval= {t0:n.height,
t1:nodes[i+1].height,
k:k,
type:(nodes[i+1].children?"coalescent":"sampling")};
interval.type==="coalescent"?k-=1:k+=1
return interval;
}
return null;
}).filter(n=>n)
}
coalescentLogLikelihood=(PopulationDemographic)=>intervals=>{
let logL =0;
for(const interval of intervals){
const intervalArea = PopulationDemographic.getIntegral(interval.t0,interval.t1);
const kChoose2 = interval.k*(interval.k-1)/2
logL += -kChoose2 * intervalArea;
if(interval.type==="coalescent"){
// logL += Math.log(kChoose2*1/PopulationDemographic.getDemographic(interval.t1));
// logL += Math.log(kChoose2)
// should we include the above term?
logL -= Math.log(PopulationDemographic.getDemographic(interval.t1))
}
}
return logL;
}
coalescentLogLikelihood(new ExponetialPopulation({N0:500,g:1}))(getIntervalList(simtree))
coalescentLogLikelihood(new ExponetialPopulation({N0:100,g:2}))(getIntervalList(simtree))
colorRange = [...d3.schemeTableau10,...d3.schemeSet3]
md`# libraries`
uuid = import('https://jspm.dev/uuid')
d3 = require('d3')
Figtree = import("https://cdn.jsdelivr.net/gh/rambaut/figtree.js@a92caa7/dist/figtree.esm.js")
// Figtree =import("https://cdn.jsdelivr.net/gh/rambaut/figtree.js@5799a65/dist/figtree.esm.js")
d3_array=require('d3-array')
Purpose-built for displays of data
Observable is your go-to platform for exploring data and creating expressive data visualizations. Use reactive JavaScript notebooks for prototyping and a collaborative canvas for visual data exploration and dashboard creation.
