Skip to content

Add selection test p-values #2190

@a-ignatieva

Description

@a-ignatieva

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 p
ts = 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())

pic1

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    enhancementNew feature or requestfutureIssues that are closed as they are not planned in the medium-term, but which are still desirable.statistics

    Type

    No type

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions