30:T409,chart = { if (n_sims > 50) { return html`Chart is not displayed if the number of simulations is greater than 50.`; } return Plot.plot({ width: width, color: { legend: false, scheme: "Set1", domain: ["A", "C", "G", "T"] }, x: { label: "length" }, marks: [ Plot.barX(data, { x: "width", y: "simulation", fill: "base", render: (i, s, v, d, c, next) => { const g = next(i, s, v, d, c); d3.select(g) .selectAll("rect") .attr("style", "cursor:pointer") .on("click", (e, d) => { show_item(e, d, g); }) .transition(500) .attr("opacity", 1); return g; } }), Plot.textX( data, Plot.stackX( { x: "sum" }, { x: "width", y: "simulation", text: "base", fill: "white", fontWeight: 600 } ) ), Plot.ruleX([0]) ] }); }31:T6c6,item_explanations = (i) => { const item = simulation_results.flat()[i]; let out = [ `We start at length ${item.start_length.toFixed( 3 )} with base ${item.base}.` ]; out.push( `The time until next mutation event for this base follows an exponential distribution of parameter λ=${item.mutation_rate}.` ); out.push( `By randomly sampling a value from this distribution, we simulate that the time until next mutation event will be ${item.mutation_time.toFixed( 3 )}.` ); let mutation_length = item.start_length + item.mutation_time; out.push( `This means that the next mutation will occur at length ${mutation_length.toFixed( 3 )}.` ); if (mutation_length > branch_length) { out.push( `As ${mutation_length.toFixed( 3 )} is greater than the branch length, no mutation is added.` ); return html`${out.join("
")}`; } const next_bases = Object.keys(item.base_rates).join(", "); const next_rates = Object.values(item.base_rates).join(", "); out.push( `The mutation rates for ${next_bases} are respectively ${next_rates}` ); const next_probs = Object.values(item.base_probs) .map((d) => d.toFixed(3)) .join(", "); out.push( `This means that the probabilities for the next mutation to be ${next_bases} are respectively ${next_probs}` ); out.push( `By randomly sampling from a uniform distribution, we simulate that the base after mutation will be ${item.mutation_base}.` ); return html`${out.join("
")}`; }32:T455,run_simulation = () => { let steps = []; let i = 0; mutable selected_item = null; let mutation_base = starting_base == "Random" ? random_base() : starting_base; let current_length = 0; while (current_length <= branch_length) { const start_base = mutation_base; const rate = -rates.get(start_base, base_rank(start_base)); const mutation_time = jstat.exponential.sample(rate); const start_length = current_length; let base_rates = null, base_probs = null; current_length += mutation_time; if (current_length <= branch_length) { let probs = compute_probs(start_base, rate); base_rates = probs.base_rates; base_probs = probs.base_probs; mutation_base = compute_new_base(base_probs); } else { mutation_base = null; } steps.push({ step: i++, base: start_base, start_length: start_length, end_length: current_length, mutation_rate: rate, mutation_time: mutation_time, base_rates: base_rates, base_probs: base_probs, mutation_base: mutation_base }); } return steps; }33:[] 2e:{"title":"Site simulations","cells":[{"id":0,"value":"# Site simulations\n","pinned":false,"hidden":false,"mode":"md"},{"id":302,"value":"## Parameters","pinned":false,"hidden":false,"mode":"md"},{"id":299,"value":"Matrix of substitution rates :","pinned":false,"hidden":false,"mode":"md"},{"id":11,"value":"Inputs.table(rates, {\n //layout: \"auto\",\n width: 400,\n format: {\n From: (x) => html`From ${x}`,\n A: format_rates,\n C: format_rates,\n G: format_rates,\n T: format_rates\n },\n header: { From: \"\", A: \"To A\", C: \"To C\", G: \"To G\", T: \"To T\" }\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":70,"value":"viewof starting_base = Inputs.radio([\"A\", \"C\", \"G\", \"T\", \"Random\"], {\n label: \"Starting base\",\n value: \"A\"\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":74,"value":"viewof branch_length = Inputs.range([0.1, 10], {\n label: \"Branch length\",\n value: 0.3,\n step: 0.1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":495,"value":"viewof n_sims = Inputs.number([1, 10000], {\n value: 1,\n label: \"Number of simulations\",\n step: 1\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":80,"value":"viewof start = Inputs.button(\"Rerun simulations\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":217,"value":"## Simulations results","pinned":false,"hidden":false,"mode":"md"},{"id":396,"value":"$30","pinned":false,"hidden":false,"mode":"ojs"},{"id":569,"value":"item_info = {\n if (n_sims > 50) {\n return html``;\n }\n\n let content = `Click an item to see calculation details.`;\n if (selected_item !== null) {\n content = item_explanations(selected_item);\n }\n return html`

${content}

`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":237,"value":"html`${show_overview()}`","pinned":false,"hidden":false,"mode":"ojs"},{"id":687,"value":"bases_chart = {\n if (n_sims < 10) {\n return html``;\n }\n return Plot.plot({\n title: \"Distribution of the last base at branch length\",\n height: starting_base == \"Random\" ? 600 : 300,\n x: { domain: [\"A\", \"C\", \"G\", \"T\"], label: \"Last base\" },\n y: {\n domain: [0, 100],\n percent: true,\n grid: true,\n label: \"Pourcentage (%)\"\n },\n fy: { label: \"Start base\" },\n color: { domain: [\"A\", \"C\", \"G\", \"T\"], scheme: \"Set1\" },\n marks: [\n Plot.barY(last_bases, {\n x: \"last_base\",\n y: \"prop\",\n fill: \"last_base\",\n fy: \"start_base\"\n }),\n Plot.text(last_bases, {\n x: \"last_base\",\n y: \"prop\",\n fy: \"start_base\",\n dy: -8,\n text: (d) => (d.prop * 100).toFixed(1) + \"%\"\n }),\n Plot.ruleY([0])\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":742,"value":"transition_chart = {\n if (n_sims < 10) {\n return html``;\n }\n return Plot.plot({\n title: \"Matrix of transition probabilities\",\n height: 200,\n grid: true,\n y: { domain: [\"A\", \"C\", \"G\", \"T\"], label: \"Start base\" },\n x: { domain: [\"A\", \"C\", \"G\", \"T\"], label: \"Last base\" },\n color: { scheme: \"Blues\", domain: [0, 1] },\n marks: [\n Plot.cell(last_bases, { y: \"start_base\", x: \"last_base\", fill: \"prop\" }),\n Plot.text(last_bases, {\n y: \"start_base\",\n x: \"last_base\",\n text: (d) => (d.prop * 100).toFixed(1) + \"%\",\n fill: (d) => (d.prop > 0.6 ? \"white\" : \"black\")\n })\n ]\n });\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":305,"value":"html`

`","pinned":false,"hidden":false,"mode":"ojs"},{"id":58,"value":"# Code\n","pinned":false,"hidden":false,"mode":"md"},{"id":563,"value":"### States","pinned":false,"hidden":false,"mode":"md"},{"id":487,"value":"mutable selected_item = null","pinned":false,"hidden":false,"mode":"ojs"},{"id":266,"value":"### Values","pinned":false,"hidden":false,"mode":"md"},{"id":26,"value":"rates = aq.table({\n From: [\"A\", \"C\", \"G\", \"T\"],\n A: [-1.916, 0.148, 0.286, 0.525],\n C: [0.541, -1.069, 0.17, 0.236],\n G: [0.787, 0.415, -0.591, 0.594],\n T: [0.588, 0.506, 0.135, -1.355]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":326,"value":"height = 100","pinned":false,"hidden":false,"mode":"ojs"},{"id":329,"value":"margin = ({ top: 30, right: 20, bottom: 30, left: 40 })","pinned":false,"hidden":false,"mode":"ojs"},{"id":332,"value":"data = {\n let sim = 1;\n let results = [];\n for (let sim_data of simulation_results) {\n let res = sim_data.map((d) => ({\n simulation: sim,\n step: `Step ${d.step}`,\n width: d3.min([d.end_length, branch_length]) - d.start_length,\n base: d.base\n }));\n results.push(res);\n sim += 1;\n }\n return results.flat();\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":61,"value":"### Functions","pinned":false,"hidden":false,"mode":"md"},{"id":583,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":480,"value":"show_item = (e, d, g) => {\n let el = d3.select(e.currentTarget);\n if (el.attr(\"selected\") == 1) {\n d3.select(g).selectAll(\"rect\").attr(\"opacity\", 1).attr(\"selected\", 0);\n mutable selected_item = null;\n } else {\n d3.select(g).selectAll(\"rect\").attr(\"opacity\", 0.3).attr(\"selected\", 0);\n el.attr(\"opacity\", 1);\n el.attr(\"selected\", 1);\n mutable selected_item = d;\n }\n return d;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":244,"value":"show_overview = (d) => {\n if (n_sims == 1) {\n let simulation_mutations = simulation_results[0].length - 1;\n return html`The simulation generated ${simulation_mutations} mutation${\n simulation_mutations == 1 ? \"\" : \"s\"\n }.`;\n }\n if (n_sims > 1) {\n let simulation_mutations = simulation_results.map((d) => d.length - 1);\n return html`On average, the simulations generated ${d3\n .mean(simulation_mutations)\n .toFixed(2)} mutations.`;\n }\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":82,"value":"simulation_results = {\n start;\n let sim = 0;\n let sims = [];\n while (sim < n_sims) {\n sims.push(run_simulation());\n sim += 1\n }\n return sims;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":508,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":178,"value":"compute_new_base = (probs) => {\n const random = jstat.uniform.sample(0, 1);\n const cum_probs = d3.cumsum(Object.values(probs));\n let i = 0;\n while (i < cum_probs.length) {\n if (cum_probs[i] > random) {\n break;\n }\n i = i + 1;\n }\n return Object.keys(probs)[i];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":136,"value":"compute_probs = (base, rate) => {\n let base_rates = rates\n .params({ base: base })\n .filter((d, $) => d.From == $.base)\n .select(aq.not(\"From\"))\n .reify()\n .objects()[0];\n base_rates = filter_object(base_rates, (d) => d >= 0);\n let base_probs = structuredClone(base_rates);\n for (const key in base_probs) {\n base_probs[key] = base_probs[key] / rate;\n }\n return {base_rates: base_rates, base_probs: base_probs}\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":63,"value":"format_rates = function (x) {\n return htl.html`

${x}`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":157,"value":"filter_object = (obj, predicate) => {\n return Object.keys(obj)\n .filter((key) => predicate(obj[key]))\n .reduce((res, key) => ((res[key] = obj[key]), res), {});\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":749,"value":"last_bases.view()","pinned":true,"hidden":false,"mode":"ojs"},{"id":681,"value":"last_bases = aq\n .from(\n simulation_results.map((sim) => ({\n start_base: sim[0].base,\n last_base: sim[sim.length - 1].base\n }))\n )\n .groupby(\"start_base\", \"last_base\")\n .count()\n .derive({ prop: (d) => d.count / op.sum(d.count) })","pinned":false,"hidden":false,"mode":"ojs"},{"id":725,"value":"random_base = () => [\"A\", \"C\", \"G\", \"T\"][Math.floor(Math.random() * 4)]","pinned":false,"hidden":false,"mode":"ojs"},{"id":56,"value":"### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":92,"value":"base_rank = (x) => [\"A\", \"C\", \"G\", \"T\"].indexOf(x)","pinned":false,"hidden":false,"mode":"ojs"},{"id":39,"value":"import { aq, op } from \"@uwdata/arquero\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":85,"value":"jstat = require(\"jstat\")","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W33","stdlib":"1","titling":"auto"}