Ways to visualize diploid tree sequences that aren't sequences of trees

[petrelharp]

In particular, we need pictures of the "first-generation" individuals that trip people up.

Edit: see below; this issue has evolved from the initial intention.

[bhaller]

Excerpted from an email I sent to Peter just now:

Well, yes, the notes do already state that the first generation is retained but is not part of the sample, and that that is necessary for recapitation. But I guess I had not really thought about / realized that first-generation individuals that are not ancestral to the sample would (a) also be present (of course), and (b) would still be there even after recapitation, but without any simulated histories (again, of course, but I never thought about it), and so (c) would need to be excluded or otherwise dealt with during subsequent analysis to avoid getting confused. That’s not an issue that I remember ever coming up, but I can see it being surprising/confusing for people. It’s one of those things where, if the tree sequences were visualized in some way, it would be immediately obvious, but when just working with them as a black-box object in Python it is easy not to realize or to think through the implications.

I know tskit currently provides some visualization tools that show a given tree graphically. I think I’m thinking of something a little different. I’d maybe like to see a visualization tool that is more just about the times and relationships of individuals, nodes, and samples, without worrying about the particular tree structures that connect them through time. This representation would let me see things like:

• there are a bunch of first-generation individuals that are retained but are not part of the sample

• only some of those first-gen individuals are ancestral to the sample

• the sample is composed of diploid individuals, but every second genome has no ancestral relationship because the individuals are actually haploid

• after recapitation, some of the first-gen individuals now coalesce back to a single ancestor, but some of them – the ones that are not ancestral to the sample – do not, but are still present

I’m imagining a depiction that is a sort of summary of the entire tree sequence, rather than a picture of a specific tree. The y-axis would be time, with the more distant past at the top I suppose (that’s what would be intuitive to me, anyway). The x-axis is meaningless and used only for visual separation of elements; ideally x positions would be chosen to minimize the number of line crossings in the drawn graph, but that’s a frill. Each node would be shown as a small dot, and each individual would be a (gray) circle that encloses two dots (for diploids); nodes that do not have a recorded individual associated with them would be shown as standalone dots not enclosed by a “individual circle”. Lines would connect a later-time dot to an earlier-time dot if any portion of the genome was inherited from that earlier-time dot (without passing through an intervening dot), so it would be more of a directed graph than a tree, and would represent an overlay of all of the tree sequences onto a single diagram. Samples would be shown as red dots, while non-samples ancestral to the sample would be black dots, and nodes not ancestral to the sample would be white dots. I’ve attached a sort of mock-up of how I imagine a model like Jean’s (a SLiM haploid clonal model – i.e., diploids where every second genome is empty and has no ancestral connections – with occasional horizontal gene transfer) would look before and after recapitation; I don’t know that it makes any actual sense as a diagram of a real tree sequence, but it illustrates what I mean, I hope. Possibly the lines of the graph could be colored to show something useful, too – perhaps the weighting of that line across all of the tree sequences (how common of a path of inheritance it is for the sample), or perhaps how many mutations occurred along all of the tree-sequences branches represented by it, or something of that sort (perhaps things like this could be options chosen by the user).

With a visualization tool like this, I could write a python script that would read in a .trees file from SLiM, recapitate, change the sample in some way, simplify, overlay mutations, etc., and I could call the visualization function at each step of the process to confirm that things look the way I intend them to look. I think this would help me immensely, since I am a very visual thinker, and I bet many users would love it.

(Addendum from a later email:) Expanding upon this idea a bit more: first-generation individuals/nodes should also use a distinctive color, and remembered individuals should be designated in some orthogonal fashion so they are also visible (bigger dots?). :->

I have attached illustrations here:

[bhaller]

@jeromekelleher and @molpopgen, comments would be welcome. I'm quite excited about this idea. @petrelharp wrote back with a caveat:

As for writing a tool that does this automatically - you'd like to draw the pedigree, sort of. Whenever I've thought about this before I've decided that it would be way too messy for anything except for extremely small examples, so it wasn't worth writing the (probably quite tricky) code. But, it might be worth a try - visualizing the tree sequence somehow would be quite a useful thing.

I agree that for typical "real" models with thousands or millions of individuals these diagrams would be very large and might not be terribly useful (although I'm not sure they would never be useful; I'd like to see what the diagram for, say, a 10-population stepping-stone model with a soft selective sweep would look like :->). But I think people could still find them very useful for developing their Python tree-seq workflow, by simulating with a downscaled model so the diagrams aren't too big, writing out their Python steps and visualizing them to make sure everything in that pipeline is working as intended, and then upscaling their model to production size (and commenting out the visualizations).

I also think it would be very useful for generating pictures to use in the documentation, for explaining concepts like the retention of first-generation ancestors, simplification, and recapitation.

[jeromekelleher]

Well, there is ts.draw_text() and ts.draw_svg():

import msprime

ts = msprime.simulate(5, recombination_rate=0.3, random_seed=2)
print(ts.draw_text())
ts.draw_svg("tmp.svg")

This is the text (might get mangled by some browsers):

5.77┊    10     ┊           ┊           ┊
    ┊  ┏━━┻━┓   ┊           ┊           ┊
4.41┊  ┃    ┃   ┊           ┊   9       ┊
    ┊  ┃    ┃   ┊           ┊ ┏━┻━━┓    ┊
2.18┊  ┃    ┃   ┊     8     ┊ ┃    8    ┊
    ┊  ┃    ┃   ┊  ┏━━┻━┓   ┊ ┃  ┏━┻━┓  ┊
0.89┊  ┃    7   ┊  ┃    7   ┊ ┃  ┃   7  ┊
    ┊  ┃  ┏━┻┓  ┊  ┃  ┏━┻┓  ┊ ┃  ┃  ┏┻┓ ┊
0.81┊  ┃  ┃  6  ┊  ┃  ┃  6  ┊ ┃  ┃  ┃ ┃ ┊
    ┊  ┃  ┃ ┏┻┓ ┊  ┃  ┃ ┏┻┓ ┊ ┃  ┃  ┃ ┃ ┊
0.04┊  5  ┃ ┃ ┃ ┊  5  ┃ ┃ ┃ ┊ ┃  5  ┃ ┃ ┊
    ┊ ┏┻┓ ┃ ┃ ┃ ┊ ┏┻┓ ┃ ┃ ┃ ┊ ┃ ┏┻┓ ┃ ┃ ┊
0.00┊ 0 4 1 2 3 ┊ 0 4 1 2 3 ┊ 3 0 4 1 2 ┊
  0.00        0.10        0.74        1.00

and SVG:

I guess this isn't what you're looking for though?

[bhaller]

Right; what I'm looking for is a diagram that combines/summarizes all of those tree diagrams into a single diagram. As Peter wrote, more like a pedigree, except that (a) only the nodes present in the tree sequence are shown, so it is not a full pedigree in the genealogical sense, and (b) lines are shown between every pair of nodes that are connected by a branch in any tree in the tree sequence, not just between parents/offspring. And the other big thing I'm looking for is color-coding and other visualization of things like which nodes are samples, which nodes are first-generation ancestors, which nodes are remembered individuals, which pairs of nodes are grouped together to form diploid individuals in the individuals table, etc., which is not visible in those existing visualizations. So this is a much higher-level view of what is going on than is provided by ts.draw_text() and ts_draw_svg(). See the diagrams I posted above, which hopefully illustrate what I'm talking about.

[molpopgen]

I don't have a concrete technical suggestion, but this seems like something that a Python package capable of drawing a network might handle?

[bhaller]

I don't have a concrete technical suggestion, but this seems like something that a Python package capable of drawing a network might handle?

Yes; if there were a package that handled choosing x-positions so as to minimize line-crossing that would be particularly useful. The rest of it is probably pretty easy to do oneself, really.

[molpopgen]

Yes; if there were a package that handled choosing x-positions so as to minimize line-crossing that would be particularly useful. The rest of it is probably pretty easy to do oneself, really.

Python (and R) have a seemingly limitless supply of options. I've been (casually) interested in this same issue, so I'll be looking for something soon.

[molpopgen]

This gets pretty close, but isn't perfect depending on how picky one is, but it took some effort to figure out.

It needs a lot of the graphviz stack on your system, including the C runtime library, as pygraphviz is a dependency in here somewhere.

import msprime
import pandas as pd
import networkx
from collections import namedtuple

Datumz = namedtuple('Datumz', ['parent', 'child'])
ts = msprime.simulate(5, recombination_rate=6, random_seed=666)

D = []
for t in ts.trees():
    for n in t.nodes():
        for c in t.get_children(n):
            D.append(Datumz(n, c))


df = pd.DataFrame(D, columns=Datumz._fields)
df['time'] = ts.tables.nodes.time[df.parent]

G = networkx.from_pandas_edgelist(
    df, 'parent', 'child', create_using=networkx.DiGraph())

A = networkx.nx_agraph.to_agraph(G)

times = df.groupby(['time'])
all_samples = [i for i in ts.samples()]
for n, g in times:
    up = [i for i in g.parent.unique() if i not in all_samples]
    A.add_subgraph(up, rank='same')

A.add_subgraph(all_samples, rank='same')
A.draw('A.png', prog='dot')

[molpopgen]

To extend this to remembered individuals, etc., you need to manipulate the horizontal groupings, too, add color, etc.. dot can do all of that.

[molpopgen]

This will change the border color and shape of nodes 0 and 1. The UI semantics are the same as the low-level program dot. In one of the loops above, I am excluding sample nodes from the sub-graph. For this simulation, that has no effect, as there is only 1 parent possible in a tree at exactly one time. For remembered nodes, you'd need to process both parents and children, ask if they are samples, and then do something fun with their color/shape/fill.

for i in [0, 1]:
    n = A.get_node(i)
    n.attr['color'] = 'red'
    n.attr['shape'] = 'box'

[molpopgen]

Seems likely that everything can be done entirely in pygraphviz: https://pygraphviz.github.io/documentation/latest/tutorial.html

I used networkx to start simply b/c I'd heard of it recently in an msprime PR.

[molpopgen]

This version respects the temporal order of the root nodes on each tree.

import msprime
import pandas as pd
import networkx
from collections import namedtuple

Datumz = namedtuple('Datumz', ['parent', 'child'])
ts = msprime.simulate(5, recombination_rate=6, random_seed=666)

D = []
for t in ts.trees():
    for n in t.nodes():
        for c in t.get_children(n):
            D.append(Datumz(n, c))


df = pd.DataFrame(D, columns=Datumz._fields)
df['time'] = ts.tables.nodes.time[df.parent]
df.sort_values(['time'], inplace=True, ascending=False)
G = networkx.from_pandas_edgelist(
    df, 'parent', 'child', create_using=networkx.DiGraph())

A = networkx.nx_agraph.to_agraph(G)


# Make a subgraph of roots, where our A->B edges
# are A is older than B.
roots = []
for t in ts.trees():
    for r in t.roots:
        if r not in roots:
            roots.append(r)

roots = sorted(roots, key=lambda x: -ts.tables.nodes.time[x])


pg = A.add_subgraph([], level="same", name="roots")
for i in range(1, len(roots)):
    pg.add_edge(roots[i-1], roots[i], style="invis")

A.add_subgraph([i for i in ts.samples()], rank='same', name="samples")

A.draw('A.png', prog='dot')

[petrelharp]

wow, love it!

[molpopgen]

With the root order preserved, at this point it should be relatively straightforward to create subgraphs representing actual samples taken at
different time points. The hard thing, maybe, is to have nested subgraphs of pairs of nodes representing diploids.

[petrelharp]

Ping everyone here - please have a look at my pull request with figures. Let me know if you can't build the docs locally (if you have sphinx set up, you just need to run make in the docs/ directory, and browse to docs/_build to look at the result. Are the figures useful? Are they in the right place? Do we need more/less? I opted to do them by hand (inkscape) rather than figure out how to do these things in python (especially the labeling!).