Extracting genotypes and haplotypes for ancestors or subsets of samples

[hyanwong]

Several people (e.g. @saurabhbelsare ) have wanted to be able to extract ancestral genotypes or haplotypes for nodes in a tree sequence. At the moment there's no built-in way to do this in tskit. The simplest way is to flag all the required nodes as samples, and use the methods directly there. For instance:

   def genotype_matrix2(ts, nodes=None, **kwargs):
        if nodes is not None:
            tables = ts.dump_tables()
            flags = tables.nodes.flags
            flags &= np.invert(np.array([tskit.NODE_IS_SAMPLE], dtype=flags.dtype))
            flags[nodes] |= tskit.NODE_IS_SAMPLE
            tables.nodes.flags = flags
            ts = tables.tree_sequence()
    return ts.genotype_matrix(**kwargs)

or

   def haplotypes2(ts, nodes=None, **kwargs):
        if nodes is not None:
            tables = ts.dump_tables()
            flags = tables.nodes.flags
            flags &= np.invert(np.array([tskit.NODE_IS_SAMPLE], dtype=flags.dtype))
            flags[nodes] |= tskit.NODE_IS_SAMPLE
            tables.nodes.flags = flags
            ts = tables.tree_sequence()
    return ts.haplotypes(**kwargs)

[hyanwong]

This should be rendered obsolete if we solve #1207 in the underlying C code.

[saurabhbelsare]

I used the haplotypes2 function, and it is a more concise version of the one I had running, but I'm having the same problem when trying to print the haplotypes without the impute_missing_data=True option. The genotype_matrix2 works, and it puts -1 in place of the missing data, but haplotypes2 gives the error

Generated haplotypes contain missing data, which cannot be represented. Either generate variants (which support missing data) or use the impute missing data option.

and no haplotypes are generated.

[hyanwong]

That's funny, because my version of haplotypes() has an option for missing_data_character, so it should be possible. Does it work with the latest tskit version (e.g. from GitHub: python3 -m pip install git+https://github.com/tskit-dev/tskit#subdirectory=python)

[saurabhbelsare]

That fixed it. I was on tskit 0.2.3. With the update (0.3.5dev1), both missing_data_character and isolated_as_missing are now working. Thanks!

[jodyhey]

I'm trying to understand if the context of these functions fits my problem. I need to generate the sequence of an immediate ancestral node to a sampled node.
The ancestral node changes as we iterate through the treesequence, so building the ancestral sequence of this type involves concatenating haplotypes drawn for each of the ancestral nodes sampled along the treesequence.

I can I use either of these functions to build this kind of sequence? (so far haplotypes2 is not doing working for this, but I don't quite follow how it works).

[andrewkern]

hi @jodyhey! so as you point out, the ancestral node to a sample will change as you move along the chromosome (or through the tree sequence) due to recombination. are you looking for the complete sequence state for one of those ancestral nodes? or something different?

if the former, you could assume the ancestral state at all sites outside of the ancestral material held in the ancestor node

[hyanwong]

Hmm, interesting question. As @andrewkern says, you can fill out missing sections with the ancestral state, or (perhaps better) with tskit.NULL, the missing data character. Alternatively (but I can't quite imagine why you would want to do this) you can mark all the immediate ancestors of a single node as samples, grab their haplotypes, and jigsaw together a new haplotype. On first reading, I thought this was what you wanted to do, but note that no individual genome that we know of actually has this haplotype.

Until #605 is solved, you can't get the variant data from the tree, so the easiest way to do this second thing is probably to go through the variants() iterator at the same time as the sites in each tree. Also the variants() iterator can output information for ancestral nodes, so I think this should work:

import numpy as np
import tskit

focal_sample = 10
tables = ts.tables
parents_of_focal = np.unique(tables.edges.parent[tables.edges.child == focal_sample])

variant_iter = ts.variants(samples=parents_of_focal, isolated_as_missing=False)
# the variant iterator will return the variants for each node defined in parents_of_focal
# so we need to map the index in that array back to the true node number
focal_parent_haplotype = {}  # store the result in a dict, keyed by position - maybe overkill
node_map = {node_id: index for index, node_id in enumerate(parents_of_focal)}
for tree in ts.trees():
    parent_node_id = tree.parent(focal_sample)
    for site in tree.sites():
        v = next(variant_iter)  # need to iterate because we can't yet get variants from a tree, see #605
        assert v.site == site
        focal_parent_haplotype[site.position] = np.array(v.alleles)[v.genotypes[node_map[parent_node_id]]]

print(focal_parent_haplotype)

[jodyhey]

thanks andrewkern (any chance you'll be at PEQG?) and thanks much hyanwong.
Indeed, such a sequence is artificial in that it never will have existed. But it is nonetheless an estimate of the immediate ancestral sequence,
(with respect to a specific sample) of a focal branch. Useful for thinking about some ancestor estimation issues.

[hyanwong]

Useful for thinking about some ancestor estimation issues.

Yes, I can see that it gives rise to interesting trains of thought. Does the code I posted make sense (and does it work for you?!)

[jodyhey]

Thanks very much.   I'm only getting parts of the logic at the moment, so I can't confirm that it does what is intended.  The number of sits is correct, and the base values look quite plausible.  e.g. in a case where 0 is focal and  there are 6 samples and the topology is ((((0,(1,2)),3),4),5)    (the root is on the edge to 5)  0parent  GCTTTTTGGGATGTCGTCGTGAGGCGGGGGCGCCTTTACGGTCATTGTCCTCTACTGTACACAAAGTAGTGTAGTTGACT 0        GCTTTTTGGGATGTCGTCGTGCGGCGGGGGCGCCTTGACGGTCATTGTCCTCTACTGTACAAAAAGTAGTGTAGGCGTCT 1        GCTTTTTGGGAGGTCGTCGTGGGGCGGGGGCGCCCTTACGGTTATTGTCCTCTACTGTGCACAAAGTAGTGGAATTGACT 2        GCTTTTTGGGAGGTCGTCGTGGGGCGGGGGCGCCCTTACGGTTATTGTCCTGTACTGAGCACAGAGTAGTGGAATTGACT 3        GATTTGTAGGACGTCGGCGCGAAGTGGGGGCGCCTTTACGGTCATCGTCCACTACTGTGCACAAAGGGGTGTAGTTGACT 4        GCCTATCGGGCGGTGATGATTAGACGAGGTCGGACGTGGAAGCCTTGTGGTCTGAGGTGCCCCACGTAGTGTAGTTGACT 5        ACCCCTTGTTCTCCCATCGTGAGGCCGACTTTGCTTGGCAACACCAACCCTCCACAATATCCATCGTACACTCGTTTAAA On Tuesday, March 29, 2022, 09:45:10 AM EDT, Yan Wong ***@***.***> wrote:

   Useful for thinking about some ancestor estimation issues.

Yes, I can see that it gives rise to interesting trains of thought. Does the code I posted make sense (and does it work for you?!) — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you were mentioned.Message ID: ***@***.***>