Efficient whole-genome association mapping using local phylogenies for unphased genotype data

2.1 Building local phylogenies

We infer local trees in two different ways, depending on whether a perfect phylogeny—i.e. a genealogy consistent with the infinite-sites mutation model without recombination—exists for a sufficiently wide region around the local site we wish to score. When building trees, we require at least m markers be used in the inference, where m is an option to the Blossoc program. In general, we recommend using larger (respectively, smaller) values of m in regions with high (respectively, low) LD. When m markers are compatible with the infinite-sites and no recombination assumption, we say that a perfect phylogeny exists, and we can construct trees directly from unphased data. When the m markers are incompatible with the assumptions, we perform a local phase inference for the m markers only and construct a tree using heuristics from Mailund et al. (2006).

Given a set G of n genotypes of the same length, the Perfect Phylogeny Haplotyping (PPH) problem is to find n pairs of haplotypes that explain G and fit a perfect phylogeny. Vijayasatya and Mukherjee (2005 and Ding et al. (2005, 2006) have independently developed linear-time algorithms for this problem. In this article, we utilize the latter algorithm, called LPPH.

Our algorithm constructs a local tree for each marker, incorporating neighboring markers as follows. To construct a local tree for a marker x, initialize X to be the set containing only x. Then, alternate the following two steps until neither is possible:

1. If X and the next marker immediately to the left together admit a PPH solution, then add that marker to X.

2. If X and the next marker immediately to the right together admit a PPH solution, then add that marker to X.

Selecting a set of markers in this way keeps the focal marker x near the center of the region.

If the resulting X contains at least m markers, where m is specified by the user, we score the tree as described in the next section. However, for a large sample size and high recombination rate, the size of regions for which PPH solutions exist tends to be very small. Such a region will contain only few markers and consequently the corresponding inferred tree will contain only few edges, making it difficult to infer reliably the clustering of cases and controls. A region with a set of m markers that does not admit a PPH solution is called an incompatible region. For such a region, we use an entropy minimization algorithm—previously considered by Halperin and Karp (2005) and Gusev et al. (2007)—to locally infer the phase of input genotypes and use the tree building algorithm described in Mailund et al. (2006) to build local phylogenies from the inferred haplotypes. The entropy minimization algorithm has phasing accuracy slightly worse than that of other widely-used methods such as fastPHASE (Scheet and Stephens, 2006) and HAP (Halperin and Eskin, 2004), while being several orders of magnitude faster (Gusev et al., 2007).

Given a haplotype h and a phasing solution ϕ for a set of n genotypes in G, define the coverage of h under ϕ, denoted by COV(h,ϕ), as the number of genotypes in G that are phased by h and some other haplotype in ϕ, plus twice the number of genotypes in G that are phased by the haplotype pair (h,h). For a fixed phasing solution ϕ, the sum of COV(h,ϕ) over all haplotypes in ϕ is equal to 2n.

As in Halperin and Karp (2005) and Gusev et al. (2007), we define the entropy of a phasing solution ϕ as

Halperin and Karp (2005) first introduced the problem of finding a phasing solution ϕ of G with the minimum entropy, and Gusev et al. (2007) later developed an accurate and highly efficient algorithm to solve the problem. Intuitively, the entropy of a phasing solution ϕ is a measure of haplotype diversity in ϕ. As is shown in Gusev et al. (2007), if the probability of a haplotype h is estimated by counting the number of times that h appears in a phasing ϕ, then maximizing the log-likelihood of phasing ϕ is equivalent to minimizing the entropy of phasing ϕ.

The phasing algorithm in Gusev et al. (2007) employs a set of short overlapping windows to phase long genotypes. In our work, we typically deal with short genotypes and hence use the following simple version of that algorithm:

1. Generate a random phasing solution ϕ for genotypes G.

2. Repeat the following:

   1. Find the pair (g, (h₁, h₂)) such that H(ϕ′) is minimized, where ϕ′ is obtained from ϕ by re-explaining g∈G with (h₁, h₂).
   2. If H(ϕ′)<H(ϕ), then let ϕ=ϕ′, else exit the loop.

3. Output phasing solution ϕ.

Missing data are handled as follows. When using our algorithms to obtain a phasing solution, we first phase those genotypes without any missing entries and calculate the frequency of each distinct haplotype in a solution ϕ′ to that subproblem. Then, for each genotype g with one or more missing entries, we test whether the non-missing entries in g can be phased using a haplotype in ϕ′ and some other haplotype. If more than one haplotype in ϕ′ satisfies this criterion, we choose the one with a higher frequency. (If there is a tie, we choose an arbitrary one among them.) Missing entries in g are then imputed according to that phasing. If none of the haplotypes in ϕ′ can be used to phase g in this way, then g is resolved arbitrarily.

2.2 Scoring local phylogenies

Once a local phylogeny is constructed, it is scored according to how well the tree helps explain the phenotype. We consider the tree a hierarchical clustering of chromosomes: each partition of the tree into subtrees defines a clustering. To each cluster c, we assign a disease risk θ_c—assumed to be independent of other clusters—and the disease status of a chromosome in cluster c is modeled as being affected with probability θ_c and unaffected with probability 1−θ_c.

Given a clustering 𝒞={c₁,…,c_n} and corresponding disease risks Θ={θ₁,…,θ_n}, the likelihood of the observed disease status is given by

where A_i denotes the number of affected leaves in cluster c_i and U_i the number of unaffected leaves in cluster c_i. Assigning a risk per cluster ignores that phenotype risk may be a function of genotypes, and thus a function of two clusters rather than one. It is a straightforward extension to model phenotype as a function of genotypes, although the implementation of the score functions would then get somewhat more involved.

When scoring a tree, the clustering and risks are nuisance parameters that we integrate out in a Bayesian approach. For the cluster specific risks, we follow the approach in Waldron et al. (2006) and choose independent uninformative β-priors π(θ)=1 and obtain

where B is the β function. To integrate out 𝒞, we again choose a uniform prior on clusters, obtaining the following final score:

where |𝒞| denotes the total number of clusterings possible in the tree. Summing over all clusterings is computationally prohibitive, so in our implementation we only sum over all sub-trees and cluster the chromosomes into only two clusters, those that are leaves in the sub-tree and those that are outside. When there are more than one perfect phylogeny consistent with a region, we average the scores over all trees for that region. This corresponds to integrating over the unknown tree with a flat prior over topologies.

As a Bayes factor, any number greater than 1 can be taken as evidence for association, while any number smaller than 1 can be taken as evidence against. However, some scores higher than 1 can still occur by chance, so we recommend doing permutation tests to judge the true significance of a score.

3.1 Ranking experiments

In genome-wide association studies, false positives are a major problem, meaning that potentially many leads may need to be followed before a true hit can be found. Hence, the quality with which a method ranks the true hits compared to spurious hits is one of the most important measures of performance.

In this section, we describe how well Blossoc performs in ranking. We simulated 100 case/control datasets, each with 1000 cases and 1000 controls, under an additive disease model with varying genetic relative risk. With A denoting the disease-predisposing allele and a the wild-type, the genotype relative risk of the heterozygote Aa is denoted by GRR. Since we assumed an additive model in our simulations, the GRR of the mutant homozygote AA was 2×GRR−1. Experiments with dominant or recessive (instead of additive) disease models produce similar results (results not shown). The wild-type risk (denoted WR in Tables 1 and 2) of being diseased was varied between 5% and 10%, while GRR =1.2, 1.4, 1.6, 1.8 and 2.0 were used. The SNP density was also varied; we used 100 SNP markers with the population-scaled recombination rates ρ=4, N_e c=100 and 400, where N_e denotes the effective population size and c the recombination rate per generation per sequence. Assuming N_e=10 000 and a recombination rate of 10⁻⁸ per adjacent pair of sites per generation, ρ=100 and 400 correspond to 250 kb and 1 Mb, respectively. Sequences were simulated using the coalescence based simulator CoaSim (Mailund et al., 2005) with the infinite-sites mutation model, and were then paired up randomly to construct diploid individuals. After assigning disease status, the disease-predisposing SNP was removed from the data.

We performed analysis both using the original phase-known version of Blossoc (Mailund et al., 2006) on the true phased data, and using our new algorithm on unphased data. We first studied how the performance of Blossoc depends on the minimum number m of markers to include when building a tree. In general, larger values of m should be used for regions with high LD and smaller values for regions with low LD. In our experiment, we tried using m=3, 5 and 10, and found m=5 to be slightly superior to m=3 and m=10, for both ρ=100 and ρ=400 (results not shown). In what follows, we therefore use only m=5.

As a way of summarizing the ranking results, we considered the mean fraction of top-10 ranked markers found within 1, 10 or 100 kb of the disease-predisposing marker location. See Table 1 for results. The results for an 2×2 allelic χ²-test, applicable to an additive disease model, are also shown there. For each parameter setting, we have highlighted the best performing method in bold print. Table 2 shows the fraction of datasets with at least one top-10 ranked marker within 1, 10 or 100 kb of the disease-predisposing marker. Measured this way, we see that, close to the true disease-predisposing marker, the phase-known Blossoc method performs the best, followed by our new phase-unknown Blossoc method, and then the χ²-method. Farther, away from the disease-predisposing marker (the 100 kb columns), however, the χ²-method generally outperforms both Blossoc methods. An explanation for this is the higher correlation between neighboring markers in the Blossoc method. One spurious signal will tend to give spurious signals at neighboring markers, and a high scoring marker far from the true locus can therefore move most (or all) of top-10 away from the true locus. Consequently, considering a raw ranking—ignoring the correlation between markers—is probably not optimal for Blossoc. Future work will include ways of ranking Blossoc scores in a more meaningful way.

3.2 Localization experiments

Another important criterion for assessing the performance of a fine-mapping method is localization; that is, how accurately the position of an untyped disease-predisposing SNP can be estimated. In our localization study, we simulated case-control datasets using CoaSim for a region corresponding to 100 kb physical distance (or ρ=40). We assumed that the region contains exactly one disease-predisposing SNP. We used GRR=1.4,1.6,1.8 and 2.0, and assumed that the wild−type risk is 5%. We simulated 100 datasets with 500 case and 500 control individuals, as well as 100 datasets with 2000 case and 2000 control individuals.

If a unique marker had the highest score, then it was chosen as our estimate of causal SNP position. Otherwise, if k > 1 consecutive markers i₁,…,i_k had the highest score, we took the midpoint between markers i₁ and i_k as our estimate. We did not encounter any case in which two markers with the highest score were separated by a marker with a lower score.

Empirical cumulative distributions of localization error ∊ (in kb) are shown in Table 3 for 500 case and 500 control individuals, and 2000 case and 2000 control individuals. Three different methods were used to analyze the data: single-marker χ²-test, Blossoc on phased data and Blossoc on unphased data; we used m=5 when using Blossoc. The faster the cumulative distribution approaches 1 as ∊ increases, the better the method. Several things are worth-while noting. First, the localization results of Blossoc on phased data and that of the method developed here for unphased data are comparable. Second, both versions of Blossoc (for phased and unphased data) perform better than χ² as GRR increases. Third, the amount of improvement increases as the number of case/control individuals increases.

We also compared Blossoc with Margarita, the method developed by Minichiello and Durbin (2006). Figure 2E of Minichiello and Durbin (2006) shows a plot of empirical cumulative distribution of localization error for 1000 case and 1000 control individuals, under an additive model with GRR=2.0 and the frequency of the disease-predisposing allele equal to 0.04. In their simulation, a region of size 1 Mb with ρ=440 was used and 300 tagging SNPs were selected. We used Blossoc to analyze the same simulated datasets. Table 4 compares the localization results of the χ²-test, Blossoc and Margarita. Clearly, both Blossoc and Margarita are more accurate than the χ²-test in terms of localization. Further, Blossoc is more accurate on unphased data than on phased data; the same behavior was observed in Margarita. In general, Margarita is more accurate than Blossoc, but this increase in accuracy is obtained at the expense of significantly longer running time; for the simulation study shown in Table 4, Blossoc took about 3 s per phased dataset and 905 s per unphased dataset, while Margarita took 118 620 s per phased dataset and 300 512 s per unphased dataset. That is, Margarita is slower than Blossoc by a factor of 40 000 for phased data and a factor 300 for unphased data. Such increases in running time can make the difference between a feasible and an infeasible analysis. More results on running time are discussed in Section 3.4.

3.3 Coriell Parkinson's disease genome-wide dataset

Parkinson's disease is a progressive neurodegenerative disorder, affecting more than one per thousand individuals Kuopio et al. (1999). Fung et al. (2006) carried out a genome-wide SNP genotyping assay of publicly available samples from a cohort of 267 Parkinson's disease patients and 270 neurologically normal controls. A total of 408 803 unique SNPs were used from the Illumina Infinium I and HumanHap300 assays. We only considered markers with less than 5% missing data, but applied no other filters. We analyzed the unphased version of this dataset using Blossoc with option m=5. The entire analysis took about 4 h on an Intel(R) Pentium(R) 4 CPU 3.00 GHz, 512 Mb RAM. In comparison, it took fastPHASE 2 days to phase only chromosome 21 of the dataset (containing 6612 SNPs). Beagle (Browning and Browning, 2007), a much faster phase inference tool, is capable of phasing the entire dataset in a reasonable time, but still takes about 8 h, while Blossoc takes 4 h including the association test.

The top-10 highest scoring SNPs found by Blossoc are shown in Table 5. Fung et al. (2006) performed a single-marker genotypic χ² association test and found 26 SNPs with uncorrected p values less than 10⁻⁴. However, none of the SNPs showed significant association after Bonferroni correction. One of the 26 SNPs was in the set of top-100 SNPs ranked by Blossoc, and other 13 of the 26 SNPs were within 100 kb from at least one SNP in Blossoc's top-100.

We also compared the SNPs in Blossoc's top-100 with the set of top-100 loci ranked by single-marker χ²-test. One third of Blossoc's top-100 were within 20 kb of at least one top-100 χ²-SNP while 58 where farther than 1 Mb away from any χ² top 100. Such an overlap seems reasonable between multi-marker and single-marker association methods.

3.4 Running time

We examined the running time of Blossoc on datasets simulated using FREGENE (Hoggart et al., 2007), software capable of simulating large-scale sequence data. Table 6 shows a summary of running time results. In light of the fact that fastPHASE takes 3 days to phase a dataset with 500 individuals and 10 000 SNPs, speed clearly is a notable advantageous feature of Blossoc. Note that the running time of Blossoc is roughly linear in the number of SNPs, and it depends more on the number of individuals than on the number of SNPs.