HMM Toolkit for Inferring IBD Segments from SNP Data

Implements continuous-time hidden Markov models (HMMs) to infer identity-by-descent (IBD) segments shared by two individuals. Supports two- and three-state models using single-nucleotide polymorphism (SNP) genotypes or genotype likelihoods. Provides posterior probabilities at each marker (forward-backward algorithm), prediction of IBD segments (Viterbi algorithm), and functions for visualising results. Supports both autosomal data and X-chromosomal data. The methodology and package are described in Vigeland et al. (2026) .

ibdfindr

The goal of ibdfindr is to detect genomic regions shared identical by descent (IBD) between two individuals, using SNP genotypes or genotype likelihoods. It fits continuous-time hidden Markov models (HMMs) with two- and three-state models for different relationship types. The main function findIBD() provides a unified workflow.

ibdfindr was introduced in Vigeland et al. (2026), where it was applied to identify a victim of the 1956 Marcinelle mining disaster in Belgium.

Installation

To get the latest official release of ibdfindr, install from CRAN as follows:

install.packages("ibdfindr")

For the development version, install from GitHub:

remotes::install_github("magnusdv/ibdfindr")

Example

library(ibdfindr)

As an example we consider the built-in dataset cousinsDemo, which contains almost 4000 SNP genotypes for two related individuals. The data were simulated assuming a relationship of first cousins. (For details on how the data was generated, the source code is available in the repository’s folder data-raw.)

head(cousinsDemo)

#>   CHROM     MARKER       MB        CM A1 A2  FREQ1 ID1 ID2
#> 1     1  rs9442372 1.083324 0.0000000  A  G 0.3890 A/G G/G
#> 2     1  rs4648727 1.844830 0.1803812  A  C 0.3628 A/C A/C
#> 3     1 rs10910082 2.487496 1.2670505  C  T 0.3740 C/T C/T
#> 4     1  rs6695131 3.084360 1.9782286  C  T 0.5719 C/T C/T
#> 5     1  rs3765703 3.675872 4.8403127  G  T 0.4058 G/G T/T
#> 6     1  rs7367066 3.887191 5.7748481  C  T 0.6739 C/C C/C

The function findIBD() conveniently wraps the key steps of the package:

• Fit a continuous-time HMM to the data (fitHMM())
• Find the most likely set of IBD segments (findSegments())
• Calculate the posterior IBD probability at each marker (ibdPosteriors())

ibd = findIBD(cousinsDemo)
#> Individuals: ID1, ID2 
#> Chromosome type: autosomal 
#> Input data: genotypes
#> Dropout probs (supplied): ID1 = 0, ID2 = 0
#> Model: 2-state HMM, unilineal relationships
#> Markers: 3915 input, 0 removed
#> ----
#> Fitting HMM parameters...
#>   Optimising `k1` and `a` jointly; method: L-BFGS-B 
#>   k1 = 0.237, a = 7.148, loglik = -7704.393
#> ----
#> Finding IBD segments...
#>   14 segments (total length: 597.21 cM)
#> Calculating IBD posteriors...
#> Analysis complete in 1.4 secs

For details of the different steps, see the documentation of the individual functions: fitHMM(), findSegments(), and ibdPosteriors().

To visualise the results we pass the output to plotIBD(). This plots the posterior probabilities on a background (grey points) showing the identity-by-state (IBS) status at each marker, i.e. whether the individuals have 0, 1 or 2 alleles in common. Inferred IBD regions are shown as red segments at the bottom of each chromosome panel.

plotIBD(ibd)

We may also inspect the identified segments as a data frame, showing the start and end positions, and the number of markers in each segment:

ibd$segments
#>    chrom    startCM     endCM   n
#> 1      1 113.000492 153.08623  54
#> 2      1 186.630319 223.39542  46
#> 3      2  67.118524 108.35905  50
#> 4      2 156.120361 227.62362  90
#> 5      4  59.474215  89.85522  41
#> 6      4 106.931772 131.69280  35
#> 7      6 102.041153 136.17734  43
#> 8      7 118.702115 149.11283  38
#> 9      8  72.976430 106.82971  55
#> 10    13  72.049460 127.22724  63
#> 11    14   1.808552  58.07117  67
#> 12    14  73.201009 116.00254  49
#> 13    18  10.592411  83.28443 101
#> 14    19  72.014020  99.15558  34

The full-sibling model

ibdfindr also supports a HMM tailored for full siblings, with three IBD states 0, 1 and 2. If we apply this to cousinsDemo it (correctly) reports no IBD2 segments, and in fact recovers most of the IBD1 segments.

ibdFS = findIBD(cousinsDemo, model = "fullsib")
#> Individuals: ID1, ID2 
#> Chromosome type: autosomal 
#> Input data: genotypes
#> Dropout probs (supplied): ID1 = 0, ID2 = 0
#> Model: 3-state HMM for full siblings
#> Markers: 3915 input, 0 removed
#> ----
#> Fitting full-sib HMM ...
#>   Optimising `a` conditional on kappa
#>   a = 1.402, loglik = -7750.889
#> ----
#> Finding IBD segments...
#>   IBD0: 31 segments (2798.14 cM, 83.2%)
#>   IBD1: 12 segments (563.93 cM, 16.8%)
#>   IBD2: 0 segments (0.00 cM, 0.0%)
#> Calculating IBD posteriors...
#> Analysis complete in 1.26 secs

Note that the likelihood reported for this model (loglik = -7750.889) is lower here than in the previous analysis (loglik = -7704.393), indicating that the full-sibling model is a worse fit to the data.

X-chromosome example

The brothersX dataset contains genotypes for two brothers typed with 246 X-chromosomal SNPs. The analysis below indicates that they share 3 IBD segments on the X chromosome.

ibdX = findIBD(brothersX)
#> Individuals: ID1, ID2 
#> Chromosome type: X (male, male) 
#> Input data: genotypes
#> Dropout probs (supplied): ID1 = 0, ID2 = 0
#> Model: 2-state HMM, unilineal relationships
#> Markers: 246 input, 0 removed
#> ----
#> Fitting HMM parameters...
#>   Optimising `k1` and `a` jointly; method: L-BFGS-B 
#>   k1 = 0.502, a = 6.711, loglik = -252.077
#> ----
#> Finding IBD segments...
#>   3 segments (total length: 88.05 cM)
#> Calculating IBD posteriors...
#> Analysis complete in 0.257 secs

plotIBD(ibdX)

Genotype likelihood data

In addition to standard genotype data, ibdfindr also supports genotype likelihoods (GLs), as typically obtained from low-pass sequencing data. To exemplify this, we use the built-in dataset sibsGL containing simulated genotype likelihoods for two full siblings. The data for each sibling is given in a data frame, containing base read counts, coverage and genotype likelihoods (AA, AC, …, TT). Here are the first three markers for the first sample:

head(sibsGL$data[[1]], 3)
#>       MARKER CHROM     POS A C G T Coverage           AA           AC
#> 1  rs9442372     1 1083324 4 0 7 0       11 9.682262e-15 6.133327e-16
#> 2  rs4648727     1 1844830 1 1 0 0        2 1.337778e-02 1.000000e+00
#> 3 rs10910082     1 2487496 0 0 0 6        6 1.457006e-15 1.457006e-15
#>             AG           AT           CC           CG           CT           GG
#> 1 1.000000e+00 6.133327e-16 1.244374e-24 2.028873e-09 1.244374e-24 2.536566e-07
#> 2 6.711409e-03 6.711409e-03 1.337778e-02 6.711409e-03 6.711409e-03 4.504302e-05
#> 3 1.457006e-15 1.594333e-02 1.457006e-15 1.457006e-15 1.594333e-02 1.457006e-15
#>             GT           TT
#> 1 2.028873e-09 1.244374e-24
#> 2 4.504302e-05 4.504302e-05
#> 3 1.594333e-02 1.000000e+00

(Note that the GLs for each marker are scaled so that the maximum value is 1. This is not a requisite; the values are rescaled internally as part of the processing.)

We first prepare the data with readGL(). This function accepts either a list of data frames, or a vector of filenames. The SNP annotation is provided in sibsGL$annot.

gl = readGL(sibsGL$data, sibsGL$annot)

We then call findIBD() with the "GL" input type and the "fullsib" model. The output can be visualised with plotIBD().

ibdSibs = findIBD(gl, input = "GL", model = "fullsib")
#> Individuals: S1, S2 
#> Chromosome type: autosomal 
#> Input data: genotype likelihoods
#> Model: 3-state HMM for full siblings
#> Markers: 3930 input, 66 removed, 3864 used
#> Removed markers (categories may overlap):
#>   66 with missing data (S1: 37, S2: 29)
#> ----
#> Fitting full-sib HMM ...
#>   Optimising `a` conditional on kappa
#>   a = 4.148, loglik = -6630.766
#> ----
#> Finding IBD segments...
#>   IBD0: 15 segments (841.79 cM, 25.0%)
#>   IBD1: 43 segments (1521.14 cM, 45.2%)
#>   IBD2: 26 segments (999.14 cM, 29.7%)
#> Calculating IBD posteriors...
#> Analysis complete in 0.462 secs

plotIBD(ibdSibs)