Estimating proportions of explained variance: a comparison of whole genome subsets

Following the publication of the ENCODE project results, there has been increasing interest in investigating different areas of the chromosome and evaluating the relative contribution of each area to expressed phenotypes. This study aims to evaluate the contribution of variants, classified by minor allele frequency and gene annotation, to the observed interindividual differences. In this study, we fitted Bayesian linear regression models to data from Genetic Analysis Workshop 18 (n = 395) to estimate the variance of standardized and log-transformed systolic blood pressure that can be explained by subsets of genetic markers. Rare and very rare variants explained an overall higher proportion of the variance, as did markers located within a gene rather than flanking regions. The proportion of variance explained by rare and very rare variants decreased when we controlled for the number of markers, suggesting that the number of contributing rare alleles plays an important role in the genetic architecture of chronic disease traits. Our findings lend support to the "common disease, rare variant" hypothesis for systolic blood pressure and highlight allele frequency and functional annotation of a polymorphism as potentially crucial considerations in whole genome study designs.


Background
The proportion of phenotypic variance explained by genetic factors is influenced by multiple variant attributes. First, an analysis of several complex traits by Yang et al showed that genic regions explain more variation than intergenic because causal variants are more likely to be located in or near the genes, particularly the protein-coding regions [1]. However, insights from the ENCODE project suggest that a number of such regions remain unidentified within the intergenic space and warrant further study [2]. Second, for traits that experience strong selection pressures, rare and very rare mutations (minor allele frequency [MAF] <5%) have been shown to contribute more variance than common mutations [3]. Genome-wide association studies, which rely on linkage disequilibrium between typed and untyped variants, are unlikely to detect effects of rare variants, which can partially explain the majority of variance, yet remain undetected [4]. Although there is clear empirical evidence that variants across the MAF spectrum are important to the genetic architecture of complex traits, there is no consensus on the relative contributions of rare and common polymorphisms to explained variance.
To quantify the relative contribution of subsets of variants (defined by either MAF or functional annotation) to genetic variance of systolic blood pressure (SBP) we fitted Bayesian linear regression models to the Genetic Analysis Workshop 18 (GAW18) data. To control for the effects of the number of variants in each of the subsets, we performed sensitivity analyses with a fixed number of markers in each category to determine the subset with the highest variance contribution per variant.

Phenotypes and covariates
SBP was first log transformed (logSBP) to control skewness and kurtosis. Subsequently, for ease of interpretation, logSBP was standardized to unit variance. To use the available longitudinal data on SBP and covariates and to maximize sample size, a cross-sectional data set was formed by selecting the first available data of 4 visits. Covariates recorded at the time of the visit used for SBP included age, gender, use of tobacco products, and use of medication for blood pressure.

Estimating variance explained by sets of genetic markers
The whole genome regression models [5] used for analysis were mixed-effects models of the following form: a regression on nongenetic covariates (sex, age, smoking, and blood pressure medication) whose effects are regarded as fixed; L l=1 x il β l is a regression on marker genotypes {x ij } whose effects {β l } are regarded as random; and ε i are independent, identically distributed normal residuals with mean equal to zero and variance σ 2 ε . Marker genotypes were expressed as deviation from the average genotype, given by 2 times the frequency of the allele coded as 1 at the corresponding marker. Intercept and fixed effects were assigned flat priors; marker effects were assigned independent, identically distributed normal priors with null mean and variance σ 2 β ; and, finally, variance parameters {σ 2 β σ 2 ε } were assigned weakly informative independent scaledinverted chi-square densities. The algorithms used to implement this model were fully described in prior publications from our group [6]. The software used for the analysis is available upon request. Inferences were based on 25,000 samples obtained after discarding 15,000 as burnin. Convergence was evaluated by visual inspection of trace plots.

Variant selection
With the goals of maximizing the biological plausibility of the analysis and reducing the requirements on computational power, we selected a set of 54,309 variants in regions of the genome enriched for blood pressure. The set of variants was identified as follows: We initially selected 84 genes implicated in pathways regulating blood pressure that were included on the SA Biosciences Human Hypertension PCR Array (Qiagen, Venlo, Netherlands). Of the 84 genes reported, only 31 were located on the odd chromosomes and thus available through the GAW18 data release; for our analysis, we selected all variants in or near these 31 genes ( ± 50 kilobases [kb] from upstream or downstream of the start or end sites of transcription, respectively). Table 1 summarizes the distribution of variants near each gene.
Subsequently, we excluded variants that did not pass quality control or had more than 2 alleles (n = 4470). Variants were further classified on the basis of functional annotation (genic, including intronic, vs. flanking) and allele frequency in the GAW18 data set (very rare if MAF <1%, rare if 1% ≥ MAF ≥5%, common if MAF >5%). To determine functional annotation, we used VCFtools [7] to concatenate all 11 odd-numbered chromosomes' variants and then converted the VCF file to ANNOVAR [8] input format. Using ANNOVAR, we then annotated the variants using Human Genome version 19 (hg19) as a reference, providing gene names (for genic markers) or flanking gene names plus distance to these genes (for markers located in the vicinity of a gene). Variants already reported in dbSNP 131 were also linked with their corresponding rsID numbers. shows the number of variants out of our total enriched set (p = 49,839) categorized by MAF and functionality. Based on these classifications, we identified 12 sets of markers and fitted whole genome regression models as described above to each set. Because of variation in the number of variants across the categories in Table 2, it was imperative to assess whether differences in the proportions of variance explained by regression on each of the marker sets were because of the nature of the marker set or the number of markers included in it. To circumvent this problem, we fitted models for each of the marker categories using a fixed number of markers (500) chosen at random from the original set. We fit each of these models 500 times, each time resampling the set of markers used. For each of the models, we quantified the proportion of variance accounted by the model using R 2 m = 1 −σ 2 εm whereσ 2 ε is the estimated residual variance of model m. We reported the average of the R 2 m values over 500 replicate runs of the model. Because the response was standardized to unit variance, R 2 m is interpretable as an R-squared statistic.  Tables 2 and 3). The distribution of variants by MAF was consistent with published observations from other populations [9]. Table 3 shows the estimates of proportion of variance explained when marker sets were defined by functional annotation and MAF. Percentage of variance explained decreased with MAF, regardless of functional annotation. Variants in genic regions generally explained more variance than those in the flanking regions, despite the smaller number of genic markers. Table 4 summarizes the results from the analysis controlling for the number of markers in each category (p = 500). Overall, estimates of explained variance slightly decreased compared to those obtained from models including all variants in each category, but the relative contributions of each class of markers remained similar, with rare and very rare variants located in genic regions explaining slightly more outcome variability than common flanking variants.

Discussion
Using sequence data from a set of 31 biologically relevant genomic regions, we established that the proportion of logSBP variance explained by genetic markers in hypertension-related pathways is roughly 11% of the phenotypic variation that remains after accounting for systematic, nongenetic effects. The proportion of variance explained by regression on markers decreases with the MAF, regardless of whether the variant is located within or near known genes. However, the observed trends across MAF categories were not statistically significant.
To interpret our findings, it is important to distinguish between estimating the proportion of variance explained by the markers, which is the focus of our article, and estimating heritability. Indeed, our estimates of proportion of variance explained by genetic factors (roughly 11% after accounting for differences as a result of nongenetic effects) are likely to be smaller than the true heritability of the trait because, as a result of imperfect linkage disequilibrium between alleles at markers and those at causal loci, some proportion of genetic variance is likely to have remained unaccounted. Our findings are consistent with evidence from several recent studies, which suggest that the genetic architecture of blood pressure is likely to involve a large number of modestly associated, and many yet undiscovered, variants [10][11][12].
The level of linkage disequilibrium between alleles at markers and those at causal loci depends on linkage disequilibrium decay and on marker density. Consequently, the proportion of variance captured by a maker set is related not only to the nature of the marker set and the genetic architecture of the trait, but also to the number of genetic markers in the set [13]. To account for the effects of the size of different marker sets, we performed an additional analysis, controlling for the number of predictors in each category, enabling direct comparisons between the contributions of very rare, rare, and common variants, as well as between those of flanking and genic markers. We found that the functional annotation of the variants has implications for the proportion of variance explained, especially at higher allele frequencies. However, the relatively small differences in variance explained suggests that flanking regions may harbor yet unidentified genes or regulatory elements that affect polygenic traits. Future studies, particularly of the ENCODE project data, will be informative in testing that hypothesis.
Although our variants were included in the models based on biological relevance, the family structure of the GAW18 data implicates identity by descent as an important contributor to shared genetic variance. To explore other contributions to SBP variance, we conducted sensitivity analyses using only unrelated individuals (data not shown) and found that the basic pattern held true, with rare variants explaining the highest proportion of outcome variance regardless of functional annotation.
As whole genome sequence data become widely available and computational software evolves in its ability to handle a large number of genetic variants, future studies may consider repeating our analyses with a finer classification of markers. Specifically, genic regions could be further subdivided into coding or noncoding regions, or into exonic, intronic, 5′ and 3′ untranslated region, downstream, upstream, and splicing variants. For both genic and flanking regions, another approach would distinguish between variants that encode microRNAs and/ or other regulatory molecules, or consider variance explained by differential DNA methylation. Additionally, future studies may evaluate whether the tradeoff between the proportion of outcome variance explained and the multiple testing burden is more favorable for randomly selected (eg, evenly spaced) variants across the genome or for variants located within a priori defined biological pathways. Finally, in our implementation we used shrinkage estimation procedures; however our proposed hypothesis could also be tested using statistical methods that perform variable selection and shrinkage simultaneously. In future, our findings may be extended beyond estimating the proportion of genetic variation to whole genome prediction, fully realizing the potential for clinical and public health applications of deep sequence data.

Conclusions
We have investigated the tradeoff between proportion of blood pressure variance explained using subsets of the whole genome sequence, and found that rare and very rare variants contribute more outcome variance regardless of their functional annotation. Table 4 Proportion of phenotypic variance of log(SBP) (averaged over 500 replicates, ± SD) explained by simultaneous regression on nongenetic covariates and on sets of equal size (500 markers), defined according to functional annotation and MAF