Rare genetic variant analysis on blood pressure in related samples

Abstract

The genetic variants associated with blood pressure identified so far explain only a small proportion of the total heritability of this trait. With recent advances in sequencing technology and statistical methodology, it becomes feasible to study the association between blood pressure and rare genetic variants. Using real baseline phenotype data and imputed dosage data from Genetic Analysis Workshop 18, we performed a candidate gene association analysis. We focused on 8 genes shown to be associated with either systolic or diastolic blood pressure to identify the association with both common and rare genetic variants, and then did a genome-wide rare-variant analysis on blood pressure. We performed association analysis for rare coding and splicing variants within each gene region and all rare variants in each sliding window, using either burden tests or sequence kernel association tests accounting for familial correlation. With a sample size of only 747, we failed to find any novel associated genetic loci. Consequently, we performed analyses on simulated data, with knowledge of the underlying simulating model, to evaluate the type I error rate and power for the methods used in real data analysis.

Background

Despite the tremendous success of genome-wide association studies (GWAS) to uncover genetic variants influencing complex traits and diseases, only a fraction of the total heritability of these traits is explained by the loci identified so far. Because GWAS focuses on common variants, a possible source of the missing heritability might be rare variants that were not included in the earlier genotyping platforms. The next logical step is to investigate rare variants, an endeavor that is now possible because of the ever-decreasing cost of sequencing.

Whole genome sequencing has the ability to uncover rare variants, but brings its own challenges. Despite a low error rate, the sheer number of base pairs sequenced makes it hard to distinguish very rare mutations from sequencing errors. Moreover, detecting association with rare variants requires very large sample sizes. Several methods to jointly analyze rare variants within a genomic region have been developed, however, and these methods have the potential to pinpoint additional variants contributing to the overall heritability of traits.

Blood pressure (BP) and hypertension are prime examples of the limitations of GWAS. Meta-analysis of GWAS from a large number of cohorts has identified multiple genetic loci over the genome that affect systolic blood pressure (SBP), diastolic blood pressure (DBP), hypertension, or a combination of these traits [1, 2]. However, the loci identified to date explain only a small portion of the total heritability in BP.

In this article, we investigate the association of rare variants in genomic regions that have been previously implicated by GWAS to identify the source of the original GWAS signal and to discover additional genetic loci influencing BP using either burden tests adjusting for familial correlation (famBT) or sequence kernel association tests (SKAT) [3] for family samples (famSKAT) [4, 5]. We also analyze rare variants genome-wide to uncover additional genomic regions harboring susceptibility variants. Finally, we use the simulated data sets with knowledge of the answer to evaluate type I error and power for famBT and famSKAT in family samples.

Methods

We used imputed single-nucleotide polymorphism (SNP) dosage files from odd-numbered chromosomes provided by the Genetic Analysis Workshop 18 (GAW18) as our genotypes in all analyses. For real data analysis, we took baseline measurements of the covariates and traits for each participant, defined as the first exam with nonmissing values for age, SBP, DBP, current use of hypertension medications (BPmeds), and current smoking (smoke). We removed participants with at least 1 missing value of these variables in all 4 exams. We also excluded participants on antihypertensive medication at the baseline we defined, resulting in a sample size of 747 participants. Table 1 provides the descriptive information for our subset of participants. Because the distribution of SBP values is highly skewed, rank-normalized SBP (rSBP) values were used in all analyses. DBP values were untransformed. We adjusted for sex, age, and smoking in all our analyses.

Table 1 Characteristics of variables of interest at the baseline

For the BP candidate gene study, we performed both common and rare-variant analysis. Common variants were defined as any variants with minor allele frequency (MAF) >5% in our subset of participants, and rare variants were variants with MAF between 0% and 5%. We performed common variant analysis as single-marker association tests using linear mixed-effect models [6] to account for familial correlation and reported the most significant SNP in each region. We performed rare variant analysis for all rare variants within each gene region with famBT and famSKAT [4, 5], using Wu weights, which is a beta distribution probability density function of the MAF with parameters 1 and 25 [3]. Both rare-variant approaches are described below.

Burden tests adjusting for familial correlation (famBT)

Assuming that the sample size is n, Y is a vector of the trait of interest; X is an n by p matrix of covariates; α is a vector of the covariate effects; G is an n by q matrix of rare variants with columns G _j ; w _j is the weight for variant j; and this is the combined genotype score:

$$g={\displaystyle \sum _{j=1}^q}{w}_j{G}_j$$

(1)

The model is

$$Y=X\alpha +g\beta +\gamma +ϵ$$

(2)

where β is the effect size for the combined genotype score, γ is the random effect vector for familial correlation, and ε is the normally distributed error. We assume $\gamma ~N\left(0,{\sigma }_G^2\Phi \right)$, $ϵ~N\left(0,{\sigma }_E^2{I}_n\right)$, where ${\sigma }_G^2$ and ${\sigma }_G^2$ are variance component parameters, and Φ is twice the kinship matrix. The model can be fitted as a linear mixed-effect model and the genotype effect can be tested as H₀ : β = 0 versus H₁ : β ≠ 0 in this framework.

SKAT for family (famSKAT)

We use the same notation as above, except that β is now a vector of length q. The model is

$$Y=X\alpha +GW\beta +\gamma +ϵ$$

(3)

where W is a diagonal matrix of weights w _j , and we assume β ~N (0, τI _q ). The genotype effects can be tested as H₀ : τ= 0 versus H₁ : τ> 0.

For the genome-wide rare-variant analysis on real data, we performed famBT and famSKAT, using both a gene-based coding and splicing variants analysis (GB) and sliding-window analysis (SW). GB was performed for each gene, using only nonsynonymous rare variants and rare variants at the splicing sites. SW was performed for all rare variants in each genomic region of 4000 base pairs (bp) length, with 2000 bp each overlapping with the previous and subsequent windows, regardless of the gene annotation.

Simulations

In addition to real data analysis, we also performed rare-variant analysis on simulated data sets, with knowledge of the underlying simulating model. To be consistent with the real data analyses, we adjusted for sex, age, and smoking in all analyses, even though simulated smoking is not associated with simulated SBP or DBP. Because we did not have missing data in the simulated data sets, we took the first exam as the baseline and excluded individuals taking antihypertensive medication at baseline. Therefore, the sample size varies slightly in different simulation replicates. We used both famBT and famSKAT for GB and SW, but we analyzed only chromosome 3 because of limited computing resources. We evaluated the type I errors of these approaches using quantitative trait Q1, which was a simulated trait not associated with any genetic variants. We also calculated the empirical power for MAP4 on DBP and rSBP.

Results

Candidate gene analysis

We chose 8 gene regions (CASZ1, MTHFR, ULK4, PLEKHA7, CSK, CSK-ULK3, PLCD3, ZNF652) that were previously reported to be associated with either SBP or DBP [1, 2]. Table 2 shows the common variant analysis and rare-variant analysis results. We report the SNPs with the lowest p values within each gene. However, the gene-based approach for SKAT and burden test for all associated regions did not reach our threshold for statistical significance (p values >0.05/8 = 0.00625). The tests did not identify evidence of association using real data. Candidate gene association results for common-variant analysis on 2 traits, rank normalized SBP and DBP, are also presented in Table 2. There were 3364 common SNPs among these gene regions. None of them was statistically significant after adjusting for multiple testing.

Table 2 Candidate gene analysis results

Genome-wide rare-variant analysis

Table 3 summarizes the genome-wide rare-variant analysis results on the real data. Because there are approximately 20,000 genes in the human genome, we used 2.5 × 10⁻⁶ as the genome-wide significance threshold for GB. Given that the human genome has approximately 3 billion bp, we tested approximately 1.5 million sliding windows, each with 4000 bp length and 2000 bp overlap. We thus used 3.3 × 10⁻⁸ as the genome-wide significance threshold for SW. However, for both GB and SW, none of the genes or sliding windows was found to be associated with the traits at genome-wide level.

Table 3 Genome-wide rare-variant analysis top findings

Simulations

We analyzed all 200 replicates for both GB and SW approaches. Table 4 summarizes the empirical type I errors. For both methods, famBT and famSKAT have correct type I error rates at α levels of 0.05 and 0.001. Table 5 shows empirical power of famBT and famSKAT for the 2 SNP selection approaches. For gene-based analysis, the sample size ranges from 742 to 783, but the number of rare coding variants in MAP4 is 6 in all 200 replicates. For each replicate, we performed both famBT and famSKAT on baseline untransformed DBP and rSBP and computed empirical power as the proportion of replicates with p values less than the corresponding thresholds. At an α level of 2.5 × 10⁻⁶, both methods have 100% power to detect association with MAP4. However, famBT has a median p value of 1.2 × 10⁻¹⁴ for DBP and 3.6 × 10^-16 for SBP, whereas famSKAT has a median p value of 9.7 × 10⁻¹⁴ for DBP and 1.6 × 10^-15 for SBP, suggesting that famBT would be more powerful than famSKAT if a more stringent α level were used. For sliding-window analysis, the sample size also ranges from 742 to 783, and there are 119 windows overlapping with the MAP4 gene. The number of rare variants in these windows ranges from 4 to 21. Because the MAP4 gene spans a region approximately 239 kb, 60 consecutive windows out of 119 fully cover the gene. For each replicate, we performed both famBT and famSKAT for all 119 windows, selected the smallest p value, and multiplied it by 60 for the purpose of adjusting for multiple testing. This adjustment is conservative because the 60 consecutive windows are correlated. Thus, the power of GB and SW may not be directly comparable. However, it is obvious that famSKAT is much more powerful than famBT in SW.

Table 4 Empirical type I errors from simulation data sets

Table 5 Empirical power from simulation data sets

Discussion

MAP4 encodes microtubule-associated protein 4. This gene is located in chromosome 3p21. The SNPs within the gene region have previously shown a genome-wide association with mean arterial pressure. The top-ranking SNP (rs319690) yields a p value of 2.69 × 10⁻⁸ [7]. MAP4 microtubule decoration restricts with beta-adrenergic receptor recycling, which might explain beta-adrenergic receptor downregulation in heart failure [8].

It is not surprising that in gene-based analysis, famBT is slightly more powerful than famSKAT because, among the 6 rare coding variants in MAP4--3_47894286, 3_47913455, 3_47957741, 3_47957996, 3_48040283, and 3_48040284--5 were simulated to be negatively associated with both SBP and DBP. The last SNP, 3_47894286, has perfect linkage disequilibrium (r = 1) with 3_47913455. In such a simulation setting, with SNPs all having the same direction of effect, the burden test should outperform most statistical approaches.

In sliding-window analysis, however, even though MAP4 is the gene most significantly associated with both SBP and DBP, some rare regulatory variants were simulated to be positively associated with the traits. As a result, famBT has almost no power to detect the association in this gene region. In contrast, famSKAT performs very well because SKAT allows effects to be in different directions. After adjusting for multiple testing, famSKAT still attains good power even at low α levels.

Conclusions

The SW method is more computationally intensive than GB because more tests are performed. However, by using SW we can generally test all possible rare variants associated with the trait, no matter where they are located. In many scenarios, intergenic variants, especially those within regulatory regions, also may be associated with quantitative traits. Thus, for rare-variant analysis on real data, unless we have strong a priori knowledge that the associated variants are nonsynonymous, we would recommend running a sliding-window analysis. By using famSKAT, we can perform rare-variant analysis on family data and have much better power than simple burden tests when there are variants with both positive and negative effects.