Generating a genome-wide tree topology count

[hyanwong]

In tskit-dev/tsinfer#989, @MatteoSebastianelli asked:

I have been generating tree sequences with tsinfer and would like to get a tree topology count across the entire genome. Would that be possible? I was thinking I should loop through my .trees files with the union function to merge them (it seems it cannot merge more than 2 files) and then use that output to count tree topologies. Does this sounds like a reasonable approach? I may be missing it, perhaps there is a functions that does this already?

[hyanwong]

If you mean that you want to summarize across multiple chromosomes, each inferred separately by tsinfer, then I think it may not be necessary to union the tree sequences: you can probably summarise the topology results across chromosomes by just adding the means (weighted by the chromosome length) and dividing by the total genome length. Note that you might want to remove centromeres and other "blank" regions from the length calculations.

Within a chromosome, however, there are fast methods for counting topologies: see #671 and the link from there to

tskit/python/tskit/combinatorics.py

Line 37 in 829fd1f

def treeseq_count_topologies(ts, sample_sets):

. We should probably make this available as an API, if people are going to use it extensively.