Is it possible to add a method to change ancestral allele after ARG was build.

[silvewheat]

Hi,
After finished tskit+tsinfer+tsdate in a large dataset, if I found several ancestral allele assignment was incorrect. Is there anyway to change the ancestral allele without re-run tskit+tsinfer+tsdate in the full dataset?

[benjeffery]

Hi @silvewheat! Thanks for your question. To do this you would need to use the tables API to both rewrite the ancestral allele at the site, and correct the set of mutations, either by simply adding one mutation at the root, or running map_mutations. If you need more pointers let me know.

[jeromekelleher]

The Tables and editing tutorial is a good place to start for help on this sort of thing @silvewheat.

[hyanwong]

Say you originally thought that the ancestral allele at site s was A, but not realise it was T. In the tree sequence you can set the ancestral allele easily enough by

tables = ts.dump_tables()
tables.sites[s] = tables.sites[s].replace(ancestral_allele="T")

but that will change the genotype of all the samples that directly inherit that ancestral state, so as Ben says, you probably either want to flip immediately back to an 'A':

tree = ts.at(ts.site(s).position)
tables.mutations.add_row(site=s, derived_state="A", node=tree.root)
# you may need to add a time for the mutation e.g. by running compute_mutation_times()
tables.compute_mutation_parents()
tables.sort()
tables.tree_sequence()

Or you want to remap all the mutations at that site to the most parsimonious distribution:

tables = ts.dump_tables()
tables.sites[s] = tables.sites[s].replace(ancestral_state='T')

tree = ts.at(ts.site(s).position)
v = tskit.Variant(ts)
v.decode(s)  # inefficient if you are doing this for lots of sites
ancestral_state, mutations = tree.map_mutations(v.genotypes, v.alleles, ancestral_state="T")
assert ancestral_state == "T"

# quick hack to delete only the mutations at this site
tables.mutations.replace_with(tables.mutations[tables.mutations.site != s])

id_map = {tskit.NULL: tskit.NULL}
for old_id, m in enumerate(mutations):
    id_map[old_id] = tables.mutations.append(m.replace(site=s, parent=id_map[m.parent]))
tables.sort()
tables.tree_sequence()

But of course, you probably have got the inference slightly wrong at that position, so you may well want to run the inference again. You can either edit the sample data file to flip the states (ask on the tsinfer GitHub discussion if you need a hand) and/or mark that site as not to be used for inference, if you aren't sure how trustworthy it is.