hyanwong
commented
8 months ago And it checks out! With ba90588f0d41c30859a00a9baa7b2fd5e381c012 I have added a test. Alternatively you can create the plot below with the following code:
import numpy as np
from matplotlib import pyplot as plt
from tests.test_gigutil import TestDTWF_one_break_no_rec_inversions_slow
from tests.test_gigutil import INVERTED_CHILD_FLAG
test = TestDTWF_one_break_no_rec_inversions_slow()
test.default_gens = 10
#test.seq_len=321 # seed 321 bugs out: need to investigate
test.test_inversion()
ts = test.ts.simplify()
inversion_above_node = np.where(ts.nodes_flags == INVERTED_CHILD_FLAG)[0]
inversion_samples = list(ts.at(150).samples(inversion_above_node[0]))
non_inversion_samples = np.zeros(ts.num_nodes, dtype=bool)
non_inversion_samples[ts.samples()] = True
non_inversion_samples[inversion_samples] = False
non_inversion_samples = np.where(non_inversion_samples)[0]
max_divergence = ts.max_time * 2
breaks = list(ts.breakpoints())
d = ts.divergence([inversion_samples, non_inversion_samples], mode="branch", windows="trees")
plt.stairs(d, breaks)
plt.ylim(max_divergence*0.95, max_divergence*1.01)
plt.xlabel("Genome pos")
plt.ylabel(f"Av genetic divergence between samples with &\nwithout the inversion (max={max_divergence})")
plt.axvline(test.inversion.left, c="green")
plt.axvline(test.inversion.right, c="red")
plt.text((test.inversion.right + test.inversion.left)/2, max_divergence, "inversion", ha="center")
plt.savefig("div.png")
77 added a simple inversion simulator, here:
https://github.com/hyanwong/GeneticInheritanceGraphLibrary/blob/77d3b1db818c97e2ddd10f2a6948435d1340a7a8/tests/test_gigutil.py#L199
We should test this is behaving as expected (e.g. like a slim simulation with inversions). For example, we can output the decapitated TS and see if there are any recombinations in the inverted region between individuals with and without the inversion (there shouldn't be)
Note that this framework should also allow simulation of nested inversions, which isn't possible in SLiM, I think?