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')
One platform to build and deploy the best data apps
Experiment and prototype by building visualizations in live JavaScript notebooks. Collaborate with your team and decide which concepts to build out.
Use Observable Framework to build data apps locally. Use data loaders to build in any language or library, including Python, SQL, and R.
Seamlessly deploy to Observable. Test before you ship, use automatic deploy-on-commit, and ensure your projects are always up-to-date.