A permutation method for detecting trend correlations in rare variant association studies

Abstract

In recent years, there has been an increasing interest in detecting disease-related rare variants in sequencing studies. Numerous studies have shown that common variants can only explain a small proportion of the phenotypic variance for complex diseases. More and more evidence suggests that some of this missing heritability can be explained by rare variants. Considering the importance of rare variants, researchers have proposed a considerable number of methods for identifying the rare variants associated with complex diseases. Extensive research has been carried out on testing the association between rare variants and dichotomous, continuous or ordinal traits. So far, however, there has been little discussion about the case in which both genotypes and phenotypes are ordinal variables. This paper introduces a method based on the γ-statistic, called OV-RV, for examining disease-related rare variants when both genotypes and phenotypes are ordinal. At present, little is known about the asymptotic distribution of the γ-statistic when conducting association analyses for rare variants. One advantage of OV-RV is that it provides a robust estimation of the distribution of the γ-statistic by employing the permutation approach proposed by Fisher. We also perform extensive simulations to investigate the numerical performance of OV-RV under various model settings. The simulation results reveal that OV-RV is valid and efficient; namely, it controls the type I error approximately at the pre-specified significance level and achieves greater power at the same significance level. We also apply OV-RV for rare variant association studies of diastolic blood pressure.

1. Introduction

For the past decade, genome-wide association studies (GWAS) have identified thousands of common variants associated with complex diseases or traits. However, recent evidence suggests that only a small proportion of the phenotypic variance can be explained by common variants (Maher, 2008; Manolio et al., 2009; Eichler et al., 2010; Gibson, 2012). Finding the sources of missing heritability has received considerable critical attention. With the advent of the next-generation of high-throughput DNA sequencing technology, an increasing number of rare variants have been detected. Recent studies have shown that rare variants have the potential to explain part of the missing heritability and may play a key role in the development of complex diseases (Bodmer & Bonilla, 2008; Nelson et al., 2012; Tennessen et al., 2012). Due to the importance of rare variants in sequencing studies, rare variant association analysis has become an increasingly important area in GWAS. To date, a large number of statistical approaches have been proposed for common variant association analysis. However, due to the low mutation rate of rare variants, traditional methods used to test single common variants usually lead to substantial bias and low power in rare variant association analysis (Li & Leal, 2008). To address the above issue, a series of burden tests have been put forward for rare variant association analysis by collapsing a group of rare variants into a specific region. Morgenthaler and Thilly (2007) collapsed the information of the rare variants in a region into a dichotomous variable and provided an approach, called cohort allelic sum test (CAST), for detecting associated rare variants. Some other burden tests for rare variant association studies include the combined multivariate and collapsing method (CMC; Bingshan & Leal, 2008), the sum test (SUM; Pan, 2009) and the weighted sum test (WSS; Madsen & Browning, 2009), among others.

It should be pointed out that all of these burden tests implicitly assume that the effects of rare variants on the phenotype are in the same direction and magnitude (after incorporating known weights), which is obviously unreasonable in GWAS. Recent studies have shown that ignoring the different directions and magnitudes of rare variant effects may lead to loss of testing efficiency (Wu et al., 2011). Hence, there remains a need for developing an efficient rare variant association test, especially when the effects of rare variants on the phenotype are in the different direction and of the same magnitude. In a seminal paper, Wu (2011) proposed a statistical method, termed the ‘sequence kernel association test’ (SKAT), for rare variant association studies. They showed that SKAT allows for different directions and magnitudes of rare variant effects and achieves greater efficiency compared with burden tests. Some extensions of SKAT can be found in the literature (SKAT-O, Lee et al., 2012; HKAT, Lin et al., 2013; W2WK, Broadaway, 2015).

All of the above methods focus on dichotomous or continuous phenotypes. However, in practice, we usually encounter situations in which the genotype or the phenotype is ordinal. For example, it is reasonable to treat the number of risk alleles or the severity of the disease as an ordinal variable. To date, a handful of methods have been proposed for the ordinal phenotype. Diao (2010) developed the variance-components methods for linkage and association analysis of ordinal traits in general pedigrees. Zhou (2016) presented a study that tested the association between rare variants and multiple traits, including ordinal traits or combinations of ordinal traits and other traits. Wang (2017) proposed a method for detecting associations between ordinal traits and rare variants based on the adaptive combination of p-values. To date, few studies have investigated the trend correlations between ordinal genotypes and ordinal phenotypes. In this paper, we put forward a method based on the γ-statistic, called OV-RV, for detecting disease-related rare variants when both genotypes and phenotypes are ordinal. Due to the extremely low mutation rates for rare variants, the asymptotic distribution of the γ-statistic is no longer the normal distribution derived by Goodman (1963). Instead of deriving the asymptotic distribution of the γ-statistic for sparse contingency tables, we employ an empirical null hypothesis by utilizing the permutation approach proposed by Fisher. We carry out extensive simulations to compare the numerical performance of OV-RV with several existing approaches in a wide range of model settings. The simulation results demonstrate that OV-RV is valid and efficient.

The remainder of this paper is organized as follows: in Section 2, we provide a brief description of the cross-contingency table in categorical data analysis. Then, we introduce a measure called γ for detecting the association between two ordinal variables and show that the asymptotic distribution of the γ-statistic is no longer applicable in rare variant association studies. To address this issue, a detailed permutation approach is provided. Extensive simulations and a real data analysis are conducted in Section 3. Section 4 contains a discussion of our results and some potential extensions of our approach.

2. Method

Suppose there are n independent subjects in a population-based study. For each subject i, we let Y _i be the phenotype and (G _i1, …, G _im) be the genotype at the m loci, where G _ij is the number of mutations in variant j for subject i. In general, G _ij ∈ {0, 1, 2}. The genetic score of the genotype for subject i is defined as

$$G_i{\rm}={\rm } \sum\limits_{{\rm i =} 1}^m {w_ig(G_{ij})},$$

where w _i is a weight and g( · ) is a link function. In practice, the selection of the weight and the use of the link function can be of various types as long as they are justified. For example, one can choose the weight utilized in Madsen and Browning (2009) to ensure that all variants in a group contribute equally. In this paper, we choose the weight w _i = 1 and the link function $g(G_{ij}) = 1_{(G_{ij} \gt 0)}$, where 1_(.) is an indicator function. At last, according to the genetic score, the genotype levels are sorted from small to large. Correspondingly, the phenotypes can be sorted from small to large in terms of the degree of the disease. For i = 1, …, n, let i and j be the numbers of different Y _i and different G _i, respectively. For ease of notation, denote by Y the phenotype and denote by G the genotype score at the m loci. Let 0, 1, …, I − 1 and 0, 1, …, J − 1 be the levels of Y and G, respectively.

2.1. The cross-contingency table

The cross-contingency table is a tool that can properly display the joint distribution of categorical variables and has been widely used in categorical data analysis. In order to express the framework of the γ-statistic explicitly, we first provide a brief description of the cross-contingency table for rare variant association studies. The cross-contingency table of the genotype level at m loci by the phenotype level is listed in Table 1, where x _ij is the number of subjects that occurs in the cell in row i and column j. Denote by π _ij the joint probability of (Y, G) in the cell of row i and column j and by {π _ij}_I×J the joint distribution of (Y, G).

Table 1. Cross-contingency table of genotype at the m loci by phenotype.

2.2. The γ-statistic

When both Y and G are ordinal, one would expect to test the monotone trend association between Y and G, where the monotone trend association refers to Y trending to increase to higher levels or trending to decrease to lower levels as the level of G increases. Define that a pair of subjects is concordant if there exists a subject that ranks higher on G and Y simultaneously. Similarly, define that a pair of subjects is discordant if there exists a subject that ranks higher on G but ranks lower on Y. Consider two independent observations randomly sampled from the joint distribution {π _ij}_I×J. For this pair of subjects, we can express the probabilities of concordance and discordance as follows:

(1)$$\Pi _c = 2\sum\limits_i {\sum\limits_j {\pi _{ij} \left(\sum\limits_{h \gt i} {\sum\limits_{k \gt j} {\pi _{hk}}} \right )}},$$

and

(2)$$\Pi _d = 2\sum\limits_i {\sum\limits_j {\pi _{ij} \left( \sum\limits_{h \gt i} {\sum\limits_{k \lt j} {\pi _{hk}}} \right )}}.$$

Then, a natural association measure to describe the monotone trend association is the difference Π_c − Π_d.

Assume that a pair is untied on both Y and G; in other words, the probability of ties Y _i = Y _j or G _i = G _j is zero. Then, Π_c/(Π_c + Π_d) and Π_d/(Π_c + Π_d) are the probabilities of concordance and discordance, respectively. Goodman (1954) suggested utilizing the difference between these probabilities to measure this trend. Specifically, the measure called γ is defined as

(3)$$\gamma = \displaystyle{{\Pi _c-\Pi _d} \over {\Pi _c + \Pi _d}}.$$

Correspondingly, the sample version is

(4)$$\mathop \gamma \limits^\wedge = \displaystyle{{C-D} \over {C + D}}, $$

where $C = \sum\limits_i {\sum\limits_j {x_{ij}}} (\sum\limits_{h \gt i} {\sum\limits_{k \gt j} {x_{hk}}} )$ and $D = \sum\limits_i {\sum\limits_j {x_{ij}}} (\sum\limits_{h \gt i} {\sum\limits_{k \lt j} {x_{hk}}} )$ are the total numbers of concordant pairs and discordant pairs, respectively.

Note that testing H ₀: Y and G are independent can be reduced to testing H ₀: γ = 0, when both Y and G are ordinal. Goodman (1963) further derived the asymptotic distribution of the γ-statistic under the null hypothesis. However, in rare variant association studies, the asymptotic distribution is no longer applicable. This is because the low mutation rate of rare variants results in most of x _ij being extremely small or even equal to zero, which in turn leads to bias of the asymptotic distribution.

2.3. The permutation approach

In this section, we provide a detailed permutation approach for estimating the distribution of the γ-statistic in what follows.

Step 1. For a = 1, …, A, execute the following steps:

1. (a) Randomly permute the original phenotype (Y ₁, Y ₂, …, Y _n);

2. (b) Generate the new cross-contingency table by matching the permuted phenotype ($\tilde{\gamma} _1^a, \tilde{\gamma} _2^a,..., \tilde{\gamma} _n^a$) and the genotype (G ₁, G ₂, …, G _n);

3. (c) Calculate the γ-statistic $\tilde{\gamma} _a$ based on the new cross-contingency table.

Step 2. Estimate the distribution of the γ-statistic under the null hypothesis:

$$F(\tilde{\gamma} \le t) = \displaystyle{1 \over A}\sum\limits_{a = 1}^A {1_{({\tilde{\gamma}} ^a \le t)}},$$

where 1_(.) is an indicator function.

3. Simulation studies

In this section, we explore the numerical performance of our method (OV-RV) and five existing methods, including CAST (Morgenthaler & Thilly, 2007), SUM (Pan, 2009), WSS (Madsen & Browning, 2009), SKAT (Wu et al., 2011) and SKAT-O (Lee et al., 2012). It is necessary to note that SKAT and SKAT-O cannot be directly used for the situation with ordinal traits. To test the associations when both the trait and the genotype are ordinal variables, one potential adjustment is to dichotomize the ordinal phenotype variables (still named SKAT and SKAT-O) and the alternative is to treat the ordered variables as continuous variables (named SKAT-C and SKAT-O-C). We compare these testing methods in terms of two aspects. First, we determine whether these methods can control the type I error at the prespecified α level. Without loss of generality, the prespecified α levels are set to be 0·05 and 0·01 in the simulations. Second, we compare the power of these methods at the same significance level. According to the scheme for generating simulated data, the simulations are divided into two cases, including a designed parameter-based simulation and a real genotype-based simulation. The simulation results are based on 1000 replications.

3.1. Simulation I

In this simulation, we set the sample size n = 500 and consider a region of loci that consists of m rare variants. Without loss of generality, m is set to be 20 and 40, respectively. We first generate ordinal genotype variables and then use continuous intermediate variables to generate ordinal phenotype variables.

For each locus j, let p _j be the minor allele frequencies (MAFs) of the corresponding rare variants. Within each region, we randomly sampled p _j from the uniform distribution U(0.001, 0.01). Under the assumptions of the Hardy–Weinberg equilibrium law, the probabilities that the genotype score G _ij has a value of 0, 1 and 2 are (1 − p _j)², 2p _j(1 − p _j) and $p_j^2$, respectively. According to the number of mutant loci in each subject, the genetic scores G _i, i = 1, …, n are classified into j ordinal categories, where $J-1 = \mathop {\max} \limits_i \sum\limits_{j = 1}^m {1_{(G_{ij} \gt 0)}}$.

To focus on the main points, we select 12 and 18 rare variants from the region of 20 and 40 rare variants as disease-causal variants, respectively. The intermediate variables T _i, i = 1, …, n are generated by the following linear model:

$$T_i =G_i \beta + \varepsilon _i, i = 1,...,n,$$

where ɛ _i, i = 1, …, n are independent and $\varepsilon _i \sim N(0,1)$, and β = d · (1,1,0,1,1,0,1,0,0,1,1,0,0,1,1,0,1,1,0,1)^T if m = 20 and β = d · (1,0,1,1,0,0,1,0,1,0,1,0,0,1,0,0,1,0,1,0,1,1,0,0,1,0,1,0,1,0,0,1,0,0,1,0,1,0,1,0)^T if m = 40. In the following simulation, the values of d are set to be 0, 0·2, 0·4, 0·6 and 0·8, respectively. It is clear that T _i and G _i are independent when d = 0. Hence, examining the control of type I errors yielded by these testing methods is under the model setting d = 0. We use the 20%, 30% and 40% sample percentiles to discretize T _i, i = 1, …, n and generate ordinal phenotype variables Y _i, which take values of 0, 1, 2 and 3. The simulation results are exhibited in Tables 2 and 3.

Table 2. Estimated type I errors of the eight methods in Simulation I.

Table 3. Estimated power results of the eight methods based on the generated genotypes.

Table 2 presents the empirical sizes of the eight methods at different prespecified significance levels. From the upper part of Table 2, we can observe that all eight methods can control type I errors at the nominal level of approximately 0·01, except for the case in which the empirical type I error of SKAT-O-C is relatively conservative when m = 40 and α = 0.01. From the lower part of Table 2, similar results are obtained when α = 0.05. Although the empirical type I error of WSS is relatively large when m = 40 and α = 0.05, it is still acceptable. These simulation results confirm the validity of the eight methods in Simulation I.

Table 3 displays the power of the eight methods at different prespecified significance levels and different parameter settings. From Table 3, we can see that the power yielded by the eight methods is decreasing when d varies from 0·8 to 0·2. Note that the larger the value of d, the stronger the trend association between the ordinal phenotype variable and the ordinal genotype variable. It is easy to interpret the aforementioned simulation results. We can also observe that the power of OV-RV uniformly dominates the other competing methods. This indicates that our OV-RV is more efficient, especially when both the phenotype and the genotype are ordinal.

3.2. Simulation II

In this section, we perform simulations to evaluate the numerical performance of OV-RV and its competing methods on more realistic simulated data. In order to simulate data with realistic linkage disequilibrium patterns, we choose the real genotypes of 697 unrelated subjects from Genetic Analysis Workshop 17 (GAW17, http://www.1000genomes.org). Specifically, we select two genes, TG and COL6A3, as candidate genes. The TG gene contains 146 single-nucleotide polymorphisms (SNPs), among which 113 out of these 146 SNPs are rare (MAF <1%), whereas the COL6A3 gene consists of 187 SNPs, and 143 out of these 187 SNPs are rare. The means of action of these genes have been revealed in several studies (Baker et al., 2005; Maierhaba et al., 2008). For example, Maierhaba (2008) pointed out that the TG gene encodes thyroglobulin and may lead to hypothyroidism and autoimmune disorders.

Likewise, for each of these two genes, we randomly selected 20 and 40 rare variants to form a region of loci, respectively. Assume that the effects of rare variants on the phenotype are in the same direction. The rest of the method for generating the ordinal phenotype values is the same as in Simulation I, and we omit the details. The detailed simulation results of the empirical sizes are listed in Tables 4 and 5. To further illustrate the superiority of OV-RV in detecting trend associations, we conduct simulations to compare the power of these methods and list the results in Tables 6 and 7.

Table 4. Estimated type I errors of the TG gene of the eight methods.

Table 5. Estimated type I errors of the COL6A3 gene of the eight methods.

Table 6. Estimated power results of the TG gene of the six methods.

Table 7. Estimated power results of the COL6A3 gene of the six methods.

Table 4 displays the empirical type I errors of the eight methods for the TG gene. The upper part of Table 4 presents the results with the significance level α = 0.01, whereas the lower part of Table 4 lists the results with α = 0.05. From Table 4, it is apparent that OV-RV controls the empirical type I errors properly at the different significance levels. We can also see that SKAT is always conservative for the TG gene. This phenomenon may be largely due to the improper dichotomization for SKAT. Table 5 presents the empirical type I errors of the eight methods for the COL6A3 gene. It can be observed that the simulation results are almost wholly consistent with those in Table 4. Although the empirical type I error yielded by OV-RV is a little aggressive when α = 0.01 and m = 20, it is still acceptable. Overall, these results further indicate that the distribution of the γ-statistic under the null hypothesis can be properly estimated by exploiting the permutation method appropriately.

Tables 6 and 7 exhibit the simulation results of the power comparisons of the six methods for the TG gene and the COL6A3 gene in Simulation II, respectively. Due to the extremely low power of SKAT and SKAT-C, we do not list their simulation results in these tables. It is clear that OV-RV shows a significant improvement in power compared with the other five methods at all model settings. By employing the γ-statistic, OV-RV can achieve greater efficiency for detecting the trend associations. Similarly, we can also conclude that the power of these methods is increasing in the parameter d. We can also determine that the power when m = 40 is uniformly larger than the corresponding power when m = 20. Under the assumption that the effects of rare variants on the phenotype are in the same direction, a larger number of causal rare variants implies a stronger trend association with the same value of d. Hence, it is easy to interpret the results.

We carry out additional simulation studies for OV-RV in testing the effects with different directions. The detailed simulation results are displayed in Additional File 1. When a small proportion of effect directions are different, the simulation results are almost wholly consistent with those in the previous simulations. However, the power of OV-RV decreases as the proportion of effects in different directions increases. This indicates that OV-RV is conservative when a large proportion of effects are of different directions. A more powerful selection of the genetic score may shed light on how to extend OV-RV to these situations, and we plan to pursue this approach in our further research.

4. Application to the detection of disease-related genes

In this section, we further apply OV-RV for the detection of disease-related genes on a real dataset called Genetic Analysis Workshop 19 (GAW19). The GAW19 dataset contains whole genome and exome sequences for odd chromosomes, gene expression measures, systolic blood pressure and diastolic blood pressure (DBP), as well as related covariates in 20 large families and 1943 unrelated individuals. Here, we focus on the 1943 unrelated individuals provided by GAW19 and consider the DBP phenotype. A series of procedures for data pre-processing are performed before carrying out association studies. We eliminate individuals who have missing phenotypes, and a total of 1851 individuals are left for analysis. In addition, we complete the missing genotype by a random sample based on the MAF.

DBP is measured in millimetres of mercury (mmHg) when the heart is at rest between beats. It has been reported that genes EBF1 and NPR3 on chromosome 5, as well as gene TMEM133 on chromosome 11, are associated with DBP (Sun et al., 2016). We apply our proposed OV-RV to test associations between these genes and DBP. From the hg19 reference (see https://www.cog-genomics.org/static/bin/plink/glist-hg19), we can obtain the gene starts and gene ends of these three genes. For each gene, genotypes are generated by selecting rare variant loci with MAF <5%. The significance level is set to be 0·05. The phenotypes are divided into four levels in terms of DBP. To be specific, phenotypes with DBP <60, 60 ⩽ DBP <80, 80 ⩽ DBP <90 and DBP ⩾ 90 correspond to levels 0, 1, 2 and 3, respectively. Due to the poor performance of the CAST, SUM and WSS methods in simulations, we only compare the performance of the remaining five methods. Detailed results are shown in Table 8. It is clear that the p-values yielded by OV-RV are uniformly smaller than those of the competing methods. We can also see that OV-RV identifies all three DBP-related genes, whereas the other competing methods identify at most one related gene. This indicates that OV-RV is more efficient at detecting disease-related genes.

Table 8. The Genetic Analysis Workshop 19 (GAW19) data shown as a list of genes associated with diastolic blood pressure.

5. Discussion

In this paper, we propose a novel method, called OV-RV, for the detection of the trend associations between ordinal genotypes and ordinal phenotypes. The γ-statistic has been successfully applied to the field of searching for the trend associations. However, the asymptotic distribution of the γ-statistic derived by Goodman (1963) is no longer valid for rare variant associations. Instead of using the asymptotic distribution directly in rare variant associations, we utilize the permutation method to estimate the distribution of the γ-statistic under the null hypothesis. Both the designed parameter-based simulation and the real genotype-based simulation illustrate that OV-RV is valid and more efficient compared with its competitors. A real data analysis on the GAW19 dataset shows that OV-RV achieves greater efficiency and can detect more disease-related genes.

Our OV-RV can also be extended in several ways. First, it has been shown that different diseases or traits usually share similar genetic mechanisms. Conducting an integrative association analysis of several traits can significantly improve testing efficiency. Hence, it is desirable to develop a method for testing associations between ordinal genotypes and multiple ordinal phenotypes. Second, as illustrated in simulations, the power of OV-RV decreases as the proportion of effects in different directions increases. This means that OV-RV is conservative when a large proportion of effects are of different directions. It would be of interest to obtain a more powerful genetic score for extending OV-RV to these situations. Third, the permutation method brings large computation costs when there is a large number of rare variants. Recently, algebraic statistics has been successfully applied in testing independence from the sparse contingency table. It may give rise to a novel method for testing trend associations from the sparse contingency table.