Obtain the stitched IBD segments between all pairs of sampled individuals.

[RomainPechey]

With Iago Bonnici (@iago-lito) and Océane Eychenne (@OceaneEYCHENNE), we were looking to obtain stitched IBD segments for all pairs of sampled individuals, which is not done natively by the ibd_segments() function. In other words, we wanted to assemble, for each pair of individuals, the adjacent segments sharing the same MRCA. We developed a script that solves this problem.

Our script has the dual advantage of being fast to execute and robust in taking into account parameters that can be added to the ibd_segments() function.

Attached you will find the code and an explanatory document schematically presenting how the combination of the convolution function and the concatenation function, which is at the heart of our approach, works.

We would be delighted to receive your feedback and comments.

Best,

Iago Bonnici, Océane Eychenne, and Romain Pechey

#################### Extract stitched IBD segments from a .trees file ####################

### Contributors :
#####    - Iago Bonnici <[email protected]> ( Research engineer at the Montpellier Bioinformatics Biodiversity platform (MBB))
#####    - Océane Eychenne <[email protected]> ( PhD student at the Institute for Evolutionary Sciences of Montpellier (ISEM) and the Center for Functional and Evolutionary Ecology (CEFE) 
#####    - Romain PECHEY <[email protected]> ( PhD student at the Institute for Evolutionary Sciences of Montpellier (ISEM) and Marine Biodiversity Exploitation and Conservation (MARBEC))
### This snippet is distibuted under the same license terms as tskit: https://github.com/tskit-dev/tskit.
# Date : 16/04/2025
# Please don't hesitate to contact us with any questions or suggestions.

import tskit
import numpy as np
import pandas as pd

# Load the .trees file
ts = tskit.load(tree_file)

# Extract unstitched IBD segments
segments = ts.ibd_segments(store_pairs=True, store_segments=True, min_span=min_span, max_time=max_time)

# Initialize empty vectors to store stitch IBD informations
Pair = []
Left = np.array([])
Right = np.array([])
Node = np.array([])
Time = np.array([])

# For each pair of individuals, extract the start and end positions of segments and the node ID of the MRCA.
# Sort these values to align them along the sequence.
for k, v in segments.items():
    kLeft = v.left
    kRight = v.right
    kNode = v.node

    order = kLeft.argsort()
    l = kLeft[order]
    n = kNode[order]
    r = kRight[order]

    # Use convolution and concatenation to merge adjacent segments with the same MRCA
    node = (np.convolve(n + (n[0] == 0), [-1, 1]) != 0)[:-1]
    discontinuities = np.concatenate(([True], l[1:] != r[:-1]))
    i = np.where(node | discontinuities)[0]

    # Recalculate left and right coordinates of IBD segments after stitching
    kPair = [k] * len(i)
    kLeft = l[i]
    kRight = np.append(r[i[1:] - 1], r[-1])
    kNode = n[i]

    # Store the values into previously defined vectors
    Pair += kPair
    Left = np.append(Left, kLeft)
    Right = np.append(Right, kRight)
    Node = np.append(Node, kNode)
    Time = np.append(Time, ts.tables.nodes.time[kNode])

# Compute segment lengths
Length = Right - Left

# Create a dataframe with stitch IBD informations
dfWF = pd.DataFrame(columns=['Pair', 'Left', 'Right', 'Length', 'Node', 'Time'])
dfWF['Pair'] = Pair
dfWF['Left'] = Left
dfWF['Right'] = Right
dfWF['Length'] = Length
dfWF['Node'] = Node
dfWF['Time'] = Time

[iago-lito]

(shall we specify licensing for the snippet @jeromekelleher or does it automatically inherits Tskit's MIT?)

[jeromekelleher]

It's up to you @iago-lito - I think code snippets like this are implicitly public domain, so go ahead and add a license if you like.

[iago-lito]

Well, it possibly depends on the country ^ ^" I think it defaults to 'all right reserved' for French law, so better add something to at least allow people to copy-paste it around. I'll have a look into it. Thank you for welcoming the snippet :)