The local ancestral values for the AmericanAdmixture_4B11 model

[jikhashkya]

Hi there,

I simulated the AmericanAdmixture_4B11 model and wanted to look at the ancestral values (i.e. local ancestral values) for the contemporary admixed population but the code I used to generate the local ancestral values showed unexpected results. I wonder if the results are messed up because of the coalescent nature of the simulation. If so, is it possible to obtain the local ancestral values for these predefined demographic models?

Simulation code:

stdpopsim HomSap \
    -c chr20 \
    -o admixed.ts \
    -g HapMapII_GRCh38 \
    -d AmericanAdmixture_4B11 \
    -s 108 AFR:100 EUR:100 ASIA:100 ADMIX:5000

Ancestral value retrieval code:

import msprime, tskit

input_ts = "./data/admixed.ts"
ts = tskit.load(input_ts)
ts = ts.simplify()
print(f"Total # of populations = {ts.num_populations}")
target_pop = ts.num_populations - 1
haplotype_ids = ts.samples(population=target_pop)
print(f"Total # of haplotypes in pop{target_pop} = {len(haplotype_ids)}")

fw = open('GT_ancestral_codes', 'w')
for n_idx, node in enumerate(haplotype_ids):
    print(f"{n_idx}. Processing node: {node}", flush=True)
    fw.write(f"ADMIX_{n_idx//2}_{n_idx % 2}")
    for idx, t in enumerate(ts.trees()):
        path = []
        v = node

        while v != tskit.NULL:
            path.append(v)
            v = t.parent(v)

        if idx == 0:
            start = int(t.interval[0])
            ancestral_code = t.population(path[-1])

        if t.population(path[-1]) != ancestral_code:
            fw.write(f",{start},{int(t.interval[0])},{ancestral_code}")
            start = int(t.interval[0])
            start_tree_idx = idx
            ancestral_code = t.population(path[-1])
    fw.write(f",{start},{int(ts.last().interval[1])},{ancestral_code}")
    fw.write('\n')
    #assert t.num_roots == 1  #, ("not coalesced! on segment {} to {}".format(
    #t.interval[0], t.interval[1]))
fw.close()

Any help would be greatly appreciated !!!

[petrelharp]

If I understand what you mean by "local ancestral values", tspop may be just what you're looking for?

[jikhashkya]

Thank you for the prompt response @petrelharp. Tspop does seem to be what I'm looking for, but more importantly, the add_census() method from msprime seems to be the answer I was looking for (pointed out in the Tspop documentation as well). For other users that might be looking for an answer, this method adds ancestral nodes at the specified time point so that for every branch of the tree there's an ancestral node. After that, the population that this ancestral node belonged to can be extracted.