The development of high-throughput biotechnologies allows the collection of omics data to study the biological mechanisms underlying complex diseases at different levels, such as genomics, epigenomics, and transcriptomics. However, each technology is designed to collect a specific type of omics data. Thus, the association between a disease and one type of omics data is usually tested individually, but this strategy is suboptimal. To better articulate biological processes and increase the consistency of variant identification, omics data from various platforms need to be integrated. In this report, we introduce an approach that uses a modified Fisher’s method (denoted as Omnibus-Fisher) to combine separate p-values of association testing for a trait and SNPs, DNA methylation markers, and RNA sequencing, calculated by kernel machine regression into an overall gene-level p-value to account for correlation between omics data. To consider all possible disease models, we extend Omnibus-Fisher to an optimal test by using perturbations. In our simulations, a usual Fisher’s method has inflated type I error rates when directly applied to correlated omics data. In contrast, Omnibus-Fisher preserves the expected type I error rates. Moreover, Omnibus-Fisher has increased power compared to its optimal version when the true disease model involves all types of omics data. On the other hand, the optimal Omnibus-Fisher is more powerful than its regular version when only one type of data is causal. Finally, we illustrate our proposed method by analyzing whole-genome genotyping, DNA methylation data, and RNA sequencing data from a study of childhood asthma in Puerto Ricans.
Because of major advances in high-throughput biotechnologies, large amounts of omics data have been collected to study the biological mechanisms underlying complex diseases at different levels, such as genomics, epigenomics, and transcriptomics. Such different types of omics data can help us understand a disease from several perspectives. However, each of the arrays or sequencing technologies is designed to collect a specific type of omics data, such as SNPs, DNA methylation markers, and RNA sequencing. Thus, the association between a complex disease and one type of omics data is usually tested individually, but this strategy is suboptimal and has some disadvantages. Researchers often find that only a small proportion of disease variation can be explained by one type of omics (e. g. , genetic) data, leading to “missing heritability” [1]. Moreover, molecular variants identified by different studies usually suffer from poor reproducibility [2,3]. Most importantly, only partial information is used for each individual analysis. Therefore, in order to better characterize biological processes and increase the consistency of variant identification, omics data from separate platforms need to be integrated and analyzed. Integrating information from different biological datasets has the potential to yield better insight into causal mechanisms of complex diseases than that from individual omics datasets. Although integrative analysis of omics data is clearly needed, the complexity of disease mechanisms, the large number of collected molecular variables, and relatively small datasets can make such analysis quite challenging. Bersanelli et al. [4] summarized a list of existing statistical approaches for integrative analysis. When developing integrative analysis methods, two inevitable issues are: 1. handling a large number of variables and 2. dealing with data from relatively small studies. In genetic studies, to handle a large number of genetic variants in a gene, gene-based approaches [5–15] have been developed to evaluate the joint effects of genetic variants in the same gene on the disease of interest. Of the existing methods, the sequence kernel-machine-based associations test (SKAT) [16,17] is a powerful, flexible, and computationally efficient test. In this kernel machine (KM) approach, the test statistic follows a mixture of chi-square distributions, and thus p-values can be computed analytically and quickly without using resampling techniques. Although gene-based tests were originally developed for genetic studies, the same concept can be applied to studies of multi-omics data. Another issue is small sample size, especially for epigenomic or transcriptomic data. For example, large-scale genome wide association studies (GWASs) have been widely conducted for genetic studies for many years, so researchers usually have hundreds or thousands of genotyped samples. However, genome- wide methylation studies are more recent, and thus researchers often have a small number of samples. Moreover, incomplete samples may be wasted when using methods requiring complete samples (e. g. , methods incorporating multi-omics data variables into one regression model). In this scenario, methods combining multiple p-values can be applied to make full use of data. For example, the p-values for association testing of a disease and SNPs, DNA methylation markers, and RNA sequencing data are be calculated separately, and then these separate p-values can be appropriately combined into one final p-value. In order to test for overall gene-level significance, we here present an approach to use a modified Fisher’s method (denoted as Omnibus-Fisher) to combine separate p-values for association testing of a disease or trait and SNPs, methylation markers, and RNA sequencing data calculated by KM regression into an overall gene-level p-value accounting for correlation between omics data. This method can be applied to either samples with all three types of omics data or samples with one or two types. To account for all possible disease models, we further extend the modified Fisher’s method to an optimal test by using perturbations. In our simulation studies, we show that a usual Fisher’s method has inflated type I error rates when directly applied to correlated omics data. In contrast, our Omnibus-Fisher test preserves the expected type I error rates when employed in correlated omics data. Moreover, the Omnibus-Fisher method has increased power compared to its optimal version when the true disease model involves all of SNPs, methylation markers, and RNA sequencing data. On the other hand, the optimal Omnibus-Fisher method is more powerful than its regular version when only one type of data is causal. Finally, we illustrate our proposed methodology by analyzing whole-genome genotyping, DNA methylation, and RNA sequencing data from a study of childhood asthma in Puerto Ricans.
When applied to samples with independent SNPs (G), methylation markers (M), and RNA sequencing data (E), all of the methods used (i. e. , the Fisher’s methods with and without considering p-value covariance [Omnibus-Fisher and usual Fisher], and the optimal test with p-values from Omnibus-Fisher as inputs [optimal Omnibus-Fisher]) had empirical Type I error rates close to the nominal level (Fig 1A and Table 1). When the usual Fisher’s method without considering covariance was applied to G and E correlated data, the Type I error rate was inflated (Fig 1B and Table 1). In contrast, the optimal and regular Omnibus-Fisher methods with considering covariance retained the desired Type I error rates as evidenced by the patterns observed in the QQ plots shown in Fig 1B and Table 1. Similar results were observed when extending to 100,000 datasets for evaluation (S1 Table). When we compared the power of the statistics on the samples with independent G, M and E (Fig 2), the power of optimal Omnibus-Fisher was consistent higher than that of the regular Omnibus-Fisher method when G was the only causal factor, but when G, M and E were all causal factors, the optimal methods had lower power. This was expected because the optimal methods automatically searched for the appropriate disease model; in contrast, the regular Omnibus-Fisher assumed that G, M and E were all causal factors. Thus, when the simulation matched the assumption of the regular version method, they performed better than the optimal version and vice versa. However, when G and M were causal factors, no methods were consistently better than another. Furthermore, similar patterns were observed, when evaluated using the samples with G and E correlated (Fig 3). Since the causal SNPs in G were correlated with E, GM causal was equivalent to G, M and E causal. Note that the usual Fisher’s method is not included in Fig 3 because of its inflated Type I error rate with correlated data. We used the proposed optimal Omnibus-Fisher statistic and its regular version to analyze the Puerto Rican childhood asthma data from WBCs for associations between asthma status and 14,808 genes with all SNPs, DNA methylation markers, and gene expression, with adjustment for age, gender and first two principal components calculated based genotypes. In addition, batch effect and cell type composition were also adjusted for DNA methylation and RNA sequencing data. We found that ZPBP2 was the most significant gene from both optimal (P = 1. 40×10−5) and regular (P = 3. 39×10−5) Omnibus-Fisher tests, although it didn’t reach a Bonferroni corrected significance level (P = 3. 38×10−6) (Fig 4). The ZPBP2 region from chromosome 17q21 has been consistently replicated as an asthma-susceptibility locus across diverse ethnic groups [18–28] including Puerto Ricans [29] and this region regulates its gene expression in Puerto Ricans [30]. In a meta-analysis of GWAS in Puerto Ricans [29], the only region associated with asthma was the ZPBP2 locus and the current genotypic dataset was analyzed as a part of the data. This gene could be served as a positive control in asthma genetic studies. In the optimal Omnibus-Fisher test, the significance of ZPBP2 as well as GSDMB was mainly driven by their genetic effect (P = 2. 89×10−6 for ZPBP2 and 2. 36×10−6 for GSDMB). Moreover, five additional genes (KAT2A, HIST1H1C, NFRKB, C14orf178 and ZNF213-AS1) were suggestively associated with asthma (P < 0. 0001 [Table 2]) from the regular Omnibus-Fisher test. Of these genes, KAT2A had moderate effects for SNPs, DNA methylation, and RNA expression separately, which could be overlooked by a single type of data analysis. The results also indicate that the optimal Omnibus-Fisher test was more powerful than its regular version when the significance was driven by one type of data. Conversely, when statistical significance was driven by two or three types of data, the regular Omnibus-Fisher had overall better power than the optimal version. These observations were generally consistent with the simulation results. Here, the optimal Omnibus-Fisher test does not outperform its regular version that assumes all types of omics data are in the disease model. Since only three types of omics data were analyzed in this study, it was still fine to assume they all were in the disease model. However, when more omics data are analyzed, the optimal test could be more useful than simply assuming all types of data are causal. We additionally output the p-value correlation for each gene across the whole genome (S1 Fig): 1. between SNPs and DNA methylation markers, 2. between SNPs and expression genes, and 3. between DNA methylation markers and expression genes. Analysis of the WBC genome-wide data with 1,116 samples and 14,808 genes took ~108. 8 hours on a single computing node with a 3 GHz CPU and 4 GB memory. Using a computer cluster with multiple nodes, we anticipate that genome-wide data analysis should be finished within hours using our proposed methods.
In this work, we developed an Omnibus-Fisher statistic using a kernel machine (KM) regression framework, which can be employed to test overall gene-level significance by combining separate p-values of association testing for a disease and SNPs, methylation markers, and expression genes, accounting for correlation between omics data. The separate p-values are calculated by gene-based KM regression. The gene-based analysis methods can improve power by testing a set of variants jointly and by reducing the multiple testing penalty. In addition, the method using a gene as the unit can easily combine different types of omics data that are mapped to the same gene and thus easily interpret the results. Since we do not know the exact disease model in reality, the extended optimal Omnibus-Fisher test can account for all possible disease models. Moreover, our proposed tests can be applied to either samples with all three types of omics data or with one or two types. In other words, samples with incomplete data can still contribute to the test statistic. The information about whether the different types of omics data are from the sample can also be accounted. In the simulation studies, we showed that using a usual Fisher’s method on correlated omics data results in an inflated Type I error rate, while the modified Fisher’s method, Omnibus-Fisher, had the correct Type I error rate because it considered the omics data correlation in the model. The Omnibus-Fisher method achieves better power performance compared to its optimal version when the true disease model involves all of SNPs, RNA expressions and DNA methylations. On the other hand, the optimal Omnibus-Fisher method has better power than its regular version assuming all types of data are causal when only one type of data is actually causal. Our real data study also shows that the regular Omnibus-Fisher test has better power than the optimal test, when two or three types of data contribute to the combined p-value. Because we only consider three types of omics data in this study, assuming they all are causal could be still acceptable. However, when more omics data are analyzed, we believe that the optimal test would be more powerful for most genes than simply assuming all types of data are causal. Nevertheless, both the optimal and regular Omnibus-Fisher tests are able to detect genes with moderate separate effects, which could be overlooked by single type of data analyses. Although the optimal Omnibus-Fisher test uses perturbation to consider the correlation between omics data and search for the optimal disease model, the genome-wide data analysis could be completed within hours using multiple CPUs (e. g. , one CPU for each chromosome). We adapt a stepwise manner to implement perturbation (e. g. , more iterations for smaller p-values) so as to save computation times when calculating large p-values. Thus, the majority of the computation time is used by genes with small p-values. However, if a large number of genes are highly associated with the phenotype, the optimal test may be infeasible due to computational intensity. In such case, the regular test is recommended. Although the regular test also involves resampling technique to calculate covariances between different types of omics data, it only requires a small number of resampling (e. g. , the default setting is 200 times). Our method framework is general and flexible. Both continuous and binary traits for independent samples can be analyzed. Covariates can be easily incorporated into the model and different covariates can be used for different omics data. The regular and optimal version of Omnibus-Fisher algorithms were implemented in R (http: //www. r-project. org) and the R package (https: //cran. r-project. org/web/packages/OmnibusFisher/index. html) is available.
We used KM regression to calculate the gene-level p-values for association testing of a disease and SNPs, methylation markers, and expression genes. First, we test the effect of SNPs. Let there be n subjects with q genetic variants. The n × 1 vector of the continuous trait y follows a linear model: y=Xβ+Gγ+ε, when the phenotypes are binary, y follows a logistic model: logitP (y=1) =Xβ+Gγ where X is an n × p covariate matrix, β is a p × 1 vector containing parameters for the fixed effects (an intercept and p– 1 covariates), G is an n × q genotype matrix for the q genetic variants of interest where an additive genetic model is assumed (i. e. , coded as 0,1, or 2 representing the copies of minor alleles) for illustration, γ is a q × 1 vector for the random effects of the q genetic variants, and ε is an n × 1 vector for the random error. The random effect γj for variant j is assumed to be normally distributed with mean zero and variance τwj; thus, the null hypothesis H0: γ = 0 is equivalent to H0: τ = 0, which can be tested with a variance component score test [17] in the mixed model. The random variable ε is assumed to be normally distributed, and is uncorrelated with γ: γ∼N (0, τW) ε∼N (0, σE2I), where W is a predefined q × q diagonal weight matrix for each variant and may use W = I when lacking of prior information, and σE2 is the error variance. Following the same rationale as in the derivation of the SKAT score statistic [31–33], the test statistic is: Q= (y−Xβ^) ′GWG′ (y−Xβ^) /σ^E2, when phenotypes are continuous, and Q= (y−μ^) ′GWG′ (y−μ^) when phenotypes are binary, where β^ is the vector of estimated fixed effects of covariates under H0 and μ^=logit−1 (Xβ^). Under the null hypothesis, the linear model is y = Xβ + ε, and the estimates are Σ^=σ^E2I=var (y−Xβ^) I β^= (X′X) −1X′y P0=I−X (X′X) −1X′; the logistic model is logit P (y = 1) = Xβ, and the estimates are Σ^=diag (μ^∙ (1−μ^) ) β^= (X′Σ^−1X) −1X′Σ^−1y P0=Σ^−Σ^X (X′Σ^X) −1X′Σ^. The statistic Q is a quadratic form and follows a mixture of chi-square distributions under H0. Thus, Q∼∑i=1qλiχ1, i2, where λi are the eigenvalues of the matrix P012GWG′P012 [34] for both continuous and binary traits. The p-values can be calculated by numerical algorithms, such as Davies’ method [35] and Kuonen’s saddlepoint method [36], which are both included in the R package. Analogously, the gene-level effects of DNA methylation markers and expression genes can be tested by replacing Gγ with Mρ and Eη, M is an n × k matrix for the k methylated loci, ρ is a k × 1 vector for the random effects of the k methylated loci, E is an n × g matrix for the RNA expression, and η is a g × 1 vector for the random effects of the RNA expression. When using microarray platform, multiple probes could map to the same gene and each probe has an expression value, which result in more than one expression value for one gene. Here, g is the number of probes for one gene. When using RNA sequencing platform, one gene can always have one expression value (i. e. , E is an n × 1 vector and η is a scalar), although it is also possible to obtain the transcript (i. e. , isoform) level expression values. The null hypothesis is ρ = 0 for testing DNA methylation markers and η = 0 for testing expression genes. It is worth to note that all three models can have the same or different null models. In order to have one single p-value to represent the significance of a gene, we propose an approach to test if the trait is associated with any SNP, DNA methylation marker, and RNA sequencing variant. This could help researchers to screen out potentially interesting genes. Thus, after obtaining the three p-values for SNPs, DNA methylation markers, and expression genes, respectively, we used a modified Fisher’s method [37] to combine the three p-values to one. In Fisher’s method, let pi (i = 1,2, …, w) be independent p-values obtained from n hypothesis tests. Under the null hypothesis that p-values follow a Uniform (0,1) distribution, the combined test statistic is equal to T=−2∑i=1wln (pi) that follows χ2w2. However, within a gene, these p-values are correlated, thus the generalized Fisher’s method cannot be used directly. To address this issue, we consider a Satterthwaite approximation by approximating a scaled T statistic with a new chi-square distribution [38]. cT≈χv2, wherec=vE (T), v=2[E (T) ]2Var (T), E (T) =E (−2∑i=1wln (pi) ) =2wand Var (T) =var (−2∑i=1wln (pi) ) =4w+2∑i<jcov (−2ln (pi), −2ln (pj) ) where w = 3 for SNPs, DNA methylation markers, and expression genes. The covariance part takes the correlations of p-values into account and can be empirically estimated by perturbations. The perturbation details are described in the following section. If the disease risk only depends on SNPs and the model with SNPs, DNA methylation markers, and expression genes is used, then the testing power will lose. Since in reality we do not know the underlying true disease model (e. g. , only SNP effect, both SNP and RNA variant effects, or all SNP, RNA variant, and DNA methylation marker effects; totally 7 combinations), it is difficult to choose the correct model. Thus, it is desirable to develop a method accommodating all possible disease models to maximize power. This can be achieved by using the minimum p-value of all possible models (7 combinations) as a new test statistic. Then, perturbation can be used to calculate the final p-value. The perturbation-based approach was described in Wu et al. [39]. For continuous phenotypes, with large n, under H0 the (y−Xβ^) /σ^E are approximately standard normal. Then each Q= (y−Xβ^) ′GWG′ (y−Xβ^) /σ^E2 is essentially comprised of a vector of standard normal variables sandwiching a square matrix. Thus, we can perturb each Q by replacing (y−Xβ^) /σ^E with a new, common vector of normal values to generate new score statistics. Following a similar procedure as described in Urrutia et al. [40]: 1. Calculate the p-values for SNPs (G), DNA methylation (M) and RNA expression (E) separately (i. e. , pG (0), pM (0), and pE (0) ) by KM regression. 2. For l ∈ {G, M and E}, compute Λl = diag (λl, 1, ⋯, λl, ml), and Vl = [vl, 1, ⋯, vl, ml] where λl, 1 ≥ λl, 2 ≥⋯≥ λl, ml are the ml positive eigenvalues of P0l12DlWlDl′P0l12 with corresponding eigenvectors vl, 1, ⋯, vl, ml, where Dl ∈ {omics data matrices G, M and E}. For example, the aforementioned P012GWG′P012 is for G. 3. Generate r (b) =[r1 (b), ⋯, rn (b) ]′ with each rj (b) ∼N (0,1). This indicates that one subject has one rj (b). If the subject has all G, M and E, the same rj (b) will be used for G, M and E, respectively. Thus, whether G, M and E come from the same subjects or different subjects are considered. 4. For l ∈ {G, M and E}, rotate r (b) using the eigenvectors to generate rl (b) =Vl′r (b). 5. Compute Ql (b) =rl (b) ′Λlrl (b) for each l and obtain a corresponding p-value, pl (b). 6. Repeat (3) - (5) B times to obtain pG (1), pG (2), ⋯, pG (B), pM (1), pM (2), ⋯, pM (B) and pE (1), pE (2), ⋯, pE (B) for some large number B. 7. Calculate the covariance between pG, pM and pE by using pG (b), pM (b), and pE (b) for b ∈ {0,1, …, B}. 8. Calculate the joint p-values of SNPs, DNA methylation and RNA expression (i. e. , for b ∈ {0,1, …, B}, pGM (b), pGE (b), pME (b), and pGME (b) ) by Omnibus-Fisher considering p-values covariance. 9. For l* ∈ {G, M, E, GM, GE, ME, and GME}; b ∈ {0,1, …, B}, set p (b) =min1≤l*≤L*pl* (b). 10. The final p-value for significance is estimated as p=B−1∑b=1BI (p (b) ≤p (0) )
In this research, we developed a statistical approach using a modified Fisher’s method (denoted as Omnibus-Fisher) to combine separate p-values of association testing for a trait and SNPs, DNA methylation markers, and RNA sequencing, calculated by kernel machine regression into an overall gene-level p-value to account for correlation between omics data. We further extended the method to an optimal version in order to consider all possible disease models.
sequencing techniques
medicine and health sciences
pulmonology
simulation and modeling
methylation
epigenetics
dna
molecular biology techniques
rna sequencing
dna methylation
chromatin
asthma
research and analysis methods
chromosome biology
gene expression
chemistry
chromatin modification
dna modification
molecular biology
research assessment
biochemistry
biochemical simulations
cell biology
nucleic acids
genetics
biology and life sciences
chemical reactions
physical sciences
computational biology
research errors
An integrative association method for omics data based on a modified Fisher’s method with application to childhood asthma