-
Notifications
You must be signed in to change notification settings - Fork 79
Closed as not planned
Labels
enhancementNew feature or requestNew feature or requestfutureIssues that are closed as they are not planned in the medium-term, but which are still desirable.Issues that are closed as they are not planned in the medium-term, but which are still desirable.statistics
Description
It would be nice to have a function to calculate the selection test described in the Relate paper (p.16 of supplement) for mutations in a given tree. Or perhaps this is something that could be added to the documentation as an example, it's easy to code up using existing functions. This is the code I have (and an example).
import msprime
import math
import numpy as np
# Get number of lineages at time just above node
def get_num_lineages(tree, node):
if tree.is_sample(node):
count = tree.num_samples()
else:
nodes_ordered = tree.timedesc()
count = 1 + np.where(nodes_ordered == node)[0]
return count
# Calculate Relate selection test p-value
def selection_pvalue(tree, mutation):
node = mutation.node
if tree.is_sample(node):
# Can only compute this for mutations with frequency >= 2
p = None
else:
N = tree.num_samples()
fn = tree.num_samples(node)
# k is the number of lineages just after the edge carrying the mutation splits into 2
k = 1 + get_num_lineages(tree, node)[0]
p = 0.0
for f in range(fn, N - k + 3):
p += (f - 1) * math.comb(N - f - 1, k - 3) / math.comb(N - 1, k - 1)
return pts = msprime.sim_ancestry(
10,
population_size = 1,
sequence_length = 1,
ploidy = 1,
recombination_rate = 0,
discrete_genome = False,
random_seed = 1)
ts = msprime.sim_mutations(ts,
rate=3,
discrete_genome=False,
random_seed=1,
model=msprime.BinaryMutationModel())for m in ts.mutations():
print(m.id, m.site, ts.site(m.site).position, selection_pvalue(ts.first(), m))0 0 0.027387595968320966 None
1 1 0.09280081023462117 0.41666666666666663
2 2 0.2044522489886731 0.5952380952380951
3 3 0.23608897626399994 None
4 4 0.4306985700968653 1.0
5 5 0.5181525482330471 0.41666666666666663
6 6 0.7783892354927957 1.0
7 7 0.8463109182193875 0.7222222222222222
8 8 0.8650202453136444 0.41666666666666663
9 9 0.9325573612004519 None
10 10 0.939127794932574 1.0
silvewheat
Metadata
Metadata
Assignees
Labels
enhancementNew feature or requestNew feature or requestfutureIssues that are closed as they are not planned in the medium-term, but which are still desirable.Issues that are closed as they are not planned in the medium-term, but which are still desirable.statistics