Computing GNN for each sample with each sample

[janaobsteter]

Hi, all! I have a question about computing the GNN and would appreciate some help. I am looking at your figure from the tsinfer paper. How does one compute the GNN of each sample with each sample (is this is what you did in this figure)? Thanks!

[percyfal]

Hi @janaobsteter!

I have been working on code to compute gnns and have successfully
managed to generate results (and plots) based on accompanying code to
the tsinfer paper (https://github.com/mcveanlab/treeseq-inference, in
particular the code in
src/plot.py
and
src/analyse_human_data.py -
maybe you have seen these already?).

Anyhow, I distilled what I believe are the relevant parts and I have
one script to compute the gnn. The script is part of a snakemake
workflow and cannot currently be run standalone, but would be easy to
modify if need be. Anyhow, the script is here
(https://github.com/percyfal/wg-genealogy-smk/blob/main/workflow/scripts/tsinfer-gnn.py),
in case it helps.

I'd be happy to discuss further as I'm currently implementing
functions to compute and plot chromosome gnns (fig d above). Here, I'm not sure I
completely understand the rationale for the underlying function. @hyanwong would this merit a separate discussion or would it fit here?

Cheers,

Per

[hyanwong]

Thanks @percyfal. I think "local" (or, preferably, "windowed") GNNs are useful, but should probably be under a different discussion thread. There is a a PR to get this in to the standard stats framework, which would mean you wouldn't have to roll your own. See #665 and #1237. It would be great to close these issues, if you have the time to get involved @percyfal!

[percyfal]

Thanks @hyanwong I'll have a look and see if I can be of any help. Meanwhile I'll open a new discussion thread related to the windowed GNNs. If I'm going to chip in I'd better understand the rationale of the algorithm :)

[hyanwong]

Sure, or continue the thread in issue #665. Either is fine.

[hyanwong]

Hi @janaobsteter . The plotting code is at https://github.com/mcveanlab/treeseq-inference/blob/0cbbb062c96ad4433d8b4d0f120f93ac2d985345/src/plot.py#L1219 and the code to generate the data to plot is at https://github.com/mcveanlab/treeseq-inference/blob/0cbbb062c96ad4433d8b4d0f120f93ac2d985345/human-data/tsutil.py#L333. We actually didn't calculate the GNN of each sample with each, we binned them into populations.

However, it would be helpful to have a recipe to do this, so I'll post a snippet below

[hyanwong]

Here's some example code for figures like (a) and (b) above. It should work in a Jupyter notebook, if you want to try it.

import demes
import msprime
import pandas as pd
import scipy
import seaborn as sns

demes_yml = """
description: Stepping stone model illustrating default value propagation.
time_units: "years"
generation_time: 25
defaults:
  epoch:
    start_size: 1000
  migration:
    rate: 0.1
demes:
  - name: deme1
  - name: deme2
  - name: deme3
migrations:
  - demes: [deme1, deme2]
  - demes: [deme2, deme3]
"""

# Generate the tree sequence
graph = demes.loads(demes_yml)
demography = msprime.Demography.from_demes(graph)
ts = msprime.sim_ancestry(samples={"deme1": 10, "deme2": 10, "deme3": 8}, demography=demography, random_seed=12)

# Generate the GNN data and put in a pandas dataframe
num_threads = 1  # could increase this for better peerformance on multicore machines
sample_sets = [ts.samples(pop) for pop in range(ts.num_populations)]
A = ts.genealogical_nearest_neighbours(
        ts.samples(), sample_sets, num_threads=1)
cols = {pop.metadata["name"]: A[:, pop.id] for pop in ts.populations()}
cols["population"] = [ts.population(ts.node(u).population).metadata["name"] for u in ts.samples()]
cols["individual"] = [ts.node(u).individual for u in ts.samples()]
df = pd.DataFrame(cols)

# Plot the data
figsize = (10, 10)
dfg = df.groupby("population").mean()
# Zscore normalise
for col in list(dfg):
    dfg[col] = scipy.stats.zscore(dfg[col])
row_linkage = scipy.cluster.hierarchy.linkage(dfg, method="average")
order = scipy.cluster.hierarchy.leaves_list(row_linkage)
x_pop = dfg.index.values[order]
cg = sns.clustermap(
    dfg[x_pop], row_linkage=row_linkage, col_cluster=False,
    figsize=figsize, rasterized=True)

[janaobsteter]

@hyanwong and @percyfal , thanks for the help! I think I can make it work now

[janaobsteter]

@yan Wong, I have a couple more question. I managed to draw the gnn heatmaps for my honeybees - on individual, lineage and country level (different fragmentation). However, it is not as nice as yours (attached the ind plot) and I would like to change the colour scale to not be linear. I've tried a couple of solutions but non of them worked. Do you have any suggestion for this?
Also, another question. This figure is only from chromosome 1. How do you do this for multiple chromosomes? Do you somehow "merge" the trees or you you average the GNN across the chromosomes? Thanks!!!

[hyanwong]

I would like to change the colour scale to not be linear

I'm not sure how you are getting values of 0...25 for your z scores? Normally they should be 0..1, or if you normalise them (as in the code), they will range from about -2 to 2. Perhaps plot them out and work out a transformation, or truncate them somehow? Since they are 0..1 I presume an arcsine transformation would be the most obvious, but YMMV.

This figure is only from chromosome 1. How do you do this for multiple chromosomes?

We just take the final dataframes and average them (weighting by the length of the chromosome). There's no reason to stick them all into one humungous TS unless you really want to. You might want to snip out the centromeres or variant-poor regions, though (see https://tsdate.readthedocs.io/en/latest/python-api.html#tsdate.preprocess_ts)