count of ibd segments between all sampled pairs

[hancockzb]

Hi all,

Quick and simple question that me, a python layman, am struggling with: I'd like to count the number of genomic segments in which individual i and individual j in range(number of individuals) have the same MRCA. Goal here is to compare number of ibd segments between modern-modern individuals sampled across space and modern-ancient individuals across space and time.

Thanks so much for your time and all the hard work you all have put into this package!

Best,
Zach

[jeromekelleher]

Hi Zach! 👋

This is actually not a simple question at all - computing IBD segments efficiently is not easy, and we don't have any public APIs for doing this at the moment. Here's a simple brute-force approach though, that you might be able to adapt for your purposes:

import itertools
import dataclasses
import msprime


@dataclasses.dataclass
class IbdSegment:
    pair: (int, int)
    mrca: int
    left: float = 0
    right: float = -1


ts = msprime.simulate(5, recombination_rate=0.01, length=40, random_seed=42)
print(ts.draw_text())

pairs = list(itertools.combinations(ts.samples(), 2))
tree = ts.first()
mrca_map = {pair: IbdSegment(pair, tree.mrca(*pair)) for pair in pairs}
ibd_segs = []

while tree.next():
    for pair in pairs:
        ibd = mrca_map[pair]
        mrca = tree.mrca(*pair)
        if mrca != ibd.mrca:
            ibd.right = tree.interval[0]
            ibd_segs.append(ibd)
            # Start a new record
            mrca_map[pair] = IbdSegment(pair, mrca, tree.interval[0])

# Flush the remaining IBD segments
for ibd in mrca_map.values():
    ibd.right = ts.sequence_length
    ibd_segs.append(ibd)

for seg in ibd_segs:
    print(seg)

This gives:

1.72┊     9     ┊       9   ┊  
    ┊  ┏━━┻━┓   ┊     ┏━┻━┓ ┊  
1.72┊  8    ┃   ┊     ┃   ┃ ┊  
    ┊ ┏┻┓   ┃   ┊     ┃   ┃ ┊  
0.87┊ ┃ ┃   7   ┊     7   ┃ ┊  
    ┊ ┃ ┃  ┏┻━┓ ┊   ┏━┻━┓ ┃ ┊  
0.44┊ ┃ ┃  6  ┃ ┊   6   ┃ ┃ ┊  
    ┊ ┃ ┃ ┏┻┓ ┃ ┊  ┏┻━┓ ┃ ┃ ┊  
0.42┊ ┃ ┃ ┃ ┃ ┃ ┊  5  ┃ ┃ ┃ ┊  
    ┊ ┃ ┃ ┃ ┃ ┃ ┊ ┏┻┓ ┃ ┃ ┃ ┊  
0.00┊ 0 1 2 3 4 ┊ 0 3 2 4 1 ┊  
  0.00        30.14       40.00

IbdSegment(pair=(0, 1), mrca=8, left=0, right=30.14286231249571)
IbdSegment(pair=(0, 2), mrca=9, left=0, right=30.14286231249571)
IbdSegment(pair=(0, 3), mrca=9, left=0, right=30.14286231249571)
IbdSegment(pair=(0, 4), mrca=9, left=0, right=30.14286231249571)
IbdSegment(pair=(0, 1), mrca=9, left=30.14286231249571, right=40.0)
IbdSegment(pair=(0, 2), mrca=6, left=30.14286231249571, right=40.0)
IbdSegment(pair=(0, 3), mrca=5, left=30.14286231249571, right=40.0)
IbdSegment(pair=(0, 4), mrca=7, left=30.14286231249571, right=40.0)
IbdSegment(pair=(1, 2), mrca=9, left=0, right=40.0)
IbdSegment(pair=(1, 3), mrca=9, left=0, right=40.0)
IbdSegment(pair=(1, 4), mrca=9, left=0, right=40.0)
IbdSegment(pair=(2, 3), mrca=6, left=0, right=40.0)
IbdSegment(pair=(2, 4), mrca=7, left=0, right=40.0)
IbdSegment(pair=(3, 4), mrca=7, left=0, right=40.0)

Is this what you were looking for? Note that this is quite a strong definition of IBD, in that the MRCAs must be identical, and this may not be what you want.