what individuals to include in tree building

[percyfal]

Hi,

I don't know if this is the right forum to post as this is not so much an issue as a question on best practice usage. Is there a mailing list or similar for usage questions?

Anyway, I'm trying out tsinfer on a data set that includes some 20 populations (species / subspecies) and three outgroups. To begin with we are interested in analyzing a population subset (~12 populations). My question is if one should build trees using only the subset or just as well use all of the populations to then subset the trees for the analyses. Also, we have derived the ancestral states based on the outgroups. Should they be excluded from tree building? Is there a best practice way forward?

Cheers,

Per

[hyanwong]

Hi Per - good to hear you are trying out tsinfer. @benjeffery has kindly turned on "GitHub discussions" and moved this question there, so hopefully others will also use this as a forum for tsinfer questions.

In answer to your question, it's early days for working out "best practice", but normally we would advise putting in as much data as is available. I think this would include the outgrips too, as long as you are reasonably sure that the ancestral state for the entire clade (including the outgroups) is correct.

It would be sensible to exclude groups that are poorly sequenced, but I suspect that it will be individual sites that will be poorly known (or perhaps where the ancestral state is uncertain). These can be excluded from the entire inference process by setting inference=False during add_site(), or excluded on a per-sample basis by setting some of the data to missing.

[jeromekelleher]

No, I'm not sure what's happened there @hyanwong - I wouldn't have expected that step to consume that much memory (This is the function it's calling into).

@percyfal, would you mind rerunning this with the debug flags turned on? (Something like -vvv should turn it all the way up)

[percyfal]

@jeromekelleher I make the tsinfer infer call from within the script (based on https://tsinfer.readthedocs.io/en/latest/tutorial.html#reading-a-vcf) but presumably this is what you meant (CLI)?
tsinfer infer -vvv data.samples --num-threads 150 --progress

2021-04-07 19:47:10,599 [1659] INFO     tsinfer.formats: Loaded SampleData(num_samples=202, num_sites=2197683)
2021-04-07 19:47:10,600 [1659] INFO     tsinfer.formats: Loaded SampleData(num_samples=202, num_sites=2197683)
2021-04-07 19:47:13,895 [1659] DEBUG    tsinfer.inference: Using C AncestorBuilder implementation
2021-04-07 19:47:13,895 [1659] INFO     tsinfer.inference: Starting addition of 2197683 sites
ga-add   (1/6)100%|████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 2.20M/2.20M [01:39, 22.1kit/s]
2021-04-07 19:48:57,625 [1659] INFO     tsinfer.inference: Finished adding sites
2021-04-07 19:48:57,910 [1659] DEBUG    tsinfer.formats: Syncronously flushing buffer
2021-04-07 19:48:57,911 [1659] DEBUG    tsinfer.formats: Flushing buffer 0: start=0 n=0
2021-04-07 19:48:57,911 [1659] DEBUG    tsinfer.formats: Buffer 0 flush done
2021-04-07 19:48:57,913 [1659] CRITICAL root: Traceback (most recent call last):
  File "/mnt/nbis-scratch/perun/valkyrian/bin/tsinfer", line 8, in <module>
    sys.exit(main())
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/__main__.py", line 5, in main
    cli.tsinfer_main()
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/cli.py", line 499, in tsinfer_main
    args.runner(args)
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/cli.py", line 143, in run_infer
    ts = tsinfer.infer(
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/inference.py", line 274, in infer
    ancestor_data = generate_ancestors(
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/inference.py", line 366, in generate_ancestors
    generator.run()
  File "/mnt/nbis-scratch/perun/valkyrian/lib/python3.9/site-packages/tsinfer/inference.py", line 939, in run
    self.descriptors = self.ancestor_builder.ancestor_descriptors()
MemoryError

EDIT: I monitored the actual memory usage and it barely registered 0.1% usage so it doesn't seem to be a matter of exhausting memory resources

[jeromekelleher]

This is very strange, thanks for your help in debugging it @percyfal. Do you think you could make the data.samples file available to us somehow? I think there must a bug of some sort here.

[percyfal]

@jeromekelleher Ok I've figured out where the code breaks, but I cannot identify what's causing it. Ideally I would like to extract the region but I would need a pointer from you as how to find it.

In the function ancestor_builder_finalise the variable focal_sites is set to NULL at some point in the iteration through the avl. I printed some information when it failed:

number of samples: 202
focal_sites is NULL
pattern_map->num_sites: 374662
pattern_map->num_samples: 202
tsk_id_t size: 4
pattern_map->num_sites * sizeof(tsk_id_t)): 1498648

1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 0|0 0|0 0|0 0|0 0|0 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1
 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 1|1 2|0 

I presume num_sites is the size of the block that is used to build a partial genealogy? Printing it for all regions typically returned single-digit values, although occasionally would jump to a couple of thousand, but never anything this big. At first I thought maybe it was caused by a dangling pointer, but having traversed the sites I don't think so.

What strikes me as odd though is that the last genotype is coded 2|0 - am I wrong to assume that we only should see 0 (ancestral) / 1 (derived)? I double-checked my samples generation script, but I had correctly made sure only to include biallelic snps coded 0 or 1, so this conversion happens later. Also, this only seems to appear for the failing chromosomes (I checked for the existence of genoypes equalling 2 in a few cases).

I've tried as best I can to narrow down the region where this happens to provide a smaller reproducible example for an issue but I don't directly see how to get the actual chromosome positions from the pattern_map struct.

[hyanwong]

A fix to the memory issue has just been pushed to main.