hyanwong
commented
8 months ago I imagined that the effect would be made worse by mispolarization of ancestral states. However, a brief test indicates that this is probably not the case:
import tsinfer
import tsdate
import numpy as np
import stdpopsim
species = stdpopsim.get_species("HomSap")
model = species.get_demographic_model("AmericanAdmixture_4B11")
contig = species.get_contig("chr22", mutation_rate=model.mutation_rate, length_multiplier=0.05)
engine = stdpopsim.get_engine("msprime")
true_data = {}
inferred_data = {}
truedate_data = {}
infdated_data = {}
n_reps = 4
s = 12
mispolarise_proportions = [0, 0.01, 0.1]
for mispolarise_proportion in mispolarise_proportions:
samples = {"AFR": s, "EUR": s, "ASIA": s, "ADMIX": s}
num_samples = sum(samples.values()) * 2
true_data[mispolarise_proportion] = np.zeros(num_samples)
inferred_data[mispolarise_proportion] = np.zeros(num_samples)
infdated_data[mispolarise_proportion] = np.zeros(num_samples)
for rep in range(n_reps):
sim_ts = engine.simulate(model, contig, samples, seed=rep+123)
sd = tsinfer.SampleData.from_tree_sequence(sim_ts, use_sites_time=False)
aa = sd.sites_ancestral_allele[:] # perfect polarisation has these all 0
to_switch = np.random.choice(np.arange(len(aa)), round(mispolarise_proportion * len(aa)), replace=False)
aa[to_switch] = 1 # assume all sites have 2 or more alleles: just switch to the first alt
sd_mispol = sd.copy() # make editable
sd_mispol.sites_ancestral_allele[:] = aa
sd_mispol.finalise()
print(sim_ts.num_samples, "samples,", sim_ts.num_trees, "trees,", sim_ts.num_sites, "sites", mispolarise_proportion, "mp")
iTS = tsinfer.infer(sd_mispol).simplify()
dTS = tsdate.variational_gamma(
tsdate.util.split_disjoint_nodes(iTS.simplify()),
mutation_rate=model.mutation_rate,
population_size=sim_ts.diversity() / (4 * model.mutation_rate),
)
dTS = tsinfer.infer(tsinfer.SampleData.from_tree_sequence(dTS, use_sites_time=True)).simplify()
for ts, data in zip(
(sim_ts, iTS, dTS),
(true_data, inferred_data, infdated_data),
):
sample_set = set([u for u in ts.samples()])
d = np.zeros(num_samples)
tot = 0
for tree in ts.trees():
if tree.num_edges > 0:
tot += tree.interval.span
n_root_attached = len(set(tree.children(tree.root)) & sample_set)
d[n_root_attached] += tree.interval.span
data[mispolarise_proportion] += d / tot
for data in (true_data, inferred_data, infdated_data):
data[mispolarise_proportion] /= n_reps
from matplotlib import pyplot as plt
for data, label, linetype in zip(
(true_data, inferred_data, infdated_data),
("Actual", "Plain inferred", "Inf+date"),
("solid", "dotted", "dashed"),
):
for i, (k, d) in enumerate(data.items()):
plt.plot(np.arange(len(d)), d, ls=linetype, c=f"C{i}", label=f"{label} ({k*100}% mispolarised)")
plt.xlabel("Number of root samples")
plt.xscale("log")
plt.xlim(0.9, 100)
plt.ylabel("proportion of genome")
plt.yscale("log")
# For tidyness, change the order of the legend items
plt.legend(*(
[ x[i] for i in [r*len(data)+i for i in range(len(data)) for r in range(0, 3)] ]
for x in plt.gca().get_legend_handles_labels()
), labelspacing=0.25)
jeromekelleher
nspope
petrelharp
In trees from real data, we often find sample nodes attached via long branches to deep in the trees, often to the root. We should investigate why this is happening.
I wondered at first if these deep-rooted samples would be partially removed if we get the sites in the correct order, which seems to produce much better dating results. But some experimentation seems to suggest this doesn't help much (and unsurprisingly, neither does a round of dating + reinference). The plot below shows that regardless of the sample size, very few sample nodes attach to the root (solid lines). But as soon as any inference is involved, even with the true site dates (dot-dashed lines) we get roughly the same distribution of samples attached to their local root as other inference attempts (dotted, dashed lines)