Skip to content

Commit

Permalink
Add interactive applet on diversity along a recombining region (#12)
Browse files Browse the repository at this point in the history
  • Loading branch information
molpopgen authored Nov 29, 2021
1 parent e026eb4 commit 4d9ae23
Show file tree
Hide file tree
Showing 2 changed files with 38 additions and 0 deletions.
1 change: 1 addition & 0 deletions book/_toc.yml
Original file line number Diff line number Diff line change
Expand Up @@ -12,6 +12,7 @@ parts:
- file: markdown/kingman_coalsecent/model_for_data
- file: markdown/kingman_coalsecent/trees
- file: markdown/kingman_coalsecent/trees_with_recombination
- file: markdown/kingman_coalsecent/recombination_diversity
- file: markdown/kingman_coalsecent/mutation_frequencies
- file: markdown/kingman_coalsecent/expectations
- file: markdown/kingman_coalsecent/distributions
Expand Down
37 changes: 37 additions & 0 deletions book/markdown/kingman_coalsecent/recombination_diversity.md
Original file line number Diff line number Diff line change
@@ -0,0 +1,37 @@
---
jupytext:
formats: md:myst
text_representation:
extension: .md
format_name: myst
kernelspec:
display_name: Python 3
language: python
name: python3
---

# Variation in diversity along a recombining genome

```{code-cell} python
import msprime
import numpy as np
import panel as pn
import holoviews as hv
hv.extension('bokeh', 'matplotlib')
def branch_diversity_along_tree(rho, n, sequence_length=1e6, num_windows=500):
wins = np.linspace(0, sequence_length, num_windows)
mids = (wins[1:] - wins[:1])/2
ts = msprime.sim_ancestry(n, ploidy=1,
sequence_length=sequence_length,
recombination_rate=rho/4/sequence_length)
assert ts.num_samples == n
diversity = ts.diversity(windows=wins, mode="branch").flatten()
return hv.Curve((mids, diversity)).opts(tools=["box_select"]).redim(y=hv.Dimension('y',
range=(0, 1.1*diversity.max()))).opts(xlabel="Window midpoint", ylabel="Diversity (branch stat)")
divplot = hv.DynamicMap(branch_diversity_along_tree,
kdims=['rho', 'n']).redim.range(rho=(0, 1000), n=(10, 100)).opts(framewise=True, width=500)
divplot
```

0 comments on commit 4d9ae23

Please sign in to comment.