Adjusting heterogeneous ascertainment bias for genetic association analysis with extended families

Background In family-based association analysis, each family is typically ascertained from a single proband, which renders the effects of ascertainment bias heterogeneous among family members. This is contrary to case–control studies, and may introduce sample or ascertainment bias. Statistical efficiency is affected by ascertainment bias, and careful adjustment can lead to substantial improvements in statistical power. However, genetic association analysis has often been conducted using family-based designs, without addressing the fact that each proband in a family has had a great influence on the probability for each family member to be affected. Method We propose a powerful and efficient statistic for genetic association analysis that considered the heterogeneity of ascertainment bias among family members, under the assumption that both prevalence and heritability of disease are available. With extensive simulation studies, we showed that the proposed method performed better than the existing methods, particularly for diseases with large heritability. Results We applied the proposed method to the genome-wide association analysis of Alzheimer’s disease. Four significant associations with the proposed method were found. Conclusion Our significant findings illustrated the practical importance of this new analysis method. Electronic supplementary material The online version of this article (doi:10.1186/s12881-015-0198-6) contains supplementary material, which is available to authorized users.


Background
Genome-wide association studies (GWASs) have been used to identify many genes involved in human diseases, and during the last decade, many disease-susceptibility variants have been identified. However, despite these successes, we have found that variants discovered from GWASs often explain only a small proportion of the heritability of diseases [1,2]. For example, SNPs significantly associated with human height explain only about 5 % of phenotypic variance, despite studies of tens of thousands of people [3]. Many reasons, such as rare causal variants and gene/gene interactions, have been attributed to this so-called "missing heritability". However, the low power induced by the multiple-testing problem is still an intractable issue in GWASs, and further investigations of the most efficient strategies for genetic association analysis are necessary.
Careful selection of samples based on phenotypes can lead to improved power for the discovery of risk variants [4][5][6][7][8][9][10][11]. One such example is the extreme discordant sibpair design in linkage analysis, which may result in a substantial increase in statistical power when compared to other sib-pair designs [11,12]. Similarly, ascertaining the extremes of quantitative phenotypes from large population cohorts has also been shown to increase the power to identify associated variants [13][14][15]. In such a design, the effect of ascertainment conditions are homogeneous between individuals, and existing methods, such as the Cochran-Armitage(CA) trend test [16], can be an efficient choice. However, in association analysis using extended families, the effects of ascertainment bias are often heterogeneous among family members, and depending on their relationships with probands, different magnitudes of ascertainment bias may be generated. In particular, the probability of each individual being affected when his or her relatives are affected is similar to the prevalence, if the heritability is small, which indicates that the heterogeneous effect of the ascertainment bias depends on the magnitude of heritability. However, the heterogeneous effects of ascertainment conditions and the influence of heritability on it have not yet been investigated, and should therefore be taken into account for association analysis.
Recently, the CA trend test was extended for association analysis of dichotomous phenotypes with family-based samples [17,18]. These statistics compares the genotype frequencies between affected and unaffected individuals, and the genetic association with family-based samples is tested by building a genotype correlation matrix with either kinship coefficients or an empirical correlation matrix estimated from large-scale genetic data. This approach has been extended to include family members with known phenotypes and missing genotypes or vice versa. By the nature of these statistics, it performs well for ascertained family-based samples and it can be an efficient choice, even for a case-control design, if the relatives' phenotype information is available. However, their statistical efficiency is affected by the heterogeneous effect of the ascertainment bias on family members, and for extended families, its effects on statistical efficiency can be substantial.
In this report, we consider the heterogeneous effects of the ascertainment bias on family members for dichotomous phenotypes. By the nature of the proposed methods, individuals with missing genotypes and non-missing phenotypes can be utilized, and incorporation of the estimated kinship matrix to the proposed statistic provided robustness against the population substructure. The proposed method consists of two steps; the probability for each family member to be affected was calculated using a latent continuous liability [19], and then this probability is incorporated into a quasi-likelihood score test. With an extensive simulation, we showed that the proposed method performed better than the existing methods, particularly for a disease with large heritability. Application of our method to Alzheimer's disease (AD) demonstrated its practical use in the detection of genetic associations in ascertained family-based samples.

Notations and statistic
We assumed that there were n families and n i family members in each family. We considered the situation where the family of size n i was ascertained because it contained a particular set of p i members, and we let q i = n ip i . We called the members of the set of p i family members "probands", and the remaining q i individuals "non-probands". To provide a clearer motivation on this concept, we randomly selected two families, family 1 and 2, from our AD data (see Fig. 1). In family 1 ( Fig. 1-(a)), individual 9 was diagnosed as AD and individuals 3-8 were selected as her relatives for genetic analysis. In family 2 ( Fig. 1-(b)), individual 3 was diagnosed as AD, and individuals 4-6 were selected. Therefore p 1 = p 2 = 1, q 1 = 6 and q 2 = 3 in this example. In real data analysis, p i is often 1 and q i = n i -1. We assumed that N individuals were available and thus N = ∑ i n i . The genotypes were coded as 0, 1, or 2, according to the number of disease alleles. x ij P and x i ' j ' N were defined as the genotypes of proband j and non-proband j' in family i and family i', respectively. Phenotypes were coded as 0 for an unaffected individual and 1 for an affected individual. If we let the prevalence of the disease be p, a missing phenotype was coded as p. We denoted the phenotypes of a proband and non-proband by y ij P and y i ' j ' N , respectively, and the vectors for genotypes and phenotypes in family i were defined by (b) Sample pedigree 2 Fig. 1 Two sample pedigree structures from AD data; (a) individual 9 of family 1 was selected as a "proband" and individual 3-8 of family 1 were selcted as "non-probands" (b) individual 3 of family 2 was selected as a "proband" and individual 4-6 of family 2 were selcted as "non-probands" We also denoted the w × w identity matrix by I w , and the w × 1 column vector 1 w indicated a vector in which all elements were 1. Let π ijj ' P and π ijj ' N be the kinship coefficient between probands j and j' in family i, and nonproband j and j' in family i, respectively. In addition, we let π ijj ' PN be the kinship coefficient between proband j and non-proband j' in family i, and let d ij P and d ij' N be the inbreeding coefficient for proband j and non-proband j' in family i, respectively. The inbreeding coefficient is the parameter that quantifies the departure from Hardy-Weinberg equilibrium (HWE) and ranges from 0 to 1. Several approaches [20,21] that can estimate d ij have been proposed. We let and R i is defined by If we let q A be the disease allele frequency, E(X i ) was 2q A 1 n i , and q A is estimated with the best linear unbiased estimator (BLUE). var(X i ) is expressed by σ 2 R i , and σ 2 is equal to 2q A (1 -q A ) under HWE.
When we analyzes the distribution of genotypes as in the FBAT approach, the statistical efficiency of the test statistic could be improved by adjustments of the phenotype with the so-called offset [22]. If we let μ ij P and μ i'j' N be offsets for proband j and non-proband j' in family i and family i', respectively, the offset vector for family i is defined as We denoted a minor allele frequency (MAF) of a variant in unaffected individuals by q. We assumed [18] that for a constant γ, where 0 < 2p + γ < 1. Then, the score for a variant [18,23] can be defined by The variance of S is and we considered the following statistic [17,18]: This statistic will be denoted by WL in the remainder of this report.

Adjusting the heterogeneous ascertainment bias
Families are often selected based on some probands, and the probability for family members to be affected depends on their relationship with the probands. Additional file 1 shows that the incorporation of conditional probability of each individual being affected to WL as offset lead to asymptotically smaller variance and therefore the adjustment of heterogeneous ascertainment bias is required to improve the statistical power of WL. This probability could be estimated with the liability model if the heritabilities, h 2 , and prevalence, p, were available. We let l ij P and l i ' j ' N be the liability of proband j and non-proband j' in family i and family i', respectively, and let L P i ¼ ðl P i1 ; …; l P ip i Þ and We assumed that each liability followed the standard normal distribution, and their joint distributions were Benchek and Morris [24] reported that significant asymptotic biases are likely to arise when the multivariate normal (MVN) liability assumption is not met and in such a case, different assumptions should be considered. We assume that M i P * and V i P * are the expectation and variances of L i P when their disease statuses are conditioned. If all probands are affected, they becomes They can be calculated with the numerical algorithms [25]. If p i is 1, both can be simply calculated. We denote the cumulative and probability density function of standard normal distribution by Ф(·) and ϕ(·). If we let c be the (1-p)th quantile of the standard normal distribution, M i P * and V i P * becomes With Pearson-Aitken formula [26,27], we could obtain the conditional mean and variance-covariance matrix of L i N given L P i > 1 p i ⋅c as follows: We denoted the jth element in M i N * by m j N * and the jth diagonal element in V i N * by v j N * . Then the probability of being affected for a non-proband under multivariate normality of the liabilities could be calculated as and this will be incorporated into the proposed statistic as offset. Thus far, we have assumed that there was a well-designed set of p i individuals who were "probands", and for this situation, we calculated the statistic as indicated and denoted FQLS 1 . However in practice, different ascertainment condition such as sequential sampling frame [28] are often utilized, and the set of p i individuals will not be well defined. For this situation, we calculated the probability for each individual to be affected under the assumption that all the other family members were "probands", and thus p i = n i -1 and q i = 1. The statistic calculated this way was denoted by FQLS 2 .

The simulation model
In our simulation studies, we considered two types of family structures; nuclear families with five offspring and the extended families that consist of 13 individuals along 3 generations (see Fig. 2). The latter will be called extended families in the remainder of this report. The disease allele frequency, p, was assumed to be 0.2. If we denoted the disease allele frequency by q A , the genotype frequencies for AA, Aa, and aa became q A 2 , 2q A (1q A ), and (1q A ) 2 under HWE, respectively, and founders' genotypes were generated under the corresponding multinomial distribution. The genotypes for non-founders were generated with randomly generated Mendelian transmission. The disease status was generated with the liability  . P-values were calculated based on 5000 replicates when the number of families was 900. The genetic effect β was assumed to be 0, and the minor allele frequency was 0.2 threshold model. Once continuous liabilities that consisted of polygenic effects and random errors were generated, they were transformed to being affected if they were larger than the threshold; and otherwise, they were considered to be unaffected. The threshold was chosen to preserve the prevalence, and prevalence was assumed to be 0.2. Continuous liability was determined by combining the phenotypic mean, polygenic effect, main genetic effect, and random error. The main genetic effect for each individual was the product of β and the number of disease alleles. If we denoted the relative proportion of the phenotypic variance attributable to the main disease gene by h a 2 , and h 2 was a heritability for continuous liability, β was calculated by For the evaluation of type-1 errors and power, h a 2 was assumed to be 0 and 0.005, respectively. Phenotypic correlations between family-members were explained by the polygenic effects. Parental polygenic effects were generated from N(0, h 2 ), and h 2 was assumed to be 0.2, 0.5, or 0.8. For non-founders, the average of maternal and paternal polygenic effects was combined with the values independently sampled from N(0, 0.5 h 2 ) for the polygenic effects of offspring. Random errors were generated from N(0, σ e 2 = 1-h 2 ). For each replicate, sampling was repeated until a given number of ascertained families was generated. Type-1 error estimates were calculated with 5000 replicates, and empirical power estimates were calculated with 1000 replicates.

Evaluation of the proposed methods with simulated data
The empirical type-1 errors for FQLS 1 and FQLS 2 were evaluated from 5000 replicates under the situation of no association (h a 2 = 0), and 900 nuclear families with five offspring in Fig. 2 were generated for each replicate. Fig. 3 shows the quantile quantile (QQ) plots from 5000 replicates, and the nominal significance levels for both methods were preserved for various significance levels. We also estimated the empirical type-1 error rates at the 0.01 and 0.05 significance levels; the empirical type-1 error estimates of FQLS 1 and FQLS 2 preserved these nominal significance levels ( Table 1). These results verified that the use of the approximation to the standard normal distribution resulted in an accurate assessment of significance for the proposed methods. Table 2 Empirical power estimates for scenario 1 when h 2 is 0.2. The empirical power estimates for scenario 1 were calculated with 1000 replicates at the both 0.01 and 0.001 significance levels. The disease allele frequency was assumed to be 0.2, and the prevalence was assumed to 0.2. The relative phenotypic variance attributable to the main disease gene was assumed to be 0.005  The empirical powers at the various significance levels were measured based on 1000 replicates at the 0.01 and 0.001 significance levels. The relative proportion, h a 2 , of phenotypic variance attributable to the main disease gene, 2p A (1p A )β 2 , was assumed to be 0.005, and nuclear and extended families in Fig. 2 were considered for the power comparison. In the first simulation setting, the numbers of nuclear families were assumed to be 100, 300, 600, 900, 1200, and 1400, and half of the families were ascertained if the number of affected family members was larger than or equal to n proband , and the other half of the families were ascertained if the number of unaffected family members was larger than or equal to n proband . Therefore, if 100 nuclear families were generated, half of nuclear families should have more than or equal to n proband affected family members, and the other half should have at least n proband unaffected family members. We assumed that the heritabilities were 0.2, 0.5, and 0.8, and results are shown in Tables 2, 3, 4, respectively. In the second simulation setting, the numbers of extended families were assumed to be 100, 300, 600, and 900, and all families were ascertained if the number of affected amily members was larger than or equal to n proband . Empirical power estimates for scenario 2 were calculated when h 2 = 0.2, 0.5, and 0.8, and the data are shown in Tables 5, 6, 7, respectively. Our results showed that either FQLS 1 or FQLS 2 was usually the most efficient statistic, and the least efficiency was provided from WL. In particular, the power gap between the proposed methods and WL was largest if h 2 was 0.8, which indicates that power improvement may be proportional to the heritability. Table 4 Empirical power estimates for scenario 1 when h 2 is 0.8. The empirical power estimates for scenario 1 were calculated with 1000 replicates at the both 0.01 and 0.001 significance levels. The disease allele frequency was assumed to be 0.2, and the prevalence was assumed to 0.2. The relative phenotypic variance attributable to the main disease gene was assumed to be 0.005 The bold text indicates the highest empirical estimate of the power for each situation Table 3 Empirical power estimates for scenario 1 when h 2 is 0.5. The empirical power estimates for scenario 1 were calculated with 1000 replicates at the both 0.1 and 0.001 significance levels. The disease allele frequency was assumed to be 0.2, and the prevalence was assumed to 0.2. The relative phenotypic variance attributable to the main disease gene was assumed to be 0.005 The bold text indicates the highest empirical estimate of the power for each situation If h 2 was 0.2, the proposed methods were only slightly better than WL. While all methods in our power comparison focused on the distribution of genotypes to calculate statistics, the proposed methods uniquely considered the heterogeneous effects of ascertainment bias among family members which were proportional to the magnitude of heritability; this explained the power improvement of the proposed methods. Furthermore the differences of empirical power estimates from WL and the proposed methods are larger for Tables 5, 6, 7 than Tables 2, 3, 4, which indicates that the heterogeneity of ascertainment condition may be positively related with family size and the proposed methods become more efficient for large families. Last our simulation results show that FQLS 2 was slightly better than FQLS 1 , and this may be induced by the uncertainty of probands in our simulation studies. Therefore, we concluded that the incorporation of a sampling scheme to the offset could make a substantial difference, and test statistic should be carefully selected depending on type of sampling scheme.The bold text indicates the highest empirical estimate of the power for each situationThe bold text indicates the highest empirical estimate of the power for each situationThe bold text indicates the highest empirical estimate of the power for each situationThe bold text indicates the highest empirical estimate of the power for each situationThe bold text indicates the highest empirical estimate of the power for each situationThe bold text indicates the highest empirical estimate of the power for each situation

Robustness of the proposed methods against the misspecification of prevalence
The statistical powers of the proposed methods may depend on the accuracy of the prevalence and we evaluated the sensitivity of the proposed method to the misspecified prevalence with simulated data. h a 2 and h 2 were assumed to be 0.05 and 0.8, and nuclear families (Fig. 2-(a)) were considered in this simulation. The number of nuclear families was assumed to be 900 and n proband was assumed to be 1, 2, or 3. Prevalence was assumed to be 0.2 for phenotype generation, and the offset for which we recommended prevalence was set to be 0.1, 0.2 or 0.3 for calculation of FQLS 1 and FQLS 2 . In particular, y ij for individuals with missing phenotypes are coded by the Table 6 Empirical power estimates for scenario 2 when h 2 is 0.5. The empirical power estimates for scenario 2 were calculated with 1000 replicates at the both 0.1 and 0.001 significance levels. The disease allele frequency was assumed to be 0.2, and the prevalence was assumed to 0.2. The relative phenotypic variance attributable to the main disease gene was assumed to be 0.005 The bold text indicates the highest empirical estimate of the power for each situation The bold text indicates the highest empirical estimate of the power for each situation assumed prevalence, and sensitivity of the proposed methods can be substantial when there are individuals with missing phenotypes. Therefore, individuals were randomly selected from non-probands, and their phenotypes were assumed to be unknown for calculation of the proposed statistics. The number of family members with missing phenotypes in each family was denoted by n missing . The empirical powers were calculated at the 0.01 significance level with 1000 replicates. Table 8 shows that the results obtained by setting prevalence to be 0.1 and 0.3 are similar to the results when the prevalence was set to be 0.2, which indicates that the power loss attributable to the misspecified prevalence is not substantial. Furthermore the empirical power estimates are positively related with n proband and inversely related with n missing . If n missing is larger than 3, the power loss may be more substantial.The bold text indicates the reference for the proposed statistics to compare with the misspecified prevalence.

Application of the proposed method to AD
AD is an irreversible, progressive brain disorder characterized by genetic heterogeneity. However, the genetic variations that contribute to AD still remain elusive. Thus, we applied the proposed method for identification of the disease susceptibility loci for AD. The heritability and prevalence of AD are approximately 0.8 [29,30] and 0.1, respectively; therefore, we chose heritabilities of 0.8, and a prevalence of 0.1 for the calculation of proposed methods. Samples were collected as part of the National Institute of Mental Health Genetics Initiative (NIMH). The NIMH Alzheimer's Disease Genetics Family Sample was used along with the information about the genotype platform (Affy 6.0) [20,31], and ethical approach and participant approval were obtained through the NIMH IRB panel. Families were selected based on the disease status of a certain family member. However the proband for each family was not clear, and FQLS 2 was uniquely applied for the proposed method. 1376 individuals from 410 families were available, and all families were nuclear. All individuals were of self-reported European ancestry. HWE for each single nucleotide polymorphism (SNP) was tested, and MAFs were estimated. SNPs for which p-values for HWE were less than 10 −6 or MAFs less than  The bold text indicates the highest empirical estimate of the power for each situation 0.05 were excluded, and therefore, 417,680 SNPs were analyzed for genetic association analysis. The choice of kinship coefficient matrix for the proposed method depends on the presence of population substructure. To confirm the presence of population substructure, multidimensional-scaling [32] analysis was performed with PLINK1.07 [20], and we constructed a multidimensional-scaling plot (Fig. 4) to provide evidence concerning the presence of population substructure. Therefore, the genomic control method [31,33,34] was used for explicit detection and correction of population stratification with common SNPs, and they were incorporated into both FQLS 2 and WL. Fig. 5 shows QQ plots for FQLS 2 and WL; these plots revealed that the presence of population substructure was appropriately adjusted. Results showed four genome-wide significant results for FQLS 2 and one significant result for WL.
Detailed information for these significant results is provided in Table 9. We also considered FBAT statistics [35]. FQLS 2 identified four genome-wide significant SNPs, while WL and FBAT identified one genome-wide significant SNP. In addition, one SNP that was significant according to WL was more significant according to FQLS 2 . The most significant result acquired using FQLS 2 was SNP4 (p = 2.10× 10 −10 ). The other three SNPs, i.e., SNP1, SNP2 and SNP3, reached the genome-wide significance level.

Discussion
Although major advances in high-density genome scans have enabled the genetic association analysis of more than 10,000 individuals, disappointing results in the mapping of many common diseases have illustrated the need for more powerful methods for detecting disease susceptibility loci. Statistical efficiency is known to being affected by the ascertainment bias, and its careful adjustment can   lead to substantial improvement of statistical power [36][37][38]. In particular, genetic association analysis has often been conducted using family-based designs, but without addressing the fact that the probability for each family member to be affected is inversely related with the familial relationship with affected probands. In this report, we proposed new methods to adjust this heterogeneity with known prevalence and heritability. Our simulation studies showed that the proposed methods provided substantial power improvement. In particular, the mis-specified heritability and prevalence can lead to the statistical power loss for the proposed methods, but it was found to be not substantial at least in our simulation studies. FQLS 1 and FQLS 2 were suggested, and FQLS 1 is an efficient choice if probands for each family are clearly defined and all remaining family members are incorporated to the genetic analysis. However, these conditions are often not satisfied, and different methods such as sequential sampling frame [28] are usually utilized. Simulation studies showed that FQLS 2 is usually better than FQLS 1 if the ascertaining condition is not clearly defined and thus we recommend FQLS 2 unless probands are clearly defined. However we considered the limited ascertainment conditions and comprehensive simulation studies are still necessary. Furthermore, the proposed method was conceptually simple and can be applied to the large families. Our methods require only a single calculation of offset for all markers, and the real data analysis could be completed with a single CPU in a few hours. For M markers and N individuals, the time complexity is O(N 3 + MN 2 ) for the proposed method. The proposed method was implemented with C++, and can be downloaded from http://healthstat. snu.ac.kr/mfqls/.
Heterogeneity between samples is an important issue in large-scale genetic analysis, and the proposed method can likely be applied to various additional scenarios with some modifications. For instance, the disease status of relatives reveals the importance of genetic components for each individual, and for this reason, such information has been used, albeit only on occasion, in genetic association analysis. The effect of relatives' disease statuses is dependent on prevalence and heritability, and the probability for each individual to be affected could be calculated with the proposed method. This probability can be used to improve the statistical efficiency of genetic association analysis. In addition, the heterogeneity of the ascertainment bias is often an important issue for genome-wide meta-analysis because samples are collected from multiple medical centers [39,40], and different sampling schemes among studies need to be adjusted to improve statistical efficiency. Therefore, we believe that proposed method can be extended to provide a statistical framework that adjusts the heterogeneity between samples.

Conclusions
We proposed FQLS method to adjust this heterogeneity with known prevalence and heritability and the software was implemented with C++. We identified several significant associations between AD and SNPs, and their potential functional information will provide the better understanding of the pathogenesis of AD. Although this study has some limitations, our proposed methods illustrated important features required for genetic analysis with family-based samples, and an extension of the proposed method to rare variant association analysis such as FARVAT [41] will be investigated in future studies.