Use of stochastic simulations to investigate the power and design of a whole genome association study using single nucleotide polymorphism arrays in farm animals

Use of stochastic simulations to investigate the power and design of a whole genome association study using single nucleotide polymorphism arrays in farm animals first-author †E-mail:benoit.auvray@agresearch.co.nz †E-mail:benoit.auvray@agresearch.co.nz benoit.auvray@agresearch.co.nz Applied Biotechnologies Group, AgResearch Limited, Invermay Research Centre, Private Bag 50034, Mosgiel 9053, New Zealand Applied Biotechnologies Group, AgResearch Limited, Invermay Research Centre, Private Bag 50034, Mosgiel 9053, New Zealand This paper presents a quick, easy to implement and versatile way of using stochastic simulations to investigate the power and design of using single nucleotide polymorphism (SNP) arrays for genome-wide association studies in farm animals. It illustrates the methodology by discussing a small example where 6 experimental designs are considered to analyse the same resource consisting of 6 006 animals with pedigree and phenotypic records: (1) genotyping the 30 most widely used sires in the population and all of their progeny (515 animals in total), (2) genotyping the 100 most widely used sires in the population and all of their progeny (1 102 animals in total), genotyping respectively (3) 515 and (4) 1 102 animals selected randomly or genotyping respectively (5) 515 and (6) 1 102 animals from the tails of the phenotypic distribution. Given the resource at hand, designs where the extreme animals are genotyped perform the best, followed by designs selecting animals at random. Designs where sires and their progeny are genotyped perform the worst, as even genotyping the 100 most widely used sires and their progeny is not as powerful of genotyping 515 extreme animals. Simulation Association study Single nucleotide polymorphism (SNP) Power Quantitative trait loci (QTL) INTRODUCTION 1 Genome-wide single nucleotide polymorphism (SNP) array technologies offer opportunities for gene discovery in farm animals (Barendse et al., 2007 2007 ). However, the use of such arrays requires careful planning due to technological or cost-related constraints or limitations inherent to association studies and other linkage disequilibrium (LD)-based methods (Lawrence et al., 2005 2005 ; Clark et al., 2005 2005 ). These limitations depend on many factors specific to each experiment. It is therefore very desirable to optimize the design of an experiment in order to maximize power given the resource at hand, before conducting any experiment using SNP arrays. Several studies have focused on computing power for association studies (Kang et al., 2004 2004 ; Guedj et al., 2007 2007 ), but these do not address the problem of potential confounding factors such as population relatedness. The use of stochastic simulations is a simple, yet powerful and versatile way to handle this problem. Indeed, replicating in silico an experiment using varying assumptions about the data and different methods of analysis permits a reasonable estimation of the power of the experiment given these factors and thus the identification of a suitable experimental design. As an example, this paper compares the power of six different designs for a genome-wide SNP association study. MATERIALS AND METHODS 1 Description of the example 2 A real animal resource composed of 6 006 animals born between 1980 and 2005 has been chosen for this example. Polygenic breeding value estimates (EBV) for a trait of interest with a heritability of ~0.2 are available for all animals. A genome-wide association study using a 60K SNP array is to be designed, and the 6 following designs are envisaged: (1) genotyping the 30 most widely used sires and all of their progeny, which represents 515 animals in total (“sires-515”), (2) genotyping the 100 most widely used sires in the population and all of their progeny, which represents 1 102 animals in total (“sires-1102”), (3) genotyping 515 animals randomly selected (“rand-515”), (4) genotyping 1 102 animals randomly selected (“rand-1102”), (5) genotyping 515 animals from the tails (i.e. the most extreme) of the phenotypic distribution (“extr-515”) or (6) genotyping 1 102 animals from the tails of this distribution (“extr-1102”). Simulation procedure 2 The methodology employed to simulate populations is as follows: 1. Twenty ancestral chromosomal segments with 41 evenly spaced SNP loci are created randomly, with all 41 loci being polymorphic (this amounts to ensuring that all the SNP minor allele frequencies (MAFs) are ≥0.05 at the beginning of the simulation). These 20 segments are paired to create 10 ancestral animals. Chromosomal segments are chosen to be 2 cM long so that SNPs are 0.05 cM apart (and 60 000 SNPs cover a region 30 M long). The middle SNP (SNP 21) is chosen to be a quantitative trait locus (QTL). The 2 cM length has been determined as being sufficient to include most SNPs in LD with the QTL. 2. The 10 ancestral animals are mated to create the next generation of 100 animals, which is in turn used to create the following generation of 100 animals and so on, with no generation overlapping. This process is carried out for 10 generations. A slight selection pressure for the QTL is applied during the whole process, in order to ensure that the QTL frequency after the 10 generations is close to a target QTL frequency. 3. The 100 animals simulated in the final generation of this process are used as founders (unknown parents) of the real pedigree, and gene-dropping (MacCluer et al., 1986 1986 ) is used to simulate the genotypes of the 6 006 animals that can potentially be genotyped. The phenotype of an animal i is simulated as: y i = x i q + u i + e i , where y i is the simulated phenotype, x i is the number of copies of allele 1 at the QTL, q is the QTL (allelic substitution) effect used in the simulation, u i is the EBV of animal i and e i is a random residual drawn from N (0,(1− h 2 ) σ u 2 / h 2 ), with the additive genetic variance ( σ u 2 ) equal to 0.543 and the heritability ( h 2 ) equal to 0.2. Simulation implementation 2 For each group, the following factors are made variable: 1. The target QTL MAF is set to 0.1 and 0.4. 2. The QTL effect is set to 0, 0.25 phenotypic standard deviations ( σ p ) and 0.50 σ p . 3. Different animals are “genotyped”, i.e., are actually used in the analysis of the simulated array, depending on the design chosen. Two hundred replicates of each combination of factors were carried out. Due to the stochastic nature of the simulation, the actual QTL MAF for each replicate differs from the targeted MAF, and at some times the QTL may even be lost due to random drift (especially when targeting the low QTL MAF). Analysis of the simulated array and power estimation 2 For the animals genotyped and for each of the 40 SNPs other than the QTL, the following mixed model is fitted: y = x b + Zu + e , with $$\operatorname{var} \left( {\begin{array}{*{20}c}{\mathbf{u}} \\{\mathbf{e}} \\\end{array} } \right) = \left( \begin{gathered}{\mathbf{\sigma }}_{\mathbf{u}}^2 {\mathbf{A}}{\text{ }}0 \\{\text{ }}0{\text{ }}{\mathbf{\sigma }}_{\mathbf{e}}^2 {\mathbf{I}} \\\end{gathered} \right)$$ , where y is a vector of simulated phenotypes, b is the fixed regression coefficient of the phenotype on the number of copies of the first SNP allele, u is a vector of random genetic additive effects, e is a vector of random residuals, x is a vector of number of copies of the first SNP allele, Z is an identity incidence matrix, I is an identity matrix and A is the additive relationship matrix computed from the pedigree. The statistic $$W = |b|/\sqrt {{\mathbf{C}}^{bb} }$$ is then calculated, where C bb is the element relating to b in the inverse of the left hand side of the mixed model equations (Henderson, 1984 1984 ) used to solve the model above. For each design, empirical “genome-wise” 5% α significance thresholds are computed as follows: 1. Assuming that genomes consist of 1 500 independent chromosomal segments of 2 cM, each with 40 SNPs, for a total of 30 M and 60 000 SNPs, a “genome-wise” Bonferroni corrected 5% α ( α genome-wise ) is computed as (Shaffer, 1995 1995 ): α genome-wise =1−(1−5%) (1/1500) =3.42×10 −5 . 2. The null distribution of the maximum W across the 40 SNPs is obtained by using each replicate of the design where a null QTL effect has been simulated, and the “genome-wise” threshold T is taken to be the 1− α genome-wise quantile of that distribution. As only 400 replicates have been carried out at most, the threshold is interpolated, and numerically very close to the maximum W . Then, for each design, the function loess from the base package of the R software (R Development Core Team, 2007 2007 ) is used (with the default parameters) to model the power by fitting locally second degree polynomials by least squares as: p = f ( MAF , q ), where p is a vector of indicators equal to 1 when W ≥ T and 0 otherwise, MAF is a vector of minor allele frequencies and q is a vector of QTL effects. RESULTS AND DISCUSSION 1 The results can be observed from various angles. Fig. 1 1 displays for each of the 6 experimental designs the power surfaces as a function of the QTL MAF and effect. Fig. 2 2 shows for each design the power curves as a function of the QTL effect size, for 3 QTL MAF (0.1, 0.25 and 0.4), while Fig. 3 3 presents the power curves as a function of the QTL MAF, for 3 QTL effect sizes (0.1 σ p , 0.25 σ p and 0.4 σ p ). The rows represent the 3 sampling methods (“sires”, “rand” or “extr”) and the columns the number of animals genotyped (515 or 1 102). The colour black corresponds to power=0 and the colour white to power=1. Contour lines of equal power for power=0.5, 0.8 and 0.9 are plotted. The white arrows exemplify the difference between designs by highlighting respectively the QTL effect size needed to reach a statistical power of 0.5 when chasing a QTL with a MAF of 0.1 and the QTL MAF needed to reach a statistical power of 0.5 when chasing a QTL with an effect of 0.2 σ p . The behaviour of the contour lines going back up for large MAF is an artefact of the power modelling procedure Power values >1 or decreasing as the QTL effect sizes increase (observed for large QTL effects and MAF) are artefacts of the power modelling procedure Power values >1 or decreasing as the QTL MAF increases (observed for large QTL effects and MAF) as well as power values <0 are artefacts of the power modelling procedure Obviously, the power is lower when fewer animals are genotyped. Maybe not surprisingly as the increased power of selective genotyping designs has been largely documented [see for instance in (van Gestel et al., 2000 2000 )], the power is also lower for the designs genotyping the sires compared to those genotyping the extreme animals, even when comparing designs “sires-1102” and “extr-515”, with the first having more than twice the number of animals genotyped. This difference can be very large when the QTL MAF is small ( < 0.2). The power curves obtained for the two “extr” designs are very similar. The designs “rand” perform slightly better than the designs “sires”. The reasons for this are unclear. The reduction of power does not seem to be due to a difference between the phenotypic distributions of both designs as these are very similar. The LD patterns are also comparable between the 2 designs, with for instance the distribution of the r 2 between adjacent SNPs having similar distributions with mean ~0.08 and variance ~0.03 for both designs. Despite having similar LD, the 2 designs have very different average additive relationship between animals, much higher (~4 times higher when genotyping 1 102 and up to ~7 times higher when genotyping 515) in the “sires” groups than in the “rand” groups. This could lead to more confounding between polygenic and QTL effects in the “sires”, thus causing a reduction of the power by decreasing the significance of the SNP effects. Overall, the power levels obtained with design “sires-515” are too low to be practically useful when the QTL MAF is small (≤0.2) or the QTL effect is small to moderate (≤0.25 σ p ), while the power levels for design “sires-1102” are not acceptable when both QTL MAF and QTL effect are small to moderate. For a QTL effect of ~0.1 σ p or under, all 6 designs have low power. Given the results, it seems that the best strategy would be to concentrate all resources into genotyping the tails of the phenotypic distribution. It also seems that more than doubling the number of animals from 515 to 1 102 is unnecessary if chasing QTL with effect sizes over 0.25 σ p . Nevertheless, power levels obtained for low QTL MAF (~0.1) and QTL effect size of 0.25 σ p or under are relatively low even for the design “extr-1102” and suggest that not many new small or rare QTL would be found using such an experimental design. CONCLUSION 1 This example presents an easy and quick way to implement stochastic simulations to investigate the power and design of using SNP arrays for genome-wide association studies. A number of general conclusions can be derived from this small study. Firstly, genotyping the extreme animals is more powerful a design than selecting animals based on relative contribution in a pedigree or at random. However, the increase of power observed with selective genotyping is of course mitigated by the fact that genotypes can be used to map QTL for different traits, thus reducing the applicability of these designs, but also by the potential greater value (for an animal industry, for instance) of having the genotypes of more widely used animals, compared to rarely used ones or random animals. Secondly, besides the obvious effects of the QTL effect size and the number of animals genotyped, the large effect that the QTL MAF has on the power has been clearly demonstrated. In conclusion, this study reaffirms the need for careful planning of genome-wide association studies using SNP arrays. References W Barendse , A Reverter , RJ Bunch , BE Harrison , W Barris , MB Thomas . A validated whole-genome association study of efficient food conversion in cattle . Genetics , 2007 . 176 ( 3 ): 1893 - 1905 . DOI: 10.1534/genetics.107.072637 http://doi.org/10.1534/genetics.107.072637 . A validated whole-genome association study of efficient food conversion in cattle AG Clark , E Boerwinkle , J Hixson , CF Sing . Determinants of the success of whole-genome association testing . Genome Res , 2005 . 15 ( 11 ): 1463 - 1467 . DOI: 10.1101/gr.4244005 http://doi.org/10.1101/gr.4244005 . Determinants of the success of whole-genome association testing M Guedj , E Della-Chiesa , F Picard , G Nuel . Computing power in case-control association studies through the use of quadratic approximations: application to meta-statistics . Ann Hum Genet , 2007 . 71 ( 2 ): 262 - 270 . DOI: 10.1111/j.1469-1809.2006.00316.x http://doi.org/10.1111/j.1469-1809.2006.00316.x . Computing power in case-control association studies through the use of quadratic approximations: application to meta-statistics CR Henderson . Applications of Linear Models in Animal Breeding , 1984 . : Guelph, Canada University of Guelph , . Applications of Linear Models in Animal Breeding SJ Kang , D Gordon , SJ Finch . What SNP genotyping errors are most costly for genetic association studies? . Genet Epidemiol , 2004 . 26 ( 2 ): 132 - 141 . DOI: 10.1002/gepi.10301 http://doi.org/10.1002/gepi.10301 . What SNP genotyping errors are most costly for genetic association studies? RW Lawrence , DM Evans , LR Cardon . Prospects and pitfalls in whole genome association studies . Philosophical Transactions of the Royal Society B Biological Sciences , 2005 . 360 ( 1460 ): 1589 - 1595 . DOI: 10.1098/rstb.2005.1689 http://doi.org/10.1098/rstb.2005.1689 . Prospects and pitfalls in whole genome association studies JW MacCluer , JL Vandeberg , B Read , OA Ryder . Pedigree analysis by computer simulation . Zoo Biol , 1986 . 5 ( 2 ): 147 - 160 . DOI: 10.1002/zoo.1430050209 http://doi.org/10.1002/zoo.1430050209 . Pedigree analysis by computer simulation R Development Core Team R: A Language and Environment for Statistical Computing , 2007 . : Vienna, Austria R Foundation for Statistical Computing , (Available from: Http://www.R-project.org Http://www.R-project.org ) . R: A Language and Environment for Statistical Computing JP Shaffer . Multiple hypothesis testing . Ann Rev Psychol , 1995 . 46 561 - 584 . DOI: 10.1146/annurev.ps.46.020195.003021 http://doi.org/10.1146/annurev.ps.46.020195.003021 . Multiple hypothesis testing S van Gestel , JJ Houwing-Duistermaat , R Adolfsson , CM van Duijn , C van Broekhoven . Power of selective genotyping in genetic association analyses of quantitative traits . Behav Genet , 2000 . 30 ( 2 ): 141 - 146 . DOI: 10.1023/A:1001907321955 http://doi.org/10.1023/A:1001907321955 . Power of selective genotyping in genetic association analyses of quantitative traits 10.1631/jzus.2007.B0802 Project supported by Ovita Limited, Dunedin, New Zealand 2007-09-19 2007-09-20 2022-01-05 Copyright © 2007, Journal of Zhejiang University Science 2007 Journal of Zhejiang University-SCIENCE B (Biomedicine and Biotechnology) 11 08