hmmIBD: software to infer pairwise identity by descent between haploid genotypes

1Segments of DNA that are shared between individuals and inherited without recombination from a common ancestor are said to be identical by descent (IBD). IBD is a fundamental concept in genetics, linking ancestry to variation due to recombination, which acts on a shorter timescale than mutation. In the field of human genetics, IBD-based analyses are used in many different applications: to map disease loci and quantitative traits, to phase and impute genotypes, and to infer demographic histories.

2Increasingly, IBD-based analyses are also being used to study haploid organisms such as the malaria parasite. Examples include studies of malaria disease transmission, malaria parasites within multiple-genotype infections, signals of parasite selection, anti-malarial resistance, relatedness across proximal parasite populations, and the relatedness of co-transmitted parasites strains. However, most existing IBD detection software (recently reviewed in) is designed for humans and other diploids. Accordingly, malaria studies have typically used one of two programs designed specifically for haploids: isoRelate, used and described in and based on, and hmmIBD, which was used in, but has not been fully described in the literature to date.

3This paper provides a full description of hmmIBD, its verification using simulated data, and validation and benchmarking of its output and performance against isoRelate. Also described is a novel feature of the program, one that allows inference of IBD between samples from distinct populations.

4hmmIBD is based on a discrete, heterogeneous, first-order HMM with two hidden states, IBD and not-IBD; mathematical details can be found in Additional file 1. It is designed to infer IBD segments shared between pairs of haploid genomes and to estimate two quantities: (1) the marginal posterior probability of the IBD state (which can be interpreted as the expected fraction of a pair of genomes that is IBD); and (2) the rate at which the genomes transition between states, parameterized by the number of outcrossed meioses since their most recent common ancestor (MRCA), which we refer to as the number of generations. These parameters are estimated iteratively using the Baum-Welch estimation-maximization algorithm; the Viterbi algorithm is then used to calculate the most probable sequence of IBD segments.

5Model specification requires three probability measures. First, initial state probabilities (IBD or not at the first position on a chromosome). These are initially set to 0.5, and then updated as the expected fraction IBD is recalculated for the entire genome under successive fits of the model. Second, the probabilities of changing state between two genomic positions. These state transition probabilities are functions of the distances between positions (in base pairs), the recombination rate and the number of generations since the MRCA (both assumed to be uniform across the genome), and the expected fraction IBD. Third, emission probabilities, which are the probabilities of the observed allelic types given IBD or not. These are calculated as follows. If two genomes are IBD at a given genomic position, they are of the same allelic type, meaning no mutations are assumed to have occurred since the MRCA. If they are not IBD, the alleles are modelled as independent draws based on the allele frequencies in the population. The probability of the observed alleles is then calculated from these probabilities by allowing for genotyping error.

6Parameters include the recombination rate (the default value in the code is for Plasmodium falciparum) and the genotyping error rate, as well as estimates of the allele frequencies for all variable sites in the input dataset. Allele frequencies can be supplied by the user; if not supplied, the program estimates them from the data. The model accommodates positions with missing data by omitting emission probabilities at those sites. The program is implemented in C; it complies with the C11 standard and requires no additional libraries. It can accept genotype data for any variant that can be represented by an integer at a single position on the chromosome (e.g. all SNPs and most small indels).

7One important assumption of the model is that all IBD regions present between two samples result from common ancestry on a similar timescale. Clearly, this need not be the case: very recent inbreeding can be present along with much older background sharing, the kind that generates linkage disequilibrium in the population. Because IBD deriving from recent common ancestry is of primary interest for many applications, hmmIBD provides an option of capping the number of generations to the MRCA in the model; its effect is to bias against breaking up segments of either state.

8The correctness of hmmIBD’s algorithm and code was verified using data simulated under the HMM on which it is based; details and full results can be found in Additional file 2. For a typical situation, the root-mean-square (RMS) error on the fraction of the genome called as IBD was 0.25 percentage points and on the number of generations was 2 generations. CPU time was linear in the number of variants and quadratic in the number of samples. The performance of hmmIBD was then validated and benchmarked against isoRelate using data created from artificially recombined field samples; details and full results can be found in Additional file 3. Both programs performed well, with accuracies, sensitivities and specificities greater than 98% (Table 1), and performance remained high when genotyping error equal to 0.5% in the artificially recombined genomes was misspecified under the model at 0.1% (Additional file 3). On average, hmmIBD was approximately 24 times faster than isoRelate, although both perform adequately in real time. (We note that isoRelate has a unique capability: by modelling the hidden state space as a set of IBD allele counts in 0, 1 or 2, it can accommodate comparisons across samples containing two genotypes.). The recombination rate under hmmIBD is assumed uniform. Given deviations within a biologically informed range, the assumption has little impact upon inference on IBD under hmmIBD using data within the accessible genome.

9An unusual feature of hmmIBD is that it can accommodate samples from different populations, even ones with very different allele frequencies. This feature should have multiple applications in studying malaria, including detection of selective sweeps that spread between populations and determining the source of imported cases in elimination settings. The effectiveness of hmmIBD can be seen by using it to examine IBD in the region of the well-known selective sweep for chloroquine resistance in P. falciparum around the gene pfcrt. Figure 1 shows the amount of IBD between field isolates from Ghana and those from Cambodia. The increase in IBD around pfcrt is clear, reflecting the fact that the resistance haplotype emerged and spread into Africa from South-East Asia. An alternative approach using hmmIBD but treating the samples as coming from a single population as in, also shown, is much less effective at detecting the cross-population IBD.

10This study provides a stand-alone description of hmmIBD, software to infer IBD between haploid genomes, and its comparison to isoRelate. Although both programs performed well, they differ in speed (hmmIBD being approximately 24 times faster than isoRelate) and in practical application. Most notably, isoRelate can handle comparisons across samples containing two genotypes where hmmIBD cannot, while hmmIBD can handle comparisons across populations where isoRelate cannot. The latter application is potentially useful for WHO malaria elimination certification, which requires no local cases for at least 3 years, but allows for imported ones.

11hmmIBD continues to perform well even when the genotyping error and recombination rate is misspecified, and appears to be relatively robust to violation of the assumed uniformity of the recombination rate. These results are unlikely to hold over regions where rates deviate greatly from their average. Exclusion of such regions in data analysed under hmmIBD is therefore recommended. It is also advisable to change the default error rate in the code if a more accurate rate is available. Similarly, if used with a species other than P. falciparum, the appropriate average recombination rate should be used. In both cases, specifying an incorrect parameter is likely to affect primarily the number of generations rather than the fraction IBD, providing regions that deviate greatly from the average parameters are excluded. Finally, large public datasets (e.g. from the Pf3k project) may provide better estimates of allele frequencies in a population than can be calculated from the data immediately available to a user; it is advisable to make use of them where available.

12There has been a recent proliferation in the number of malaria studies using IBD. As interest in IBD grows, the need to provide comprehensive details of software used to infer IBD increases. hmmIBD is the only program known to the authors that is designed specifically for haploid malaria genomes and is capable of comparing samples across populations with different allele frequencies. This will likely be a useful feature for malaria elimination efforts, since it could facilitate identification of imported malaria cases.