--- jupytext: text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.10.3 kernelspec: display_name: Python 3 language: python name: python3 --- ```{currentmodule} tskit ``` (sec_tskit_viz)= ```{code-cell} ipython3 :tags: [remove-cell] import io import string import msprime import tskit def viz_ts(): dem = msprime.Demography.from_species_tree("((A:900,B:900)ab:100,C:1000)abc;", initial_size=1e3) samples = {"A": 20, "B": 20, "C":20} # take 20 diploids from terminal populations A, B, C ts_full = msprime.sim_ancestry( samples, demography=dem, sequence_length=5e4, recombination_rate=1e-8, random_seed=1234 ) ts_full.dump("data/viz_ts_full.trees") first_4_nodes = (0, 1, 2, 3) # ids of the first 4 sample nodes (here, 2 individuals from A) ts_tiny = ts_full.simplify(first_4_nodes) # a tiny 4-tip TS tip_orders = { tuple(u for u in t.nodes(order="minlex_postorder") if t.is_sample(u)) for t in ts_tiny.trees() } topologies = {tree.rank() for tree in ts_tiny.trees()} # Check we have picked a random seed that gives a nice plot of 7 trees assert tip_orders == {first_4_nodes} and len(topologies) > 1 and ts_tiny.num_trees == 8 ts_tiny.dump("data/viz_ts_tiny.trees") eight_nodes = first_4_nodes + (40, 41, 80, 81) # Add nodes from individuals in B & C ts_small = ts_full.simplify(eight_nodes) # a small 8-tip TS ts_small.dump("data/viz_ts_small.trees") ts_small_mutated = msprime.sim_mutations(ts_small, rate=1e-7, random_seed=342) # 3rd tree should have first site with 2 muts first_site_tree_2 = next(ts_small_mutated.at_index(2).sites()) assert len(first_site_tree_2.mutations) == 2 # mutation 8 should be above node 16 in the 1st tree assert ts_small_mutated.site(8).mutations[0].id == 8 assert ts_small_mutated.site(8).mutations[0].node == 16 ts_small_mutated.dump("data/viz_ts_small_mutated.trees") def viz_selection(): sequence_length = 5e4 sweep_model = msprime.SweepGenicSelection(position=sequence_length/2, s=0.01, start_frequency=0.5e-4, end_frequency=0.99, dt=1e-6) ts_selection = msprime.sim_ancestry(12, ploidy=1, model=[sweep_model, msprime.StandardCoalescent()], population_size=1e4, recombination_rate=2e-8, sequence_length=sequence_length, random_seed=222, ) tables = ts_selection.dump_tables() tables.nodes.clear() tables.nodes.metadata_schema = tskit.MetadataSchema.permissive_json() for node in ts_selection.nodes(): metadata = {} if node.is_sample(): metadata["name"] = f"Sample {string.ascii_uppercase[node.id]}" tables.nodes.add_row(node.replace(metadata=metadata)) ts_selection = tables.tree_sequence() ts_selection.dump("data/viz_ts_selection.trees") def create_notebook_data(): viz_ts() viz_selection() # create_notebook_data() # uncomment to recreate the tree seqs used in this notebook ``` # Visualization **Yan Wong** It is often helpful to visualize a single tree --- or multiple trees along a tree sequence --- together with sites and mutations. {ref}`Tskit ` provides methods to do this, outputting either plain ascii or unicode text, or the more flexible [Scalable Vector Graphics](https://www.w3.org/TR/SVG11/) (SVG) format. The first two sections of this tutorial give details of these methods, using a few 50kb tree sequences generated by {ref}`msprime ` as examples: one in which selection has occurred, and others which involve subdivision into 3 populations (labelled _A_, _B_, and _C_). ```{code-cell} ipython3 :"tags": ["hide-input"] import tskit # To see how these tree sequences were made, inspect the removed cell at the top of this # notebook, visible in e.g. the `.md` (markdown) version downloadable from the menubar ts_tiny = tskit.load("data/viz_ts_tiny.trees") ts_small = tskit.load("data/viz_ts_small.trees") ts_full = tskit.load("data/viz_ts_full.trees") ``` If you just want a quick look at visualization possibilities, you might want to skip the explanations and just browse some {ref}`sec_tskit_viz_examples`, which contain fully reproducible code. :::{note} This tutorial is primarily focussed on showing a tree sequence as a set of marginal trees along a genome. The section titled {ref}`sec_tskit_viz_other` provides examples of other representations of tree sequences, or the processes that can create them. Additionally, other software tools exist that can plot genealogies (ARGs) in tree sequence format. For example, in these tutorials we use the [tskit_arg_visualizer](https://github.com/kitchensjn/tskit_arg_visualizer) tool to provide simple graph-centric visualizations. You might also want to explore tools such as [ARGscape](https://www.argscape.com), [Lorax](https://lorax.ucsc.edu), [tsbrowse](https://tskit.dev/tsbrowse/docs/stable/intro.html), or even high-level analysis tools like [TwisstNTern](https://github.com/HilaLifchitz/twisstntern_v2). ::: ## Text format The {meth}`TreeSequence.draw_text` and {meth}`Tree.draw_text` methods provide a quick way to print out a tree sequence, or an individual tree within it. They are primarily useful for looking at topologies in small datasets (e.g. fewer than 20 sampled genomes), and do not display mutations. ```{code-cell} ipython3 # Print a tree sequence print(ts_tiny.draw_text()) print("The first tree in the tree sequence above, but replacing some node ids with names:") print(ts_tiny.first().draw_text( # An example of how to change or omit node labels: unspecified nodes are omitted # The same convention applies to SVG graphics node_labels={0: "Alice", 1: "Bob", 2:"Chris", 3: "Dora", 6: "MRCA"} )) ``` ## SVG format Most users will want to use the SVG drawing functions {meth}`TreeSequence.draw_svg` and {meth}`Tree.draw_svg` for visualization. Being a vectorised format, SVG files are suitable for presentations, publication, and {ref}`editing or converting ` to other graphic formats; some basic forms of animation are also possible. Both functions produce an SVG string which is automatically drawn if the string is the result of the last call in a Jupyter notebook: ```{code-cell} ipython3 from IPython.display import display svg_size = (800, 250) # Height and width for the SVG: optional but useful for this notebook svg_string = ts_tiny.draw_svg( size=svg_size, y_axis=True, y_label=" ", # optional: show a time scale on the left time_scale="rank", x_scale="treewise", # Force same axis settings as the text view title="A basic SVG plot", ) display(svg_string) # If the last line in a cell, wrapping this in display() is not needed ``` By default, sample nodes are drawn as black squares, and non-sample nodes are drawn as black circles (but see below for ways to e.g. hide or colour these node symbols - NB. apologies to US readers: the British spelling of "colour" will be used in the rest of this tutorial). ### Axes and scales For ease of drawing, the text representation and the SVG image above use unconventional non-linear X and Y coordinate systems. By default, the SVG output uses a more conventional linear scale for both time (Y) and genome position (X), and indicates the position of each tree along the genome by an alternating shaded background. Although more intuitive, linear scales can obscure some features of the trees, for example causing labels to overlap: ```{code-cell} ipython3 ts_tiny.draw_svg(size=svg_size, y_axis=True) ``` One way to avoid overlapping labels on the Y axis is to use the `y_ticks` parameter, which will be used in most subsequent examples in this tutorial. (sec_tskit_viz_larger_tree_sequences)= ### Larger tree sequences So far, we have plotted only very small tree sequences. To visualize larger tree sequences it is sometimes advisable to focus on a small region of the genome, possibly even a single tree. The `x_lim` parameter allows you to plot the part of a tree sequence that spans a particular genomic region: here's a slightly larger tree sequence with 8 samples, but where we've restricted the amount of the tree sequence we plot: ```{code-cell} ipython3 x_limits = [5000, 15000] # Create evenly-spaced y tick positions to avoid overlap y_tick_pos = [0, 1000, 2000, 3000] display(ts_small.draw_svg( size=svg_size, y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits, title="The tree sequence between positions {} and {}".format(*x_limits) )) third_tree = ts_small.at_index(2) display(third_tree.draw_svg(title="A plot of the third tree")) ``` As the number of sample nodes increases, internal nodes often bunch up at recent time points, obscuring relationships. Setting `time_scale="rank"`, as in the first SVG plot, is one way to solve this. Another is to use a log-scale on the time axis, which can be done by specifying `time_scale="log_time"`, as below. To compare node times across the plot, this example also uses the `y_gridlines` option, which puts a very faint grid line at each y tick (if you are finding the lines difficult to see, note that the line intensity, along with many other plot features, can be modified through {ref}`styling `, which we also use in this example to avoid overlapping text by shrinking the node labels and rotating those associated with leaves; styling is detailed {ref}`later in this tutorial `). (sec_tskit_viz_large_tree_sequence)= ```{code-cell} ipython3 wide_fmt = (1200, 280) # Create a stylesheet that shrinks labels and rotates leaf labels, to avoid overlap node_label_style = ( ".node > .lab {font-size: 80%}" ".leaf > .lab {text-anchor: start; transform: rotate(90deg) translate(6px)}" ) ts_full.first().draw_svg( size=wide_fmt, time_scale="log_time", y_gridlines=True, y_axis=True, y_ticks=[1, 10, 100, 1000], style=node_label_style + "svg > .title text {font-size: 150%}", title="A larger tree, on a log timescale", ) ``` For even larger numbers of samples, you can plot the trees for a subset of samples by applying the {meth}`TreeSequence.simplify` method (make sure to specify `filter_nodes=False` to retain the same node IDs). For a fancier solution, see {ref}`sec_tskit_viz_SVG_examples_larger_plots` below. ### Plotting mutations By default the SVG visualization also plots sites and mutations on the tree or tree sequence (this can be disabled using the `omit_sites` parameter). For example, adding mutations to the 8-sample tree sequence above gives the following plot: each mutation is marked by a red cross on the branch where it occurs. Symbols are either placed at the mutation's known time, or (if the mutation time is unknown or the `time_scale` parameter has been set to `"rank"`) spaced evenly along the branch. By default, each mutation is also labelled with its {class}`mutation ID`. If the X axis is shown (which it is by default when drawing a tree sequence, but not when drawing an individual tree) then the sites are plotted using a black tickmark above the axis line. Each plotted mutation at the site is then overlaid on top of this as a red downwards-pointing chevron. ```{code-cell} ipython3 ts_mutated = msprime.sim_mutations(ts_small, rate=1e-7, random_seed=342) ts_mutated.draw_svg(y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits) ``` Note that, unusually, the rightmost site on the axis has more than one stacked chevron, indicating that multiple mutations in the tree occur at the same site. These could be mutations to different allelic states, or recurrent/back mutations. In this case the mutations, 14 and 15 (above nodes 1 and 6) are recurrent mutations from T to G. ```{code-cell} ipython3 :"tags": ["hide-input"] ts_mutated = tskit.load("data/viz_ts_small_mutated.trees") site_descr = str(next(ts_mutated.at_index(2).sites())) print(site_descr.replace("[", "[\n ").replace("),", "),\n ").replace("],", "],\n")) ``` #### Which mutations are shown? When using the `x_lim` parameter, only the mutations in the plotted region are shown. For the third tree in the tree sequence visualization above, we thus haven't plotted mutations above position 15000. We can see *all* the mutations in the tree by changing the plot region, or simply plotting the tree itself: ```{code-cell} ipython3 tree3 = ts_mutated.at_index(2) tree3.draw_svg( size=(300, 300), title=f"Third tree, from {int(tree3.interval.left)} bp to {int(tree3.interval.right)} bp", ) ``` (sec_tskit_viz_extra_mutations)= However, when plotting a single tree it may not be evident that identical branches may exist in several adjacent trees, indicating an {ref}`edge ` that persists across adjacent trees. For instance the rightmost branch in the tree above, from node 10 down to 7, exists in the previous two trees too. Indeed, this edge has a mutation on it at position 6295, in the first tree. This mutation is not plotted in the tree above, but if you want *all* the mutations on each edge to be plotted, you can set the `all_edge_mutations` parameter to `True`. This adds any extra mutations that are associated with an edge in the tree but which fall outside the interval of that tree; by default these mutations are drawn in a slightly different shade (e.g. mutation 64 below). ```{code-cell} ipython3 tree3.draw_svg( size=(300, 300), all_edge_mutations=True, title=( 'Third tree as above, but with' 'visible edges showing all mutations' ) ) ``` (sec_tskit_viz_labelling)= ### Labelling and annotation Although the default node and mutation labels show unique identifiers, they aren't terribly intuituive. The `node_labels` and `mutation_labels` parameters can be used to set more meaningful labels (for example from the tree sequence {ref}`sec_metadata`). See {ref}`sec_tskit_viz_dynamic_effects` if you want to dynamically hide and show such labels. ```{code-cell} ipython3 nd_labels = { # An array of labels for the nodes # Set sample node labels from metadata. Here we use the population name, but you might want # to use the *individual* name instead, if the individuals in your tree sequence have names n.id: ts_mutated.population(n.population).metadata["name"] for n in ts_mutated.nodes() if n.is_sample() } mut_labels = { # An array of labels for the mutations mut.id: f"{mut.inherited_state}→{mut.derived_state}" for mut in ts_mutated.mutations() } ts_mutated.draw_svg( size=(1000, 300), y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits, node_labels=nd_labels, mutation_labels=mut_labels, ) ``` #### Annotating genome regions To annotate genomic regions along the X axis of the tree sequence plot, you can pass a dictionary of `x_regions`, mapping `(start, end)` tuples to labels. Below we have also used css to hide all but the first and last x axis tick label to avoid visual clashing, and hidden all the mutation and node labels by setting them to `{}`: ```{code-cell} ipython3 x_regions = {(5_000, 12_000): "Gene 1", (21_123, 33_321): "Gene 2"} hide_internal_x_tick_labels = ( ".x-axis .tick .lab {visibility: hidden}" ".x-axis .tick:first-of-type .lab, .x-axis .tick:last-of-type .lab {visibility: visible}" ".x-regions .r1 rect {fill: cyan}" # 2nd region (ID 1) tagged with class="r1" ) ts_mutated.draw_svg( size=(1000, 300), y_axis=True, y_ticks=range(0, int(ts_mutated.max_time), 1000), node_labels={}, mutation_labels={}, style=hide_internal_x_tick_labels, x_regions=x_regions, ) ``` #### Arbitrary annotation It is also possible to add arbitrary annotations to the SVG plot, as detailed in the next section. To locate positions for annotation, see the examples in {ref}`sec_tskit_viz_svg_plot_internals`. (sec_tskit_viz_adding_bespoke)= ### Adding bespoke SVG The `preamble` option allows arbitrary SVG text to be added at the start of the plot. This can be useful to annotate plots, produce legends, etc. although it requires some knowledge of the SVG graphics language (see the {ref}`sec_tskit_viz_legend_example` later in this tutorial for a node colour legend). ```{code-cell} ipython3 mut_labels = { mut.id: f"{mut.inherited_state}{int(site.position)}{mut.derived_state}" for site in ts_mutated.sites() for mut in site.mutations } svg_text=( '' '' 'Mutation labels:' '<inherited_state><POSITION><derived_state>' '' ) ts_mutated.draw_svg( size=(800, 300), title="A plot with a rudimentary legend to describe the mutation label format", y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits, node_labels=nd_labels, mutation_labels=mut_labels, preamble=svg_text, ) ``` #### Plotting side-by-side As any conceivable SVG commands can be added (including nesting one SVG inside another), using `preamble` is extremely flexible, but can be fiddly. You can, for example, plot one tree next to a completely different one by adding one into the preamble of the other. This is demonstrated in the helper function below. It adds space into the first plot without rescaling the image itself by using the `canvas_size` parameter. Further plots are then added into the preamble, moved to the left using `root_svg_attributes={"x", x_offset}`. ```{code-cell} ipython3 def draw_svg_side_by_side( drawables, *, size=(200, 200), sizes=None, padding=40, canvas_size=None, per_svg_kwargs=None, **kwargs, ): """ Plot multiple Tree or TreeSequence objects side-by-side by embedding the SVG for each later object in the preamble of the first. :param list drawables: A list of Tree or TreeSequence objects. :param tuple size: The size of each individual SVG plot, if `sizes` is not given. :param list sizes: An optional list of sizes, one per plot. :param int padding: The horizontal gap between adjacent plots. :param tuple canvas_size: The overall canvas size. If None, infer it from `size`. :param list per_svg_kwargs: An optional list of dicts of keyword arguments to pass to each individual `draw_svg` call. These are merged on top of `**kwargs`. :param kwargs: Common keyword arguments passed to each `draw_svg` call. """ if len(drawables) == 0: raise ValueError("Need at least one drawable") if per_svg_kwargs is None: per_svg_kwargs = [{} for _ in drawables] if len(per_svg_kwargs) != len(drawables): raise ValueError("per_svg_kwargs must have the same length as drawables") if sizes is None: sizes = [size for _ in drawables] if len(sizes) != len(drawables): raise ValueError("sizes must have the same length as drawables") if canvas_size is None: canvas_size = ( sum(s[0] for s in sizes) + (len(drawables) - 1) * padding, max(s[1] for s in sizes), ) preamble = [] x_offset = sizes[0][0] + padding for j, drawable in enumerate(drawables[1:], start=1): svg_kwargs = dict(kwargs) svg_kwargs.update(per_svg_kwargs[j]) svg_kwargs["size"] = sizes[j] svg_kwargs["root_svg_attributes"] = { **svg_kwargs.get("root_svg_attributes", {}), "x": x_offset, } preamble.append(drawable.draw_svg(**svg_kwargs)) x_offset += sizes[j][0] + padding first_kwargs = dict(kwargs) first_kwargs.update(per_svg_kwargs[0]) first_kwargs["size"] = sizes[0] first_kwargs["canvas_size"] = canvas_size first_kwargs["preamble"] = first_kwargs.get("preamble", "") + "".join(preamble) return drawables[0].draw_svg(**first_kwargs) ``` ```{code-cell} ipython3 draw_svg_side_by_side( [tree3, tree3], # these could be different trees from different tree sequences size=(300, 300), per_svg_kwargs=[{"title": "3rd tree"}, {"title": "3rd tree, no sites", "omit_sites": True}], ) ``` (sec_tskit_viz_reordering nodes)= ### Visually reordering nodes As _tskit_ is not primarily a visualization library, there are no methods for rotating branches when visualizing specific trees. Moreover, a tidy arrangement for one tree could be a poor arrangement for an adjacent tree in the tree sequence. However, because the default ordering is to put lower numbered leaves to the left, it is possible to obtain an arbitrary ordering by changing the node IDs as required (e.g. using the `subset` method). If you are labelling nodes using metadata, this will all be fine, but if you are using the default ID labelling scheme, you'll should provide labels that map back to the original IDs, to avoid confusion. Here's an example, using a generic reordering function to reverse the visual order of nodes 0...4 () ```{code-cell} ipython3 import numpy as np def reorder_leaves(tree, leaf_order_ids): # Return a new tree in an identical tree sequence but with the node IDs # swapped to draw leaf order as close possible to the provided leaf_order_ids # You will need to plot the node labels using metadata or the returned node_map ts = tree.tree_sequence all_nodes = np.arange(ts.num_nodes) leaves = [u for u in tree.nodes(order="minlex_postorder") if tree.is_leaf(u)] all_nodes[leaves] = leaf_order_ids reorder_ts = ts.subset(all_nodes, reorder_populations=False, remove_unreferenced=False) return reorder_ts.at_index(tree.index), all_nodes # swap the order of 0, 1, 2, 3, 4 orig_tree = ts_mutated.first() new_tree, node_map = reorder_leaves(orig_tree, [4, 3, 2, 1, 0, 5, 7, 6]) draw_svg_side_by_side( [orig_tree, new_tree], per_svg_kwargs=[ {"title": "Original tree"}, {"title": "Leaf nodes reordered", "node_labels": {u: v for u, v in enumerate(node_map)}}, ], ) ``` This approach can be used to define functions to plot ladderized trees, see e.g. [this GitHub discussion](https://github.com/tskit-dev/tskit/discussions/3160). (sec_tskit_viz_styling)= ### Styling The SVG output produced by tskit contains a large number of [classes](https://www.w3.org/TR/SVG2/styling.html#ClassAttribute) which can be used to target different elements of the drawing, allowing them to be hidden, styled, or otherwise manipulated. This is done by passing a [cascading style sheet (CSS)](https://www.w3.org/TR/SVG/styling.html) string to `draw_svg`. A common use of styles is to colour nodes by their population: ```{code-cell} ipython3 styles = [] # Create a style for each population, programmatically (or just type the string by hand) for colour, p in zip(['red', 'green', 'blue'], ts_full.populations()): # target the symbols only (class "sym") s = f".node.p{p.id} > .sym " + "{" + f"fill: {colour}" + "}" styles.append(s) print(f'"{s}" applies to nodes from population {p.metadata["name"]} (id {p.id})') css_string = " ".join(styles) print(f'CSS string applied:\n "{css_string}"') ts_full.first().draw_svg( size=wide_fmt, node_labels={}, # Remove all node labels for a clearer viz style=css_string, # Apply the stylesheet ) ``` Colouring nodes by population makes it immediately clear that, while the tree structure does not exactly reflect the population divisions, there's still considerable population substructure present in this larger tree. :::{todo} The (older) {meth}`Tree.draw` function also has a `node_colour` argument that can be used to colour tree nodes, which is used in some of the other tskit tutorials. Under the hood, this function simply sets appropriate SVG styles on nodes. We intend to make it easier to set colours in a similar way: see https://github.com/tskit-dev/tskit/issues/579. ::: The CSS string used to style the tree above takes advantage of the general classes defined in a tskit SVG file: a node symbol always has a class named `sym`, which is contained within a [grouping element](https://www.w3.org/TR/SVG2/struct.html#Groups) of class `node`. Moreover, elements such as `node` have *additional* classes, such as `p1`, indicating that the node in this case belongs to the population with ID 1. #### Available css classes Here are the css classes in a tskit SVG which can be used to style specific elements. ##### Within the plotting area * `tree`: a grouping element containing each tree * `node`: a grouping element within a tree, containing a node and its descendant elements such as a node symbol, an edge, mutations, and other nodes. * `mut`: a grouping element containing a mutation symbol and label * `extra`: an extra class for mutations {ref}`outside the tree ` * `lab`: a label element (for a node, mutation, axis, tick number, etc.) * `sym`: a symbol element (e.g. a node, mutation, or site symbol) * `edge`: an edge element (i.e. a branch in a tree) * `root`, `leaf` and `sample`: additional classes applied to a node group if the node is a root node, a leaf node or a sample node * `unknown_time`: a class added to `mut` groups if the time of the mutation is {data}`tskit.UNKNOWN_TIME`. * `rgt` and `lft`: additional classes applied to labels for left- or right-justification ##### Outside the plotting area * `axes`: a grouping element containing the X and Y axes, if either are present * `x-axis`, `y-axis`: more specific grouping elements contained within `axes` * `tick`: a single tick on an axis, containing a tickmark line and a label * `site`: a grouping element representing a site (plotted on the X axis), containing a site symbol (a tick line) and zero or more `mut` groups, each containing a chevron-shaped mutation symbol * `background`: the shaded background of a tree sequence plot * `grid`: a gridline ##### ID-based classes Elements have additional classes based on the IDs of trees, edges, nodes, parent (ancestor) nodes, individuals, populations, mutations, and sites. These class names start with a single letter (respectively `t`, `e`, `n`, `a`, `i`, `p`, `m`, and `s`) followed by a numerical ID. For example, here's a typical node in a tskit SVG plot: ``` ... ``` This corresponds to node 7, the rightmost leaf in the third tree in the mutated tree sequence (plotted in the previous section but one). The classes indicate that it has an immediate ancestor (parent) node with ID 10 (`a10`), and that the node belongs to an {ref}`individual ` with ID 3 (`i3`). The classes `n7` and `p2` tell us that the node ID is 7 and is from the population with ID 2 (`p2`). Other ID classes on the node tell us about the mutations above that node, of which there are two in this case, with IDs 16 and 17 (`m16`, `m17`); those mutations are associated with site IDs 15 and 16 (`s16`, `s17`). Other grouping elements apart from nodes can also contain ID-based classes. For example the `tree` group contains the ID of the tree (e.g. `t0`), the `site` group on the X axis contains the site ID (e.g. `s15`) the `mut` class contains the mutation ID (e.g. `m16`), and so on. #### CSS selector quick reference If you don't do this all the time it's not easy to remember what the various separators mean, so here's a quick reference (for more, see [these docs](https://developer.mozilla.org/en-US/docs/Web/CSS/CSS_selectors)): - `abc` (Type `abc`, like `g` for a `...` tag) - `.xyz` (Class `xyz`) - `#uvw` (ID `uvw`) - `,` (Selector list, means "or") - `>` (Child combinator) - "` `" (Descendant combinator, a space) - `+` (Next-sibling combinator) - `~` (Subsequent sibling combinator) - `|` (Namespace separator) #### Styling graphical elements The classes above make it easy to target specific nodes or edges in one or multiple trees. For example, we can colour branches that are shared between trees: ```{code-cell} ipython3 css_string = ".a15.n9 > .edge {stroke: cyan; stroke-width: 2px}" # branches from 15->9 ts_small.draw_svg(time_scale="rank", size=wide_fmt, style=css_string) ``` Or rather than identifying shared branches using the same parent and child, you can use the edge ID, which allows effects like colouring the edges by their span, to emphasize the routes through which more genomic information has been inherited. Below we do this for a single tree, which avoids making a styling rule for every edge in the tree sequence: ```{code-cell} ipython3 tree = ts_small.at_index(2) edges = tree.edge_array[tree.edge_array != tskit.NULL] spans = tree.tree_sequence.edges_right[edges] - tree.tree_sequence.edges_left[edges] values = ((1-spans/spans.max()) * 255).astype(int) style = "".join([f".edge.e{e} {{stroke:#{v:02X}{v:02X}{v:02X}}}" for e, v in zip(edges, values)]) tree.draw_svg(style=style + ".edge {stroke-width: 2px}") ``` By generating the css string programatically, you can target all the edges present in a particular tree, and see how they gradually disappear from adjacent trees. Below, for example the branches in the central tree have been coloured red, as have the identical branches in adjacent trees. The central tree represents a location in the genome that has seen a selective sweep, and therefore has short branch lengths: adjacent trees are not under direct selection and thus the black branches tend to be longer. These (red) shared branches extending far on either side represent shared haplotypes, and this shows how long, shared haplotypes can extend much further away from a sweep than the region of reduced diversity (which is the region spanned by the short tree in the middle). For visual clarity, node symbols and labels have been turned off. ```{code-cell} ipython3 :"tags": ["hide-input"] ts_selection = tskit.load("data/viz_ts_selection.trees") ``` ```{code-cell} ipython3 css_edge_targets = [] # collect the css targets of all the edges in the selected tree sweep_location = ts_selection.sequence_length / 2 # NB: sweep is in the middle of the ts focal_tree = ts_selection.at(sweep_location) for node_id in focal_tree.nodes(): parent_id = focal_tree.parent(node_id) if parent_id != tskit.NULL: css_edge_targets.append(f".a{parent_id}.n{node_id}>.edge") css_string = ",".join(css_edge_targets) + "{stroke: red} .sym {display: none}" css_string += ( # Rotate the position labels etc ".x-axis .ticks .lab {text-anchor: start; transform: translate(6px) rotate(90deg)}" ".x-axis .title .lab {text-anchor: start}" ) wide_tall_fmt = (1200, 400) ts_selection.draw_svg( style=css_string, size=wide_tall_fmt, canvas_size=(wide_tall_fmt[0], wide_tall_fmt[1] + 30), x_lim=[1e4, 4e4], node_labels={}, ) ``` :::{note} Branches in multiple trees that have the same parent and child do not always correspond to a single {ref}`edge ` in a tree sequence: for example, edges have the additional constraint that they must belong to _adjacent_ trees. ::: (sec_tskit_viz_dynamic_effects)= #### Dynamic effects In the previous example, the large size of the plotted trees meant that, for clarity, the node and mutation labels were turned off (in that case by passing empty mappings to the `node_labels` and `mutation_labels` parameters). Nevertheless, it can be useful to identify nodes and mutations, and this can be done dynamically (on "mouseover") by setting the CSS [display](https://www.w3.org/TR/SVG2/render.html#VisibilityControl) property to `none` vs `initial`, and combining it with the CSS `:hover` pseudoclass. Here's an example using a region from within the previous example: ```{code-cell} ipython3 from IPython.display import HTML # add some mutations ts = msprime.sim_mutations(ts_selection, rate=2e-8, random_seed=1) node_label_css = ( # hide node labels by default "#hover_example .node > .sym ~ .lab {display: none}" # Unless the adjacent node or the label is hovered over "#hover_example .node > .sym:hover ~ .lab {display: inherit}" "#hover_example .node > .sym ~ .lab:hover {display: inherit}" ) mut_label_css = ( # hide mutation labels by default "#hover_example .mut .sym ~ .lab {display: none}" # Unless the adjacent node or the label is hovered over "#hover_example .mut .sym:hover ~ .lab {display: inherit}" "#hover_example .mut .sym ~ .lab:hover {display: inherit}" ) optional_css = ( # These are optional, but setting e.g. the node label text to bold with grey stroke # and a black fill, serves to make black text readable against a black tree "svg#hover_example {background-color: white}" "#hover_example .tree .plotbox .lab {stroke: #CCC; fill: black; font-weight: bold}" "#hover_example .tree .mut .lab {stroke: #FCC; fill: red; font-weight: bold}" ) HTML(ts.draw_svg( style=optional_css + node_label_css + mut_label_css, y_axis=True, y_ticks={0: "0", 500: "", 1000: "1000"}, x_lim=[2.3e4, 2.7e4], root_svg_attributes={"id": "hover_example"}, # Label node by name in metadata, if it exists, else node ID node_labels={u.id: u.metadata.get("name", f"NodeID={u.id}") for u in ts.nodes()}, mutation_labels={ # Label mutation by site position, prev state, and new state m.id: ( f"pos {s.position:g}: " + (s.ancestral_state if m.parent<0 else ts.mutation(m.parent).derived_state) + f"→{m.derived_state}" ) for s in ts.sites() for m in s.mutations }, )) ``` :::{note} Above we have wrapped the svg in an IPython {class}`~ipython:IPython.display.HTML` class, and given the SVG a unique id as described below in {ref}`sec_tskit_viz_styling_more_about`. This forces the SVG plot to be rendered inline (rather than inside an `` tag), allowing the hover functionality to work in all supported Jupyter notebook implementations. However, depending on your Jupyter setup, the `HTML()` wrapper may not be necessary. ::: Using the transformations discussed in the next section, it is also possible to animate SVG images, as shown in the {ref}`sec_tskit_viz_SVG_examples_animation` code within the {ref}`sec_tskit_viz_SVG_examples` section near the end of this tutorial. (sec_tskit_viz_styling_transforming_and_masking)= #### Transforming and masking elements We can also use styles to transform elements of the drawing, shifting them into different locations or changing their orientation. For example, {ref}`earlier in this tutorial ` we used the following CSS string to rotate leaf labels: ```css .leaf > .lab {text-anchor: start; transform: rotate(90deg) translate(6px)} ``` Transformations not only allow us to shift elements about, but also resize and skew them. When applied to both symbols and labels this can create rather different formatting styles: ```{code-cell} ipython3 css_string = ( # Draw large yellow circles for nodes ... ".node > .sym {transform: scale(2.2); fill: yellow; stroke: black; stroke-width: 0.5px}" # ...but for leaf nodes, override the yellow circle using a more specific CSS target ".node.leaf > .sym {transform: scale(1); fill:black}" # Override default node text position to be based at (0, 0) relative to the node pos # Note that the .tree specifier is needed to make this more specific than the default # positioning which is targeted at ".lab.lft" and ".lab.rgt" ".tree .node > .lab {transform: translate(0, 0); text-anchor: middle; font-size: 7pt}" # For leaf nodes, override the above positioning using a subsequent CSS style ".node.leaf > .lab {transform: translate(0, 12px); font-size: 10pt}" ) ts_small.first().draw_svg(style=css_string) ``` Note that when transforming elements, parts of the drawing may be plotted outsize of of the standard canvas, so the `canvas_size` option is particularly useful when performing more radical CSS transformations, for example to create {ref}`sec_tskit_viz_SVG_examples_3D`: ```{code-cell} ipython3 skew = 0.8 # How skewed the trees are, in radians # CSS transforms used to skew the trees style = f".tree .plotbox {{transform: skewY({skew}rad)}}" # Shift the x axis to make room for the skewed trees shift_axis = (10, 50) style += f".x-axis {{transform: translate({shift_axis[0]}px, {shift_axis[1]}px)}}" # Must define a bigger canvas size so we don't crop the axis off size = (800, 200) # width, height of svg canvas_size = (size[0] + shift_axis[0], size[1] + shift_axis[1]) ts_tiny.draw_svg(size=size, x_scale="treewise", style=style, canvas_size=canvas_size) ``` :::{note} Using `transform` in styles is an SVG2 feature, and has not yet been implemented in the software programs Inkscape or librsvg. Therefore if you are {ref}`converting or editing ` the examples above, the transformed elements may be positioned incorrectly. For changing symbols, the `symbol_size` option can be used to simply change the size of all symbols in the plot, but otherwise you may need to use the `chromium` workaround documented {ref}`here `. ::: Although it is hard to change the style of a node symbol, the visible area of the symbol can be modified using the `clip-path` CSS property. This can be useful to show, for instance, a triangle to summarise the descendants of a MRCA. ```{code-cell} ipython3 # Check that MRCA of 2 & 3 is node 4 in all trees, assumed later assert all([4 == tree.mrca(2, 3) for tree in ts_tiny.trees()]) styles = [ # Set all node labels to be rotated and small ".node > .lab {text-anchor: start; transform: rotate(90deg) translate(6px); font-size: 8px}", # Hide all nodes descending from node 4. We then treat node 4 as a summary node ".n4 > .node {display: none}", # Use clipping & scaling to change the symbol for node 4 into a summary triangle ".n4 > .sym {clip-path: polygon(50% 50%, 75% 75%, 25% 75%); transform: scale(8.0, 8.0)}", # Make the font bigger for this summary node label ".n4 > .lab {transform: rotate(90deg) translate(14px); font-size: 16px}" ] node_labels = {0: "Nd. 0", 1: "Nd. 1", 4: "Two samples"} ts_tiny.draw_svg( size=(800, 300), x_scale="treewise", time_scale="log_time", style="".join(styles), node_labels=node_labels, ) ``` In the example above we simply hid the descendant topology for each "summary MRCA", meaning more horizontal space was taken up than expected. For a more sophisticated example, see {ref}`sec_tskit_viz_SVG_examples_larger_plots` in which some descendant samples are actually removed from the tree sequence entirely, and their MRCA is changed into a sample node instead. #### Styling and SVG structure To take full advantage of the SVG styling capabilities in tskit, it is worth knowing how the SVG file is structured. In particular tskit SVGs use a hierarchical grouping structure that reflects the tree topology. This allows easy styling and manipulation of both individual elements and entire subtrees. Currently, the hierarchical structure of a simple 2-tip SVG tree produced by tskit looks something like this: ``` Mutation 0 Node 1 Node 0 Root (Node 2) ``` And in a tree sequence plot, the SVG simply consists of a set of such trees, together with groups containing the background and axes, if required. ``` ... ... ... ... ``` #### Styling subtrees The nested grouping structure makes it easy to target a node and all its descendants. For instance, here's how to draw all the edges of node 13 and its descendants using a thicker blue line: ```{code-cell} ipython3 edge_style = ".n13 .edge {stroke: blue; stroke-width: 2px}" nd_labs = {n: n for n in [0, 1, 2, 3, 4, 5, 6, 7, 13, 17]} ts_mutated.draw_svg(x_lim=x_limits, node_labels=nd_labs, style=edge_style) ``` This might not be quite what you expected: the branch leading from node 13 to its parent (node 17) has also been coloured. That's because the SVG node group deliberately contains the branch that leads to the parent (this can be helpful, for example, for hiding the entire subtree leading to node 13, using e.g. `.n13 {display: none}`). To colour the branches *descending* from node 13, you therefore need to target the nodes nested at least one level deep within the `n13` group. One way to do that is to add an extra `.node` class to the style, e.g. ```{code-cell} ipython3 edge_style = ".n13 .node .edge {stroke: blue; stroke-width: 2px}" # NB to target the edges in only (say) the 1st tree you could use ".t0 .n13 .node .edge ..." ts_mutated.draw_svg(x_lim=x_limits, node_labels=nd_labs, style=edge_style) ``` If you want to colour the branches descending from a particular mutation (say mutation 7) then you need to colour not only the edges, but also *part* of an edge (i.e. the line that connects a mutation downwards to its associated node). The tskit SVG format provides a special `` element to enable this, which is normally made invisible using `fill: none stroke: none` in the default stylesheet. Here's an example of activating this normally-hidden line: ```{code-cell} ipython3 default_muts = ".mut .lab {fill: gray} .mut .sym {stroke: gray}" # all other muts in gray m8_mut = ( ".m8 .node .edge, " # the descendant edges ".mut.m8 line, " # activate the hidden line between the mutation and the node ".mut.m8 .sym " # the mutation symbols on the tree and the axis "{stroke: red; stroke-width: 2px}" ".mut.m8 .lab {fill: red}" # colour the label "8" in red too ) css_string = default_muts + m8_mut ts_mutated.draw_svg(x_lim=x_limits, node_labels=nd_labs, style=css_string) ``` If you want to colour the branches *above* a node, you can use the :meth:`~tskit.Tree.ancestors` method, and chain the css specifiers together using `,`: ```{code-cell} ipython3 focal_node = 9 tree = ts_mutated.first() target_nodes = [focal_node] + list(tree.ancestors(focal_node)) css_string = ",".join([f".node.n{u} > .edge" for u in target_nodes]) + "{stroke: red}" tree.draw_svg(style=css_string, omit_sites=True) ``` #### Restricting styling Sometimes the hierarchical nesting leads to styles being applied too widely. For example, since style selectors include all the descendants of a target, to target just the node itself (and not its descendants) a slightly different specification is required, involving, the "`>`" symbol, or [child combinator](https://www.w3.org/TR/selectors-3/#child-combinators) (we have, in fact, used it in several previous examples). The following plot shows the difference when all decendant symbols are targeted, versus just the immediate child symbol: ```{code-cell} ipython3 node_style1 = ".n13 .sym {fill: yellow}" # All symbols under node 13 node_style2 = ".n15 > .sym {fill: cyan}" # Only symbols that are an immediate child of node 15 css_string = node_style1 + node_style2 ts_small.draw_svg(y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits, style=css_string) ``` Another example of modifying the style target is *negation*. This is needed, for example, to target nodes that are *not* leaves (i.e. internal nodes). One way to do this is to target *all* the node symbols first, then replace the style with a more specific targeting of the leaf symbols only: ```{code-cell} ipython3 hide_internal_symlabs = ".node > .sym, .node > .lab {display: none}" show_leaf_symlabs = ".node.leaf > .sym, .node.leaf > .lab {display: initial}" css_string = hide_internal_symlabs + show_leaf_symlabs ts_small.draw_svg(y_axis=True, y_ticks=y_tick_pos, x_lim=x_limits, style=css_string) ``` Alternatively, the `:not` selector can be used to target nodes that are *not* leaves, so the following style specification should produce the same effect in SVG viewers that support it (note, however, as of v1.2 Inkscape does not appear to support this selector). ``` style_string = ".node:not(.leaf) > .sym, .node:not(.leaf) > .lab {display: none}" ``` (sec_tskit_viz_styling_more_about)= #### More about styling NOTE: if your SVG is embedded directly into an HTML page (a common way for jupyter notebooks to render SVGs), then according to the HTML specifications, any styles applied to one SVG will apply to all SVGs in the document. To avoid this confusing state of affairs, we recommend that you tag the SVG with a unique ID using the `root_svg_attributes` parameter, then prepend this ID to the style string: ```{code-cell} ipython3 ts_small.draw_svg( x_lim=x_limits, root_svg_attributes={'id': "myUID"}, style="#myUID .background * {fill: #00FF00}", # apply any old style to this specific SVG ) ``` SVG styles allow a huge amount of flexibility in formatting your plot, even extending to animations. Feel free to browse the {ref}`examples ` for inspiration. (sec_tskit_viz_converting)= ### Converting and editing SVG #### Converting [Inkscape](https://inkscape.org) is an open source SVG editor that can also be scripted to output bitmap files. [Imagemagick](https://imagemagick.org/) is a common piece of software used to convert between image formats. It can be configured to delegate to one of several different SVG libraries when converting SVGs to bitmap formats. Currently, both the librsvg library or the Inkscape library produce reasonable output, although librsvg currently misaligns some labels due to [ignoring certain SVG properties](https://gitlab.gnome.org/GNOME/librsvg/-/issues/414). (sec_tskit_viz_converting_note)= :::{note} A few stylesheet specifications, such as the `transform` property, are SVG2 features, and have not yet been implemented in Inkscape or librsvg. Therefore if you use these in your own custom SVG stylesheet (such as the example above where we rotated leaf labels), they will not be applied properly when converted with those tools. For custom stylesheets like this, a workaround is to convert the SVG to PDF first, using e.g. the programmable chromium engine: ``` chromium --headless --print-to-pdf=out.pdf in.svg ``` The resulting PDF file can then be converted by Inkscape, retaining the correct transformations. ::: #### Editing the SVG - Editing can be done in [Inkscape](https://inkscape.org) (subject to the note above) :::{todo} Tips on how to cope with the hierarchical grouping when editing (e.g. in Inkscape using `Extensions menu > Arrange > Deep Ungroup`, but note that this will mess with the styles!) ::: (sec_tskit_viz_examples)= ## Examples ### Text examples #### Tree orientation In the text format, trees (but not tree sequences) can be displayed in different orientations ```{code-cell} ipython3 :"tags": ["hide-input"] from IPython.display import HTML import msprime ts = msprime.sim_ancestry(4) orient = "top", "left", "bottom", "right" html = [] for o in orient: tree_string = ts.first().draw_text(orientation=o) html.append(f"
{tree_string}
") HTML((" "*10).join(html)) ``` (sec_tskit_viz_SVG_examples)= ### SVG examples #### A standard ts plot Note that this tree sequence also illustrates a few features which are not normally produced e.g. by `msprime` simulations, in particular a "empty" site (with no associated mutations) at position 50, and some mutations that occur above root nodes in the trees. Graphically, root mutations necessitate a line above the root node on which to place them, so each tree in this SVG has a nominal "root branch" at the top. Normally, root branches are not drawn, unless the `force_root_branch` parameter is specified. ```{code-cell} ipython3 :"tags": ["hide-input"] # Make a tree sequence with multiple mutations including some above the root import io import tskit def make_unusual_ts(): nodes = io.StringIO( """\ id is_sample population individual time metadata 0 1 0 -1 0 1 1 0 -1 0 2 1 0 -1 0 3 1 0 -1 0 4 0 0 -1 0.1145014598813 5 0 0 -1 1.11067965364865 6 0 0 -1 1.75005250750382 7 0 0 -1 5.31067154311640 8 0 0 -1 6.57331354884652 9 0 0 -1 9.08308317451295 """ ) edges = io.StringIO( """\ id left right parent child 0 0 100 4 0 1 0 100 4 1 2 0 100 5 2 3 0 100 5 3 4 80 85 6 4 5 80 85 6 5 6 6 80 7 4 7 85 91 7 4 8 6 80 7 5 9 85 91 7 5 10 91 100 8 4 11 91 100 8 5 12 0 6 9 4 13 0 6 9 5 """ ) sites = io.StringIO( """\ position ancestral_state 4 A 6 0 30 Empty 50 XXX 91 T """ ) muts = io.StringIO( """\ site node derived_state parent time 0 9 T -1 15 0 9 G 0 9.1 0 5 1 1 9 1 4 C -1 1.6 1 4 G 3 1.5 2 7 G -1 10 2 3 C 5 1 4 3 G -1 1 """ ) return tskit.load_text(nodes, edges, sites=sites, mutations=muts, strict=False) ts = make_unusual_ts() ts.draw_svg() ``` #### Highlighted mutations Specific mutations can be given a different colour. Moreover, the descendant lineages of specific mutations can be coloured and the branch colours overlay each other as expected. Note that in this example, internal node labels and symbols have been hidden for clarity. ```{code-cell} ipython3 :"tags": ["hide-input"] ts = make_unusual_ts() # Defined in the first SVG example css_string = ( ".edge {stroke: grey}" ".mut .sym{stroke:pink} .mut text{fill:pink}" ".mut.m2 .sym, .m2>line, .m2>.node .edge{stroke:blue} .mut.m2 .lab{fill:blue}" ".mut.m3 .sym, .m3>line, .m3>.node .edge{stroke:cyan} .mut.m3 .lab{fill:cyan}" ".mut.m4 .sym, .m4>line, .m4>.node .edge{stroke:red} .mut.m4 .lab{fill:red}" # Hide internal node labels & symbols ".node:not(.leaf) > .sym, .node:not(.leaf) > .lab {display: none}" ) ts.draw_svg(style=css_string, time_scale="rank", x_lim=[0, 30]) ``` #### Leaf, sample & isolated nodes By default, sample nodes are square and non-sample nodes circular (at the moment this can't easily be changed). However, neither need to be at specific times: sample nodes can be at times other than 0, and nonsample nodes can be at time 0. Moreover, leaves need not be samples, and samples need not be leaves. Here we change the previous tree sequence to make some leaves non-samples and some samples internal nodes. To highlight the change, we have plotted sample nodes in green, and leaf nodes (if not samples) in blue. ```{code-cell} ipython3 :"tags": ["hide-input"] def swap_samples_ts(): tables = make_unusual_ts().dump_tables() # Defined in the first SVG example tables.mutations.clear() tables.sites.clear() flags = tables.nodes.flags flags[:] = [0, 0, 1, 1, 0, 0, 0, 1, 1, 0] tables.nodes.flags = flags return tables.tree_sequence() ts = swap_samples_ts() css_string=".leaf .sym {fill: blue} .sample > .sym {fill: green}" ts.draw_svg( style=css_string, x_scale="treewise", time_scale="rank", y_axis=True, y_gridlines=True, x_lim=[0, 10], ) ``` :::{note} By definition, if a node is a sample, it must be present in every tree. This means that there can be sample nodes which are "isolated" in a tree. These are drawn unconnected to the main topology in one or more trees (e.g. nodes 7 and 8 above). ::: #### A fancy formatted plot Here we have activated the Y axis, and changed the node style. In particular, we have coloured nodes by time, and increased the internal node symbol size while moving the internal node labels into the symbol; node labels have also been plotted in a sans-serif font. Axis tick labels have been changed to avoid potential overlapping (some Y tick labels have been removed, and the X tick labels rotated). ```{code-cell} ipython3 :"tags": ["hide-input"] import msprime import numpy as np def make_standard_ts(): seed=370009 return msprime.sim_ancestry( 7, ploidy=1, sequence_length=1000, random_seed=seed, recombination_rate=0.001) ts = make_standard_ts() # Thin the tick values so we don't get labels within 0.01 of each other y_ticks = ts.tables.nodes.time y_ticks = np.delete(y_ticks, np.argwhere(np.ediff1d(y_ticks) <= 0.01)) css_string = ( ".tree .lab {font-family: sans-serif}" # Normal X axis tick labels have dominant baseline: hanging, but it needs centring when rotated + ".x-axis .tick .lab {text-anchor: start; dominant-baseline: central; transform: rotate(90deg)}" + ".y-axis .grid {stroke: #DDDDDD}" + ".tree :not(.leaf).node > .lab {transform: translate(0,0); text-anchor:middle; fill: white}" + ".tree :not(.leaf).node > .sym {transform: scale(3.5)}" + "".join( f".tree .n{n.id} > .sym {{fill: hsl({int((1-n.time/ts.max_root_time)*260)}, 50%, 50%)}}" for n in ts.nodes() ) ) ts.draw_svg(size=(1000, 350), y_axis=True, y_gridlines=True, y_ticks=y_ticks, style=css_string) ``` #### Ticks, labels, and gridlines Y tick labels can be specified explicitly, which allows time scales to be plotted e.g. in years even if the tree sequence ticks in generations. The ``title`` class allows axis titles to be moved out of the way of tick labels. Finally, grid lines associated with each y tick can also be changed or even hidden individually using the CSS [nth-child pseudo-selector](https://www.w3.org/TR/2018/REC-selectors-3-20181106/#nth-child-pseudo), where tickmarks are indexed from the bottom. Below is an example of all 3 techniques, drawing on an example from the {ref}`sec_msprime_introgression` tutorial: ```{code-cell} ipython3 :"tags": ["hide-input"] import msprime time_units = 1000 / 25 # Conversion factor for kya to generations def make_introgression_ts(sequence_length, random_seed=None): """ Function from the introgression tutorial - see there for justification """ demography = msprime.Demography() # The same size for all populations; highly unrealistic! Ne = 10**4 demography.add_population(name="Africa", initial_size=Ne) demography.add_population(name="Eurasia", initial_size=Ne) demography.add_population(name="Neanderthal", initial_size=Ne) # 2% introgression 50 kya demography.add_mass_migration( time=50 * time_units, source='Eurasia', dest='Neanderthal', proportion=0.02) # Eurasian 'merges' backwards in time into Africa population, 70 kya demography.add_mass_migration( time=70 * time_units, source='Eurasia', dest='Africa', proportion=1) # Neanderthal 'merges' backwards in time into African population, 300 kya demography.add_mass_migration( time=300 * time_units, source='Neanderthal', dest='Africa', proportion=1) return msprime.sim_ancestry( recombination_rate=1e-8, sequence_length=sequence_length, samples=[ msprime.SampleSet(1, ploidy=1, population='Africa'), msprime.SampleSet(1, ploidy=1, population='Eurasia'), # Neanderthal sample taken 30 kya msprime.SampleSet(1, ploidy=1, time=30 * time_units, population='Neanderthal'), ], demography = demography, random_seed=random_seed, ) ts = make_introgression_ts(20 * 10**6, random_seed=1) base_size = (1200, 500) x_shift = 60 css = f".tree-sequence {{transform: translateX({x_shift}px)}}" # Move rightwards css += ".sample .lab {text-anchor: start; transform: rotate(90deg) translate(6px); font-size: 80%}" css += ".x-axis .tick .lab {font-size: 85%}" css += ".y-axis .title {transform: translateY(250px)}" # Move out of the way of y_tick labels css += ".y-axis .tick .grid {stroke: lightgrey}" # Default gridline type css += ".y-axis .ticks .tick:nth-child(3) .grid {stroke-dasharray: 4}" # 3rd line from bottom css += ".y-axis .ticks .tick:nth-child(3) .grid {stroke: magenta}" # also 3rd line from bottom css += ".y-axis .ticks .tick:nth-child(4) .grid {stroke: blue}" # 4th line from bottom css += ".y-axis .ticks .tick:nth-child(5) .grid {stroke: darkgrey}" # 5th line from bottom y_ticks = { 0: "0", 30: "30", 50: "Introgression event", 70: "European origin", 300: "Neanderthal origin", 1000: "1000", } ts.draw_svg( size=base_size, x_lim=(0, 25_000), time_scale="log_time", node_labels = {0: "Africa", 1: "Europe", 2: "Neanderthal"}, y_axis=True, y_ticks={y * time_units: lab for y, lab in y_ticks.items()}, y_gridlines=True, style=css, canvas_size=(base_size[0] + x_shift, base_size[1]), ) ``` (sec_tskit_viz_SVG_examples_larger_plots)= #### Simplifying larger plots It is common to want to visualise a tree sequence with many samples and trees. If there are many trees, the `max_num_trees` parameter can be used to just show those at the start and end of the genome. To reduce the size of each tree, multiple samples can be clustered into a single representative clade. If that clade has the same set of descendant samples throughout the tree sequence, {ref}`sec_simplification` can be used to turn the MRCA of these samples into a sample node itself, while removing the original descendants. By using the scaling and masking method described in {ref}`sec_tskit_viz_styling_transforming_and_masking` this summary MRCA can be shown as a large triangle, of size proportional to the number of samples underneath it. The example below shows how a tree sequence of 40 sample nodes can be visualised relatively compactly using these techniques: ```{code-cell} ipython3 :"tags": ["hide-input"] import msprime import numpy as np def clonal_mrcas(ts, most_recent=True): """ Identify the nodes in a tree sequence which define clonal subtrees (i.e. in which the samples descending from that node are identical and show identical relationships to each other over the entire tree sequence). This includes, at its limit, nodes with only a single descendant sample. :param bool most_recent: If True, and the clonal node is a unary node, return IDs of the most recent node that defines the clonal subtree. In this case, the returned IDs represent cases where the node is either a tip or a coalescent point. :return: a list of nodes defining constant subtrees over the entire tree sequence :rtype: list """ for interval, edges_out, edges_in in ts.edge_diffs(): if interval.left==0: is_full_length_clonal = np.ones(ts.num_nodes, dtype=bool) # nodes start clonal else: for e in edges_in: is_full_length_clonal[e.parent] = False for e in edges_out: is_full_length_clonal[e.parent] = False clonal_nodes = np.where(is_full_length_clonal)[0] tables = ts.dump_tables() edges = tables.edges # only keep edges where both the child and the parent are full-length clonal nodes keep_edge = np.logical_and( np.isin(edges.child, clonal_nodes), np.isin(edges.parent, clonal_nodes), ) tables.edges.set_columns( left = tables.edges.left[keep_edge], right=tables.edges.right[keep_edge], parent=tables.edges.parent[keep_edge], child=tables.edges.child[keep_edge], ) clonal_ts = tables.tree_sequence() assert clonal_ts.num_trees == 1 # Also remove all the edges ascending from removed edges tree = clonal_ts.first() non_clonal_ancestors = set() deleted_edges = np.logical_not(keep_edge) for u in np.unique(ts.edges_parent[deleted_edges]): while u != tskit.NULL and u not in non_clonal_ancestors: non_clonal_ancestors.add(u) u = tree.parent(u) non_clonal_ancestors = np.array(list(non_clonal_ancestors)) tables = clonal_ts.dump_tables() remove_edge = np.isin(tables.edges.parent, non_clonal_ancestors) tables.edges.replace_with(tables.edges[np.logical_not(remove_edge)]) clonal_ts = tables.tree_sequence() tree = ts.first(sample_lists=True) clonal_tree = clonal_ts.first(sample_lists=True) clonal_nodes = [] for root in clonal_tree.roots: # Clonal trees should subtend the same set of samples as in the original tree assert set(tree.samples(root)) == set(clonal_tree.samples(root)) u = root if most_recent: # decend to the first coalescent node (i.e. MRCA) while clonal_tree.num_children(u) == 1: u = clonal_tree.children(u)[0] clonal_nodes.append(u) return clonal_nodes ts = msprime.sim_ancestry( 20, population_size=1e2, sequence_length=1e4, recombination_rate=1e-6, random_seed=83) clones = clonal_mrcas(ts) # Simplify but keep the same node IDs using filter_nodes=False ts_simp = ts.simplify(clones, filter_nodes=False) styles = [ ".node > .lab {font-size: 9px}", ".leaf > .lab {text-anchor: start; transform: rotate(90deg) translate(5px); font-size: 12px}", ] node_labels = {u: u for u in range(ts.num_nodes)} for u in ts.samples(): node_labels[u] = f"S{u}" tree = ts.first() for u in clones: num_samples = tree.num_samples(u) if num_samples > 1: node_labels[u] = f"{num_samples} samples" styles.append( f".n{u} > .sym {{clip-path: polygon(50% 50%, 100% 100%, 0% 100%);"+ f"transform: scale({(num_samples-1)/5 + 1}, 4.0)}}" + f".n{u} > .lab {{transform: rotate(90deg) translate(15px); font-size: 13px}}" ) ts_simp.draw_svg( size=(1000, 400), style="".join(styles), node_labels=node_labels, time_scale="log_time", y_axis=True, y_ticks=[0, 1, 10, 100], max_num_trees=4) ``` (sec_tskit_viz_SVG_examples_3D)= #### 3D effects We can use various CSS transforms, as {ref}`discussed previously`, to skew the trees and stagger them. With a bit of trigonometry, this can create flexible and tolerably good 3D effects for presentations, etc. ```{code-cell} ipython3 :"tags": ["hide-input"] import math import msprime def make_7_tree_4_tip_ts(): ts = msprime.sim_ancestry( 4, ploidy=1, random_seed=889, sequence_length=1000, recombination_rate=0.001) ts = msprime.sim_mutations(ts, rate=2e-3, random_seed=123) # Check we have picked a random seed that gives a nice plot of 7 trees tip_orders = { tuple(u for u in t.nodes(order="minlex_postorder") if t.is_sample(u)) for t in ts.trees() } topologies = {tree.rank() for tree in ts.trees()} assert tip_orders == {(0, 1, 2, 3)} and len(topologies) > 1 and ts.num_trees == 7 return ts ts = make_7_tree_4_tip_ts() # Set some parameters: these can be adjusted to your liking tree_width = 100 height = 200 # Normal height for tree + x-axis y_step = 40 # Stagger between trees (i.e. 0 for all trees in a horizontal line) skew = 0.6 # How skewed the trees are, in radians width = tree_width * ts.num_trees + 20 + 20 # L & R margins in draw_svg = 20px angle = math.atan(y_step/tree_width) ax_mv = y_step, (ts.num_trees - 1) * y_step + math.tan(skew) * (tree_width * .9) # CSS transforms used to skew the axis and stagger + skew the trees style = f".x-axis {{transform: translate({ax_mv[0]}px, {ax_mv[1]}px) skewY(-{angle}rad)}}" for i in range(ts.num_trees): # Stagger each tree vertically by y_step, transforming the "plotbox" tree container style += ( f".tree.t{i} > .plotbox " + "{transform:" + f"translateY({(ts.num_trees - i - 1) * y_step}px) skewY({skew}rad)" + "}" ) # Define a bigger canvas size so we don't crop the moved trees from the drawing size = (width, height) canvas_size = (width + y_step, height + ts.num_trees*y_step + math.tan(skew)*tree_width) ts.draw_svg(size=size, x_scale="treewise", style=style, canvas_size=canvas_size) ``` (sec_tskit_viz_legend_example)= #### Legend example A classic case for a legend is to explain node formatting, e.g. when nodes are coloured by their population. Here's an example of how it can be done: ```{code-cell} ipython3 :"tags": ["hide-input"] ts = tskit.load("data/viz_ts_full.trees") styles = [ f".node.p{p.id} > .sym " + "{" + f"fill: {colour}" + "}" for colour, p in zip(['red', 'green', 'blue'], ts_full.populations()) ] # styles created as in previous examples # Make a legend: first create a surrounding box and a title legend = '' legend += 'Key' # Now make the legend lines, one for each population. Setting classes that match those # used for normal nodes means that styled colours are automatically picked-up. legend += "".join([ f'' # an SVG group f'' # Square symbol f'Population {p.metadata["name"]} (id={p.id})' # Label for p in ts_full.populations() if p.id < 3 ]) ts.first().draw_svg( size=(1200, 250), node_labels={}, # Remove all node labels for a clearer viz style="".join(styles), # Apply the stylesheet preamble=legend, ) ``` (sec_tskit_viz_svg_plot_internals)= #### SVG plot internals Internally, SVG plots are produced by creating an object of class `.drawing.SvgTree` or `.drawing.SvgTreeSequence`, then calling the `.draw()` method on the object. These internal classes and methods are deliberately undocumented, as they may be subject to change at any time. Nevertheless, accessing them can be useful e.g. to extract internal settings, such as the X and Y position of point in SVG space, for more complex programmatic annotations. ```{code-cell} ipython3 :"tags": ["hide-input"] ## WARNING - this code uses internal functions that may change in future tskit versions def node_positions(svgtree): # the internal `node_x_coord` dict maps node IDs to horizontal pos x = svgtree.node_x_coord # the internal `timescaling.transform` function returns vertical positions node_times = svgtree.ts.nodes_time # grab node times from the associated ts if svgtree.time_scale == "rank": # NB: the "rank" scale doesn't use raw ages within_tree_node_times = np.unique(node_times[svgtree.tree.preorder()]) node_times = np.searchsorted(within_tree_node_times, node_times) y = svgtree.timescaling.transform(node_times) return x, y # Create an SvgTree object using the same parameters as for draw_svg() internal_obj = tskit.drawing.SvgTree( ts_small.first(), canvas_size=(300, 200), time_scale="log_time", y_axis=True, y_ticks=[0, 10, 100, 1000], node_labels={}, ) horiz_node_pos, vert_node_pos = node_positions(internal_obj) x10, y10 = horiz_node_pos[10], vert_node_pos[10] # Node 10 x8, y8 = horiz_node_pos[8], vert_node_pos[8] # Node 8 internal_obj.preamble = ( f'' f'This is node 10!' f'' f'This is node 8!' ) internal_obj.draw() ``` Coupled with rotations etc, this allows some quite sophisticated viz possibilities, such as creating a "tanglegram" to compare two trees in a tree sequence: ```{code-cell} ipython3 :"tags": ["hide-input"] ## WARNING: uses internal tskit apis which could change def tanglegram( ts, tree_indexes=None, titles=None, *, size=None, order=None, separation=None, line_gap=None, node_labels=None, style=None, x_axis=None, x_label=None, x_ticks=None, x_gridlines=None, y_axis=None, y_label=None, y_ticks=None, y_gridlines=None, return_node_maps=None, **kwargs ): r""" Create an SvgTree object describing a "tanglegram" that compares the topology leading to leaf nodes on two trees in a tree sequence, by drawing them facing each other, and plotting lines between the same leaves in each tree. The object can be turned into an SVG string or drawn using `returned_obj.draw()`. The separate plots can be styled as the left and right SVGs are given `lft_tree` and ``rgt_tree` classes. Node IDs in any stylesheet will need to be remapped by specifying `return_mappings=True` and using the returned mappings to change the indexes. By default, plot titles based on the tree indexes and compare the first and last non-empty tree. If want to plot two trees from separate tree sequences, you can concatenate them together using :meth:`tskit.TreeSequence.concatenate()`, providing a node_mapping to match leaf or sample nodes between the two trees, and then pass the resulting tree sequence to this function. .. note:: This does not "untangle" the trees to minimise the number of line crossings, but simply plots them using the default "minlex" order. If you do wish to untangle them, you will need to use an external program to calculate the leaf orders, and pass them in through the `order` parameter. :param TreeSequence ts: The tree sequence from which to take trees :param tuple tree_indexes: A tuple of two integers, the indexes of the two trees to compare. By default take the first and last non-empty trees in the tree sequence. :param tuple titles: A tuple of two strings, the titles to use for the left and right trees. By default show "Tree X". To show no titles, provide titles=(None, None). :param tuple size: A tuple of two integers, the width and height of the SVG image. By default, use the default (square) size per tree, giving (400, 200) in total. :param tuple order: A tuple of two lists of integers, the order in which to plot the leaves. Either list can be `None`, meaning default to ``minlex_postorder``. If either is a list, it must contain unique integers corresponding to the IDs of the leaves in the corresponding tree, with the length of each list matching the number of leaves in the tree. Default: ``None`` treated as ``(None, None)``. :param float separation: The distance between the base of each tree. Default: ``None`` treated as a standard distance of 64px. :param float line_gap: The distance between tangle_lines and each leaf on the tree. If None, draw tangle lines of equal length in the middle of the plot :param str path: A path to which the SVG will be saved. :param dict node_labels: A dictionary mapping node IDs to label to plot. See :meth:`Tree.draw_svg()` for details. :param str style: a string of CSS styles. See :meth:`Tree.draw_svg()` for details. :param bool x_axis: Should we plot two horizontal axes showing time underneath the trees. Note that this corresponds to the y-axis in a conventional SVG plot, as tanglegrams are rotated 90°. :param str x_label: X axis label (equivalent of ``y_label`` in ``draw_svg()``). :param tuple[Union[list, dict]] x_ticks: Location of the tick marks on the two time axes. This is a tuple of two values, each of which is equivalent the the ``y_ticks`` value in ``draw_svg()``. :param bool x_gridlines: Whether to plot vertical lines behind the tree at each y tickmark (equivalent of ``y_gridlines`` in ``draw_svg()``). :param bool y_axis: Equivalent of `x_axis` parameter in ``draw_svg()``. Probably not what you want. :param bool y_label: Equivalent of `y_label` parameter in ``draw_svg()`` Probably not what you want. :param bool y_ticks: Dummy option: has no effect :param bool y_gridlines: Dummy option: has no effect :param bool return_node_maps: Instead of just returning an ``SvgTree``, return a tuple of ``(SvgTree, left_node_map, right_node_map)``. :param \**kwargs: Additional keyword arguments to pass to :meth:`Tree.draw_svg()`, such as `time_scale` ("log_time", "rank", etc) :returns: A tuple of an SvgTree object (that can be plotted by calling obj.draw()) and a left and a right node mapping. """ def node_positions(svgtree): x = svgtree.node_x_coord node_times = svgtree.ts.nodes_time if svgtree.time_scale == "rank": within_tree_node_times = np.unique(node_times[svgtree.tree.preorder()]) node_times = np.searchsorted(within_tree_node_times, node_times) y = svgtree.timescaling.transform(node_times) return x, y def make_reverse_map(node_map): reverse_map = np.zeros_like(node_map) kept = node_map != tskit.NULL reverse_map[node_map[kept]] = np.arange(len(node_map))[kept] return reverse_map def reorder_tree_nodes(tree, node_order): # given a node order and a tree, make a new tree and return that and the order node_map = np.arange(tree.tree_sequence.num_nodes) node_map[np.sort(node_order)] = node_order ts = tree.tree_sequence.subset(node_map, False, False, False) return node_map, ts.at_index(tree.index) def get_valid_leaf_order(tree, node_order): # take a node ordering and return a leaf ordering on the reordered-node tree if len(np.unique(node_order)) != len(node_order): raise ValueError("Order must contain unique integers") node_map, tree = reorder_tree_nodes(tree, node_order) leaves = np.array([u for u in tree.nodes(order="minlex_postorder") if tree.is_leaf(u)]) return node_map[leaves] if y_ticks is not None: raise ValueError("Invalid option") if y_gridlines is not None: raise ValueError("Invalid option") if order is None: order = (None, None) if tree_indexes is None: tree_indexes = ( 1 if ts.first().num_edges == 0 else 0, -2 if ts.last().num_edges == 0 else -1, ) lft = ts.at_index(tree_indexes[0]) rgt = ts.at_index(tree_indexes[1]) if titles is None: titles = (f"Tree {lft.index}", f"Tree {rgt.index}") if separation is None: extra_sep = 0 else: extra_sep = separation - 64 if size is None: # Note that the width is twice the default tree height, as these are rotated size = (200 * 2, 200) w = size[0] / 2 # width (tree height) of one of the plotted trees, after 90° rotation height = size[1] style = ( ".lft_tree > g.tree, .lft_tree > g.tangle_lines {transform: translate(0, " + str(height) + "px) rotate(-90deg);}" + ".lft_tree > g.tree .node > .lab {text-anchor: start; transform: rotate(90deg) translate(4px);}" ".lft_tree > .title {transform: translate(" + str(w/2) + "px);}" ".rgt_tree > g.tree {transform: translate(" + str(w) + "px, 0) rotate(90deg);}" ".rgt_tree > g.tree .node > .lab {text-anchor: end; transform: rotate(-90deg) translate(-4px);}" ".rgt_tree > .title {transform: translate(" + str(w/2) + "px);}" ".lft_tree .axes .y-axis .title text {transform: translate(11px) rotate(90deg);}" ".rgt_tree .axes .y-axis .title text {transform: translate(-11px) rotate(-90deg);}" ".rgt_tree .axes .y-axis .ticks .lab {text-anchor: end; transform: rotate(180deg);}" ) + (style or "") # For tree 1 we need to reverse the plotting order of leaves, so the leftmost # tip appears at the top when the tree is rotated 90° anticlockwise. We do this # by reordering using `subset()`, so minlex order will reverse the current order # This also means we need to re-adjust the node labels, as the sample IDs will change if node_labels is None: node_labels = {u: str(u) for u in np.arange(ts.num_nodes)} if order[0] is None: leaves = np.array([u for u in lft.nodes(order="minlex_postorder") if lft.is_leaf(u)])[::-1] else: leaves = get_valid_leaf_order(lft, order[0][::-1]) lft_node_map, lft = reorder_tree_nodes(lft, leaves) lft_rev_map = make_reverse_map(lft_node_map) # Have to change the node labels, because even provided ones will be targeting the wrong IDs lft_node_labels = {u: node_labels[v] for u, v in enumerate(lft_node_map) if v in node_labels} if order[1] is None: # We do not reorder the RH tree, so the node IDs should stay as-is # TODO - we could check the leaf IDs match here rgt_node_labels = node_labels rgt_node_map = rgt_rev_map = np.arange(ts.num_nodes) else: rleaves = get_valid_leaf_order(rgt, order[1]) rgt_node_map, rgt = reorder_tree_nodes(rgt, rleaves) if set(rleaves) != set(leaves): raise ValueError("Leaf IDs in the two trees are not the same") rgt_node_labels = {u: node_labels[v] for u, v in enumerate(rgt_node_map) if v in node_labels} rgt_rev_map = make_reverse_map(rgt_node_map) kwargs["size"] = (height, w) kwargs["order"] = "minlex_postorder" kwargs["y_label"] = x_label # Swapped because of 90° rotation kwargs["y_gridlines"] = x_gridlines # Swapped because of 90° rotation kwargs["x_axis"] = y_axis # Swapped because of 90° rotation kwargs["x_label"] = y_label # Swapped because of 90° rotation # We'll embed the right tree and the tangle lines within the left tree later, via the preamble svgtree_lft = tskit.drawing.SvgTree( lft, title=titles[0], canvas_size=(w * 2 + extra_sep, height), node_labels=lft_node_labels, root_svg_attributes={'class': 'lft_tree'}, y_axis=x_axis, y_ticks = None if x_ticks is None else x_ticks[0], **kwargs, ) svgtree_rgt = tskit.drawing.SvgTree( rgt, title=titles[1], canvas_size=(w, height), node_labels=rgt_node_labels, root_svg_attributes={'class': 'rgt_tree', 'x': w + extra_sep}, y_axis='right' if x_axis else x_axis, # Swapped because of 90° rotation y_ticks = None if x_ticks is None else x_ticks[1], style=style, **kwargs, ) x_lft, y_lft = node_positions(svgtree_lft) # Here we just need any list of leaves (order doesn't matter, as long as it's consistent) tip_h_lft = [x_lft[u] for u in lft_rev_map[leaves]] tip_w_lft = y_lft[lft_rev_map[leaves]] tip_w_lft = np.full_like(tip_w_lft, w - 10) if line_gap is None else (tip_w_lft + line_gap) x_rgt, y_rgt = node_positions(svgtree_rgt) tip_h_rgt = [x_rgt[u] for u in rgt_rev_map[leaves]] tip_w_rgt = y_rgt[rgt_rev_map[leaves]] tip_w_rgt = np.full_like(tip_w_rgt, w - 10) if line_gap is None else (tip_w_rgt + line_gap) lines = [ f'' for x1, y1, x2, y2 in zip(tip_h_lft, tip_w_lft, tip_h_rgt, tip_w_rgt) ] svgtree_lft.preamble = ( # Add the RH tree plus lines as the preamble '' + ''.join(lines) + '' + svgtree_rgt.draw() ) return (svgtree_lft, lft_rev_map, rgt_rev_map) if return_node_maps else svgtree_lft # Now run the function on an edited tree sequence (make nodes 1, 4, 5 a bit older, to demo) tables = ts_mutated.dump_tables() tables.mutations.time = np.full_like(tables.mutations.time, tskit.UNKNOWN_TIME) tables.nodes[1] = tables.nodes[1].replace(time=1000) tables.nodes[4] = tables.nodes[4].replace(time=500) tables.nodes[5] = tables.nodes[5].replace(time=2000) non_ultrametric_ts = tables.tree_sequence() tanglegram(non_ultrametric_ts, node_labels={u:u for u in ts_mutated.samples()}, mutation_labels={}, line_gap=15).draw() ``` Here's a larger tanglegram example (Note that if you want to compare two trees in their own separate tree sequences, you can concatenate them together using {meth}`.TreeSequence.concatenate`). ```{code-cell} ipython3 ts = msprime.sim_ancestry( samples={0: 20, 1:5}, sequence_length=1e5, demography=msprime.Demography.island_model([1000, 1000], migration_rate=0.01), recombination_rate=1e-8, random_seed=123, ) css = ".node.sample.p0 .sym {fill: green} .node.sample.p1 .sym {fill: darkred}" tanglegram( ts, size=(600, 500), line_gap=3, node_labels={}, x_axis=True, time_scale="log_time", x_ticks=[[0, 1, 10, 100, 1000]] * 2, # same x_ticks on lft and rgt style=css).draw() ``` With a bit more effort, you can even identify identical clades in the same trees (although styling can be complex because of having to remap node IDs): ```{code-cell} ipython3 :"tags": ["hide-input"] from hashlib import blake2b # Use the blake2b hash to identify nodes with the same samples ts = ts.simplify(list(range(20))) tg, lft_map, rgt_map = tanglegram( ts, size=(600, 300), line_gap=3, node_labels={}, x_axis=True, time_scale="log_time", x_ticks=[[0, 1, 10, 100, 1000]] * 2, # same x_ticks on lft and rgt return_node_maps=True) ltree = ts.first() rtree = ts.last() def hashfunc(tree, u): return blake2b(" ".join(str(v) for v in sorted(tree.samples(u))).encode(), digest_size=20).digest() hashdict1 = {hashfunc(ltree, u): u for u in ltree.nodes() if not ltree.is_sample(u)} hashdict2 = {hashfunc(rtree, u): u for u in rtree.nodes() if not rtree.is_sample(u)} shared_hashes = hashdict1.keys() & hashdict2.keys() css = "" css += ",".join([f".lft_tree > .tree .n{lft_map[hashdict1[h]]} > .sym" for h in shared_hashes]) + f"{{r: 3px; fill: magenta; stroke: black;}}" css += ",".join([f".rgt_tree > .tree .n{rgt_map[hashdict2[h]]} > .sym" for h in shared_hashes]) + f"{{r: 3px; fill: magenta; stroke: black;}}" # Add the extra styles into the preamble legend = ( '' '' '= shared clades' ) tg.preamble = f"" + legend + tg.preamble tg.draw() ``` (sec_tskit_viz_SVG_examples_animation)= #### Animation The classes attached to the SVG also allow elements to be animated. Here's a [d3.js](https://d3js.org)-based animation of sucessive subtree-prune-and-regraft (SPR) operations, using the {ref}`ARG representation ` of a tree sequence to allow identification of pruned edges. ```{code-cell} ipython3 :"tags": ["hide-input"] from IPython.display import HTML import msprime def make_full_arg_for_spr_animation(): # created with record_full_arg needed to track recombination nodes (branch positions) # random_seed chosen to produce a ts whose leaves are plotted in the same order return msprime.sim_ancestry( 5, ploidy=1, sequence_length=10000, recombination_rate=0.00005, random_seed=6787, model="smc_prime", record_full_arg=True) css_string = ( "#anim_svg {background-color: white} " "#anim_svg .node:not(.sample) > .lab, #anim_svg .node:not(.sample) > .sym {display: none}" ) html_string = r""" %s """ ts = make_full_arg_for_spr_animation() HTML(html_string % ( ts.draw_svg(root_svg_attributes={"id": "anim_svg"}, style=css_string), ts.num_trees, )) ``` (sec_tskit_viz_other)= ## Other visualizations As well as visualizing a tree sequence as, well, a sequence of local trees, or by plotting {ref}`statistical summaries `, other visualizations are possible, some of which are outlined below. (sec_tskit_viz_other_graph)= ### Graph representations A tree sequence can be treated as a specific form of (directed) [graph](https://en.wikipedia.org/wiki/Graph_(discrete_mathematics)) consisting of nodes connected by edges. Standard graph visualization software, such as [graphviz](https://graphviz.org) can therefore be used to depict tree sequence topologies. Alternatively, the [tskit_arg_visualizer](https://github.com/kitchensjn/tskit_arg_visualizer) project will draw a interactive `tskit` graph directly, in which nodes can be dragged horizontally, and embedded (local) trees highlighted by hovering over the "genome bar" underneath the graph. Below is an example, showing an `msprime` "full ARG" tree sequence. In this case, nodes have only 2 children in the graph, so `edge_type="ortho"` can be used to draw a traditional "[Ancestral Recombination Graph](sec_args)" style plot: ```{code-cell} ipython3 :"tags": ["hide-input"] %%javascript require.config({paths: {d3: 'https://d3js.org/d3.v7.min'}}); require(["d3"], function(d3) {window.d3 = d3;}); ``` ```{code-cell} ipython3 import msprime import tskit_arg_visualizer tip_order = [3, 0, 1, 2, 6, 7, 4, 5] # Found by trial and error for this seed full_arg_ts = msprime.sim_ancestry( 4, sequence_length=1000, recombination_rate=0.001, record_full_arg=True, random_seed=3) d3arg = tskit_arg_visualizer.D3ARG.from_ts(ts=full_arg_ts) d3arg.draw(width=500, height=500, edge_type="ortho", sample_order=tip_order); ``` For tree sequences that may not be "full ARGs", the default `edge_type="line"` is preferable. The example below also uses uses the `variable_edge_width` option to emphasise which edges have wider spans, and `show_mutations` to display mutations on edges in an interactive style (hovering over the mutation or the genome bar will reveal their locations along the chromosome): ```{code-cell} ipython3 import msprime import tskit_arg_visualizer ts = msprime.sim_ancestry( 4, sequence_length=1000, recombination_rate=0.001, record_full_arg=True, random_seed=3) # simplify into a standard (non "full ARG") tree sequence ts = ts.simplify() ts = msprime.sim_mutations(ts, rate=0.001, random_seed=5) d3arg = tskit_arg_visualizer.D3ARG.from_ts(ts=ts) tip_order = [3, 0, 1, 2, 6, 7, 4, 5] drawinfo = d3arg.draw( width=500, height=500, edge_type="line", variable_edge_width=True, sample_order=tip_order, show_mutations=True, label_mutations=True, ) print(f"Extra styling is possible by targetting CSS at this unique id: #{drawinfo.uid}") ``` For more general graph plots, it can be helpful convert the tree sequence to a [networkx](https://networkx.org) graph first, as described in the {ref}`sec_args_other_analysis` section of the {ref}`sec_args` tutorial. This provides interfaces to graph plotting software such as [graphviz](https://graphviz.org), which provides the `dot` layout engine for directed graphs: ```{code-cell} ipython3 :"tags": ["hide-input"] ## Networkx conversion code taken from the ARG tutorial import networkx as nx import pandas as pd import tskit def to_networkx_graph(ts, interval_lists=False): """ Make an nx graph from a tree sequence. If `intervals_lists` is True, then each graph edge will have an ``intervals`` attribute containing a *list* of tskit.Intervals per parent/child combination. Otherwise each graph edge will correspond to a tskit edge, with a ``left`` and ``right`` attribute. """ D = dict(source=ts.edges_parent, target=ts.edges_child, left=ts.edges_left, right=ts.edges_right) G = nx.from_pandas_edgelist(pd.DataFrame(D), edge_attr=True, create_using=nx.MultiDiGraph) if interval_lists: GG = nx.DiGraph() # Mave a new graph with one edge that can contai for parent, children in G.adjacency(): for child, edict in children.items(): ilist = [tskit.Interval(v['left'], v['right']) for v in edict.values()] GG.add_edge(parent, child, intervals=ilist) G = GG nx.set_node_attributes(G, {n.id: {'flags':n.flags, 'time': n.time} for n in ts.nodes()}) return G ``` ```{code-cell} ipython3 import networkx as nx from IPython.display import SVG def graphviz_svg(networkx_graph): AG = nx.drawing.nx_agraph.to_agraph(networkx_graph) # Convert to graphviz "agraph" nodes_at_time_0 = [k for k, v in networkx_graph.nodes(data=True) if v['time'] == 0] AG.add_subgraph(nodes_at_time_0, rank='same') # put time=0 at same rank return AG.draw(prog="dot", format="svg") G = to_networkx_graph(ts, interval_lists=True) # Function from the ARG tutorial print("Converted `ts` to a networkx graph named `G`") print("Plotting using graphviz...") SVG(graphviz_svg(G)) ``` Alternatively, you can read the `graphviz` positions back into `networkx` and use the `networkx` drawing functionality, which relies upon the [matplotlib](https://matplotlib.org) library. This allows modification of node colours and symbols, labels, rotations, annotations, etc., as shown below: ```{code-cell} ipython3 :"tags": ["hide-input"] from matplotlib import pyplot as plt import string def get_graphviz_positions(networkx_graph): AG = nx.drawing.nx_agraph.to_agraph(networkx_graph) # Convert to graphviz "agraph" nodes_at_time_0 = [k for k, v in networkx_graph.nodes(data=True) if v['time'] == 0] AG.add_subgraph(nodes_at_time_0, rank='same') # put time=0 at same rank AG.layout(prog="dot") # create the layout, storing positions in the "pos" attribute return {n: [float(x) for x in AG.get_node(n).attr["pos"].split(",")] for n in G.nodes()} pos=get_graphviz_positions(G) edge_labels = { edge[0:2]: "\n".join([f"[{int(i.left)},{int(i.right)})" for i in edge[2]["intervals"]]) for edge in G.edges(data=True) } samples = set(ts.samples()) nonsamples = set(range(ts.num_nodes)) - samples plt.figure(figsize=(10, 4)) # Sample nodes as dark green squares with white text nx.draw(G, pos, node_color="#007700", font_size=9, node_size=250, node_shape="s", nodelist=samples, edgelist=[]) nx.draw_networkx_labels(G, pos, font_size=9, labels={u: u for u in samples}, font_color="white") # Others as blue circles with alphabetic labels nx.draw(G, pos, node_color="#22CCFF", font_size=9, edgecolors="black", nodelist=nonsamples, edgelist=[]) nx.draw_networkx_labels(G, pos, font_size=9, labels={u: string.ascii_lowercase[u] for u in nonsamples}) nx.draw_networkx_edges(G, pos, edge_color="lightgrey", arrows=False, width=2); nx.draw_networkx_edge_labels(G, pos, edge_labels, font_size=7, rotate=False); ``` Note, however, that finding node and edge layout positions that avoid too much overlap can be tricky, even for the graphviz layout engine, and there is no easy functionality to place nodes at specific vertical (time) positions. (sec_tskit_viz_other_demographic)= ### Demographic processes If you are generating a tree sequence via a {ref}`Demes ` model, then you can visualize a schematic of the demography itself (rather than the resulting tree sequence) using the [DemesDraw](https://grahamgower.github.io/demesdraw/) software. For example, here's the plotting code to generate the {ref}`demography plot` from the "{ref}`sec_what_is`" tutorial: ```{code-cell} ipython3 :"tags": ["hide-input"] import matplotlib_inline import matplotlib.pyplot as plt %matplotlib inline matplotlib_inline.backend_inline.set_matplotlib_formats('svg') import demes import demesdraw def size_max(graph): return max( max(epoch.start_size, epoch.end_size) for deme in graph.demes for epoch in deme.epochs ) # See https://popsim-consortium.github.io/demes-docs/ for the yml spec for the file below graph = demes.load("data/whatis_example.yml") w = 1.5 * size_max(graph) positions = dict(Ancestral_population=0, P=-w, Q=w) ax = demesdraw.tubes(graph, positions=positions, seed=1) plt.show(ax.figure) ``` ### Geography :::{todo} How to get lat/long information out of a tree sequence and plot ancestors (or a tree) on a geographical landscape. :::