Mutations vs recombinations to resolve recent ancestry

[hyanwong]

A quick check to see if we have more information in shared recombinations than in mutations for e.g. doubletons in a large simulation. The result for 40, 000 samples is

Num doubleton nodes with a mutation 2055 / 5728
Num doubleton nodes with > 1 parent 2968 / 5728

Note that this is not quite reasonable, because even if the nodes have > 2 parents, we might not be able to spot it if the parents have identical sequences on either the L or the R of the breakpoint. I suspect this is quite unlikely, though (I have not come up with an efficient way of checking this).

It should, I think, be possible to infer this sort of information from real data. I'm not quite sure how, though....

import msprime
import numpy as np
import tqdm

gts = msprime.sim_ancestry(
    20_000,  # A reasonably large human dataset
    population_size=1e4,
    sequence_length=5e7,  # 50Mb
    recombination_rate=1e-8,
    model=[
        msprime.DiscreteTimeWrightFisher(duration=500),
        msprime.StandardCoalescent(),
    ]
)
ts = msprime.sim_mutations(gts, rate=1e-8)
uniq_child_parent = np.unique(np.column_stack((ts.edges_child, ts.edges_parent)), axis=0)
nd, count = np.unique(uniq_child_parent[:, 0], return_counts=True)
multipar = set(np.where(count > 1)[0])

more_than_doubleton = np.zeros(ts.num_nodes, dtype=bool)
sometimes_doubleton = np.zeros(ts.num_nodes, dtype=bool)
# find nodes which are only doubletons
for t in tqdm.tqdm(ts.trees(), total=ts.num_trees):
    for u in t.parent_array[ts.samples()]:
        ns = t.num_samples(u)
        if ns > 2:
            more_than_doubleton[u] = True
        if ns == 2:
            sometimes_doubleton[u] = True

doubletons = np.logical_and(
    np.logical_not(more_than_doubleton), sometimes_doubleton)
print(
    "Num doubleton nodes with a mutation",
    np.sum(np.isin(np.where(doubletons)[0], ts.mutations_node)),
    "/",
    np.sum(doubletons)
)
print(
    "Num doubleton nodes with > 1 parent",
    np.sum(np.isin(np.where(doubletons)[0], np.where(count > 1)[0])),
    "/",
    np.sum(doubletons)
)

[hyanwong]

Presumably this is essentially a feature of whether a branch above a node is involved in a recombination or a mutation, both of which are proportional to branch length under the coalescent / neutral mutation model

[hyanwong]

Oops - there's a bug in the code. I use num_children_array, whereas what I really wanted was num_samples_array (which doesn't exist)