Density-based clustering in haplotype analysis for association mapping

General methods

All analyses were carried out with knowledge of the true location of susceptibility loci.

Marker names from the chromosome 6 dense SNP scan are abbreviated here such that "denseSNP6_N" will be denoted as "SNP N". We tested the markers flanking the HLA-DRB1 locus (DRB1, coincident with SNP 3437, 49.5 cM) and locus D (between SNPs 3916 and 3917, 54.6 cM) for redundancy using BEST [6], and removed one redundant marker, SNP 3434, from the region near DRB1. We explored patterns of linkage disequilibrium (LD) and assessed the significance of nonzero LD by the likelihood ratio test in HaploView version 3.32 [7].

Testing for dHWE was carried out using the exact test in HaploView, and, for confirmation, the exact and χ² tests as implemented in the R genetics library package, version 1.2.0. Analyses were performed on sets of 1500 cases – one affected sib chosen at random from each affected-sib pair (ASP) family – and 2000 unrelated controls. All cases and controls in each set were from the same replicate.

We used the ASSOC program in S.A.G.E. [8] for case-control tests of association. The transmission disequilibrium test (TDT) was performed on trios of parents and affected child as implemented in the S.A.G.E. program TDTEX. Trios were selected with one random offspring from all 1500 ASP families in each replicate. In addition, the generalized family-based approach implemented in FBAT version 1.7.2 [9, 10] was carried out in parallel on sets of 1500 complete ASP families, with a null hypothesis of no linkage or association.

Association mapping by linear regression with clustering of haplotypes

We have extended the regression-based approach for association testing of Schaid et al. [3] to incorporate haplotype groups via a density-based clustering algorithm [1]. The primary goal was to reduce the dimensionality of the regression by clumping together haplotypes that are likely to have diverged recently, whether through mutation or recombination. Posterior haplotype probabilities from unphased data are obtained from the Decipher program in S.A.G.E. [8] and are imported into a modified version of the HapMiner program [1] as haplotype weights for clustering. Each pair of haplotypes is assigned a similarity score, a generalization of several scores previously described [11], which is converted to a distance metric on the interval [0,1] [1]. Clusters are formed in regions of high density (haplotype weight). A haplotype is designated a "core" haplotype if enough density, determined by the density threshold MinPts, is located within a given distance ε from it. Haplotypes within this ε neighborhood are clustered together. We modified HapMiner such that very common haplotypes, defined by the parameter p_min, are never clustered together. This prevents improper grouping of ancient haplotypes. For the analyses presented here, we selected a value for p_min of 1/(2k), where k is the number of haplotypes present with a frequency large enough to include in the GLM (see below).

Cluster assignments for all possible haplotypes are imported into the haplo.score function in HaploStats [2, 3]. This method first estimates haplotype frequencies by the expectation maximization algorithm, and uses the frequencies to calculate posterior probabilities of haplotype pairs for each individual, assuming HWE. The posterior probabilities, in turn, are incorporated into a score test for association based on the likelihood function of a particular GLM – for case-control data, a logistic model – in which each haplotype is assigned a model coefficient [3]. A global score test for association is asymptotically distributed, under the null model of no association (all coefficients equal to 0), as a χ² random variable with degrees of freedom equal to the rank of the variance matrix for the score statistic. In our cluster-based approach, parameter estimates are obtained for clusters, rather than for haplotypes. We calculated the variance of the score statistic as per the generalized score test of Boos [12] as implemented by Tzeng et al. [13] because we found that variance calculation in Schaid et al. [3] based on the Louis information [14] inflated the type I error of the test when covariates were included in the analysis (data not shown). To prevent numerical instability and loss of power resulting from estimation of rare haplotypes [3], only haplotypes or clusters with frequencies above 0.002 were included.

Power and type I error analyses on simulated datasets

We carried out studies of power and type I error of association mapping methods on subsets of all 100 replicates. Cases were randomly selected from the offspring of ASP families, such that no sample contained both sibs from any family; controls were chosen at random from the unrelated controls. All individuals within a sample were taken from the same replicate. We adjusted the sample size for each analysis to yield moderate (40–60%) power from the haplotype-only test. Where necessary, we included sex and the number of DRB1*04 alleles in the model as covariates, to reduce the signal strength to a level useful for power comparisons; this adjustment was not part of the analysis. Type I error in association-positive regions was assessed by randomizing case/control status (and covariates, if any) relative to genotype data by permutation.

A second null-model analysis was performed at locations far removed from the HLA region and locus D. The mean haplotype diversity, measured as the number of unique haplotypes present in the phased data, for all sets of six and eight consecutive markers was estimated over replicates 1–50 in two large regions on chromosome 6q comprising SNPs 11001–12000 and 15001–16000. Four levels of diversity were defined: low (10th percentile of all marker sets), medium (50th), high (90th) and very high (99th). Samples were extracted at selected locations at each diversity level, and score tests for association were carried out as above (without permutation).

Deviation from HWE

We observed extensive dHWE in the neighborhood of HLA-DRB1 (Fig. 1). Forty-three SNPs within a region of 170 markers spanning from 31.1 Mb to 33.0 Mb on chromosome 6 showed significant dHWE, with p < 0.001 in at least five of the 100 replicates. Six consecutive markers flanking DRB1 gave p-values below 0.001 in every replicate and median p-values of less than 10^-6 from HaploView over all replicates. The three methods for testing HWE produced nearly identical results (data not shown). As expected, there was marked loss of heterozygosity in these markers. Because the score test for association assumes HWE, we took care to assess type I error for the two methods within this region (see below).

Single-marker association

As a preliminary step toward haplotype association analysis, we performed exploratory single-marker scans on case-control samples and trios. An initial scan with ASSOC on case-control samples gave evidence of an extremely strong genetic signal in the HLA region (data not shown). Evidence for association of most SNPs near DRB1 with RA was, overall, highly significant under both TDTEX and FBAT (Fig. 2). Nevertheless, a few SNPs interspersed among the markers in this region yielded p-values of unusually low significance by one or both methods. These results, which were consistent across replicates, could not be explained entirely by the informativity of the markers: the correlation between p-value and the number of families fully informative for the TDT (i.e., with both parents heterozygous) was low when the TDT had unusually low power. Specifically, of the eight markers with median p-values greater than 0.1, six were within the top three quintiles of the percentage of informative families (i.e., ≥36%). We expected a uniformly strong correlation between the available information on transmission and power of the TDT, and we were unable to explain this surprising result.

Power analyses at HLA-DRB1 and locus D

Ability to detect association at locus D was markedly reduced with adjustment for sex and number of DRB1*04 alleles, with power at the 0.01 significance level falling almost to background (Table 2). This observation strongly suggests that DRB1 and locus C are providing most of the genetic signal at locus D. Low but significant LD (|D'| < 0.1; LOD score > 2 for H₀: D' = 0) was observed between SNP 3437 at DRB1 and SNP 3917, 1.6 kb from locus D (data not shown). Given the overwhelming effect of DRB1 on RA, this small level of LD may explain the association between RA and haplotypes at locus D.

Assessment of type I error

False-positive rates measured at HLA-DRB1 in the absence of covariates matched expectation when trait values were permuted relative to genotypes, both with and without haplotype clustering, as did those measured at locus D. With adjustment for sex and DRB1*04 alleles, cluster-based score tests again yielded proper type I error, whereas the test was somewhat conservative without clustering, returning a false-positive rate of about 0.025 at a significance level of 0.05. Results were similar for all six- and eight-marker haplotypes examined in power analyses (data not shown). Thus, our cluster-based score test appears to be robust to the severe dHWE encountered within the HLA region.