Identifying mutation regions for closely related individuals without a known pedigree

The problem

Suppose that there is a (hidden) pedigree containing many (e.g., 5 to 7) generations, where we only have the genotype data on a chromosome for the individuals in the latest generation (or the latest two generations). Here those individuals with given genotype data are referred to as the input individuals. For each input individual, we also know if such an individual is diseased or normal. An example is given in Figure 1, where there are five generations in the pedigree and we only have the genotype data for the individuals in the dashed rectangle at the bottom of the pedigree. In such a figure, a square represents a male, while a circle represents a female. Moreover, a filled square (respectively, circle) represents a diseased male (respectively, female), while an unfilled square (respectively, circle) represents a normal male (respectively, female). Furthermore, if two squares (respectively, circles) enclose the same number in the figure, then they correspond to the same male (respectively, female) and their sides (respectively, circumferences) are dashed.

We also say that the pair of haplotype segments (h,h ^′ )can explain g.

Throughout this paper, we study the dominant inheritance situation, where sharing of one mutation allele can cause a disease phenotype. The general problem is as follows: we are given two sets of genotypes on the whole chromosome D = {G₁,G₂,…,G _k } and N = {G_{k + 1},G_{k + 2},…,G _n }, where the k genotypes in D are from diseased individuals and the n − kgenotypes in N are from normal individuals. The n individuals in D and N are closely related (in the same hidden pedigree). The objective is to detect the mutation regions on the chromosome, where all the diseased individuals share a common haplotype segment on the mutation region and none of the normal individuals has such a common haplotype segment on the mutation region. Note that each individual has two haplotype segments on each region. If we know the haplotypes of each input individual over the chromosome, the shared mutation regions can be computed by finding the haplotype segments which are shared by all the diseased individuals but by none of the normal ones. The true mutation region is a shared mutation region containing the disease gene. Therefore, to solve the problem, the key issue is to infer the haplotype segments.

The task of inferring the haplotypes of each individual over the whole chromosome is extremely hard. For our purpose, we divide the whole chromosome into a set of disjointed length L segments, where L is a parameter to be determined later. For each length L segment, we try to infer the two haplotype segments of each individual based on the following mathematical model.

The algorithm for MRHIP

An example of a type II case is the following: D ^R  = {g₁ = 111,g₂ = 121}, and N ^R  = {g₃ = 112,g₄ = 102}. Based on g₁ and g₂, the shared center haplotype h ^R must be 111. However, g₃indicates that normal individuals also have a haplotype segment 111 on R which is identical to the shared center h ^R . Thus, condition (2) in MRHIP does not hold.

In our algorithm, we first compute the center haplotype h ^R based on the diseased genotype segments in D ^R . We look at the positions in R one by one. Based on C1 and C2, if one of the diseased individuals has genotype value 0, then the haplotype value of h ^R at this position should be 0; if one of the diseased individuals has genotype value 1 at a position, then the haplotype value of h ^R at this position should be 1. If there exists a position p at which one diseased individual has genotype value 0 and the other diseased individual has genotype value 1, then a conflict occurs and position p is called a conflicting position. Once a conflict occurs, we simply conclude that there is no MRHIP solution on this segment R. We say the Type I False occurs in this case. If all the diseased individuals have genotype value 2 at a position p in R, then the haplotype value of h ^R cannot be determined at this step. We call such a position the wild card position and put a ∗ at the wild card position to indicate that the haplotype value of h ^R will be determined later. The detailed procedure (Procedure P1) for computing the center haplotype segment h ^R is given as follows:

Without loss of generality, for each diseased individual $g_i^R\in D^R$, we set $h_{i,1}^R=h^R$ for i ≤ k. Then we can set $h_{i,2}^R$ in such a way that $(h_{i,1}^R,h_{i,2}^R)$ is a haplotype pair for $g_i^R$ with i ≤ k. Note that, the values of $(h_{i,1}^R,h_{i,2}^R)$ at wild card positions in R are still not yet determined. Here we refer to $h_{i,2}^R$ as a partially inferred haplotype segment on R. Let $\mathrm{intQ}=\{h_{1,2}^R,h_{2,2}^R,\dots ,h_{k,2}^R\}$, where if two haplotype segments $h_{i,2}^R$ and $h_{i^{\prime },2}^R$ are identical, then we just keep one of them. Note that if any haplotype segment $h_{i,2}^R\in \mathrm{intQ}$ is undetermined at a position p, then all the haplotype segments in intQ are undetermined at p.

After partially determining h ^R and $(h_{i,1}^R,h_{i,2}^R)$ for every $g_i^R$ in D ^R , we use a heuristic method to infer the haplotype segments for $g_i^R\in N^R$. Since we want to minimize the total number of resulting distinct haplotype segments on R, our strategy is to let the inferred haplotype segments for $g_i^R\in N^R$ share as many haplotype segments as possible. This is actually Clark’s inference rule [22].

We can use h to solve$g_j^R$ by constructing two haplotype segments $(h_{j,1}^R,h_{j,2}^R)$ as follows:

A genotype segment $g_i^R\in N^R$ is solved if the pair of haplotype segments $(h_{i,1}^R,h_{i,2}^R)$ for $g_i^R$ are (partially) determined. In our algorithm, we use P to store the set of genotypes in N ^R that have not been solved. Initially, P = N ^R . We then use the haplotype segments in Q one by one and try to solve each of the genotypes in P. After trying to use a h ∈ Q to solve all $g_j^R$’s in P, we delete h from Q. Two haplotype segments on R in Q are compatible if there does not exist any position p in R such that one segment has value 0 and the other segment has value 1 at p. For two compatible haplotype segments h₁ and h₂ on R, we can merge them to form one haplotype segment h, where the value of h is determined as 0 or 1 if at least one of h₁ and h₂ is 0 or 1 and the value of h remains undetermined if both h₁and h₂ are undetermined. Again, for any position p, if the value of h₁or h₂ is not identical to that of h, then the value of h₁or h₂ is changed from undetermined to 0 or 1. Thus, we have to update some of the previously inferred haplotype segments accordingly.

When we use a h ∈ Q to solve a $g_j^R$ in P, we can obtain another new haplotype segment h ^′ on R. If h ^′ is compatible with a haplotype segment in Q, we then merge them. Note that, h ^′ might be compatible with more than one haplotype segment in Q. In this case, we arbitrarily choose a compatible haplotype segment in Q and merge the two haplotype segments. If h ^′ is not compatible with any haplotype segment in Q, we add h ^′ into Q.

We give an example to illustrate the above process.

Algorithm 1

Mutation Region Haplotype Inference

Input: Two sets of genotype segments $D^R=\{g_1^R,g_2^R,\dots ,g_k^R\}$ and $N^R=\{g_{k+1}^R,g_{k+2}^R,\dots ,g_n^R\}$ on R.

Output: True if there is a solution. Type I false if two diseased haplotype segments are conflict at a position in R; Type II false otherwise.

1: Compute the center haplotype h ^R as in Procedure P1.

2:if h ^R does not exist then

3:return Type I False

4:else

5:

set $h_{i,1}^R=h^R$ for $g_i^R\in D^R$;

6: compute $h_{i,2}^R$ according to $h_{i,1}^R$ and C1 and C2 for each $g_i^R\in D^R$; Set $Q=\{h_{1,2}^R,h_{2,2}^R,\dots ,h_{k,2}^R\}$ (removing identical segments) and P = N ^R .

7:end if

8:while Q ≠ ∅and P ≠ ∅ do

9: Delete a haplotype segment $h_{i,2}^R$ from Q;

10:if$h_{i,2}^R$ can solve $g_j^R$then

11: use $h_{i,2}^R$ to solve $g_j^R$. Add the newly obtained haplotype segments (after merging compatible segments) into Q and delete $g_j^R$ from P.

12:end if

13:end while

14:if there are at least 2 genotype segments in P then

15:if there exists a pair of genotype segments $g_j^R$ and $g_{j^{\prime }}^R$ that can share a haplotype segment on R then

16:

infer the haplotype segments of $g_j^R$ and $g_{j^{\prime }}^R$ and insert them (after merging) into Q, goto line 8.

17:end if

18:if any inferred haplotype segments for some $g_j^R\in N^R$ on R is identical to h ^R then

19:return Type II False.

20:end if

21:

fix $h_{s,1}^R$ and $h_{s,2}^R$ for each $g_s^R\in P$ so that $h_{s,1}^R\ne h^R$ and $h_{s,2}^R\ne h^R$.

22:if Line line 21 fails then

23:return Type II False

24:else

25:return True.

26:end if

27:end if

The following is an example to illustrate the case when Q becomes empty.

The algorithm for the whole chromosome

For an input segment R on a chromosome, if Algorithm 1 returns true, then R is a valid segment. In order to get the mutation regions, we decompose the whole chromosome into a set of disjointed length L segments. (In this paper, we performed experiments on chromosomes with about 110,000 SNP sites. In this case, we set L = 500.) For each segment, we run Algorithm 1 to test if the segment is valid. After finding all the valid segments, we repeatedly merge two valid segments into a long valid segment if the two segments are within 2L SNPs and Algorithm 1 returns Type II false on all the segments in the gap.After the above merging process, we obtain several long valid segments. For each such long valid segment [sb,se), we run Algorithm 1 on the three segments [sb−0.5L,se + 0.5L), [sb−0.2L,se + 0.2L) and [sb,se) and select the longest one (denoted as [b,e)) which returns true. Since we impose that Algorithm 1 returns Type II false for the segments in gaps in the merging process, we can always ensure that Algorithm 1 returns true for [sb,se). Extending [b,e) to the left and right:

After we obtain R = [b,e) as discussed above, we try to extend the segment [b,e) to the left and right. On the segment R=[b,e), we have inferred $h_{i,1}^R$ and $h_{i,2}^R$ for each $g_i^R\in D^R\cup N^R$. These $h_{i,1}^R$ and $h_{i,2}^R$ form a collection of disjointed sets H₁,H₂,…,H _m , where each H _k (1 ≤ k ≤ m) is a set of identical haplotypes in $\{h_{i,1}^R,h_{i,2}^R|i=1,2,\dots ,n\}$ on R.

The extension process stops when we find conflicts in both directions. The extended region obtained from [b,e) is denoted as [rb,re]. After the extension process, our program reports all the mutation regions obtained in the algorithm. The complete algorithm to find the mutation regions on the whole chromosome is shown as Algorithm 2.

Algorithm 2

The algorithm for the whole chromosome

Input: Two sets of genotype on the whole chromosome D = {G₁,G₂,…,G _k }and N = {G_{k + 1},G_{k + 2},…,G _n }.

Output: The detected mutation regions

1:Decompose the whole chromosome into segments of length L = 500.

2:for each segment do

3:

use Algorithm 1 to test if the segment is valid.

4:end for

5:Merge the valid segments (see The algorithm for the whole chromosome in Implementation, the first paragraph) to form longer segments.

6:for each segment [sb,se)obtained in line 5 do

7:

Select the longest segment of [sb − 0.5L,se + 0.5L), [sb − 0.2L,se + 0.2L)and [sb,se)which will return true by calling Algorithm 1. Denote it as [b,e).

8:

Extend [b,e)to get a candidate mutation region.

9:end for

10:Output all the mutation regions obtained in the for loop of lines 6-9.

Experiments on simulated data

In order to evaluate the performance of our method and the feasibility of the mathematical model proposed in this paper, we write a program in C++ to produce simulated data. The program takes a pedigree (e.g., Figure 1) and the haplotype data for the whole chromosome of each founder in the pedigree as the input. It generates the haplotype data for the remaining individuals in the pedigree using the standard χ² model for recombination with m (the degree of freedom divided by 2) equal to 4 ([25]) and according to the male/female averaged genetic map for chromosome 1 downloaded from HapMap (http://hapmap.org). Also see [26]. The haplotype data of a non-founder in the pedigree are generated to randomly inherit one strand of the four-strand chromatid bundle from each parent of the non-founder. A mutation point is selected uniformly at random from the SNP sites of the chromosome. Each diseased offspring is forced to inherit (from each of its parents) the strand with the mutation point and the normal offspring are forced to inherit the strand without the mutation point. In this way, we can guarantee that there is exactly one true mutation region. Note that the true mutation region must be a shared mutation region. See Implementation for the definition. Moreover, since we know the haplotype data of all the individuals in the simulations, we can easily find the shared mutation regions. By definition, there may exist more than one shared mutation region.

To generate the simulated data, we randomly chose some of the haplotype data for chromosome 1 of 170 unrelated Japanese in Tokyo and Han Chinese in Beijing in the database of HapMap project (http://hapmap.ncbi.nlm.nih.gov/) as the haplotype data for each founder. There are about 110,000 SNP sites on this chromosome.

Recall that our program takes two sets of individuals D and N and their genotype data as input. After generating the haplotype data of each individual, we only use some of the individuals (in the dashed rectangle at the bottom of the pedigree, say, e.g., Figure 1) and their associated genotype (obtained from the simulated haplotype data) data as the input of our program.

To evaluate the performance of our method, we used different pedigrees to evaluate our algorithm. Figures 1, 2, 3 and 4 are pedigrees of 5 generations with 2 to 5 diseased individuals in the dashed rectangle at the bottom. Those individuals in the dashed rectangle at the bottom of each pedigree are the input of our program.

The correctly detected mutation regions are the intersection of the regions reported by the computer program and the true mutation region. Here, precision is defined as the number of SNPs in the correctly detected mutation regions divided by the total number of SNPs in the regions output by the program. The value of recall is defined as the number of SNPs in the correctly detected mutation regions divided by the total number of SNPs in the true mutation region. So, if the value of recall is 1, then all the SNPs in the true mutation region have been reported by the program. Similarly, if precision is 1, then all the reported SNPs are in the true mutation region.

We performed 200 experiments for each pedigree. Since there are about 110,000 SNP sites on the chromosome, we set L = 500. For each region [rb,re]reported by our algorithm, we define a score as follows: Let DH be the number of distinct haplotypes on this region and LENGTH = (re − rb + 1)the length of this region. Then the score of this region is defined as SCORE = (2∗n − DH)∗LENGTH, where n is the total number of input genotypes in D∪N. This score can balance the length of the mutation region and the number of distinct haplotype segments on the region. With the longer region and smaller DH, the score becomes higher. To illustrate the quality of our program, we report the results when our program reports the region with the highest score and the first three regions with the highest scores, respectively. In fact, our program does not need this score in the computation. The program simply reports all the mutation regions. See the Genotype data error handling in Discussion.

From Table 1, we can see that the values of recall are from 82.71% to 98.07% and the values of precision are from 23.07% to 39.94% in the four pedigrees. When the number of diseased individuals is increased, the values of both recall and precision are improved significantly. When there are 4 or 5 diseased individuals, the value of recall is more than 97%. That is, the program can report most of the SNPs in the true mutation region.

In practice, one can often get the genotype data for the individuals of the latest two generations. Thus, we study this case by looking at different input individuals based on the pedigrees in Figures 1, 2, 3 and 4.

Now, we study different sets of input individuals in the latest two generations of Pedigree 1. These different sets of input individuals in the latest two generations in the pedigree are given in Figure 5, a square/circle with a slash indicates that such individual is not included as part of the input though the individual is used in generating the simulated data. For the rest of test, we performed 200 experiments for each case and show the average values. Table 2 shows the results for the different sets of input individuals in Figure 5. The individuals in the latest and latest two generations are not distinguished in our algorithm. We just input the genotype for all the individuals without the slash. Again, the setting is similar to that of Table 1. From Table 2, we can see that the values of recall for different inputs are close to 99% except for 2d-3fam-3, 2d-2fam-2 and 2d-2fam-3, where the input contains only 2 or 3 diseased individuals. Comparing Table 1 and Table 2, we can see that more input individuals do help improve the values of precision and recall.

We also performed similar experiments for Pedigree 2-4 (see Figures 2, 3 and 4). The results are similar to that in Table 2 and are given in Additional file 1.

Similar to the case of 5 generations, for pedigrees with 6 and 7 generations, we also tested various cases when some individuals of the latest two generations are available as input individuals. The results for 6 and 7 generations are similar to that of 5 generations. The detailed results are given in the Additional file.

A real case study

To illustrate the usefulness of our program, we applied our method to a set of real data originally from the phase II HapMap database and was studied in [27]. In [27], the authors studied two CEU (Utah residents with European ancestry from the CEPH collection) families (parent-offspring trios) CEPH 1341 and CEPH 1375 (see Figure 6). They identified a segment (107M to 110M) on chromosome 9 shared by the four individuals NA06991, NA10863, NA06985 and NA12264. For this set of data, there are totally 168,321 SNP sites on the chromosome after the unknown genotypes are eliminated from the database. There are totally 6,519 SNP sites between 107M and 110M on chromosome 9 starting at the 122348-th SNP site and ending at the 128866-th SNP site.

We applied our program with the six individuals in the two families CEPH 1341 and CEPH 1375 as input and set the four individuals NA06991, NA10863, NA06985 and NA12264 as diseased individuals. We set L = 500. Our algorithm found several segments shared by the four diseased individuals. The lengths of all the reported segments are approximately 500 SNP sites except the longest ones. All these segments are shown in Table 5. Table 5 shows the starting and ending point of the segments, the number of distinct haplotypes on the segment (DN), and the score for each segment. As there is only one conflicting position 126786 between segments [124561,126785] and [126787,129451], we should consider such a conflicting position as a data error. Therefore, segment [124561,129451] should be the predicted mutation. This segment starts at the 124561-th SNP site and ends at the 129451-th SNP site. The details are shown in Figure 7. We can see that the shared segment found by PLINK in [27] (the blue line) starts at the 122348-th SNP site while the starting position of our reported segment (the red line) is 124561. For the subsegment [122348, 124560] that we did not report, we found 79 conflicting positions in this subsegment containing 2,213 SNP sites. (See the filled dots on the blue line.) However, on the segment [124561, 129451] containing 4,891 SNP sites reported by our program, there is only one conflicting position. This is strong evidence that the subsegment [122348, 124560] is not shared by all the four diseased individuals. We also looked at the segment [129452, 131664] with length of 2,213 SNP sites on the right of our reported segment [124561, 129451], and found 52 conflicting positions among the 2,213 SNP sites. We can see that the subsegments [122348, 124560] and [129452, 131664] on the left and right of our reported segment have approximately the same number of conflicting positions.

Segment	DN	Score
[20935,21543]	8	2436
[41826,42748]	8	3692
[54648,55612]	8	3860
[58385,59106]	8	2888
[59895,60644]	8	3000
[75984,76525]	8	2168
[93972,94699]	9	2184
[97475,98191]	8	2868
[110440,111044]	8	2420
[124561,126785]	8	8900
[126787,129451]	9	7995

Haplotype inference methods

As discussed before, if the haplotype data for each input individual are known, the problem of finding the true mutation region is straightforward. Currently, there are several population-based phasing methods that can give accurate haplotype segments [28–30]. However, these methods can only phase a small number of SNPs effectively and take an extremely long time to infer the haplotype for the whole chromosome. LRP in [31] can phase more than 1,000 SNPs simultaneously within a reasonable time. However, it is still very slow for phasing 110,000 SNPs of a whole chromosome (as our program does). Moreover, they cannot directly report the true mutation region for a set of input individuals. On the other hand, our program can complete the computation in less than 20 seconds for about 110,000 SNPs with about 10 to 20 individuals.

Related mutation region detection methods

To our knowledge, all the existing software packages (except PLINK in [27]) need a clearly given pedigree as part of the input. If the pedigree is not known, most of the software packages do not work. Our algorithm deals with the case where the input individuals are closely related but the pedigree is not given.

Merlin is widely used for linkage analysis, where a pedigree is required as part of the input. It works well on SNP data due to the use of sparse trees. However, it can only analyze pedigrees of moderate size. When the family size is big, a large memory space is needed and the computation cannot successfully be completed. As shown in [19], Merlin cannot report the results for some pedigrees, e.g., P14 and P16 in [19], where there are less than 16 input individuals. However, our program can deal with the cases where the number of input individuals is large. The first row in Table 3 of the Additional file shows the results for 20 (not including those without known genotype data) input individuals. We can see that our method can give very high precision and recall in this case (without taking the pedigree as part of the input). Therefore, our algorithm can handle some cases which cannot be handled by Merlin.

The rule-based algorithm in [10] uses a set of heuristics for haplotype inference with a given pedigree. It can give very accurate results when the number of family members is large enough and for each nuclear family the genotype data for both parents are available. However, it does not work well when the genotype data of one of the parents are missing in the nuclear family. If the data for both parents are missing, it does not work. LIden [19] is an extended software package of the algorithm in [10]. It focuses on handling the case where the genotype data for one of the parents in a nuclear family are missing for the entire chromosome. But it still does not work when the genotype data for both parents in a nuclear family are missing or the family pedigree is not given in the input.

PLINK in [27] and Beagle in [30, 32] can identify the shared haplotype segment between two individuals based on population-based linkage analysis. But it cannot automatically identify the mutation region taking a set of individuals as the input, which is expected to give more precise prediction. The above real case study has illustrated this.

Genotype data error handling

The real datasets often contain errors. Handling the genotype data errors is an important issue in practice. For our program, we have a pre-process step to delete all the SNPs containing missing data. That is, if the genotype data for an input individual at an SNP site are missing, we delete this SNP site from the input. Without this step, we cannot get reasonable results for the real case study. When the genotype data contain errors, it is hard to detect and correct them. The errors may affect our program’s results in two ways: (a) an SNP site in the true mutation region may become a conflicting position due to error; and (b) the number of distinct haplotype segments to explain the genotype data is increased. When (b) occurs, our score for the detected mutation regions becomes worse. When (a) occurs once, our program reports two detected regions with the conflicting position in between. See [124561, 126785] and [126787, 129451] in Table 5 of the real case study. When this kind of error occurs many times, our program reports many regions separated by a few SNPs in the middle. When the user looks at the results of our program, it is possible to realize that the few SNPs between two closely located reported regions are due to errors. This is similar to other linkage analysis programs such as Merlin, where each SNP site has a score, and the user decides a region (by ignoring fluctuations) with high scores as the true mutation region.

Our program may work for the situation where the input individuals are from multiple families. Our algorithm tries to find regions shared by all the diseased individuals. Thus, as long as the diseased individuals from multiple families share the same (or similar) haplotype segment on the true mutation region, our program should be able to find such region. Even if the haplotype segments from different families on the true mutation region are slightly different, the program should be able to report several smaller regions with a few missing SNPs in the middle. Again, it is possible for users to figure out the whole true mutation by ignoring the few missing SNPs in the middle.