There are generally three steps to isolate a disease linkage-susceptibility gene: genome-wide scan, fine mapping, and, last, positional cloning. The last step is time consuming and involves intensive laboratory work. In some cases, fine mapping cannot proceed further on a set of markers because they are tightly linked. For years, genetic statisticians have been trying different ways to narrow the fine-mapping results to provide some guidance for the next step of laboratory work. Although these methods are practical and efficient, most of them are based on IBD data, which usually can be inferred only from the genotype data with some uncertainty. The corresponding methods thus have no greater power than one using genotype data directly. Also, IBD-based methods apply only to relative pair data. Here, using genotype data, we have developed a statistical hypothesis-testing method to pinpoint a SNP, or SNPs, suspected of responsibility for a disease trait linkage among a set of SNPs tightly linked in a region. Our method uses genotype data of affected individuals or case-control studies, which are widely available in the laboratory. The testing statistic can be constructed using any genotype-based disease-marker disequilibrium measure and is asymptotically distributed as a chi-square mixture. This method can be used for singleton data, relative pair data, or general pedigree data. We have applied the method to simulated data as well as a real data set; it gives satisfactory results.

---

RECENTLY, genome-wide scans have been widely used in the study of complex genetic diseases such as cardiovascular diseases, obesity, diabetes, schizophrenia, etc., due to the advance in biological science that hundreds of markers could be genotyped quickly with reduced cost. Subsequent fine-mapping studies have been also frequently reported, which narrow the linkage region to a disease trait to about one or a few centimorgans. However, very few of the studies reach the final step of positional cloning to isolate the gene responsible for the linkage to a complex disease. Part of the reason is that the process involves genomic DNA spanning millions of base pairs at the linkage region, sequencing large amounts of the overlapped genomic DNA fragments, and genotyping tens or hundreds of markers in the region, which take intensive work in an ordinary laboratory. In some cases, fine mapping cannot proceed farther on a set of markers because they are tightly linked. For years, genetic statisticians have been trying to develop parametric and/or nonparametric methods to pinpoint the linkage to one or very few markers suspected to be truly responsible for the linkage of a disease trait and to exclude those only in linkage disequilibrium (LD) to the susceptibility markers.

Difficulty in the identification of specific disease-predisposing alleles may result due to multiple genetic factors (TAIT and HARRISON 1991; THOMSON 1991). GREENBERG (1993) and HODGE (1993) considered the analysis of "necessary" vs. "susceptibility" loci in which the associated marker allele itself increases disease susceptibility but is neither necessary nor sufficient for disease expression. The conditioning method is one of the typical statistical tools for studying such problems. FULKER et al. (1999) developed a conditioning method using the variance component model. This method tests both linkage and association at the same time, so that it provides the result whether a locus is the candidate locus to the trait or is just in LD with the candidate locus. This idea was further expanded by CARDON and ABECASIS (2000), in which a combined linkage and association method using the variance components model is proposed. VALDES and THOMSON (1997) and SIEGMUND et al. (2001) used the conditioning method to narrow down the association region. LAZZERONI and LANGE (1998) proposed such a framework in the transmission/disequilibrium test. Furthermore, SORIA et al. (2000) considered a conditioning argument to pinpoint the linkage of the G20210A mutation in the prothrombin gene to the disease gene. On the other hand, BLANGER et al. (2000) studied a Bayesian variance components method, and HORIKAWA et al. (2000) used a modified association study method, which identified a single-nucleotide polymorphism (SNP), SNP43, that showed significant association with the evidence for linkage with type 2 diabetes.

Recently, SUN et al. (2002) proposed a statistical method for this problem. They used a conditioning hypothesis-testing procedure to pinpoint, among a set of tightly trait-linked genes, a single or a few susceptible markers, using identity-by-descent (IBD) data from affected sibships. This method is based on the genome-wide scan result, which identified a region showing strong linkage with a putative trait. Often markers in such a region are tightly linked among themselves. The goal of the method is to identify which of those markers are truly responsible for the linkage and which are merely tightly linked to such putative markers. This method is practical in application and yielded good results in their simulation studies.

However, most of the existing methods for this problem use IBD data on paired family members. Usually IBD data are not fully available in practice and can be inferred only from genotype data with uncertainty and often inconsistently from different methods used. Inference based on them has no greater power than that based on genotype data, unless the IBD data are a sufficient statistic for the parameters underlying the model. Also, IBD-based methods apply only to relative pair data.

Here we present a method for this problem by formulating a set of conditional hypothesis testing, in this respect similar to that in SUN et al. (2002), but we use genotype data instead and the testing statistic is different in nature from theirs. Using any genotype-based trait-marker disequilibrium measure, the testing statistics are constructed by successively conditioning on each of the tightly linked SNP sites. Our method is nonparametric: it does not require model specification or phase information in the data. It applies to family data of arbitrary structure, including singleton data, in which each individual comes from a different and independent family. Under the null hypothesis of being the sole susceptible site, each of these statistics follows asymptotically a chi-square mixture distribution. The corresponding P values are easily obtained via simulation.

The data:

Let A be the unknown disease allele, for which we want to infer its position in the human genome. Assume that there are J identified SNP markers, M_j (j = 1, ... , J), with alleles M_jk (k = 1, 2), which are brought to our attention because of their tight linkage to the disease allele. A natural question is whether all of them are susceptible genes for the linkage or some of them show disease linkage just because of their strong linkage with the true susceptible gene(s). Our goal here is to identify the true susceptible SNP(s), if any, among them.

For ease of explanation we first describe our method for singleton data and then extend it to general pedigree data in a later section. We now describe a general procedure for the conditional inference of this problem; the construction of the specific testing statistic is detailed later. Let G = (G₁, ... , G_J) be a general notation for the composite SNP genotype at all the SNP loci, where G_j = (g_j1, g_j2) is its allelic notation; G_nj = (g_nj1, g_nj2) be the observed genotype of the nth individual at the jth SNP locus (n = 1, ... , N; j = 1, ... , J); and G_n = (G_n1, ... , G_nJ) be the vector notation of the observed composite genotype of individual n. The data to be used are G₁, ... , G_N, the observed composite genotypes of N individuals at J SNP loci each.

Here we assumed the common practice that at each SNP locus there are two different alleles in the population; we code them as 1 and 2, although the same value from alleles at different loci may have different allelic meaning. At each locus, we code the genotype as G_nj = 0 when g_nj1 g_nj2, G_nj = I when (g_nj1, g_nj2) = (1, 1), and G_nj = II when (g_nj1, g_nj2) = (2, 2). Note that we have two representations of a SNP genotype, one allelic and one numerical. Which one(s) will be used, even in the same expression, depends on convenience.

The disequilibrium measure and the conditioning principle:

The proposed conditional testing procedure uses testing statistics, which are constructed via the conditional version of any trait-marker LD measure using genotype data. We first state the conditional testing principle and then give the specific forms of the testing statistic for some particular data designs.

Now we describe the trait-marker LD measure. Let p_A be the population frequency of the disease allele A, q_jk be those of allele k of marker j, and P_A,jk be those of the haplotype (A, M_jk). Let D_A,jk = P_A,jk – p_Aq_jk be the LD measure between the disease allele A and allele k of marker j. Since the position of A is unknown, p_A, P_A,jk, and thus D_A,jk cannot be directly estimated from the observed data; instead various quantities are constructed to infer it.

When D_A,jk is positive, the marker allele M_jk is more likely to be associated with the disease-susceptible allele A than would be expected by chance. The disequilibrium measures D_A,jk are among the main tools for finding the association between a marker locus (loci) and the disease locus. There are numerous ways to construct inference statistic from the D_A,jk's, some using relative pair IBD data at markers, and some using marker genotype data (BENGTSSON and THOMSON 1981; LEHESJOKI et al. 1993; FEDER et al. 1996; NIELSEN et al. 1998). Here we develop the conditional version of the genotype-based method.

Let q_jk|i be the population frequency of allele k of the jth SNP conditional on the ith SNP genotype (k = 1, 2). Let P_A,jk|i be the population frequency of the disease-SNP haplotype (A, M_jk) at the jth marker locus conditional on the ith SNP genotype and P_jr|i be that of the homozygote SNP genotype r at the jth SNP locus conditional on the ith SNP genotype (r = I, II). We choose the conditional LD measure at the jth locus, given the ith SNP genotype, as

(1)

Note that P_Aj2|i – p_Aq_j2|i = –(P_Aj1|i – p_Aq_j1|i), so only one of the marker alleles is needed to define this disequilibrium. Our motivation to use the conditional LD measure is that if marker i is the sole susceptible site of linkage to the disease allele, then the genotype data from this site constitute a sufficient statistic for this measure, or, in other words, it will explain all the disequilibria in the region. Thus, conditioning on the site of interest, the disequilibria parameters D_j|i vanish from the conditional distribution of the data, for all j i.

In the following we explain what the conditioning actually means in practice. Suppose we have genotype data for 502 individuals on two SNP loci, each locus has two alleles, and each allele takes one of the two forms that we coded as 1 and 2. The genotype at each locus is thus represented as (1, 1), (1, 2) = (2, 1) and (2, 2). The supposed observed genotype frequencies for the two loci are given in Table 1.

View this table: In this window In a new window   TABLE 1 Genotype frequencies at two loci

 
By conditioning on the genotype at the first locus being (1, 1), we mean the subgroup of 156 individuals whose genotype at the first locus is (1, 1). Within this subgroup, the genotype at the second locus is denoted as locus 2/(1, 1) and similarly for conditioning on the first locus genotype being (1, 2) and (2, 2). Thus conditioning on the first locus genotype (1, 1), (1, 2), and (2, 2) separately, the data are divided into three nonoverlap subdata sets, and we obtain the genotype frequency of the second locus as shown in Table 2. Likewise, conditioning on the second locus genotypes separately, we get the genotype frequency of the first locus as shown in Table 3.

View this table: In this window In a new window   TABLE 2 Genotype frequencies at locus 2 conditioning on locus 1

 

View this table: In this window In a new window   TABLE 3 Genotype frequencies at locus 1 conditioning on locus 2

 
In conditional testing, test statistics are constructed with the data in Tables 2 and 3. For example, to test the hypothesis that locus 1 is the only susceptible site, then conditioning on it we obtain three subtables. If the hypothesis is true, the LD vanish on each of the subtables, and the test statistics constructed from them should manifest nonsignificance.

The hypotheses and testing statistics:

We are interested in testing the null hypothesis H_i: among the set of markers, SNP marker i is the sole cause of the linkage to the disease locus. Here we assume background effects on the linkage are negligible; see the DISCUSSION for more details on this. Under H_i, D_j|i = 0 for all j i. For each fixed i, testing statistics S_j|i are constructed, usually a function of the empirical version _j|i (j i), such that they tend to be small under H_i and large otherwise. H_i can be decomposed as H_i = H_i1 H_i2, where H_ik is the hypothesis: genotype k at site i is the sole susceptibility SNP reasonable for the trait LD in the region. Note that when H_ik is rejected, we can conclude only that the SNP does not contribute to the LD in the region, that the single or multiple causal polymorphism may not be among those that are typed, or that there is more than one source of such contribution. By testing the sequence of {H_ik}, we can find a confidence set, which may consist of a single SNP or several SNPs, or it may be empty. This set may be more accurately inferred by testing the hypothesis of multiple SNPs as in a later subsection. Our method can be used to detect a more detailed local relationship by testing the more detailed hypothesis H_j|i, LD at site j is completely caused by site i, or even the finer hypothesis H_j|ik, LD at site j is completely caused by genotype k of site i.

These last hypotheses are inferred using the statistics S_j|ik, which are the corresponding versions of the S_j|i's for the H_j|ik's. For recessive disease, the conditional statistic notation S_j|ik means . The S_j|ik's are constructed of the form , and the random column vector is jointly asymptotic normal under H_ik. Let _ik be the asymptotic variance matrix of X_ik and = (₁, ... , _J–1) be its eigenvalues. Usually _ik and thus can be estimated by their empirical version. The particular forms of the S_j|ik's are given later for different data designs.

Asymptotic distribution of the testing statistic:

Let us consider H_ik. Its testing statistic is given by

To get the asymptotic distribution of S_+|ik() or S_+|ik under H_ik, we first give a general result for the distribution of quadratic form of normal random variables. The proof is given in the APPENDIX.

Proposition:

Let X = (X₁, ... , X_d)' be a nondegenerate normal random vector: X N(0, ) (i.e., || 0), with eigenvalues = (₁, ... , _d); A is a d-dimensional positive definite symmetric matrix with eigenvalues = (₁, ... , _d); the _i's and the _i's keep the same order in the diagonalization. We have

1. 
   

   where the Y²_j's are independent and identically distributed (IID) ²₁ random variables.
   

2. Let = diag(₁, ... , _d) and = (₁, ... , _d); then
   

   Especially, when A = I_d, we have
   

   
   

Remark:

1. The case or A being degenerate is not of much interest and can be avoided easily in the construction of the testing statistic.
   
2. It requires that and be of the same order; this can be done using the same orthogonal matrix (matrices) in the diagonalization of and A. More conveniently, since it actually used only the _jj's, they are just the eigenvalues of A (or A).
   
3. Using i or ii is a matter of choice. i is simpler in forming the ² statistic but not in computing the quantiles or P values, while the order of the _jj's does not matter. ii involves computing A^1/2 in forming the ² statistic, and the order of _jj and that of X_j must match. In practice, this is not trivial; however, it is simpler in computing the quantiles or P values using the existing ² tables.
   
4. Given and , the density of ²_d can be derived by the multiple convolution formula, and thus its th quantile and/or the P value of the observed statistic can be obtained. But, more conveniently, for a given level , the th quantile and/or the P value of the observed statistic can be consistently estimated by their empirical versions.
   

To sample from ²_d, we sample Y²₁, ... , Y²_J–1 from ²₁ independently; then ₁₁Y²₁ + ... + _ddY²_d is a sample from ²_d.

The ² linear combination is the general form of the quadratic form of normals. When the X_j's are independent, _j = Var(X_j); when the X_j's are IID and . There are some other similar results about the quadratic form of normals (GRAYBILL and MARSAGLIA 1957; GOOD 1969; KHATRI 1980, 1982; ANDERSON and STYAN 1982). Our result is independent and not of the same formulations and conditions as the others.

Let the eigenvalues (in their original order) of _ik be = (₁, ... , _J–1); by ii of the Proposition, we have (see APPENDIX):

Corollary:

Under H_ik, asymptotically

and

where the Y²_j's are IID ²₁ random variables.

Thus for given 0 < < 1, the asymptotic level test of H_ik is given by the rejection rule: the P value of the observed S_+|ik is smaller than , or S_+|ik > Q_J–1(, ), the th quantile of the ²_J–1 distribution.

Note that our method requires only the genotype information and allele counts at each locus. It does not require phase information in diploids; thus it is practical in applications.

In the following we give the specific forms of the S_+|ik's [S_+|ik()'s] under some commonly used settings; those of the S_+|ik()'s are the same and are omitted.

Multiple susceptible loci:

Our method can be extended to the case of multiple susceptible loci without conceptual difficulty, but with more involved computations. Consider the hypothesis H_{i1k1,...,irkr} (1 r < J) that the composite genotypes (k₁, ... , k_r) at loci (i₁, ... , i_r) are the true susceptible ones. The corresponding testing statistics S_{j|i1k1,...,irkr} are constructed similarly as before. The only difference is now the inference set, the conditional data set, consisting of those individuals whose alleles at loci (i₁, ... , i_r) are (k₁, ... , k_r), and

which is asymptotically ²_J–r, and = (₁, ... , _J–r) is the eigenvalue of the asymptotic variance matrix _{i1k1,...,irkr}, which is estimated the same way as the single susceptible locus case, but uses the current inference data set.

For fixed r, there are J!(J – r)!/r! of such tests across different choices of loci combinations, and 2^r of such tests for each choice of loci combination. So the total number of tests will be 2^rJ!(J – r)!/r!.

Note that the above construction of the testing statistic is general; its inference behavior depends on the particular statistic used. The general form of the testing statistic is asymptotically a chi-square mixture, which is centralized under H_ik and noncentralized otherwise. The functional form of the parameters of interest entering the noncentrality parameter in the chi-square mixture will explain the behavior of the test in terms of asymptotic power. We give more detail on this for specific tests used in the following sections.

Now we explain how to construct the S_+|ik's in this type of data. In the case J = 1, assume the two SNP alleles are M and , and let A be the disease allele. Let p_A, q_M, and P_AM be the population frequency of the alleles A and M and haplotype AM, respectively, and let D_AM = P_AM – p_Aq_M be the LD. For clarity we first assume the disease is recessive and P(Affected|AA) = 1. Under these assumptions, FEDER et al. (1996) and more specifically NIELSEN et al. (1998) discovered the relationship

where is the probability that an individual will exhibit the disease due to causes other than this locus, and is the prevalence of the disease in the population. This equality enables us to detect the marker-disease association by testing Hardy-Weinberg disequilibrium at the marker locus without using IBD information. In fact the connection between the marker allele frequencies and the marker-disease LD is kept if we use only the numerator in the above equality, and this will simplify the computation. That is,

(2)

We derive a conditional version of (2) to serve our purpose.

Since all individuals are affected in this study, we drop off the index "Affected" to simplify the notations. We want to test the hypothesis H_ik: SNP type k at locus i is the sole cause of the LD in the region. Let P_jr|ik be the population frequency of genotype r (r = I, II) of locus j given one's genotype being k at locus i, q_jr|ik be that of allele r (r = 1, 2) at locus j given one's genotype being k at locus i, _j be the probability that an individual will exhibit the disease due to causes other than locus j, and D_j|ik be the disequilibrium corresponding to the conditional LD measure. Now the same derivation of (2) leads to

(3)

Under H_ik, all association of SNP j is completely explained by genotype k of locus i; thus D_j|ik = 0 and hence T_j|ik = 0 (j i).

We comment that our method works for a general disease model; in this case T_j|ik is still a function of D_j|ik but the expression is more involved (see NIELSEN et al. 1998, pp. 1533–1534), and under H_ik we still have T_j|ik = 0 (j i); hence the test is still valid. In this case, the power and error rate computation will be more involved. The same comment applies to the case-control section also.

Now we construct testing statistics for H_ik (i = 1, ... , J). The consistent estimates _jr|ik of P_jr|ik and _jr|ik of q_jr|ik are given by

where is the total number of individuals with the ith SNP genotype being k, and we rearrange them as the first, second, ... , and the N_ikth individual. I_n,jr|ik,(= 0, 1) is the indicator that the nth individual among this set has genotype type r on the jth locus given he (she) has genotype k at locus i, and

where J_n,j1|ik(= 0, 1, 2) is, for the nth individual, the number of times allele 1 occurs at locus j, given one's genotype being k at locus i. The estimate of T_j|ik is

Let _ik = (_j|ik : j i) be the (J – 1) dimensional column vector. Under H_ik, _ik is asymptotically N(0, _ik) for some matrix _ik to be identified later. Let = (₁, ... , _J–1) be all the eigenvalues of _ik, and . By the Corollary, under H_ik asymptotically

and _ik is estimated by

(4)

(APPENDIX), where

and

Here means matrix direct summation, which results in a (J – 1) x 2(J – 1) dimensional matrix. From _ik, we obtain the estimated eigenvalues .

Similarly, for H_i, let N_i = N_i1 + N_i2, _r = N_ir/N_i,

and _i = (_j|i : j i). Let _i be the asymptotic matrix of _i and = (₁, ... , _2(J–1)) be all the eigenvalues of _i. Note that _i = _i1 + _i2, and _i1 and _i2 are independent, so under H_i, _i is asymptotically N(0, _i), with . Its estimate is obtained as , and _ir is constructed as before.

Let

Under H_i, S_+|i ²_J–1.

We remark that in the above the asymptotic variance matrices _jk are estimated the same way as for the IID data. In general the familial data are not IID, and the above variance matrices are dealt with differently. Usually, in the positive dependent case, the asymptotic variance matrix will be larger, in the sense of generalized variance—the determinant of the variance matrix—and consequently will tend to have larger eigenvalues than the IID case, such as the singleton data case. In the case of homogeneous familial structure, more accurate estimates can be obtained. We study the above methods for general pedigree data in the extension section later.

In some of the existing methods for this problem, e.g., SUN et al. (2002), the conditional IBD sharing statistics are computed at each site given the genotype at that site. In this way the statistic can test whether each of the sites is the sole susceptible site, but will not be able to find the more detailed relationship between sites when the null hypothesis of only one susceptible site is rejected, while our test statistic can be used to reveal more detailed relationship. If H_ik is accepted, it is reasonable to say that the connection between site j and the disease locus is due to genotype k of site i.

By the asymptotic normality of _ik and (3), when H_ik is false, S_+|ik will be asymptotically a noncentral ²_J–1 (, µ), with noncentrality parameter

It is clear that H_ik is true if and only if D_j|ik = 0 (j i). In terms of µ, the null hypothesis is rephrased as H_ik: µ = 0. For a given level (= P(reject H_ik|H_ik is true)), the parameters , _js, , and the D_j|ik's, the asymptotic power of the test is

Here Q_J–1(, ) is the th quantile of the noncentral ²_J–1 distribution, which can be simulated by the sampling method after the Remark of the Proposition, but with Y₁, ... , Y_J–1 independent, and Y_j from N(µ_j, 1) with .

For this particular test statistic, since the power is an increasing function of µ, H_ik will be more accurately rejected when the _j(1 – )_j's and the conditional D_j|ik's are large, and _j is small or the disease is relatively rare. Likewise, H_ik will be more correctly accepted when the _j(1 – _j)'s and the conditional D_j|ik's are small (i.e., mainly explained by allele k of locus i), and the disease is relatively common.

Likewise, the error rate, the probability of false acceptance, is

Let q_M|A and q_M|U denote marker M population frequencies for the affected (case) and unaffected (control) individuals. BENGTSSON and THOMSON (1981) and LEHESJOKI et al. (1993) gave the following LD measure:

The conditional version of the above is

Let N_A and N_U be the number of affected and unaffected individuals, and

where

is the total number of "affected" individuals with the ith SNP genotype being k,

where J^A_n,jr|ik and J^U_n,jr|ik are the same as the J_n,jr|ik before, but here for affected and unaffected individuals. is the total number of unaffected individuals whose ith SNP genotype is k. Let N_ik = N_A,ik + N_U,ik. Assume N_A,ik/N_{ik A,ik} and N_U,ik/N_{ik U,ik} = 1 – _A,ik. To test H_ik, let

Under H_ik, asymptotically

where is the vector of eigenvalues of the matrix _ik. Let

Then _ik is estimated by

(5)

(APPENDIX), where, for singleton data, the affected and the unaffected are independent, so ,

and

Similarly, to test H_i, let

Then under H_i, asymptotically S_+|i ²_J–1, and is the vector of eigenvalues of _i, which is estimated by

(6)

(APPENDIX), where, for singleton data, ,

where , N_i = N_i1 + N_i2, N_A,i = N_A,i1 + N_A,i2, N_U,i = N_U,i1 + N_U,i2, _A,i = N_A,i/N_i, and _U,i = N_U,i/N_i = 1 – _A,i. Similarly,

Other LD measures can also be used, for example, the trend test statistic (ARMITAGE 1955; DEVLIN and ROEDER 1999).

As in the affected individual case, when H_ik is not true, S_+|ik is asymptotically noncentral ²_J–1, where

Given , , _j, , p_A, the q_j2's, and the D_j1|ik's, the power and error rate can be computed by simulation as before, but with Y₁, ... , Y_J–1 independent, with Y_j from N(µ_j, 1), where

Here, the power and probability of correct acceptance of H_ik depend on , , p_A, the q_j2, and the D_j1|ik's. The power is maximum when the conditional D_j1|ik's are maximum, and the test is more likely to accept H_ik when the D_j1|ik's are small. Their relationships with the other parameters can be analyzed similarly.

As mentioned earlier, the only difference in our methods between general pedigree data and the singleton data is the estimations of the corresponding asymptotic variance matrices. A simple method for this purpose can be found in the work of G. E. BONNEY, V. APPREY and A. YUAN (unpublished data), without any assumption on the data and no extra parameters introduced for the dependence. We illustrate this with the affected familial data, which for the case-control family data is similar. For such data, the estimations for the genotype/allele frequencies in the previous sections are not IID averages; we rewrite them as IID versions, so that their asymptotic variance matrices can be computed easily. First we assume the data have the same familial structure. Suppose there are M families with S individuals each (N = MS). We redefine _jr|ik as

where is the total number of families in which at least one individual with SNP type k at locus i, I_ik(s, m), is the indicator that in the mth family, there are s individuals with SNP type k at locus i. I_jr|ik (s, m) is the indicator that there are s individuals in family m with SNP type r on the jth locus, given the family is in the group with SNP type k on the ith locus. Let . Then for fixed (jr, ik), {I_jr|ik(m) : m = 1, ... , M} is an IID sequence, and for different (jr, ik) and (j'r', i'k'), {I_jr|ik(m) : m = 1, ... , M} and {I_j'r'|i'k'(m) : m = 1, ... , M} are independent. Similarly, _jr|ik is redefined as

where J_jr|ik(s, m) is the count that there are s SNP allele r in family m on the jth locus, and their SNP type is k on the ith locus. Let . Then for fixed (jr, ik), {J_jr|ik(m) : m = 1, ... , M} is an IID sequence, and for different (jr, ik) and (j'r', i'k'), {J_jr|ik(m) : m = 1, ... , M} and {J_j'r'|i'k'(m) : m = 1, ... , M} are independent.

Let _j|ik and _ik be as before but with _jI|ik, _jII|ik, and _j1|ik replaced by the above versions. Let . Now it is clear that the in (4) can be replaced by the consistent estimator for this case as

where

and _ik = ', and is the same as in (4).

More generally, suppose that there are L different familial structures in the data set, with size M_l each, and the lth structure has S_l individuals per family (l = 1, ... , L). Let

where is the total number of families with the structure l in which at least one individual with SNP type k at locus i, I_ik and I_jr|ik, is the counterpart of I_ik(s, m) and jr|ik(s, m), respectively, for familial structure l. Let . Then for fixed (l, jr, ik), is an IID sequence, and for different (l, jr, ik) and (l', j'r', i'k'), and are independent. Let define the estimate of P_jr|ik as

Similarly, let

where J_jr|ik is the counterpart of J_jr|ik(s, m) for familial structure l. Let , and define

Now for this general pedigree data, let _j|ik and _ik be as before but with _j1|ik, _j2|ik, and _j1|ik replaced by the above versions. Let , and we assume , then a consistent estimate of _ik is given by

(7)

(APPENDIX), where

and

For the test of H_i, or the case of case-control data, testing statistics and the corresponding asymptotic variance matrices can be obtained in a similar way; we omit the details here.

Here we use simulated data to illustrate our method. To exhibit the applicability of our method, we use singleton data, which is out of the scope of the IBD-based methods. We simulate the data G₁, ... , G_N, where G_n = (G_n1, ... , G_nJ) (n = 1, ... , N) and G_nj = (g_nj1, g_nj2), the two alleles at SNP site j for the nth individual. The g_njk's are coded as 1, 2 for its possible two alleles. We assume phase is known to simplify the simulation process, so that for each n, the two haplotypes (g_n11, ... , g_nJ1) and (g_n12, ... , g_nJ2) are independent. In this example, we take J = 6, so all the vectors G_n = (G_n1, ... , G_n6) are random samples from the population genotype S = (S₁, ... , S₆), and S_j = (s_j1, s_j2) is the genotype at the jth site. We assume genotype (1, 1) at the third SNP site is responsible for all the LD with the disease allele A; the other first alleles, s_j1(j 3), in this region are tightly linked to s₃₁.

Now the haplotypes S⁽¹⁾ = (s₁₁, ... , s₆₁) and S⁽²⁾ = (s₁₂, ... , s₆₂) are independent and the s_j2's are independent within themselves. Denote and as the two haplotypes of the nth individual. To sample such data, for each n we need only to sample G_n from S⁽¹⁾ and G_n from S⁽²⁾ independently. Let q_A = 0.8 be the frequency of the disease allele allele A = 1 among the affected individuals, q⁽¹⁾ = (q₁₁, ... , q₆₁) be the frequencies of S⁽¹⁾ = (1, ... , 1), and q⁽²⁾ = (q₁₂, ... , q₆₂) be that of S⁽²⁾ = (1, ... , 1). To sample from S⁽²⁾ is trivial; i.e., just sample g_nj2 independently from B(q_j2), the Bernoulli distribution with probability q_j2 of getting 1 and probability 1 – q_j2 of getting 0. To sample G_n, we need to sample from a joint Bernoulli distribution with probability q⁽¹⁾. Such a joint distribution can be specified in the form

(COX 1972; FITZMAURICE and LAIRD 1993), where and are parameters and exp {–A(, )} is the normalizing constant and W is all the cross-product terms of S⁽¹⁾, including all the second- and higher-order terms. This distribution can be sampled using the Gibbs sampler (GEMAN and GEMAN 1984). But the specification of the joint Bernoulli distribution has some subjectivity and the sampling scheme is not simple. Instead, we use a normal discretization method to sample it. We use high correlation for linkage. Let be the corresponding correlation matrix of the J + 1-dimensional normal distribution for (A, S⁽¹⁾),

Note that this matrix corresponds to a strong connection between A and s₃₁, but not between A and (s₁₁, s₂₁, s₃₁, s₄₁, s₅₁, s₆₁); it also corresponds to a strong connection between s₃₁ and (s₁₁, s₂₁, s₃₁, s₄₁, s₅₁, s₆₁). Thus all the loci have apparent linkage with the disease allele A.

To sample the composite genotypes from the above distribution, let X = (x_A, x₁, ... , x₆) be a sample from the normal distribution N(0, ); if x_j < ^–1(q_j1), we assign g_nj1 = 1; otherwise g_nj1 = 0, (j = 1, ... , 6), where ^–1(q) is the qth quantile of the standard normal distribution. Since q₃₁ is the proportion of allele 1, at locus 3, which is linked to the disease allele, in the affected population, the two alleles at locus 3 are in Hardy-Weinberg disequilibrium. The disease is recessive. We make the corresponding conditional probability P(s₃₂ = 1|s₃₁ = 1) high, say 0.8, among the affected individuals. In the simulation, we used a high frequency of q_j1 = q = 0.1, 0.2, ... , 0.9, (j 3) for allele 1 at each locus, to see how this affects the results.

By the same way we simulated control data, in which the two haplotyes are sampled the same way as G_n above. Together with the previous affected data we have case-control data, and the analysis is displayed in Table 6.

View this table: In this window In a new window   TABLE 6 Case-control data: values of observed S+|j1 (P value)

 
Specifically, the sampling scheme has the following three steps:

For each n = 1, ... , N, (N = 1000):

1. Draw a sample X = (x_A, x₁, ... , x_J) from the normal distribution N(0, ); if x_j < ^–1(q_j1), we assign g_nj1 = 1; otherwise g_nj1 = 0₃ (j = 1, ... , 6). Then we get the sample G⁽¹⁾ = (g_n11, ... , g_nJ1).
   
2. If g_n31 = 1, set q₃₂ = P(s₃₂ = 1|s₃₁ = 1) = 0.8, else q₃₂ = 0.1. For each j = 1, ... , J, draw X from U(0, 1), the uniform distribution on [0, 1]; if X < q_j2 assign g_nj2 = 1; otherwise assign g_nj2 = 0. Then we get a sample G⁽²⁾ = (g_n12, ... , g_nJ2).
   
3. G_n = (G⁽¹⁾, G⁽²⁾) is a sample from S.
   

When the two alleles at each locus are in Hardy-Weinberg disequilibrium, we use a two-dimensional normal with mean (0, 0) and variance matrix = (1, r; r, 1) with r = 0.2 to model their dependence. For each n, we first get the sample G⁽¹⁾ = (g_n11, ... , g_nJ1) from (x₁, ... , x_J) as before, then for each j = 1, ... , J separately, sample y_j from the conditional distribution N(rx_j, 1 – r²). If y_j < ^–1(q_j2), assign g_nj2 = 1, otherwise 0.

To simulate the case-control data, we choose q = 0.6 for the case and q = 0.25 for the control.

Simulated data:

We constructed the test statistics S_+|ik, (i = 1, ... , J; k = 1, 2) and computed the corresponding eigenvalues = (₁, ... , _J–1), using the method described in the Remark after the Proposition to compute the ² P value under the null hypotheses. Since in the simulation the sole linkage with the disease allele comes from s₃₁, we expect H₃₁ will be accepted, and the other H_jk's will be rejected. Table 4 is a summary of the observed values of the S_+|j1's for the H_jk's, for different choices of q, with corresponding P values in parentheses. We simulated and computed data for q = 0.1, 0.2, ... , 0.9; we display only part of them to save space.

View this table: In this window In a new window   TABLE 4 Affected individual data: values of observed S+|j1 (P value) for different q

 
For each testing statistic S_+|jk, there is a set of nonnegative eigenvalues = (₁, ... , _J–1). Their magnitude plays an important role in determining the asymptotic P value of the observed S_+|jk. For a given observed value of S_+|jk and fixed number of loci J, a roughly larger eigenvalue total || (defined as ₁ + ... + _J–1) results in a larger P value, and vice versa. Although for two sets of eigenvalues ₁ = (₁₁, ... , _1,J–1) and ₂ = (₂₁, ... , _2,J–1), even if |₁| = |₂|, the corresponding distributions ²(₁) and ²(₂) may not be equal, and they are equal if and only if ⁽¹⁾ = ⁽²⁾, where ^(k) = (_k,(1), ... , _k,(J–1)) is the ordered version of _k (k = 1, 2).

We display in the following the eigenvalues _j = (_j1, ... , _j5) for the S_+|j1's, for the case q = 0.7.

and

We find that in most cases the P values of S_+|31 suggest acceptance of H₃₁ with high confidence, and those for S_+|j1 (j 3), suggest rejection of H_j1, except for the case q = 0.9, in which the P values of S_+|51 and S_+|61 are also significant, along with that of S_+|31. We regard this last case as exceptional, in which the over-high proportion of allele 1 at each locus blurred the identifiability of the problem (think of the extreme case of q 1; the corresponding locus contributes nearly no information for the problem). Thus, in all these cases, the true hypothesis H₃₁ is accepted with high confidence, and the other false ones, H_j1, are rejected; i.e., the true disease-linkage-related allele 1 at locus 3 is correctly identified among all six loci that are all in LD with the disease locus.

To investigate the influence of the deviation from Hardy-Weinberg on our method, we simulated the data for this case, in which we use the allelic correlation r 0 at each locus for the deviation from Hardy-Weinberg equilibrium (HWE). The disease allele population frequency is fixed at q = 0.7 and the results are displayed in Table 5.

View this table: In this window In a new window   TABLE 5 Affected individual data: values of observed S+|j1 (P value) for different r

 
In the non-HWE case, it seems that the true picture becomes more difficult to recover as the deviation from HWE increases. In general, significant departures from HWE are not expected, but if observed, caution should be taken in applying this method (if genotyping error is present, for example). In particular, in situations in which nonrandom mating is a known confounder because of inbreeding or population structure, care should be exercised.

For the case-control data, we used q = 0.6 for the case and q = 0.25 for the control; HWE is assumed, and again locus 3 is the only connection to the disease allele. The results are shown in Table 6. It is seen that again, for the case-control data SNP locus 3 is correctly identified, and all the other loci are rejected as sources of cause for LD in the region.

The following is a tabulation of power of the test for the above simulated data, using the above and some combinations of , _j = (j i), , and D_j|31 = d(j 3). To get a sense of the power behavior of our methods, we choose J = 6, = ₁ as shown before. The noncentrality parameter µ involves 2J – 1 parameters in it. It is impractical to investigate and tabulate the influence of each of the 2J – 1 parameters to the power. Instead, we investigate the influence of µ to the power, with the given genetic model. Each given value of µ, corresponding to a 2(J – 1)-dimensional parameter subspace, is given by the formula for µ. Table 7 shows the display of power for both the affected individual data and the case-control data, for some choices of the level and the parameter µ. We comment that for the above specification of the parameter µ, the power of the tests for both the affected individual data and that for the case-control data are the same.

View this table: In this window In a new window   TABLE 7 Power for some given parameter values

 
Since the µ in the power of the test for affected individual data and that for the case-control data have different expressions, more detailed power computation can be obtained by the specification in terms of all the parameters involved.

Application to real data:

Non-insulin-dependent diabetes mellitus-1 data:

We first apply our method to the non-insulin-dependent diabetes mellitus-1 (NIDDM1) data used in SUN et al. (2002) and list our results along with theirs in Table 8. We see that, for these data, the two methods yield quite different, although not contradictory, results. With the method of Sun et al., loci 2 and 12 are most likely responsible for the LD, while by our method, loci 2, 4, 6, 7, 8, 9, 11, 12, 13, 14, 16, 17, 18, 19, 20, and 22 all likely contribute to the LD in the region. One possibility for the difference of the two methods might be that the calpain-10 region has some patterns of LD that are not understood—violating one of the assumptions of the methods. Since the truth in the data is unknown, we do not comment on the performances of the two methods on these data. It is not uncommon in the hypothesis test context, even for methods based on the same type of data, that different methods may have different results, even contradictory ones. In principle, methods using genotype data have no less power in inference than those using IBD data. Here it is too early to comment on the pros and cons for the two types of methods. A formal assessment may involve long-term and large-scale studies. At least our method provides the user more options and a flexible tool for this problem. Also, more methods will give us more strength in the inference. If the methods give consistent results, this will strengthen our confidence in decision; if they do not or are contradictory, the problem may need further investigation. We may perform the hypothesis tests on the current confidence set and continue this way to get a final confidence set of SNPs, in which all of them are accepted as possible sources of LD in the region. We do not pursue this in detail here because of space limitation.

View this table: In this window In a new window   TABLE 8 Results from the NIDDM1 data

 

Diabetes data:

Next, in a diabetes study, 280 individuals with type 2 diabetes were genotyped at a large number of SNP sites. First we find those SNPs with strong linkage to the trait and then use our method to identify the susceptible one(s). We use the measure of NIELSEN et al. (1998) to detect the marker-disease association, which is given by

where _ij|Affected and _i|Affected are the estimated frequencies of marker genotype A_iA_j and allele A_i from the observed affected individuals and m is the total number of alleles. They showed that this marker Hardy-Weinberg disequilibrium measure is proportional to the square of the disease-marker LD measure. Under the null hypothesis that there is no disease-marker LD, ²_HW is approximately distributed as a ² variable with degrees of freedom m(m – 1)/2.

After computing the value of ²_HW at each marker and their corresponding P values, we found that 13 of the markers significantly indicate strong evidence of disease-marker disequilibrium. To apply our method, we choose a set of six SNPs, and we code them as sites 1–6 for simplicity. The ²_HW values are displayed in Table 9, along with their P values in parentheses.

View this table: In this window In a new window   TABLE 9 Values of observed 2HW

 
We see from this table that all six loci are very tightly linked to the trait. Now we use our method to identify which one of the six SNPs is the sole true cause of linkage, if any. The computed values of the conditional testing statistics and their P values are in