# Fast Genome-Wide QTL Association Mapping on Pedigree and Population Data

###### Abstract

Since most analysis software for genome-wide association studies (GWAS) currently exploit only unrelated individuals, there is a need for efficient applications that can handle general pedigree data or mixtures of both population and pedigree data. Even data sets thought to consist of only unrelated individuals may include cryptic relationships that can lead to false positives if not discovered and controlled for. In addition, family designs possess compelling advantages. They are better equipped to detect rare variants, control for population stratification, and facilitate the study of parent-of-origin effects. Pedigrees selected for extreme trait values often segregate a single gene with strong effect. Finally, many pedigrees are available as an important legacy from the era of linkage analysis. Unfortunately, pedigree likelihoods are notoriously hard to compute. In this paper we re-examine the computational bottlenecks and implement ultra-fast pedigree-based GWAS analysis. Kinship coefficients can either be based on explicitly provided pedigrees or automatically estimated from dense markers. Our strategy (a) works for random sample data, pedigree data, or a mix of both; (b) entails no loss of power; (c) allows for any number of covariate adjustments, including correction for population stratification; (d) allows for testing SNPs under additive, dominant, and recessive models; and (e) accommodates both univariate and multivariate quantitative traits. On a typical personal computer (6 CPU cores at 2.67 GHz), analyzing a univariate HDL (high-density lipoprotein) trait from the San Antonio Family Heart Study (935,392 SNPs on 1357 individuals in 124 pedigrees) takes less than 2 minutes and 1.5 GB of memory. Complete multivariate QTL analysis of the three time-points of the longitudinal HDL multivariate trait takes less than 5 minutes and 1.5 GB of memory. The algorithm is implemented as the Ped-GWAS Analysis (Option 29) in the Mendel statistical genetics software package, which is freely available for Macintosh, Linux, and Windows platforms from http://genetics.ucla.edu/software/mendel.

Key words: genome-wide association study; pedigree; kinship; score test; fixed-effects models, multivariate traits

## 1 Introduction

Genome-wide association studies (GWAS) are now at a crossroads. After the discovery of thousands of genes influencing hundreds of common traits ^{9}, much of the low-hanging fruit has been plucked ^{13}, ^{36}. Because of the enormous sample sizes of current studies, new trait genes are still being uncovered. Unfortunately, most entail small effects. Is it possible that inheritance is predominantly polygenic, and a law of diminishing returns has set in? The push to exploit rare variants is one response to this dilemma. The previous generation of geneticists relied on linkage to map rare variants. Linkage mapping fell from grace because of its poor resolution. Reducing a genome search to a one or two megabase region leaves too large an expanse of DNA to sift. The real gold of linkage mapping may well be its legacy pedigrees ^{26}. Pedigree data is particularly attractive in association studies because it permits control of population substructure and study of parent-of-origin effects. Related affecteds are also more likely to share the same disease predisposing gene than unrelated affecteds. Even in population-based association studies, taking into account estimated identity-by-descent (IBD) information is apt to reduce false positives and increases power. The recent availability of dense marker data from genotyping chips enables quick and accurate estimation of global and even local IBD ^{7}.

Mendel | FaST-LMM | GEMMA | |
---|---|---|---|

Multi-threaded operation | Yes | Yes | No |

Can estimate kinships via SNPs | Yes | Yes | Yes |

Imports & exports kinship estimates | Yes | Yes | Yes |

Allows retained co-variates | Yes | Yes | Yes |

Allows linear constraints on co-variates | Yes | No | No |

Can use either LRT or score test | Yes | No | Yes* |

Allows multivariate analysis | Yes | No | Yes |

Can perform multiple univariate analyses | Yes | No | No |

Allows variance components | Yes | No | No |

Analyzes X-linked loci | Yes | No | No |

Automatic SNP filtering on MAF | Yes | No | Yes |

Allows non-additive SNP models | Yes | No | No |

Detects outlier pedigrees | Yes | No | No |

Detects outlier individuals | Yes | No | No |

Can simulate genotype/phenotype data | Yes | No | No |

Reads in fractional genotype values | No | Yes | Yes |

Geneticists turned to random sample and case-control data because of the relative ease of collecting population data and the computational challenges posed by pedigrees. The tide of computational complexity is now beginning to turn. To handle pedigree data in association testing, statistical geneticists have proposed semiparametric methods such as the generalized linear mixed model (GLMM) ^{1}, ^{3} and generalized estimating equations (GEE) ^{5}, ^{4}. Although such methods work for both quantitative and binary traits, they are compromised by current restrictions that reduce power. The GEE approach requires input of a working correlation structure for each pedigree. The kinship coefficient matrix is a natural candidate. However, current implementations require the same working correlation matrix across all clusters, which implicitly requires all pedigrees to have the same structure ^{4}. This is a dubious and restrictive assumption. In the limited context of case-control studies, recent methods such as ^{32}, roadtrips ^{33}, and FPCA ^{46} correct for pedigree and ethnically induced correlations by exploiting dense marker data. Other authors attack the same issues more broadly from the GLMM perspective ^{10}, ^{41}, ^{23}. Korte et al. ^{12} generalizes GLMM to multivariate traits. Models based on the transmission-disequilibrium test (TDT) ^{30} and its generalization, the family-based association test (FBAT) ^{15}, ^{16}, ^{35}, ^{37}, ^{38}, are promising but ignore covariates and polygenic background. See Van Steen ^{34} for a recent overview of FBAT methods for GWAS. We treat all of these extensions in a unified framework consistent with exceptionally fast computing.

The present paper re-examines the computational bottlenecks encountered in association mapping with pedigree data. It turns out that the previous objections to pedigree GWAS can be overcome. Kinship coefficients can be based on explicitly provided pedigree structure or estimated from dense markers when genealogies are missing or dubious. Frequentist hypothesis testing usually operates by comparing maximum likelihoods under the null and alternative hypotheses. Maximization of the alternative likelihood must be conducted for each and every marker. Score tests constitute a more efficient strategy than likelihood ratio tests. This is the point of departure taken by Chen and Abecasis ^{6}, but they use approximations that we avoid. The glogs program ^{31} makes similar approximations in the case-control setting. Here we consider arbitrary pedigrees and multivariate quantitative traits. Score tests require no additional iteration under the alternative model. All that is needed is evaluation of a quadratic form combining the score vector and the expected information matrix at the maximum likelihood estimates under the null model. Although it takes work to assemble these quantities, a careful analysis of the algorithm shows that fast testing is perfectly feasible.

In our implementation of score testing, the few SNPs with the most significant score-test p-values are automatically re-analyzed by the slightly more powerful, but much slower, likelihood ratio test (LRT). Our fixed effects (mean component) model assumes Gaussian variation of the trait; the two alleles of a SNP shift trait means. There is no confounding of association and linkage. This framework carries with it several advantages. First, it applies to random sample data, pedigree data, or a mix of both. Second, it enables covariate adjustment, including correction for population stratification. Third, it accommodates additive, dominant, and recessive SNP models. Fourth, it also accommodates both univariate and multivariate traits. And fifth, as just mentioned, it fosters both likelihood ratio tests and score tests. The mean component model is now implemented in our software package Mendel for easy use by the genetics community. In addition, Mendel provides a complete suite of tools for pedigree analysis, including GWAS data preparation and manipulation, pedigree genotype simulation (gene dropping), trait simulation, genotype imputation, local and global kinship coefficient estimation, and pedigree-based GWAS (ped-GWAS) ^{21}, ^{20}.

The competing software packages EMMAX ^{11}, MMM ^{27}, FaST-LMM ^{23}, ^{24}, and GEMMA ^{44}, ^{45} already implement variance component models for quantitative trait locus (QTL) analysis. Exhaustive comparison of Mendel to each of these programs is beyond the scope of the current paper. We limit our comparisons of Mendel to the state-of-art packages FaST-LMM and GEMMA, arguably the fastest and most sophisticated of the competition. Table 1 summarizes some of the qualitative features of these packages. Our numerical examples also demonstrate an order of magnitude advantage in speed of Mendel over FaST-LMM and GEMMA. This advantage stems from our careful formulation of the score test and our exploitation of the multicore processors resident in almost all personal computers and computational clusters.

## 2 Methods

### 2.1 QTL Association Mapping with Pedigrees

QTL association mapping typically invokes the multivariate Gaussian distribution to model the trait values over a pedigree. The observed trait value of person can be either univariate or multivariate. For simplicity we first assume is univariate and later indicate the necessary changes for multivariate . The standard model ^{17}
collects the corresponding trait means into a vector and the corresponding covariances into a matrix and represents the loglikelihood of a pedigree as

(1) |

where denotes the determinant function and the covariance matrix is typically parametrized as

(2) |

Here the variance component is the global kinship coefficient matrix capturing additive polygenic effects, and is a condensed identity coefficient matrix capturing dominance genetic effects. When pedigree structure is explicitly given, these genetic identity coefficients are easily calculated ^{17}. With unknown or dubious genealogies, the global kinship coefficient can be accurately estimated from dense markers ^{7}. The household effect matrix has entries if individuals and belong to the same household and 0 otherwise. Individual environmental contributions and trait measurement errors are incorporated via the identity matrix . Mendel’s implementation of this model can include both the two standard variance classes, additive and environmental, as well as the two extra variances classes, dominance and household. Inclusion of additional variance classes has no significant effect on Mendel’s speed of computation.

Genotype | Additive | Dominant | Recessive |
---|---|---|---|

`1/1` |
–1 | –1 | –1 |

`1/2` |
0 | –1 | +1 |

`2/2` |
+1 | +1 | +1 |

In general, a mixed model for QTL association mapping captures polygenic and other random effects through and captures QTL fixed effects through . Let denote the full vector of regression coefficients parameterizing . In a linear model one postulates that for some predictor matrix incorporating relevant covariates such as age, gender, and diet. In testing association against a given SNP, is augmented by an extra column whose entries encode genotypes according to one of the models (additive, dominant, and recessive) shown in see Table 2. To accommodate imprecise imputation in an additive model, these encodings can be made fractional. The corresponding component of , , is the SNP effect size. In likelihood ratio association testing one contrasts the null hypothesis with the alternative hypothesis . In testing a univariate trait, the likelihood ratio statistic asymptotically follows a distribution. In testing a multivariate trait with components, each row of must be replicated times. The likelihood ratio statistic then asymptotically follows a distribution. To implement likelihood ratio testing, iterative maximum likelihood estimation must be undertaken for each and every SNP under the alternative hypothesis. This unfortunate requirement is the major stumbling block retarding pedigree analysis.

Score tests serve as convenient substitutes for likelihood ratio tests. The current paper describes how to implement ultra-fast score tests for screening SNPs. Only SNPs with the most significant score test p-values are further subjected to the more accurate likelihood ratio test. An advantage of the likelihood ratio method is that it estimates effect sizes. In contrast, the score test only requires parameter estimates under the null hypothesis and involves no iteration beyond fitting the null model. The score vector is the gradient of the loglikelihood , where the full parameter vector includes variance components such as the additive genetic variance in addition to the regression coefficient vector . The transpose of the score is a row vector called the first differential of . The expected information is the covariance matrix of the score vector. It is well known that the expected value of the observed information matrix (negative second differential) coincides with ^{28}. The score statistic

is evaluated at the maximum likelihood estimates under the null hypothesis with the parameter of the alternative hypothesis set to 0.

### 2.2 Fast Score Test for Individual SNPs

Under the multivariate model, the expected information matrix for a single pedigree can be written in the block diagonal form

(3) |

where denotes the vector of variance parameters ^{17}. For independent pedigrees, the loglikelihoods (1) and corresponding score vectors and expected information matrices add. Hence, the block diagonal form of is preserved. Because the inverse of a block diagonal matrix is block diagonal, the score statistic splits into a piece contributed by the variance components plus a piece contributed by the mean components. The maximum likelihood estimate under the null model is a stationary point of the loglikelihood. Thus, the variance components segment of the score vector vanishes. We therefore focus on the mean components segment of the score vector.

If the pedigrees are labeled , then the pertinent quantities for implementing the score test are

where is the residual for pedigree and the covariance matrix for pedigree is determined by equation (2). See Chapter 8 of Lange ^{17} for a detailed derivation of the score and expected information. Since the score statistic is calculated from estimated parameters under the null model, residuals do not change when we expand the null model to the alternative model keeping . Calculation of the maximum likelihood estimate under the null is accomplished by a
quasi-Newton algorithm whose initial step reduces to Fisher scoring ^{22}, ^{17}.

For pedigree under the alternative hypothesis, the design matrix can be written as , where is the design matrix under the null hypothesis and conveys the genotypes at the current SNP. In testing a univariate trait, the entries of are taken from Table 2. If allele counts are imputed under the additive model, then the entries of may be fractional numbers drawn from the interval . In testing a multivariate trait with components, each row of must be replicated times. The only exceptions to this rule occur for people missing some but not all component traits;
otherwise, the covariance matrix for pedigree decomposes into a sum of Kronecker products ^{17}. Regardless of whether the trait is univariate or multivariate, one must compute the quantities

At the maximum likelihood estimates under the null model, the partial score vector vanishes. Hence, the score statistic for testing a SNP can be expressed as

where

In the score statistic , the covariance matrices and residual vectors are evaluated at the maximum likelihood estimates under the null model. Large sample theory says that asymptotically follows a distribution.

These formulas suggest that we precompute and store the quantities , , and for each pedigree and the overall sum at the maximum likelihood estimates under the null hypothesis. From these parts, the basic elements of the score statistic can be quickly assembled. The most onerous quantity that must be computed on the fly as each new SNP is encountered is . If there are people in pedigree , then computation of the quadratic form requires arithmetic operations. This looks worse than it is in practice since the entries of are integers (–1, 0, and 1) in the absence of fractional imputation. This simplification allows one to avoid a fair amount of arithmetic. Assembling the remaining parts of the score statistic requires arithmetic operations.

Individuals missing univariate trait values are omitted from analysis. Individuals missing some but not all components of a multivariate trait are retained in analysis. The proper adjustments for missing data are made automatically in the score statistic because sections of Gaussian random vectors are Gaussian.

SNPs with minor allele counts below a user-designated threshold are also omitted from analysis. Note that if the minor allele count across a study is 0, then the given SNP is mono-allelic and worthless in association testing. Mendel’s default threshold of 3 is motivated by the rule of thumb in contingency table testing that all cells have an expected count of at least 3. For a multivariate trait, a SNP may fall below the threshold for some component traits but not for others. This situation can occur when each trait displays a different pattern of missing data across individuals. Mendel retains such anomalous SNPs only for those component traits with a sufficient number of minor alleles. Again, proper adjustments are made automatically within the score test statistic to account for partial data.

Mendel’s analysis yields a score test p-value for each SNP. For the user-designated most significant SNPs, Mendel’s subsequent likelihood ratio test outputs an estimated SNP effect size, a standard error of that estimate, and the fraction of the total variance explained by that SNP. For a multivariate trait, Mendel outputs a SNP effect size and associated standard error for each component trait. In the initial analysis under the null model with no SNPs, Mendel provides estimates with standard errors of all mean and variance components included in the model. Finally, an estimate of heritability with standard error is also provided.

The extension of the score test to the multivariate -distribution is straightforward ^{18}. Suppose equals the degrees of freedom of the -distribution and equals the number of observed person-trait combinations for pedigree . The sections of the score and expected information pertinent to the mean components for the pedigree reduce to

where is the residual and is the associated Mahalanobis distance. A sensible choice for is its estimate under the null model.

### 2.3 Kinship Estimation From SNPs

Mendel can either calculate the global kinship coefficient matrix from the provided pedigree structures or estimate it from dense genotypes. In global kinship estimation Mendel’s default uses an evenly spaced 20% of the available SNPs, and only compares pairs of individuals within defined pedigrees. Hence, is block diagonal. Users can trivially elect to exploit a larger fraction of the available SNPs or estimate kinship for all pairs of individuals. Given selected SNPs, Mendel estimates the global kinship coefficient of individuals and based on either the genetic relation matrix (GRM) method

or the method of moments (MoM) ^{7}, ^{19}

where is the minor allele frequency at SNP , is the number of minor alleles in ’s genotype at SNP , and

is the observed fraction of alleles identical by state (IBS) between and . The GRM method is Mendel’s default. In general, one can think of the GRM method centering and scaling each genotype, while the MoM method uses the raw genotypes and then centers and scales the final result.

### 2.4 Other Utilities for Handling Pedigree Data

To encourage thorough testing of new statistical methods, such as the current Ped-GWAS score test, we have implemented both genotype and trait simulation in our genetic analysis program Mendel ^{20}. Mendel does genotype simulation (gene dropping) subject to prescribed allele frequencies, a given genetic map, and Hardy-Weinberg and linkage equilibrium. If one fixes founder haplotypes and simulates conditional on these, then the unrealistic assumption of linkage equilibrium can be relaxed. Missing data patterns are respected or imposed by the user. It is also possible to set the rate for randomly deleting data and to simulate genotypes for people of mixed ethnicity by defining different ancestral populations, each with its own allele frequencies. If this feature is invoked, then each pedigree founder should be assigned to a population.

Trait simulation can be layered on top of genotype simulation. Mendel simulates either univariate traits determined by generalized linear models or multivariate Gaussian traits determined by variance component models. The biggest limitations are the restriction to a single major locus and the generalized linear model assumption that trait correlations are driven solely by this locus. Variance component models enable inclusion of environmental effects and more complicated correlations among relatives. In the variance component setting, univariate as well as multivariate Gaussian traits can be simulated. Most variance component models are built on Gaussian distributions, but Mendel allows one to replace these by multivariate -distributions. Thus, users can investigate robust statistics less prone to distortion by outliers. More theoretical and implementation details appear in the Mendel documentation ^{20}.

## 3 Results

### 3.1 The San Antonio Family Heart Study

We analyzed a real data set collected by the San Antonio Family Heart Study (SAFHS) ^{25}. The data consist of 3637 individuals in 211 Mexican American families. High-density lipoprotein (HDL) levels were measured at up to three time points for each of the 1429 phenotyped individuals. These traits are denoted HDL, HDL, and HDL, measured at corresponding ages AGE, AGE, and AGE. Some of the phenotyped individuals have HDL measurements at only one or two of the time points. Of the 1429 phenotyped individuals, 1413 were genotyped at 944,427 genome-wide SNPs. The genotyping success rate exceeded 98% in 1388 of these individuals over 124 pedigrees. The largest family contains 247 individuals; five others also contain more than 90 individuals. The smallest pedigree was a singleton. Genotyping success rates were above 98% for 935,392 SNPs.

### 3.2 Comparison with FaST-LMM and Gemma

Trait | SNP | Chr. | Base Pair | MAF | ||||
---|---|---|---|---|---|---|---|---|

Position | Mendel default | Mendel all-pairs | FaST-LMM | GEMMA | ||||

rs7303112 | 12 | 97,596,023 | 0.00455 | 10.21 | 10.71 | 7.63 | 7.24 | |

HDL | rs8040647 | 15 | 32,304,988 | 0.00147* | 7.44 | 7.56 | 7.35 | 7.45 |

rs9972594 | 15 | 32,421,102 | 0.00147* | 7.44 | 7.56 | 7.37 | 7.46 | |

rs7167103 | 15 | 32,830,477 | 0.00147* | 7.44 | 7.56 | 7.35 | 7.44 | |

HDL | rs7100957 | 10 | 28,207,332 | 0.00183* | 8.84 | 8.95 | 8.88 | 8.82 |

HDL | rs17060933 | 8 | 22,510,029 | 0.00382 | 8.23 | 8.28 | 8.61 | 8.59 |

rs7303112 | 12 | 97,596,023 | 0.00644 | 9.89 | 9.94 | |||

HDL | rs16925210 | 10 | 25,308,103 | 0.00217 | 8.15 | 8.33 | ||

with | rs7091416 | 10 | 25,318,381 | 0.00217 | 8.15 | 8.33 | Not | Not |

constrained | rs10075658 | 5 | 148,911,957 | 0.00144* | 8.16 | 8.21 | Available | Available |

covariates | rs7733139 | 5 | 145,977,990 | 0.00217 | 7.36 | 7.34 | ||

rs7100957 | 10 | 28,207,332 | 0.00870 | 7.20 | 7.30 | |||

rs7303112 | 12 | 97,596,023 | 0.00644 | 9.82 | 9.88 | 11.08 | ||

HDL | rs16925210 | 10 | 25,308,103 | 0.00217 | 8.04 | 8.23 | 3.53 | |

without | rs7091416 | 10 | 25,318,381 | 0.00217 | 8.04 | 8.23 | Not | 3.52 |

constrained | rs10075658 | 5 | 148,911,957 | 0.00144* | 8.12 | 8.17 | Available | 3.47 |

covariates | rs7733139 | 5 | 145,977,990 | 0.00217 | 7.41 | 7.40 | 3.47 | |

rs7100957 | 10 | 28,207,332 | 0.00870 | 7.19 | 7.30 | 4.48 | ||

rs10083226 | 13 | 104,434,452 | 0.00219 | 7.10 | 7.31 | 2.14 |

For fair comparisons, we directed Mendel to estimate SNP-based global kinship coefficients for all pairs of individuals ignoring the input pedigrees. This is the default in FaST-LMM and GEMMA. In addition, we ran Mendel’s default in which the coefficients are estimated only for pairs of individuals within the same input pedigree. We also slightly adjusted some of the default quality control thresholds so the programs would be analyzing roughly the same set of SNPs and individuals. For example, by default Mendel filters SNPs with fewer than three occurrences of the minor allele in the data; in contrast, FaST-LMM only filters SNPs with zero occurrences of the minor allele, and GEMMA filters SNPs with minor allele frequency (MAF) . All other defaults were observed throughout. Users can easily adjust the Mendel analysis parameters via its control file and the FaST-LMM and GEMMA analysis parameters via their command line.

We first carried out three univariate QTL analyses of HDL, HDL, and HDL, using SEX and AGE, AGE, or AGE as covariates. We then ran a multivariate QTL analysis of HDL, HDL, and HDL jointly, which we refer to as HDL. For the multivariate analysis, the most appropriate configuration is to constrain the effects of the SEX and AGE covariates to be the same on all three HDL measurements. Such linear constraints are imposed in Mendel via a few simple lines in its control file. FaST-LMM and GEMMA do not allow constraints on covariates. Therefore, we also ran a multivariate analysis with only the SEX covariate and no constraints. With no constraints, SEX will have a slightly different effect on each component phenotype in the multivariate analysis. For instance, Mendel’s default run estimated a female effect of on HDL, on HDL, and on HDL. FaST-LMM cannot perform any multivariate analyses.

Table 3 reports all SNPs with MAF that achieve genome-wide significance (p-values less than ) as reported by at least one software package. For the univariate analyses, each software package found the same set of significant SNPs, except that one of GEMMA’s p-values was slightly short of the significance threshold. Figure 1 shows a Manhattan plot and a Q-Q plot from the HDL analysis by Mendel given kinship estimates for all pairs of individuals. The results for the other analyses, both univariate and multivariate, were similar. Each Mendel all-pairs univariate analysis had genomic control in the range 1.002 to 1.006; in default mode, was in the range 0.992 to 1.022. The various Q-Q plots and associated values show there is no systematic biases in the data or analysis. In the all-pairs Mendel HDL analysis, the grand mean (intercept) was 49.0 0.8. The SEX covariate was significant in all null models. For example, in the all-pairs Mendel HDL analysis with constrained covariates, the SEX effect was 2.4 0.3 for females and, by design, the opposite for males. The AGE covariate was not significant in any run. For example, again in the all-pairs HDL analysis with parameter constraints, the AGE effect was 0.04 0.02. In the null model for the all-pairs Mendel HDL analysis, the additive variance was estimated as 78.8 9.9, and the environmental variance was estimated as 78.1 7.2. This gives an overall heritability estimate for HDL of 0.50 0.04. Similar variance estimates were seen in other null models.

For the multivariate analysis without parameter constraints, Mendel is able to include almost twice as many individuals in the analysis as GEMMA (see Table 4). GEMMA only includes individuals phenotyped at all component traits and covariates. This probably explains why Mendel finds several more SNPs with significant p-values than GEMMA.

Program | Trait Analyzed | Analyzed | RunTime | RAM | |
---|---|---|---|---|---|

Samples | SNPs | (min:sec) | (GB) | ||

Mendel default | 1357 | 935,392 | 1:51 | 1.2 | |

Mendel all-pairs | HDL | 1357 | 935,392 | 7:49 | 1.2 |

FaST-LMM | 1397 | 941,546 | 76:11 | 30.0 | |

GEMMA | 1397 | 919,050 | 206:54 | 0.4 | |

Mendel default | 818 | 935,392 | 1:33 | 1.1 | |

Mendel all-pairs | HDL | 818 | 935,392 | 3:25 | 1.1 |

FaST-LMM | 840 | 934,216 | 49:44 | 18.0 | |

GEMMA | 840 | 914,051 | 180:21 | 0.3 | |

Mendel default | 914 | 935,392 | 1:38 | 1.1 | |

Mendel all-pairs | HDL | 914 | 935,392 | 3:54 | 1.1 |

FaST-LMM | 939 | 937,208 | 54:58 | 20.0 | |

GEMMA | 939 | 918,626 | 182:26 | 0.3 | |

Mendel default | HDL | 1388 | 935,392 | 4:08 | 1.2 |

Mendel all-pairs | with | 1388 | 935,392 | 83:24 | 1.2 |

FaST-LMM | constrained | Not Available | |||

GEMMA | covariates | Not Available | |||

Mendel default | HDL | 1388 | 935,392 | 3:49 | 1.2 |

Mendel all-pairs | without | 1388 | 935,392 | 80:04 | 1.2 |

FaST-LMM | constrained | Not Available | |||

GEMMA | covariates | 712 | 912,318 | 630:37 | 0.6 |

Table 4 tallies the run times and memory footprints from each analysis on a typical personal computer with adequate RAM to accommodate FaST-LMM (6 CPU cores at 2.67 GHz, with 48 GB total RAM). Even when estimating the global kinship coefficients for all pairs of individuals, each univariate QTL run took Mendel less than 8 minutes to read, quality check, and analyze the data for kinship estimates and association tests, roughly 10% of the time required for FaST-LMM and 5% of the time required by GEMMA. (For GEMMA, the kinship estimation and association tests are run separately. The run times reported here are their total.)

The three programs use different association test strategies: Mendel performs score tests for all SNPs and LRTs for the top SNPs; FaST-LMM performs LRTs; and GEMMA by default performs Wald tests, but the user can change this to LRTs or score tests. For the univariate analyses on a six-core computer, excluding estimation of kinship coefficients, GEMMA’s run times under the Wald test and LRT options were roughly similar to FaST-LMM’s; GEMMA’s run time under the score test option was roughly double Mendel’s in all-pairs mode. This is impressive given GEMMA’s lack of multithreading. It is kinship estimation, which in practice can be done once per data set, that is substantially slower in GEMMA (running roughly 135 minutes) than in FaST-LMM or Mendel (less than 1 minute).

Each trivariate QTL run took Mendel less than 90 minutes. Mendel required roughly one-eighth the time of GEMMA while analyzing almost twice as many individuals. Mendel is also memory efficient. The univariate and multivariate runs each required less than 1.5 GB of memory, which is well below the amount of RAM in a typical computer. FaST-LMM’s memory usage is more than 15 times larger than Mendel’s. GEMMA uses even less memory than Mendel but is considerably slower.

## 4 Discussion

We have implemented an ultra-fast algorithm for QTL analysis of pedigree data or mix of population and pedigree data. In our opinion Mendel’s comprehensive environment for genetic data analysis is a decided advantage. In addition to its exceptional speed and memory efficiency, Mendel can handle multivariate quantitative traits and detect outlier trait values and pedigrees. Most competing programs ignore multivariate traits and outliers altogether.

A recent review of univariate QTL analysis packages for family data ^{8} shows that all the explored packages obtain similar results, leaving speed, features, and ease of use as the important factors in choosing between them. Once the current version of Mendel came out, the authors of the review were kind enough to add a comment (http://www.plosgenetics.org/annotation/listThread.action?root=81847www.plosgenetics.org/annotation/listThread.action?root=81847) to their article observing that Mendel was now the fastest and one of the easiest to use packages they reviewed.

In the SAFHS example data set we used with HDL phenotypes, all the significant SNPs we found had MAF . Due to these low MAFs, we do not claim these SNPs are strong candidates for further study. However, the key point here is that all four methods found the same SNPs, at least for the univariate analyses. We also note that the p-values are quite similar regardless of whether one uses kinship estimates between all individuals (Mendel’s all-pairs mode) or only between individuals within the same input pedigrees (Mendel’s default mode). This suggests that the input pedigree structures for this data set are substantially correct and complete, with few mistaken or hidden relationships. Obviously, this may not be true for other data sets. By supplying good kinship estimates ignoring pedigree structures, the currently reviewed packages make the hard fieldwork of relationship discovery superfluous.

A future version of Mendel will address its failure to read fractional genotype values. This is simply a logistical issue, as all Mendel’s internal genotype computations are already handled as floating point operations. Another imminent feature is a fourth style of kinship coefficient estimation that allows the user to force theoretical kinship coefficients for pairs of individuals within the same pedigree and estimated kinships for all other pairs.

By supplying a comprehensive, fast, and easy to use package for GWAS on quantitative traits in general pedigrees, we hope to encourage exploitation of family-based data sets for gene mapping. A gene mapping study should collect as large a sample as possible consistent with economic constraints and uniform trait phenotyping. If the sample includes pedigrees, all the better. One should not let the choice of statistical test determine the data collected; on the contrary, the data should determine the test. Here we have argued that score tests can efficiently handle unrelated individuals, pedigrees, or a mixture of both. For human studies, where controlled breeding is forbidden, nature has provided pedigrees segregating every genetic trait. Many of these pedigrees are known from earlier linkage era studies and should be treasured as valuable resources.

Let us suggest a few directions for future work. The current method works marker by marker and is ill equipped to perform model selection. Lasso penalized regression is available to handle model selection for case-control and random sample data ^{40}, ^{39}, ^{43}, ^{42} and can be generalized to variance component models. Although we have generalized the score test to distributions such as the multivariate , extending it to discrete traits may be out of reach. For likelihood based methods, there simply are no discrete analogues of the Gaussian distribution that lend themselves to graceful evaluation of pedigree likelihoods. Treating case/control data as a 0/1 quantitative variable is a possibility that has been explored by Pirinen et al. ^{27}. The GEE method is another fallback option because it does not depend on precise distributional assumptions.

In rare variant mapping, grouping related SNPs in a variance component may be a good alternative to the mean component models used here. Each variant may be too rare to achieve significance in hypothesis testing. Fortunately, aggregating genotype information within biological units such as genes or pathways offer better power than marginal testing of individual SNPs. See Asimit and Zeggini ^{2} for a recent review of aggregation strategies. Kwee et al. ^{14} have successfully applied a variance component model for association testing of SNP sets in a sample of unrelated subjects. Rönnegård et al. ^{29} consider score tests for random effects models in the context of experimental line crosses. Score tests may well be the key to implementing random effect models in pedigrees. However, the computational demands are apt to be more formidable than those encountered here with fixed effects models. In particular, if tests are based simply on local identity-by-descent (IBD) sharing, then the boundaries between pedigrees disappear, and the entire sample collapses to one large pedigree. The required local kinship coefficients can again be well estimated from dense markers, but this demands more computation than the estimation of global kinship coefficients under the mean components model advocated here ^{7}. Since inversion of a pedigree covariance matrix scales as the cube of the number of individuals in the pedigree, treating the entire sample as a single pedigree will put a practical upper limit on sample size. There are other issues in implementing variance component models such as assigning p-values and dealing with multivariate traits that are best left to a separate paper.

## Acknowledgement

#### Funding:

NIH (GM053275, HG006139, MH059490) and NSF (DMS-1310319). Burroughs Wellcome Fund Inter-school Training Program in Metabolic Diseases (KHKC).

## References

- Amin et al. 2007 Amin, N., van Duijn, C. M., and Aulchenko, Y. S. (2007). A genomic background based method for association analysis in related individuals. PLoS ONE, 2(12):e1274.
- Asimit and Zeggini 2010 Asimit, J. and Zeggini, E. (2010). Rare variant association analysis methods for complex traits. Annual Review of Genetics, 44(1):293–308.
- Aulchenko et al. 2007 Aulchenko, Y. S., de Koning, D.-J., and Haley, C. (2007). Genomewide rapid association using mixed model and regression: A fast and simple method for genomewide pedigree-based quantitative trait loci association analysis. Genetics, 177(1):577–585.
- Chen et al. 2011 Chen, M.-H., Liu, X., Wei, F., Larson, M. G., Fox, C. S., Vasan, R. S., and Yang, Q. (2011). A comparison of strategies for analyzing dichotomous outcomes in genome-wide association studies with general pedigrees. Genetic Epidemiology, 35(7):650–657.
- Chen and Yang 2010 Chen, M.-H. and Yang, Q. (2010). GWAF: an R package for genome-wide association analyses with family data. Bioinformatics, 26(4):580–581.
- Chen and Abecasis 2007 Chen, W.-M. and Abecasis, G. R. (2007). Family-based association tests for genome-wide association scans. American Journal of Human Genetics, 81(5):913–926.
- Day-Williams et al. 2011 Day-Williams, A. G., Blangero, J., Dyer, T. D., Lange, K., and Sobel, E. M. (2011). Linkage analysis without defined pedigrees. Genetic Epidemiology, 35(5):360–370.
- Eu-ahsunthornwattana et al. 2014 Eu-ahsunthornwattana, J., Miller, E. N., Fakiola, M., Wellcome Trust Case Control Consortium 2, Jeronimo, S. M. B., Blackwell, J. M., and Cordell, H. J. (2014). Comparison of methods to account for relatedness in genome-wide association studies with family-based data. PLoS Genetics, 10(7):e1004445.
- Hindorff et al. 2009 Hindorff, L. A., Sethupathy, P., Junkins, H. A., Ramos, E. M., Mehta, J. P., Collins, F. S., and Manolio, T. A. (2009). Potential etiologic and functional implications of genome-wide association loci for human diseases and traits. Proceedings of the National Academy of Sciences, 106(23):9362–9367.
- Kang et al. 2010 Kang, H. M., Sul, J. H., Service, S. K., Zaitlen, N. A., Kong, S.-Y., Freimer, N. B., Sabatti, C., and Eskin, E. (2010). Variance component model to account for sample structure in genome-wide association studies. Nature Genetics, 42(4):348–354.
- Kang et al. 2008 Kang, H. M., Zaitlen, N. A., Wade, C. M., Kirby, A., Heckerman, D., Daly, M. J., and Eskin, E. (2008). Efficient control of population structure in model organism association mapping. Genetics, 178:1709–1723.
- Korte et al. 2012 Korte, A., Vilhjalmsson, B. J., Segura, V., Platt, A., Long, Q., and Nordborg, M. (2012). A mixed-model approach for genome-wide association studies of correlated traits in structured populations. Nature Genetics, 44(9):1066–1071.
- Ku et al. 2010 Ku, C. S., Loy, E. Y., Pawitan, Y., and Chia, K. S. (2010). The pursuit of genome-wide association studies: where are we now? Journal of Human Genetics, 55(4):195–206.
- Kwee et al. 2008 Kwee, L. C., Liu, D., Lin, X., Ghosh, D., and Epstein, M. P. (2008). A powerful and flexible multilocus association test for quantitative traits. American Journal of Human Genetics, 82(2):386–397.
- Laird et al. 2000 Laird, N. M., Horvath, S., and Xu, X. (2000). Implementing a unified approach to family-based tests of association. Genetic Epidemiology, 19(Suppl 1):S36–S42.
- Lange and Laird 2002 Lange, C. and Laird, N. M. (2002). On a general class of conditional tests for family-based association studies in genetics: the asymptotic distribution, the conditional power, and optimality considerations. Genetic Epidemiology, 23(2):165–180.
- Lange 2002 Lange, K. (2002). Mathematical and Statistical Methods for Genetic Analysis. Statistics for Biology and Health. Springer-Verlag, New York, 2nd edition.
- Lange et al. 1989 Lange, K., Little, R. J. A., and Taylor, J. M. G. (1989). Robust statistical modeling using the distribution. Journal of the American Statistical Association, 84(408):881–896.
- Lange et al. 2014 Lange, K., Papp, J. C., Sinsheimer, J. S., and Sobel, E. M. (2014). Next-generation statistical genetics: Modeling, penalization, and optimization in high-dimensional data. Annual Review of Statistics and Its Application, 1(1):279–300.
- Lange et al. 2013 Lange, K., Papp, J. C., Sinsheimer, J. S., Sripracha, R., Zhou, H., and Sobel, E. M. (2013). Mendel: The Swiss army knife of genetic analysis programs. Bioinformatics, 29(12):1568–1570.
- Lange et al. 2005 Lange, K., Sinsheimer, J. S., and Sobel, E. (2005). Association testing with Mendel. Genetic Epidemiology, 29(1):36–50.
- Lange et al. 1976 Lange, K., Westlake, J., and Spence, M. A. (1976). Extensions to pedigree analysis iii. variance components by the scoring method. Annals of Human Genetics, 39(4):485–491.
- Lippert et al. 2011 Lippert, C., Listgarten, J., Liu, Y., Kadie, C. M., Davidson, R. I., and Heckerman, D. (2011). FaST linear mixed models for genome-wide association studies. Nature Methods, 8(10):833–835.
- Listgarten et al. 2012 Listgarten, J., Lippert, C., Kadie, C. M., Davidson, R. I., Eskin, E., and Heckerman, D. (2012). Improved linear mixed models for genome-wide association studies. Nature Methods, 9(6):525–526.
- Mitchell et al. 1996 Mitchell, B. D., Kammerer, C. M., Blangero, J., Mahaney, M. C., Rainwater, D. L., Dyke, B., Hixson, J. E., Henkel, R. D., Sharp, R. M., Comuzzie, A. G., VandeBerg, J. L., Stern, M. P., and MacCluer, J. W. (1996). Genetic and environmental contributions to cardiovascular risk factors in Mexican Americans: The San Antonio Family Heart Study. Circulation, 94(9):2159–2170.
- Ott et al. 2011 Ott, J., Kamatani, Y., and Lathrop, M. (2011). Family-based designs for genome-wide association studies. Nature Reviews Genetics, 12(7):465–474.
- Pirinen et al. 2013 Pirinen, M., Donnelly, P., and Spencer, C. C. A. (2013). Efficient computation with a linear mixed model on large-scale data sets with applications to genetic studies. Annals of Applied Statistics, 7(1):369–390.
- Rao 2009 Rao, C. (2009). Linear Statistical Inference And Its Applications. Wiley, 2nd edition.
- Rönnegård et al. 2008 Rönnegård, L., Besnier, F., and Carlborg, O. (2008). An improved method for quantitative trait loci detection and identification of within-line segregation in f2 intercross designs. Genetics, 178(4):2315–2326.
- Spielman and Ewens 1998 Spielman, R. S. and Ewens, W. J. (1998). A sibship test for linkage in the presence of association: the sib transmission/disequilibrium test. American Journal of Human Genetics, 62(2):450–458.
- Stanhope and Abney 2012 Stanhope, S. A. and Abney, M. (2012). GLOGS: a fast and powerful method for GWAS of binary traits with risk covariates in related populations. Bioinformatics, 28(11):1553–1554.
- Thornton and McPeek 2007 Thornton, T. and McPeek, M. S. (2007). Case-control association testing with related individuals: A more powerful quasi-likelihood score test. American Journal of Human Genetics, 81(2):321–337.
- Thornton and McPeek 2010 Thornton, T. and McPeek, M. S. (2010). ROADTRIPS: Case-control association testing with partially or completely unknown population and pedigree structure. American Journal of Human Genetics, 86(2):172–184.
- Van Steen 2011 Van Steen, K. (2011). Perspectives on genome-wide multi-stage family-based association studies. Statistics in Medicine, 30(18):2201–2221.
- Van Steen and Lange 2005 Van Steen, K. and Lange, C. (2005). PBAT: a comprehensive software package for genome-wide association analysis of complex family-based studies. Human Genomics, 2(1):67–69.
- Visscher et al. 2012 Visscher, P. M., Brown, M. A., McCarthy, M. I., and Yang, J. (2012). Five years of GWAS discovery. American Journal of Human Genetics, 90(1):7–24.
- Won et al. 2009a Won, S., Bertram, L., Becker, D., Tanzi, R., and Lange, C. (2009a). Maximizing the power of genome-wide association studies: A novel class of powerful family-based association tests. Statistics in Biosciences, 1(2):125–143.
- Won et al. 2009b Won, S., Wilk, J. B., Mathias, R. A., O’Donnell, C. J., Silverman, E. K., Barnes, K., O’Connor, G. T., Weiss, S. T., and Lange, C. (2009b). On the analysis of genome-wide association studies in family-based designs: A universal, robust analysis approach and an application to four genome-wide association studies. PLoS Genetics, 5(11):e1000741.
- Wu et al. 2009 Wu, T. T., Chen, Y., Hastie, T., Sobel, E. M., and Lange, K. (2009). Genome-wide association analysis by lasso penalized logistic regression. Bioinformatics, 25(6):714–721.
- Wu and Lange 2008 Wu, T. T. and Lange, K. (2008). Coordinate descent algorithms for lasso penalized regression. Annals of Applied Statistics, 2(1):224–244.
- Zhang et al. 2010 Zhang, Z., Ersoz, E., Lai, C.-Q., Todhunter, R. J., Tiwari, H. K., Gore, M. A., Bradbury, P. J., Yu, J., Arnett, D. K., Ordovas, J. M., and Buckler, E. S. (2010). Mixed linear model approach adapted for genome-wide association studies. Nature Genetics, 42(4):355–360.
- Zhou et al. 2011 Zhou, H., Alexander, D., Sehl, M., Sinsheimer, J. S., Sobel, E. M., and Lange, K. (2011). Penalized regression for genome-wide association screening of sequence data. Pacific Symposium on Biocomputing, 2011:106–117.
- Zhou et al. 2010 Zhou, H., Sehl, M. E., Sinsheimer, J. S., and Lange, K. (2010). Association screening of common and rare genetic variants by penalized regression. Bioinformatics, 26(19):2375–2382.
- Zhou and Stephens 2012 Zhou, X. and Stephens, M. (2012). Genome-wide efficient mixed-model analysis for association studies. Nature Genetics, 44:821–824.
- Zhou and Stephens 2014 Zhou, X. and Stephens, M. (2014). Efficient multivariate linear mixed model algorithms for genome-wide association studies. Nature Methods, 11(4):407–409.
- Zhu and Xiong 2012 Zhu, Y. and Xiong, M. (2012). Family-based association studies for next-generation sequencing. American Journal of Human Genetics, 90(6):1028–1045.