Analysis of high-density single-nucleotide polymorphism data: three novel methods that control for linkage disequilibrium between markers in a linkage analysis

The recent availability of rapid, accurate, and relatively low-cost genome-wide, high-density single-nucleotide polymorphism (SNP) panels is changing the study of many complex diseases. Linkage analysis to track inheritance of chromosomal regions in pedigrees no longer must rely on highly informative but sparsely spaced conventional microsatellite markers. Rather, SNPs, which are far more abundant than microsatellite markers, have the capacity to yield a higher information content, and hence an improved potential for localizing disease genes [1].

Classical linkage analysis assumes linkage equilibrium between markers. Applying classical linkage packages to genome-wide SNP panels may generate biased results, if parental genotypes are missing, as SNP panels will include markers that are in linkage disequilibrium (LD). LD between markers may inappropriately weight shared haplotypes in the likelihood calculations, which can lead to inflated LOD scores, resulting in false-positive evidence of linkage [2]. To overcome the potential bias in a linkage analysis by using SNP data, ideally LD must either be incorporated into the analysis or markers in high LD must be eliminated. The commonly used, freely available SNPLINK [3] program eliminates high-LD SNPs using a simple algorithm between contiguous pairs of SNPs. However, we have previously shown that a non-contiguous LD structure is more likely to model the complex pattern of recombination and mutation that have been observed to exist in candidate genes [4, 5].

Here we propose three novel methods to control for LD between markers in a linkage analysis that do not require a contiguous LD structure: one method that allows for LD between markers using graphical modeling and two methods that eliminate markers in high LD using principal components. Our three novel methods were compared to the previously published analysis generated by the SNPLINK [3] program. We used genome-wide scan data made available for the Genetic Analysis Workshop (GAW15), Problem 2. The North American Rheumatoid Arthritis Consortium (NARAC) provided 5744 SNPs from a genome-wide scan of 757 families, approximately 90% of which were Caucasian. The NARAC Caucasian families were previously analyzed by Amos et al. [6], who found their two highest linkage results on chromosome 6 (LOD = 18.53) and chromosome 21 (LOD = 11.59), not accounting for LD. However, when markers in high LD were eliminated using SNPLINK [3], Amos et al. [6] found that the chromosome 6 peak remained (LOD = 16.14) but the chromosome 21 peak disappeared (LOD = 1.11). By reanalyzing the chromosome 6 and 21 peaks with three novel LD control methods, we verify the previous results and suggest alternative methods for controlling for high LD without requiring contiguous LD structure.

Details regarding NARAC subject enrollment, phenotype designation, and genotype information are published elsewhere [6]. Using the total 404 SNPs available for chromosome 6 and the 104 SNPs available from chromosome 21, we checked for genotype errors using the software CheckErrors [7], which uses graphical modeling to calculate the posterior probability of genotype errors in pedigrees. For chromosome 6 we eliminated 16 SNPs and for chromosome 21 we eliminated 1 SNP because of errors. We performed a linkage analysis on each of these cleaned files using the multipoint Markov-chain Monte Carlo (MCMC) linkage method MCLINK [8], which computes the robust multipoint TLOD [9] statistic. For each chromosome, we identified peak regions (TLOD ≥ 2.0), not accounting for LD. To ensure that the peak could be resolved in other analyses, we included ~20 SNPs on either side of the peak and called these data our "Complete" SNP sets. These linkage results are likely to be biased because of underlying LD as 63.7% of the parents in the data set were not genotyped.

The four methods (i.e., the three novel methods and SNPLINK) were each applied to the "Complete" SNP data files. For the graphical modeling of LD combined with linkage analysis method, we used the new McLink software http://www-genepi.med.utah.edu/~alun/software/, which we will refer to as McLink-LD in this manuscript. For the three LD elimination methods, we used both the original MCLINK [8] and Merlin [10] software packages to perform the linkage analyses. Allele frequencies were estimated from observation of all genotyped individuals at each locus, rather than only the founders. Genotype data was available for only ~15% of the pedigree founders, and because the total number of individuals with genotype data was large (n = 1991 individuals) the estimated allele frequencies from all genotyped individuals provides reasonable estimates of the founder allele frequencies, even without adjustment for familial relationships. The original pedigree structure for all 757 pedigrees was used for all analyses, except for the Merlin analyses, in which 56 large pedigrees were removed because of memory limitations. For the genetic map, we assumed that 1 Mb was equivalent to 1 cM. Although this is a simplistic assumption, when inter-marker distances are relatively short and when a dense marker map is used, the assumption has been shown to produce nearly identical linkage results as a more detailed genetic map [11]. We performed a parametric analysis using a dominant model, assuming a minor allele frequency of 0.01 and a penetrance of 0.008, 0.5, and 0.5 for carriers of none, one, and two disease alleles, respectively. Our parametric model results in an overall prevalence rate of 7.5 cases per 1000 individuals, which is consistent with published prevalence rates for rheumatoid arthritis [12].

Graphical illustrations of the linkage results are provided for McLink-LD and the LD elimination methods analyzed in MCLINK in Figures 1 and 2. Because the Merlin peak profiles were similar to the MCLINK output, we do not provide a graphical illustration of those results. The maximum TLOD scores from MCLINK and McLink-LD and the HLOD scores from Merlin are displayed in Table 2.

The linkage results after controlling for LD all consistently verified the conclusions of Amos et al. [6], that is that the chromosome 6 peak is a true-positive result and the chromosome 21 peak is a false-positive result. For chromosome 21, all LD control methods behaved similarly with peak values of ~1.0. However, larger peak differences across methods were observed for chromosome 6 (range: 9.07 to 16.07). Peak position also varied on chromosomes 6 and 21. For chromosome 6, all methods that removed high-LD SNPs peaked at 28.96 cM, whereas McLink-LD peaked at 38.61 cM. We note that the HLA-DRB1 gene, the best characterized single genetic risk factor contributing to rheumatoid arthritis, maps between two SNPs included in our data sets at 28.96 cM and 33.1 cM, which is the peak location obtained from analysis of the Complete data set and all LD elimination methods. Although here we reported TLOD results as this is currently the only output of McLink-LD, rheumatoid arthritis is a complex disease, and heterogeneity LOD (HLOD) scores may more accurately reflect the true position of a peak. However, it should be noted that the HLOD scores from Merlin were reported at identical (or nearly identical) peak position values for the three LD-elimination methods using MCLINK at both the chromosome 6 and 21 regions.

While a method that incorporates the underlying LD structure should more accurately reflect the true linkage peak height and position for a chromosomal region, it may be too early to consider McLink-LD the 'gold standard'. MCMC methods can be problematic due to computational requirements and poor mixing properties and it may not be suitable for use in all cases. We note that an additional feature of the new McLink-LD program is maximization of LOD scores across recombination fractions on the (0, 1) interval. The (0, 1) interval may better be able to incorporate model errors and improve mixing properties of the graphical LD model. Maximizing over the (0, 1) interval, we observed a peak TLOD score of 13.75 at 37.04 cM, which is more similar to the results for the high-LD SNP elimination methods. However, because the high-LD SNP elimination methods all used the traditional (0, 0.5) interval, the comparison is not ideal. We also note that maximization of LOD scores over the (0, 1) interval resulted in general symmetry of the LOD function about a recombination fraction of 0.5. This is most likely due to a heavy reliance on two-generation pedigrees in the GAW15 NARAC data resource with no phase information. Any asymmetry in the chromosome 6 and 21 regions may be due to inadequate modeling of parameters [21] or failure of the sampler to mix between phase states in the few three or more generation families.

The three high-LD elimination methods all peaked at the same position, but with different peak heights. We examined residual LD between all pairs of SNPs for each of the LD elimination methods. Using a threshold for high LD of D' > 0.7, the number of pairs exceeding the threshold was less than 0.7% compared to the total number of pairwise comparisons for each method across both chromosomes, suggesting that the majority of high-LD SNPs pairs were captured. However, we do note that for chromosome 6, although the percentages were small, the number of pairs of SNPs exceeding the threshold tended to follow the same pattern as the peak TLOD results (i.e., LD threshold pattern: PCA-genotype < PCA-haplotype < SNPLINK). Thus, differences in peak height most likely are due to residual LD in each of the methods. Although not investigated here, it is also possible that the PCA-genotype method, which is the most conservative, may have diminished power compared to the other two LD elimination methods as more SNPs were eliminated. The PCA methods provided a logical approach to elimination of SNPs that were not based on a contiguous LD structure or random elimination of SNPs from a high-LD set as does SNPLINK. Because an automated PCA package for removing high-LD SNPs is yet to be developed, we recommend that PCA be used only for follow-up of regions with at least suggestive linkage evidence.

Control for LD is an essential component of analyzing high-density SNP linkage data, and inflated linkage peaks may result if some method for controlling for LD is not implemented. Because we observed differences in linkage peak height and position across the four methods studied here, further work is needed to explore these and other LD control methods.