add relatedness calculation based on tree-seq table data

[bhaller]

This is sort of a question for @petrelharp. (I apologize if I'm absent-minded and have asked it before!) Is it possible, and straightforward, to calculate true relatedness (percent of IBD segments, I guess?) between individuals A and B based on the recorded information in the tree-sequence tables, without actually constructing the tree sequence (which is very expensive)? There is a lot of demand for this sort of thing (just now, again, in the Vienna workshop), because people want to be able to model outbreeding/inbreeding, behavior that varies with kinship, etc. It kinda seems like the information is there in the tables, and maybe simplification is not needed to do the calculation since you'd just work upwards from A and B to their common ancestor(s)? But the details are beyond me. :-> Thoughts?

[gregorgorjanc]

This preprint On ARGs, pedigrees, and genetic relatedness matrices On ARGs, pedigrees, and genetic relatedness matrices describes the method that can be run on tree sequence, but I am confused by your question asking "based on the recorded information in the tree-sequence tables, without actually constructing the tree sequence (which is very expensive)?". You do need a tree sequence to calculate branch relatedness.

[bhaller]

This preprint On ARGs, pedigrees, and genetic relatedness matrices On ARGs, pedigrees, and genetic relatedness matrices describes the method that can be run on tree sequence, but I am confused by your question asking "based on the recorded information in the tree-sequence tables, without actually constructing the tree sequence (which is very expensive)?". You do need a tree sequence to calculate branch relatedness.

Hi @gregorgorjanc! (Was there supposed to be a link in your comment?) Yes, what you say is precisely my question. With a tree sequence, computing relatedness is easy (and in fact tskit provides the necessary machinery). In SLiM, we don't have a tree sequence, however; we only have the information in the (node, site, edge...) tables. Making a tree sequences out of those tables is quite slow and heavyweight, not least because simplification has to be performed. My question is whether it might be possible to write an algorithm that would compute relatedness directly from the information in the tables, by walking upward through the relevant entries in the table, without actually constructing the tree sequence. If so, that would be very useful in the context of SLiM. But maybe it is not possible (without, in effect, doing all of the work that you would have to do to construct the tree sequence to answer the question).

[gregorgorjanc]

https://pubmed.ncbi.nlm.nih.gov/40093116/

I understand your question now. Tricky and beyond me atm.

[petrelharp]

The thing you describe - percent IBD - is I think most naturally calculated in a "bottom-up" type framework (that's how the IBD stuff works), and so, like simplify, doesn't really need the tree sequence (since you only need the trees for left-right). So one could possibly use the existing IBD code on the tables to do this. However, a few major issues with this:

1. I don't think it would actually work out of the box since I think it's written to work on a tree sequence (even if it could be easily modified to work on tables).
2. It would need tables to be sorted (like simplify).
3. What, exactly, would want to be calculated? "IBD" only makes sense relative to some reference time, since after all everyone is IBD everywhere if you look back far enough. Probably people would like to know percent IBD in the last two or three generations?
4. With N individuals there are N^2 measures of relatedness; certainly we don't want all of them. How would the result be returned? I guess the interface could be something like spatial operations.

Here's a thing you could do: suppose that everyone's parents and grandparents are in the tree sequence (as individuals). We could then compute percent IBD between all siblings, cousins, parent-offspring, niblings, etc in O(N) time because everyone has a pointer to their parents (at least the ones they got genome from) in the edge table. Now, suppose you remembered everyone, but with permanent=F and retainCoalescentOnly=T. Then I think maybe you'd have all the information you need? So, you'd have to do something like: construct a list that gives for each ancestor a list of their descendants and which bits of genome they inherited (do this by looping over alive individuals and using information in (grand)parentPedigreID and the edge table); then report back on all pairs sharing overlapping bits in the ancestors.

[vsudbrack]

I confess I struggle a bit to see the point of having the expected relatedness calculated for you in a genetically-explicit simulation (where the exact IBD can be calculated for any pair of individuals, or for any pair of haplosomes within an individual). I guess it also goes back to the point (3) of @petrelharp

But I just thought of mentioning that I think Jérôme Goudet can be very helpful on this brainstorming, he works on the interface of pedigrees and metrics of kinship (IBD, Fis, Froh, Jacquard coefficients). I’m afraid he is not here on GitHub (I don’t seem to be able to tag him), but I’ll drop him an email with this thread and maybe he either joins us here or sends you an email with his thoughts.

[bhaller]

Hmm. Well, people want to make simulations where individuals tend to either mate with, or avoid mating with, close relatives. SLiM's pedigree tracking only goes back to grandparents, which isn't quite far enough for practical purposes, I gather, and adding even one more generation of pedigree tracking information would take quite a lot of additional memory and runtime; it's expensive. I'm not sure why IBD doesn't cut it, exactly; that's a good question. But perhaps relatedness calculated by mutations is a pretty imperfect and stochastic approximation of the truth (which the tree-sequence tables could give you), and perhaps organisms often know "the truth" or something close to it? Humans, for example, generally know who their relatives are, and in most cultures avoid mating with them, even if their genetics would suggest that such a mating would happen to be safe; they don't know their own genetics, after all. I imagine some other animals do similar things. I dunno. Anyway, it's a thing people ask for, and it seems perhaps possible to provide it. If there are strong opinions that it's not meaningful/useful/correct, and that we shouldn't provide it, then that's useful feedback, but I'd like to better understand why.

@petrelharp Why would the tables need to be sorted? Can't we walk backwards to parents, grandparents, etc., at whatever index those entries are at? Where does the requirement for sorting come from? (It's painful to satisfy, so that why I'm asking.) Re: (3), I guess 3-5 generations probably? Re: (4), I don't understand; isn't the calculation of consanguinity, given tree-sequence information, well-defined conceptually? Where do the "N^2 measures of relatedness" and "spatial operations" aspects come from? And finally, re: "remembered everyone", is that really necessary (because it's again painful)? If a given ancestor was simplified away, don't we get the same answer, in terms of consanguinity, by just going up to the next ancestor up? I.e., wasn't that ancestor simplified away precisely because they don't matter, and the next ancestor up gives the same information that would have been provided by the ancestor that was simplified away? Sorry for the dumb questions, I just don't really understand the simplification algorithm very deeply so all of this is a bit opaque to me. :-O

[gregorgorjanc]

@bhaller tracking pedigree should not be too hard/memory hungry and you could calculate pedigree relatedness for the currently live individuals only in a quick fashion using Colleau's algorithm https://pubmed.ncbi.nlm.nih.gov/12270102 - all operations are based on sparse matrix. The output is still O(n^2) so this won't scale memory-wise for huge sims (many individuals)! This is all pedigree, hence expected relatedness since some base population, but that should suffice for avoiding mating close relatives.

[bhaller]

@bhaller tracking pedigree should not be too hard/memory hungry and you could calculate pedigree relatedness for the currently live individuals only in a quick fashion using Colleau's algorithm https://pubmed.ncbi.nlm.nih.gov/12270102 - all operations are based on sparse matrix. The output is still O(n^2) so this won't scale memory-wise for huge sims (many individuals)! This is all pedigree, hence expected relatedness since some base population, but that should suffice for avoiding mating close relatives.

Hi @gregorgorjanc. Keeping ids for two parents is two 64-bit ints per individual; keeping grandparents is another four 64-bit ints. Keeping great-grandparents would be another eight 64-bit ints. The impact, in terms of both memory footprint and the amount of time required to set up those values, is prohibitive. Keep in mind that people run SLiM models with as many as a billion individuals – it adds up! :->

[gregorgorjanc]

You have to store it as table/vectors

id mother father
1 0 0
2 0 0
3 1 2
4 1 2
5 3 4
etc

So, it recursive and just 2 extra info per individual, which would be less than what you do now (2 for parents and 4 grandparents) ;)

If you have time point info you can subset this table/vectors by time, say keeping only 2 generations as you do now.

Alternative to pedigree would be to grab DNA loci giving you X matrix of genotypes (ind x loci) and G=ZZ^T is a measure of genomic/realised relatedness (Z column centred X).

[bhaller]

You have to store it as table/vectors

id mother father 1 0 0 2 0 0 3 1 2 4 1 2 5 3 4 etc

So, it recursive and just 2 extra info per individual, which would be less than what you do now (2 for parents and 4 grandparents) ;)

Hmm, interesting suggestion. At present this information is kept inside the individual objects, and the parental individuals don't necessarily exist any more, so the previous data would be lost; that's why I store the grandparent info even though it's redundant in some sense (everybody's grandparent is somebody else's parent). But yes, all the information could be kept in a separate data structure instead as you suggest. But then (a) finding the information for a given individual then requires searching the data structure (looking up a row for a given individual ID would require a binary search each time, it seems to me, unless additional indexing information is kept), so doing relatedness calculations might be significantly slower than now; and (b) older entries would need to be cleared out eventually to avoid the data structure growing without bound, which would be somewhat complicated, especially in a nonWF model with age structure, probably requiring a bottom-to-top painting algorithm. This suggested design is worth considering, since it would solve this problem of the pedigree information only going back to grandparents; but it is not entirely simple and would add overhead in some ways. And of course pedigree-based relatedness is different from, and less precise in some sense than, relatedness based upon actual genetic inheritance along the chromosome, so this is somewhat orthogonal to whether we can/should be able to measure relatedness based on the tree sequence data.

Alternative to pedigree would be to grab DNA loci giving you X matrix of genotypes (ind x loci) and G=ZZ^T is a measure of genomic/realised relatedness (Z column centred X).

Yeah, I think that's more or less @vsudbrack and I were discussing above, right? Using the mutations as an approximation for the truth expressed by the pedigree/trees. But it is only an approximation; it depends on mutational density; and it's much slower to calculate, probably.

[petrelharp]

I confess I struggle a bit to see the point of having the expected relatedness calculated for you in a genetically-explicit simulation

There's more to relatedness than genetic relatedness! (And, if you're doing a simulation of a single chromosome (or less), the expected relatedness might stand in for the rest of the genome.)

[petrelharp]

@petrelharp Why would the tables need to be sorted? Can't we walk backwards to parents, grandparents, etc., at whatever index those entries are at?

Yes - sorry, that was also my point. That if you wanted to do this following existing tree sequence tools, to go back to an arbitrary time, then you'd want them to be sorted. But the "just two generations" proposal doesn't require it since you just follow pointers.

Re: (3), I guess 3-5 generations probably?

Once you go beyond 1 or 2 generations, the precise number seems very arbitrary but also matters a great deal.

Re: (4), I don't understand; isn't the calculation of consanguinity, given tree-sequence information, well-defined conceptually? Where do the "N^2 measures of relatedness" and "spatial operations" aspects come from?

The issue is that N^2 is a large number; you don't want to have to calculate all N^2 numbers. So the obvious thing to do is to only (calculate and) return relatedness values for those pairs of individuals that are closely related, so there is an API question of how to actually do that. The spatial interaction code in SLiM is an example of such an API (since there are N^2 pairwise distances, but we don't want all of them...).

And finally, re: "remembered everyone", is that really necessary (because it's again painful)? If a given ancestor was simplified away, don't we get the same answer, in terms of consanguinity, by just going up to the next ancestor up?

Yes - recall that if you do permanent=F then individuals can be simplified away; that's the point. I suspect that doing permanent=F, retainCoalescentOnly=T then you'll have exactly those individuals you need for the calculation.

But: this seems like a giant pain, honestly.

[bhaller]

@petrelharp Why would the tables need to be sorted? Can't we walk backwards to parents, grandparents, etc., at whatever index those entries are at?

Yes - sorry, that was also my point. That if you wanted to do this following existing tree sequence tools, to go back to an arbitrary time, then you'd want them to be sorted. But the "just two generations" proposal doesn't require it since you just follow pointers.

OK, I think we're on the same page. For 3, 4, 5 generations you could still just follow the pointers, yes? We'd only need to sort if we wanted to generate a tree sequence object. Yes?

Re: (3), I guess 3-5 generations probably?

Once you go beyond 1 or 2 generations, the precise number seems very arbitrary but also matters a great deal.

Huh. I'm not sure I see why that is true. So where does that leave us, then? Is the very idea of genetic relatedness nonsensical and arbitrary, then, because the answer you get depends strongly on how far back you look, which is arbitrary? That doesn't seem quite right to me. I mean, kin selection etc. is a real phenomenon, no? Inbreeding vs. outbreeding is real, no? So I'm not really getting this objection.

Re: (4), I don't understand; isn't the calculation of consanguinity, given tree-sequence information, well-defined conceptually? Where do the "N^2 measures of relatedness" and "spatial operations" aspects come from?

The issue is that N^2 is a large number; you don't want to have to calculate all N^2 numbers. So the obvious thing to do is to only (calculate and) return relatedness values for those pairs of individuals that are closely related, so there is an API question of how to actually do that. The spatial interaction code in SLiM is an example of such an API (since there are N^2 pairwise distances, but we don't want all of them...).

Aha, I see. Hmm.

And finally, re: "remembered everyone", is that really necessary (because it's again painful)? If a given ancestor was simplified away, don't we get the same answer, in terms of consanguinity, by just going up to the next ancestor up?

Yes - recall that if you do permanent=F then individuals can be simplified away; that's the point. I suspect that doing permanent=F, retainCoalescentOnly=T then you'll have exactly those individuals you need for the calculation.

I see, hmm. But we need to remember everyone? Just allowing simplification to happen as it normally happens is not sufficient? That surprises me, since I thought simplification only removed things that are not informative anyway.

But: this seems like a giant pain, honestly.

Yes, agreed. At the same time, it's something people really want to be able to do – better relatedness calculations, during a simulation to affect the simulation's dynamics. I think perhaps I like @gregorgorjanc's idea of how to make the memory/runtime requirements of pedigree tracking less onerous, and then extend pedigree tracking to cover more generations back in time, and provide better pedigree-based relatedness. That won't make everyone happy, but I think it will make many people happy. But I'm still not sure how to handle the other issues here, such as how far back in time we want the relatedness calculations to go (and your statement that the result you get depends heavily upon that choice). I was thinking that if I just made it possible to go, say, 4 generations back instead of 2, that would be unambiguously good and would make everybody happy; but maybe not. Maybe I need to allow the user to specify the number of generations back that the pedigree tracking data should be kept, and use that many generations for the relatedness calculation?

Anyway, it's a can of worms, but I do think it ought to be handled at some point since the demand for it is so high; it's probably the number one feature request I get during workshops now. (Number two being simulating indels, speaking of cans of worms. :->)