PIGS: improved estimates of identity-by-descent probabilities by probabilistic IBD graph sampling

IBD graph

An IBD graph is constructed over a set of N haplotypes at a genomic locus as follows. Each haplotype is represented by a node and there exists an edge between nodes if the two haplotypes to which the nodes correspond are IBD at the locus. Valid IBD graphs obey a transitivity property such that if individuals 1 and 2 are IBD and individuals 2 and 3 are IBD, then individuals 1 and 3 are IBD [20]. An IBD graph is transitive if the edges obey the transitivity property, otherwise the graph is intransitive and can not represent the true state of IBD at the locus. Due to the transitivity property, all connected components of a valid IBD graph at a locus are cliques. We leverage the clique property of IBD graphs to improve the pairwise probabilities of IBD output by existing software packages such as Refined IBD [14].

Probabilistic IBD graph

Given the probability of IBD between all pairs of N haplotypes at a genomic locus, we construct a probabilistic IBD graph G _P = (N, P) as follows. Each haplotype i is represented by a node n _i ∈ N. For every pair of haplotypes there is an edge assigned probability p _ij ∈ P where p _ij is the probability of IBD between haplotypes i and j at that locus.

Updating the pairwise IBD graph

Our objective is to update the probability of each pair of individuals being IBD by conditioning on the probabilities of all pairs in the graph. The intuition is best understood with an example. Consider a probabilistic IBD graph with three nodes (i ∈ {1, 2, 3}) and suppose the initial pairwise probabilities have been assigned by a pairwise IBD-calling algorithm such that p₁₂ = 0.9, p₁₃ = 0.9, and p₂₃ = 0.1. The edges e₁₂ and e₁₃ have a higher probability of IBD than e₂₃, but given that true IBD graphs obey transitivity, the probability of e₁₂ and e₁₃ conditioned on edge e₂₃ will be lower. Similarly, the probability of e₂₃ will be higher when conditioned on e₁₂ and e₁₃ as shown in Figure 1b. By the transitive rule when 2 of the 3 edges in a triangle have a high probability of IBD we expect the third edge to have a high probability of IBD as well.

Computing exact conditional probabilities requires computing the probability of every transitive IBD graph, which has a sample space of size O(2^h(h−1)/2), where h is the number of haplotypes. Unfortunately, this is computationally infeasible to enumerate and so we develop a sampling method that can be used to efficiently approximate conditional edge probabilities.

Efficient computation of conditional IBD probabilities

We start by generating the probabilistic IBD graph for a given genomic location. We only consider the unique positions along the genome where the IBD graph changes, or more specifically, the points where the initial IBD segments begin or end. Analyzing other positions would be redundant because the positions provide no information about how the IBD graph changes. An initial graph is generated by adding in all edges output by Refined IBD [14] that pass a LOD score threshold. Alternative pairwise IBD probability methods may be used if desired. We identify the connected components of this graph because edges that are part of disjoint components have no effect on each other when computing updated probabilities.

For each connected component c, we construct G _P by translating the Refined IBD LOD scores to probabilities (see Section Converting LOD scores to probabilities). A connected component c may have edges that were never called by Refined IBD and thus have a probability of 0. We assign uncalled edges the probability ϵ = 0.0046 in order to ensure that P (g | G _P ) > 0 and that the edge can be sampled (discussed in detail in Section Converting LOD Scores to Probabilities). Then instead of enumerating the set of all possible valid graphs V inducible by G _P we sample from V. We define N _g as the current sum of probabilities of all sampled graphs so far and N _ij as the current sum of probabilities of all sampled graphs containing edge e _ij . At any stage in the sampling process, the estimate of the conditional probability that individuals i and j are IBD is ${\widehat{p}}_{ij}=\frac{N_{ij}}{N_g}$. If all valid graphs are sampled with equal probability, this converges to the exact conditional probability shown in equation (2). The sampling procedure is given in Algorithm 1.

Algorithm 1 Graph sampling

Input: G _P

Output: N _g , N _ij (∀i, j; i ≠ j)

   N _ij = 0, N _g = 0

   for all i, j do

      if p _ij ≥ 0.99 then e _ij = 1

      else e _ij = 0

   Add edges to make all connected components of G _p cliques

   Compute P(g | G _P )

   N _g + = P (g | G _P )

   for all i, j where e _ij = 1 do N _ij + = P (g | G _P )

   while ${\widehat{p}}_{ij}$has not converged ∀i, j do

      Sample a random e _ij and set e _ij = 1 with probability p _ij

      Ensure graph transitivity Compute P (g | G _P )

      N _g + = P (g | G _P )

      for all i, j where e _ij = 1 do N _ij + = P (g | G _P )

If σ ≈ 0, then edges with p _ij ≈ 1 or 0 will almost never be sampled since the selection weights of such edges will be infinitesimally small. Similarly if σ is too large, then all edges will be assigned similar selection weights and as a result graphs will be sampled uniformly instead of in proportion to their probability. We selected σ = 0.234 because it allows for efficient convergence times (see Convergence properties and runtime).

This weighted sampling scheme assures that edges with p _ij ≈ 1 or 0 are sampled less often than edges with p _ij ≈ 0.5. Intuitively this makes sense because we sample proposed IBD graphs in proportion to their respective P (g | G _P ). Edges with p _ij ≈ 1 induce a proposed graph with a greater P (g | G _P ) when they are present. Similarly edges with p _ij ≈ 0 induce a proposed graph with a greater P (g | G _P ) when they are missing. Thus, to sample the most probable graphs more often, edges with high values of p _ij should typically have e _ij = 1 and edges with low values of p _ij should typically have e _ij = 0. Changing the state of an edge can cause the proposed graph g to be intransitive. Therefore we add or remove edges from g to ensure transitivity. At each iteration if an edge has p ≥ 0.99 then p _ij is set to 1 so that we never sample very high probability edges.

Algorithm 2 Ensure graph transitivity

Input: G _P , e _ij that was just added or removed

Output: Transitive G _P

   if e _ij = 1 then

      Add edges to make all connected components of G _p cliques

   else

      S _i = nodes connected to i, S _j = nodes connected to j

      for all p _mk < 0.99 do e _mk = 0

      for each connected component X, where |X| > 1 do

         ${\bar{p}}_{iX}=\frac{p_{ix}}{\left|X\right|},{\bar{p}}_{jX}=\frac{p_{jx}}{\left|X\right|}$∀ nodes x ∈ X

         if ${\bar{p}}_{iX}>{\bar{p}}_{jX}$then Add all nodes of X to S _i

         else if ${\bar{p}}_{iX}<{\bar{p}}_{jX}$then Add all nodes of X to S _j

         else Flip a fair coin

         if ${\bar{p}}_{k{S}_i}>{\bar{p}}_{k{S}_j}$then Add k to S _i

         else if ${\bar{p}}_{k{S}_i}<{\bar{p}}_{k{S}_j}$then Add k to S _j

         else Flip a fair coin

      Make S _i and S _j cliques

We continue selecting edges until we reach a convergence criterion or a user-set time limit has been reached. As a convergence criterion we check if all edge probabilities change less than 1 × 10⁻¹¹ for 5000 sequential iterations. With larger graphs that may never reach the convergence criteria, we allow the user to set a runtime limit. After convergence or hitting the user runtime limit, we output all edges and their respective updated probabilities. We tried a variety of sampling schemes to explore the space of graphs and selected this one due to its performance over simulated datasets (see Application to simulated data).

Converting LOD scores to probabilities

To find the relationship between LOD scores and the true positive rate of IBD segments we ran Refined IBD on simulated data (see Application to simulated data) using a LOD score cutoff of 0.1 and a length cutoff of 0.1 centimorgans. A true positive segment is defined as a predicted segment that is at least 50% true IBD. We fit a curve to the observed relationship between LOD score and true positive rate of IBD segments (see Figure 2). The equation of our curve is of the form p = (2o + af)/(o + f ) where o = posterior odds, f = (prior * (10³)/.997)−(prior * (10³)), a = (1 − LOD)³/7 if LOD ≤ 1, and a = −0.15 otherwise. The values for f and a were chosen to maximize the fit of the curve.

Since the Refined IBD [14] LOD score is the negative base 10 log likelihood of one shared haplotype divided by the likelihood of no shared haplotypes, we use Bayes rule for odds to convert a LOD score into a posterior odds: $O\left(A|B\right)=O\left(A\right)*\frac{L\left(A|B\right)}{L\left(A^c|B\right)}$, where O(A | B) is the posterior odds, O(A) is the prior odds, and $\frac{L\left(A|B\right)}{L\left(A^c|B\right)}$ is the likelihood ratio.

For the prior, we use the probability of any two individuals in the sample being IBD at any point in the genome, ϵ = 0.0046. This is the average proportion of the genome shared IBD between all pairs of individuals estimated using results from Refined IBD over simulated data. For edges that have a probability of 0 (i.e. an edge with no pairwise call), we assign a probability equal to the prior because otherwise these edges would never be sampled and graphs would have a probability of 0.

Ideally, the relationship between LOD score and true positive rate is given by p = odds/(1 + odds). However, the relationship between LOD score and true positive rate in our sample of simulated individuals deviates from this theoretical relationship. Our function and p = odds/(1 + odds) are of the same form (i.e. g(x) = (cx + d)/(mx + n)). This served as our motivation in defining the conversion function. Lastly, given that the simulated data we generated is reflective of European population growth, the relationship between LOD score and true positive rate may differ in other populations (see Discussion).

Merging results across graphs and inferring new segments

IBD segments can span multiple regions and our method analyzes IBD at a single region. The probability of IBD between two individuals can therefore be output at multiple adjacent regions by our method. Furthermore, the IBD probability may be assigned a different value in each region due to the inexact nature of the sampling method. If the same IBD segment is assigned different probabilities across multiple loci we use the maximum value across all regions.

Once an IBD graph is analyzed using the sampling procedure, edges that were previously missing (i.e. those that were not called by Refined IBD) are output with a start and stop site that is equal to the intersection of all IBD segment boundaries in the graph. Since we do not look in the region for sequence identity between haplotypes we can only output the probability that IBD exists somewhere within the region. These new segments may also overlap with other called IBD segments. In order to reconcile overlapping IBD segments, we merge them provided that they pass a probability threshold set by the user and that they lie on the same haplotype. As the final probability, we use the maximum ${\widehat{p}}_{ij}$ of the merged segments. For all analyses presented here, we only merged segments that had a probability of 0.99 or greater.

Creating simulated IBD data

We generated simulated genotype data as previously described by [14]. To start, we use Fastsimcoal [21] to generate phase known DNA sequence data of 2000 diploid individuals. A single individual is represented as one chromosome consisting of ten independent 30 MB regions, each with a mutation rate of 2.5 × 10⁻⁸ and a recombination rate of 10⁻⁸. The population simulated begins with an effective population size of 3000 diploid individuals with a growth rate of 1.8% at time t = 300 (where t is the number of generations ago from the present). Moving forward in time, the growth rate was changed to 5% and to 25% at times t = 50 and t = 10 respectively, resulting in a final effective populations size of 24,000,000 at t = 0. The simulation is reflective of European population sizes estimated from the linkage disequilibrium of common variants [22].

Using the DNA sequence data we create genotype data by first filtering single nucleotide polymorphisms (SNPs) that were not bi-allelic with a minor allele frequency (MAF) less than 2%. Next, we choose 10,000 variants uniformly by MAF (where 2% ≤ MAF ≤ 50%) per 30 MB region. This SNP density is in line with that of a 1,000,000 SNP genotyping array. Finally, we remove all phase information and apply a genotyping error at a rate of .05% by turning heterozygous genotypes into homozygous genotypes and vice versa. Using the simulated genotype data, we use Refined IBD [14] to phase the data and call pairwise IBD. We define true IBD segments as those segments longer than or equal to 0.1 centimorgan. A potential consequence of this approach to creating simulated data is that the resulting IBD graph may not completely obey transitivity.

Convergence properties and runtime

We ran PIGS 25 times and calculated ${\delta }_l^{25}$ which is δ _l averaged over all 25 runs. From Figure 3 we see that for graphs with 3 to 7 nodes, edges are within 1% of true conditional probability after 5000 iterations. For 8 node graphs, the probabilities are within 15% of the true ${\dot{p}}_{ij}$ after 5000 iterations and within 5% within 7500 iterations. We recorded the average runtime of the 25 runs and show the results in (Table 1). While it is computationally feasible to sample until convergence for small graphs, this approach will not scale to genome-wide IBD studies of a large number of individuals. Instead PIGS takes as input a user specified time limit for sampling each region.

Application to simulated data

Ultimately, the metrics of merit are the IBD calls themselves, not IBD probabilities. IBD calls can be made from IBD probabilities using a thresholding approach in which all probabilities exceeding a threshold are output as IBD. Alternatively, methods such as DASH [12], EMI [19], and IBD-Groupon [18] leverage the clique nature of IBD graphs to output cliques over a region as opposed to IBD pairs. The choice of IBD calling method is a function of the objective of the study. For example, DASH was designed specifically for association testing in which individuals in a clique are given a psuedo-genotype of 1 and all others are given a pseudo-genotype of 0. Other testing methods examine the distribution of IBD between cases and controls [13, 9, 10] and rely on IBD calls that powerfully and accurately cover true IBD segments. For population genetics purposes such as inferring demographic history [5], the distribution of IBD segments sizes is the figure of merit.

This diversity of uses of IBD precludes any single metric as being the gold standard for assessing the quality of IBD calls. Therefore, we compare several different methods of computing IBD probabilities and calling IBD over a range of metrics. We compare a thresholding approach to calling IBD applied to PIGS probabilities as well as Refined IBD LOD scores. We also examine the behavior of the clique-calling approaches DASH [12] and EMI [19] when applied to Refined IBD output and PIGS output. We attempted to include IBD-Groupon but in its current implementation some hard-coded parameters make it unsuitable for the sample sizes we examined here. This will be addressed in a future release (personal communication with Dan He).

We created simulated genotype data on ten 30 MB regions for 2000 individuals (see Creating simulated IBD data). We generated IBD calls from Refined IBD by using a LOD threshold of 3. For PIGS, we first generated pairwise graphs from Refined IBD by using a LOD threshold of 0.1 and a segment length cutoff of 0.1 centimorgans. PIGS was then run over the pairwise graphs for a maximum of 2 minutes and IBD calls were made using a probability threshold of 0.99. IBD calls for Germline were generated using their suggested parameters "-haploid -bin out -min_m 1 -bits 32 -err_hom 1 -err_het 1" [17] after phasing genotype data using fastIBD [11]. DASH was run over the Refined IBD calls passing a LOD threshold of 3. All results were filtered to have a minimum segment length of 0.5 centimorgans.

Identification of IBD segments

For a given genomic locus, the power of tests comparing the distribution of IBD in cases or between cases and controls [13, 10, 9], is a function of the number of true IBD segments intersected by predicted segments. We therefore performed an analysis of the total number of true IBD segments intersected by IBD calls from Refined IBD, Germline, and PIGS. The results shown in Figure 4a show that PIGS substantially outperforms Refined IBD for small IBD segments. DASH was not included in this analysis because it was not designed for this purpose and the resulting error rates were 10 fold higher than PIGS and Refined IBD even at 1 centimorgan segments. For predicted segments of size 0.5, 0.6, 0.7, 0.8, 0.9, and 1 centimorgans, there was an increase of 95%, 43%, 27%, 17%, 12%, and 9% in the number of predicted segments intersecting a true segment over Refined IBD. For predicted segments of size 1.1, 1.2, and 1.3 centimorgans, Germline was able to detect 60%, 27%, and 12% more segments than PIGS but the calls were less accurate (See Accuracy of IBD segments).

In order to assess the error rate we examined the fraction of segments that did not intersect any true IBD segment. Note that this is error rate may be inflated due to the fact that true segments are required to be at least 0.1 centimorgans (see Creating simulated IBD data). The results shown in Figure 4b demonstrate that PIGS has nearly identical error rates to Refined IBD at small segments. However, at 0.5 centimorgans the error rate increases from 0.3% to 0.7%; this is a modest increase relative to the 95% increase in the number of segments identified. Germline, in general, was similar to the other methods in terms of error rates for segments between 1 and 2 centimorgans.

Accuracy of IBD segments

In population genetics settings, such as inferring demography [23, 5], methods often rely on the distribution of IBD segment lengths. The figure of merit here is related to the accuracy of predicted segments recovered. We first examined power, the average proportion of true IBD segments that were overlapped by predicted segments. For true segments between 0.5 and 2.5 centimorgans our method had modestly greater power (Figure 5a). This came at the expense of a slight increase in false discovery rate (FDR) as shown in Figure 5b. The false discovery rate is defined as the average proportion of predicted segments that does not overlap a true IBD segment. This is somewhat expected since new segments from PIGS use existing IBD segment boundaries to approximate the new start and stop sites. The greatest decrease comes at 0.5 centimorgans where PIGS predicts 95% more segments than Refined IBD. However, the difference in FDR between Refined IBD and PIGS is still less than 5% (10% versus 14%). On the other hand, for segments of size between 0.6 and 1.5 centimorgans, PIGS predicts 23% more segments while keeping the FDR within 1% of Refined IBD.

We also examined the true positive rate, defined as the percentage of predicted segments with at least 50% overlap with a true segment. Compared to Refined IBD the true positive rate for PIGS drops slightly for segments that are smaller than 1 centimorgan but the difference is less than 1% for all sizes except for at 0.5 centimorgans where it is 3% (Figure 6). The reason for this drop in performance is at least partly due to the fact that we add new segments according to the IBD graph without specifically examining the sequence. Given the results of the previous section, the most likely explanation is that the PIGS predicted segments intersect true IBD segments, but not at the 50% threshold required by definition of a true positive. Based on these results PIGS could be used for population genetics purposes, but users should take into account the slight increase in error rates for smaller segment sizes.

Application to real data

Identification of IBD segments

We applied PIGS, RefinedIBD, DASH, and EMI to 489 Latino trios from the Genetics of Asthma in Latino Americans (GALA) cohort [25]. The availability of trio genotype data allows us to phase the genotypes with high accuracy by taking into account the rules of Mendelian segregation. The increased phasing accuracy in turn boosts the power to detect IBD segments because phasing errors are a major source of difficulty in calling IBD [26]. To evaluate how well a given method is able to identify segments of IBD in real data we used IBD segment calls made using Refined IBD in a trio-aware mode (trio-IBD segments). Trio-IBD segments were thresholded at a length of 0.1 centimorgans and at a LOD score of 3. We then asked how many IBD calls made by a given method without access to the near-perfect trio phasing overlapped with trio-IBD segments. DASH and EMI were run in the same way as we described for identifying cliques (using PIGS and RefinedIBD as input), however the resulting clique edges were converted into IBD calls and merged with the original input. We include the clique calls here as IBD calls because we do not know the structure of the real IBD graph. All IBD calls were thresholded at a segment size of 0.5 centimorgans.

As shown in Figure 7, when considering PIGS and Refined IBD calls at 0.5 centimorgans, there is an increase of 10% in the number of segments identified by PIGS over Refined IBD. After applying DASH and EMI to the input of both methods we see an increase of 8% and 7%, respectively, for PIGS input. It is clear that both DASH and EMI improve the power of both main approaches to detect IBD for use in association studies regardless of the segment size. DASH and EMI seem to perform similarly in terms of boosting power when called segments are bigger than 0.8 centimorgans, but EMI appears to have the upper hand for anything smaller. For example, at 0.5 centimorgans the difference between EMI and DASH for PIGS input is 8% but at 0.8 centimorgans the difference is only 0.8%. Across all segment sizes, we see increases of 4%, 3%, and 2.5% for PIGS, P-DASH, and P-EMI over their Refined IBD counterparts. The increases are more modest than in the simulated data, most likely due to the fact that without sequencing data we are underpowered to detect small segments of IBD even when trio phased genotypes are available.

PIGS and Refined IBD called 3134591 and 2968480 segments that overlapped with at least one of trio-IBD segment, respectively, which equates to a 6% increase. Similar increases are also seen when using PIGS input to DASH and EMI, with a 4% increase (3177234 versus 3047734) for DASH and a 3% increase (3263594 versus 3158818) for EMI. 1330207 PIGS and 803527 Refined IBD calls did not overlap with any trio-IBD segments. Because we only have access to true positives in the real data, there is no perfect way to determine the false positive rate of any of these methods, and it could be argued that PIGS increases the power to detect IBD at the expense of a higher false positive rate. To determine if this was the case, we made random IBD calls along the genome. As an example consider the calls of size 0.5 centimorgans, where PIGS made 349221 calls and Refined IBD made 262064 calls, a difference of 87157 calls. Of these, 207122 PIGS and 187582 Refined IBD segments overlap a trio-IBD segment, which is an increase of 19540 segments. After making 87157 random calls (of length 0.5 centimorgans) we only identified 453 segments compared to the 19540 we observed originally. This means that even if Refined IBD were to make an additional 87157 random guesses along the genome we would not expect Refined IBD to have the same power as PIGS, showing the increase in performance is not entirely due to false positives. Furthermore, if we assumed all non-overlapping segments are false positives both Refined IBD and PIGS have an error rate over 20%, which is not reflective of the simulations where the error rate for both methods was below 1% (see Figure 4b) [14].

Given these results, we conclude that the increased performance in PIGS was not driven by the extra IBD calls and that the majority of the non-overlapping segments are indeed true as suggested by the simulation results. Assuming that the true false positive rate in real data is similar to simulation data, the large increase in predicted segments that overlap trio-IBD segments when using PIGS (with or without a clique calling method) shows the potential for substantial power increases when using PIGS for IBD mapping studies.