Any simplified instructions for a unified genealogy?

[lakishadavid]

I need help going from a tree sequence to a unified genealogy.

My use case is a small number of individuals (n = 38, with plans to expand, but will still be a relatively small data set). I think I'm following enough to create a tree sequence because my results include 9593 trees and 11030 nodes. Output: "Sample file created for 76 samples (38 individuals) with 37281 variable sites."

The tskit tutorials that I have seen don't cover a unified genealogy and the code at https://github.com/awohns/unified_genealogy_paper/tree/master/all-data offers so much that I can't see what I need for my simple case. Is there documentation or a tutorial that presents a simple use case?

From what I understand from the unified genealogy paper, the steps are:

1. infer tree sequence (done)
2. date it using tsdate
3. quality controls based on the count of ancestral haplotypes, edges, variable sites, and mutations.
4. TMRCA between each pair of haplotypes - I did not see the script that would do this as https://github.com/awohns/unified_genealogy_paper/blob/master/src/tmrcas.py gets the mean between populations.
5. Hierarchical clustering on average pairwise TMRCA values - I am clueless as to where this script was located.

I can do close to what I want on a single tree but wondered what was already available for the tree sequence as a unified genealogy:

which_tree = inferred_ts.at_index(3)
ancestor = which_tree.mrca(10, 20)
which_tree.time(ancestor)

# Identifies lineage from the sample to the oldest anecestor of all individuals
for u in which_tree.samples():
    path = []
    v = u
    while v != tskit.NULL:
        path.append(v)
        v = which_tree.parent(v)
    print(u, "->", path)

[awohns]

Hi @lakishadavid, thanks for your message and my apologies for the delayed response!

A few points:

• Have you been able to date your inferred tree sequence with tsdate?
• Are you using tsinfer with a recombination map?
• Quality control of an inferred tree sequence can be subjective and is definitely an area which is ripe for further exploration, but some good things to look at are the number of trees in your inferred tree sequence compared to the number of variant sites (do you have genotyping or sequence data?) as well as the number of edges per sample. The latter analyses can be accomplished with the print_sample_edge_stats() function here. If any samples have many more edges than average, that is a good indication that something might be wrong with them.
• Regarding calculating the TMRCA between each pair of haplotypes, what you have written is close to what you would want for this. The code I've added below is a simple way of finding the average TMRCA between each pair of haplotypes.
• Furthermore, the code below also contains simple plotting code to create a hierarchically clustered plot similar to figure 2 in our "Unified genealogy of modern and ancient genomes" paper. This is where we performed that plotting function in the github repo for the paper.

By the way, I plan to write a tutorial which will help guide users seeking to infer a tree sequence from modern and ancient genomes, please keep an eye out at https://tsdate.readthedocs.io/en/latest/ for this!

Code to get the pairwise (log) average TMRCA for each pair of haplotypes and plot the result with hierarchical clustering on the log-average times:

import itertools
import numpy as np
import seaborn as sns

samples = inferred_ts.samples()
haplotype_pairs = {(pair[0], pair[1]): [] for pair in itertools.combinations(samples, 2)}
    
for tree in inferred_ts.trees():
    for pair in itertools.combinations(samples, 2):
        haplotype_pairs[(pair[0], pair[1])].append(tree.tmrca(pair[0], pair[1]))
average_tmrcas = np.zeros((inferred_ts.num_samples, inferred_ts.num_samples))
for pair, tmrcas in haplotype_pairs.items():
    average_tmrcas[pair[0], pair[1]] = np.exp(np.mean(np.log(tmrcas)))

sns.clustermap(average_tmrcas)