Application to whole-genome, whole-population data

We used GERMLINE to detect IBD in SNP array data from 3000 individuals, essentially the entire adult population of the Island of Kosrae, Micronesia (see Methods for dataset description). We first phased the entire population using the Beagle localized haplotype clustering tool (Browning and Browning 2007). We then applied GERMLINE to scan for shared segments over 10 cM in length according to a consensus database (Duffy 2006) of standard genetic maps (Lien et al. 2000; Kong et al. 2004).

The available pedigree for the Kosrae samples brings forth multiple pairs of related individuals as positive controls. Simplified to a single summary statistic, the overall fraction of the genome shared by a related pair, their analysis (Fig. 1) provides additional validation for the GERMLINE method, agreeing with theoretical expectation (Table 1) for relatives up to four meioses apart. Distant relationships, however, show more sharing than expected, suggesting a background of hidden relatedness in these descendants of a small, isolated ancestral population.

Open in a separate window

Figure 1.

Expected and detected genome-wide sharing (Kosrae Cohort). (Blue, ▲) Expected genome-wide sharing; (green, ■) detected genome-wide sharing.

We further demonstrate the utility of segmental IBD analysis beyond the single statistic of genome-wide sharing. To this end, we compare and resolve relationships that are indiscriminate using genome-wide averages alone. Specifically, we juxtaposed pairs of individuals that are half-siblings (one shared parent) versus those that are related through a complete avuncular relationship. These two relationships are expected to have the same amount of overall sharing statistics in terms of both , the genome-wide proportion of IBD, as well as Z₁, the overall probability of sharing one allele by IBD. In contrast, avuncular and half-sib pairs are expected to differ in segment length distribution. In accordance with expectation, Figure 2 presents a significant difference between such pairs of Kosraens in the average shared segment length, as identified by GERMLINE. No such difference is observed when comparing and Z₁. With these measures alone, one cannot distinguish between avuncular and half-sib individuals, whereas shared segment analysis facilitates resolving these relationships by average segment length as a classifier.

Open in a separate window

Figure 2.

PLINK metrics and share length for equally related pairs. Comparison of PLINK and Z₁ values with GERMLINE share length for individuals of equal relationship coefficients. (Left) PLINK values; (middle) PLINK Z₁ values; (right) GERMLINE share length (cM). Error bars, 99% CI.

We note that in our attempt to confirm these results on segmental sharing with the PLINK algorithm, it was only able to identify whole-chromosome sharing. With PLINK's focus at less-related cohorts, this tool may need specific tuning for resolving relatedness in the inbred Kosraean data set. Computationally, PLINK required 556 h to complete analysis of the eight shortest chromosomes, while GERMLINE processed the same data in 40.3 h (30.8 for Beagle phasing and 9.5 for GERMLINE analysis).

Identical matching across subsets of H

When considering a large fraction of a chromosome, true IBD may not span all of the available SNPs. Therefore, in a whole-genome or whole-chromosome context we are interested in detecting partial matches, or pairs of individuals that share a common recent ancestor only along a segment of a chromosome. As such, we establish a defined threshold on the minimum length for an IBD shared segment. The choice of L_IBD corresponds to the expected segment length for the most distantly related individuals we aim at detecting (see Table 1). Formally, we define the length L(j,j′) of an interval between columns j and j′ as the genetic distance between the corresponding genes. A pair of haplotypes (i,i′) is defined to be sharing a segment in a SNP region [j .. j′] if their included SNPs are identical and L(j,j′) exceeds L_IBD. The problem is now to find, given H, a set of quartets (i,i′,j,j′) such that (i,i′) shares the segment [j .. j′]. We now propose a divide and conquer strategy for using match to discover such shared segments. Our goal is to identify pairs of long identical segments that are shorter than the length of H. We can approximate this by dividing the columns of H into equal vertical intervals, or slices, and finding pairs of individuals that match along contiguous slices. We thus distinguish between a haplotype, which is an entire row of H and a word that represents the part of a haplotype that intersects a slice of H. A match between two individuals along several slices can be considered extended if the words of these individuals also match in the succeeding slice, otherwise the match terminates. A shared segment can thus be redefined as belonging to a pair of individuals that extend across several word pairs, and will represent an IBD segment rounded to the nearest slice break.

Formally, our algorithm accepts as input H and iteratively uses match to generate a set M′ of all shared segments in H. We vertically slice H into nonoverlapping, equal width submatrices of δ columns, with each slice denoted as H_k. The algorithm is a dynamic program that scans slices along the chromosome and maintains sets of the terminated and extendable matches within the current slice. To avoid redundant conversion from SNPs to slices, we will henceforth refer to a match (i,i′; j . j′) as (i,i′; m . m′) such that j = mδ and j′ = (m′ + 1)δ − 1. At each slice k = 0→(s/δ) − 1, we compute independent M_k = match(H_k), the set of identical matches at k. As detailed in EXTEND (Algorithm 2), we extend complementary matches in neighboring slices by examining M_k and generating a corresponding set M_k′ that contains all extended matches; naturally, each extended match contains m, the start of the match, in the range of from 0 to k-1 and m′, the end of the match, equal to k. In the first slice, define M₀′ = M₀; subsequently, initiate M_k′ = M_k, and search through all matches M_k′ for extendable matches from M_k-1′. Where they exist, we updated the starting position for matches in slice k to be that of the match in slice k − 1. Similarly, terminated matches present in M_k-1′ but not M_k, are either discarded or added as IBD to M′, depending on their length. Upon completing the final iteration, all matches in M_{(s/δ) − 1} are discarded or added to M′ in the same manner.

Genotyping error

Up until this point, we have ignored the effect genotyping errors may have on identifying matches. While modern genotyping platforms achieve accuracy levels >99% (Paynter et al. 2006), across many slices the chance of an error becomes non-negligible, even in a single sample. Across thousands of samples and along densely typed, complete chromosomes, the presence of errors approaches certainty. For IBD matching, random errors are unlikely to produce false positives. We are, however concerned about error-induced false negatives. Intuitively, a true IBD match would be present in several consecutive slices and may be detected by matching in any of them. Assuming random, uniform, and independent error rate per SNP, the number of mismatching allele calls between a pair of IBD haplotype intervals of length δ SNPs is Poisson distributed with parameter λ = 2δ. Across an IBD segment of length s_IBD SNPs [where L(s_IBD) = L_IBD], the expected number of matching slices can be modeled as

and typically being flanked by nearly identical slices. Specifically, the chance of these words to be identical is (1 − )^2δ, with the length of a complete observed match in an IBD region thus geometrically distributed.

Nearly identical matching

In practice, with δ such that E(N_IBD) >1, long IBD segments will likely have identical matches therein. The entire segment is identified by merging nearly identical flanking matches. determines the allowed mismatches in otherwise identical intervals. Specifically, conservatively assuming = 0.01; L_IBD = 2,000 SNPs; δ = 100 SNPs; E(N_IBD) = 2.7, allowing one mismatching bit/slice. The haplotypes corresponding to such slices are defined as being nearly identical. In MERGE-PARTIAL (Algorithm 3), we modify EXTEND to detect nearly identical matches by post-processing L′ to include nearly identical matches.

All of the modules integrate into a complete procedure with input matrix H and output set L of all contiguous identical or nearly identical matches in H, based on predefined length and mismatch thresholds. As described in Algorithm 4, H is divided into slices that are analyzed sequentially by MATCH, MERGE, and MERGE-PARTIAL. Unextendable matches from previous slices are discarded or output as L′. This procedure deals also with missing data.

Algorithm complexity

Computationally, the algorithm has a significant gap between worst-case and typical scenarios. Specifically, the time complexity is highly dependent on expected number of matches, which is determined by the underlying population structure. In general, if the average number of matches per word is m, the complete computation requires O(sn) to build the dictionary and O(sm) to attempt extension on all matches. In the worst case, where L is formed by the Cartesian square of the 2n haplotypes, m = 4n² and GERMLINE is comparable to pairwise exhaustive search. In practice, if we make a naïve assumption of independence of sites, the expected number of matches occurring at a slice at random is

where p_s and q_s are the allele frequencies. A more realistic assumption would acknowledge local LD within each segment. If we denote a set of population haplotype frequencies f of size f_n, the expected number of matches occurring at a word is

Even in large samples sizes, this factor is low enough where overall complexity approaches O(sn).

Gap likelihood scoring

To prioritize the segmental gaps that were most likely to be representative of a deletion, we developed a scoring function based on the number of mismatching SNPs and levels of homozygosity in the gap. From per-gap data we estimated the probability p for an independent SNP to have a mismatch. We registered the number of mismatches n and the total number of SNPs k occurring in each respective pair showing a particular gap. The binomial distribution term is therefore a likelihood score for each gap being due to mismatches consistent with the local rate of errors. Furthermore, a deletion is expected to be typed as completely homozygous, and we filtered the IBD-mismatching set for nearly complete loss of heterozygosity.

Implementation

GERMLINE was implemented in C++ and made available at http://www.cs.columbia.edu/~itsik/Software.htm. All experiments were conducted on a Linux node of 2 × 2.4 GHz Xeon CPUs with 2 GB of memory.

We thank the anonymous reviewers of this manuscript for their support and insightful comments. A.G. was supported by NIH grant 5 U54 CA121852 and the NSF Graduate Research Fellowship; I.P. was supported by NIH grant 5 U54 CA121852 and NSF grant NOA CCF-0829882; J.K.L. was supported in part by an NIH/NIDDK grant (5R01 DK60089) and by a Massachusetts Biomedical Research Corporation Tosteson Postdoctoral Fellowship. This work was further supported by grants from the Starr Foundation and Howard Hughes Medical Institute.