Genome-wide QTL and eQTL analyses using Mendel

Pedigree genome-wide association studies (GWAS) (Option 29) in the current version of the Mendel software is an optimized subroutine for performing large-scale genome-wide quantitative trait locus (QTL) analysis. This analysis (a) works for random sample data, pedigree data, or a mix of both; (b) is highly efficient in both run time and memory requirement; (c) accommodates both univariate and multivariate traits; (d) works for autosomal and x-linked loci; (e) correctly deals with missing data in traits, covariates, and genotypes; (f) allows for covariate adjustment and constraints among parameters; (g) uses either theoretical or single nucleotide polymorphism (SNP)–based empirical kinship matrix for additive polygenic effects; (h) allows extra variance components such as dominant polygenic effects and household effects; (i) detects and reports outlier individuals and pedigrees; and (j) allows for robust estimation via the t-distribution. This paper assesses these capabilities on the genetics analysis workshop 19 (GAW19) sequencing data. We analyzed simulated and real phenotypes for both family and random sample data sets. For instance, when jointly testing the 8 longitudinally measured systolic blood pressure and diastolic blood pressure traits, it takes Mendel 78 min on a standard laptop computer to read, quality check, and analyze a data set with 849 individuals and 8.3 million SNPs. Genome-wide expression QTL analysis of 20,643 expression traits on 641 individuals with 8.3 million SNPs takes 30 h using 20 parallel runs on a cluster. Mendel is freely available at http://www.genetics.ucla.edu/software.


Background
The classical variance component model has been a powerful tool for mapping quantitative trait loci (QTLs) in pedigrees. Polygenic effects are effectively modeled by introducing an additive genetic variance component operating on the kinship coefficient matrix. With unknown or dubious pedigree structure, global kinship coefficients can be accurately estimated from dense markers using either the genetic relationship matrix (GRM) or the method of moments. In GWAS (genome-wide association studies), the 2 alleles of a SNP (single nucleotide polymorphism) shift trait means and can be tested as a fixed effect. However, fitting a variance component model is computationally challenging, especially when it has to be done for a large number of markers. In the newly released version of the Mendel software [1], Option 29 implements an ultrafast score test for pedigree GWAS. Score tests require no additional iteration under the alternative model. Only SNPs with the most promising score-test p values are further subject to likelihood ratio testing, thus achieving a good compromise between speed and power for large-scale QTL analysis. In this paper, we demonstrate the capabilities of Mendel on the Genetic Analysis Workshop 19 (GAW19) sequencing data.

Methods
QTL association mapping typically invokes the multivariate normal distribution to model the observed Tvariate trait Y ∈ ℝ n × T over a pedigree of n individuals. The standard model [2] collects the means of the responses vec(Y) into a vector v and the corresponding covariances into a matrix Σ and represents the loglikelihood of a pedigree as where the covariance matrix is typically parametrized as Here Φ is the global kinship matrix capturing additive polygenic effects, and Δ 7 is a condensed identity coefficient matrix capturing dominance genetic effects. For Φ, Mendel can use (a) the theoretical kinship matrix from provided pedigree structures; (b) SNP-based estimates for the kinship of pairs of people within each pedigree; or (c) SNP-based estimates for the entire global kinship matrix ignoring pedigree information. To estimate kinship coefficients from dense SNP data, Mendel employs either the GRM or the method of moments [3,4]. The household effect matrix H has entries h ij = 1 if individuals i and j are in the same household and 0 otherwise. Individual environmental contributions and trait measurement errors are incorporated via the identity matrix I. QTL fixed effects are captured through the mean component υ = Aβ for some predictor matrix A and vector of regression coefficients β. To test a SNP against a Tvariate trait, A is augmented with T extra columns holding the allele counts at the SNP, and the corresponding regression coefficients are jointly tested for association [5]. For longitudinal measurements of covariates such as smoke, age, and blood pressure medication (BPmed), we may either assume time varying effect sizes or constrain their effect sizes at different time points to be the same. The latter tactic leads to a more parsimonious and interpretable model and can be easily enforced by setting appropriate parameter constraints in Mendel's control file, which lists the user's choice of model parameters. In Mendel, SNPs with the most impressive test score p values (top 10 by default) are further tested by the more accurate, but slower, likelihood ratio method, thus achieving a good compromise between speed and power for large-scale QTL analysis. We refer readers to our companion manuscript [6] for more model and implementation details.

Results and discussion
Family data Size and power study using simulated traits (SIMPHEN.  The power to detect the 6 functional variants in the MAP4 gene on chromosome 3 is evaluated from the 200 simulation replicates of the trivariate traits systolic blood pressure (SBP) and diastolic blood pressure (DBP). Type I errors are evaluated from the univariate Q1 trait, which does not involve a major gene. Our analysis includes covariates sex, age, BPmed, smoke, and their pairwise interactions, and uses the theoretical kinship matrix as the additive polygenic variance component. We constrain the covariate effects to be equal across 3 time points. To read in all the data and run standard quality control (QC) procedures took just under 5 min. QC excluded 10,603 SNPs and 110 individuals based on genotyping success rates below 98 %. The remaining 8,338,071 SNPs and 849 individuals were analyzed. The subsequent ped-GWAS analysis ran in 73 min for all results reported in Fig. 1 and Tables 2 and 3. Because we excluded rare SNPs with low minor allele frequencies <0.03 across 849 individuals, p values were calculated for only 3,084,046 SNPs. Accordingly the genome-wide significance threshold is 1.62 × 10 −8 or 7.79 on the log 10 scale; the threshold for a false discovery rate (FDR) of 0.05 is 4.19 × 10 −8 or 7.38 on the log 10 scale. Table 2 lists the estimates for environmental effects and their interactions under the null model (no SNPs included). Figure 1 displays the Manhattan and quantilequantile (Q-Q) plots. The genomic inflation factor of 1.023 indicates no systematic bias. One SNP passes the Bonferroni-corrected genome-wide significance level, and 3 SNPs pass the FDR 0.05 threshold. They are listed in Table 3.

Unrelated data
A second data set consists of exome sequence calls, blood pressure phenotypes at a single time point, and simulated phenotypes on a large set of unrelated individuals. Like the family data set, these individuals are Mexican Americans; however, they were independently ascertained and do not overlap with the family data set.
Size and power study using simulated traits (SIMPHEN.  The data set provides 200 simulation replicates of the trait SBPs and DBP. However, GAW19 organizers did not distribute the exact simulation details, except to state that "The set of causal variants is somewhat different since this is exome data rather than the full sequence data that was provided last time, and so not all of the GAW18 variants, regulatory ones in particular, are present in the new data set." This precludes a precise size and power study. For ease of comparison, we tested 5 of the 6 variants displayed in Table 1 (for family data) against the bivariate trait (SBP, DBP) for all 200 simulation replicates and report the rejection rates in Table 4. In the model, we include covariates sex, age, BPmed, smoke, and their pairwise interactions, and use the SNP-based   Fig. 3 and Table 5. Because we exclude rare variants with a MAF equal to or less than 0.01 in all individuals, p values were calculated for 52,314 SNPs. Accordingly, the genome-wide significance threshold is 9.56 × 10 −7 or 6.02 on the log 10 scale. Estimated environmental effects and their interactions and variance components under the null model (no SNPs included) are listed in Table 5. Figure 3 displays the Manhattan and Q-Q plots. The genomic inflation factor of 1.001 indicates no systematic bias. No SNPs pass the genome-wide significance level or FDR 0.05 threshold.

Conclusions
All analyses in this article use Mendel v14.3, which is freely available at http://www.genetics.ucla.edu/software. Pedigree GWAS (Option 29) in Mendel proves to be an extremely efficient and versatile implementation for large-scale QTL analysis. Most competing programs ignore multivariate traits and outliers altogether. See Zhou et al [6] for a side-by-side comparison with the Factored Spectrally Transformed Linear Mixed Model (FaST-LMM) [7] and GEMMA (Genome-wide Efficient Mixed-Model Analysis) [8] programs. Here we have emphasized Mendel's flexibility in specifying the global kinship matrix, adjusting for confounding, and capturing interactions. These assets, plus its raw speed, make it an ideal environment for QTL mapping. Mendel continues to mature, and geneticists are advised to give it a second look for genetic analysis [1]. In rare variant mapping, each variant may be too rare to achieve significance in hypothesis testing. Grouping related SNPs in a variance component may be more powerful than the mean component models used here. Extending Mendel to test variance component is among the focuses of our current work.