Methodology Article  Open  Published:
Integrated genomic analysis of biological gene sets with applications in lung cancer prognosis
BMC Bioinformaticsvolume 18, Article number: 336 (2017)
Abstract
Background
Burgeoning interest in integrative analyses has produced a rise in studies which incorporate data from multiple genomic platforms. Literature for conducting formal hypothesis testing on an integrative gene set level is considerably sparse. This paper is biologically motivated by our interest in the joint effects of epigenetic methylation loci and their associated mRNA gene expressions on lung cancer survival status.
Results
We provide an efficient screening approach across multiplatform genomic data on the level of biologically related sets of genes, and our methods are applicable to various disease models regardless whether the underlying true model is known (iTEGS) or unknown (iNOTE). Our proposed testing procedure dominated two competing methods. Using our methods, we identified a total of 28 gene sets with significant joint epigenomic and transcriptomic effects on oneyear lung cancer survival.
Conclusions
We propose efficient variance componentbased testing procedures to facilitate the joint testing of multiplatform genomic data across an entire gene set. The testing procedure for the gene set is selfcontained, and can easily be extended to include more or different genetic platforms. iTEGS and iNOTE implemented in R are freely available through the inote package at https://cran.rproject.org//.
Background
Burgeoning interest in integrative analyses has produced a rise in studies which incorporate data from multiple genomic platforms. In general, there are two methods of integrating genomic data [1]. The first is horizontal integration, where genomic data from different studies but of the same type (e.g. multiple gene expression microarray studies) are combined, sometimes across labs, cohorts, and platforms. The second is vertical integration, where multiple levels of ’omics data (e.g. DNA variation, methylation, and gene expression) are gathered on the same subjects and analyzed. A useful distinction to be made in methods for vertical integrative approaches involves whether the multiplatform data are assessed via a “screenandclean” paradigm [2, 3], where each platform is analyzed independently to screen for and select a subset of significant candidates to use in a combined analysis (i.e. a sequential integration analysis), or whether the multiplatform data are assessed simultaneously (i.e. a joint integration analysis).
Most integrative studies employ approaches that primarily rely on dimension reduction methods to accommodate the high dimensionality of analyzing multiple platforms [4, 5]. These techniques seek to synthesize complex genetic information into summary statistics, potentially at the cost of discarding large quantities of data which might still be mechanistically informative. And while methods development for nonreductive multiplatform integrative analysis has become more common in recent years [6, 7], these methods are mainly restricted to candidate gene interrogations, and do not encapsulate the highly likely networklevel interactions between diseaseriskconferring genes. Of course, numerous tests of gene sets are available [8–10] – but few that also include the integration of additional genomic platforms.
Additionally, literature for conducting formal hypothesis testing on an integrative gene set level is considerably more sparse than that for estimation. For example, integrative methods for identifying potential risk pathways include strategies that employ Bayesian mixture modeling [11–14], Bayesian graphical models [13], Bayesian network models [15], nonnegative matrix factorization [16–18], and weighted gene correlation network approaches [3]. To our knowledge, methods for joint integrative testing of any kind are small in number; for gene sets, there is a variant of GSEA [4, 5], and for candidate gene approaches there are a few multivariate and mediation methods [6, 7, 19]. Although effect estimation is informative when candidate gene sets/networks are already identified or hypotheses are welldefined, an efficient screening approach across multiplatform genomic data is critical for hypothesis generation. Therefore, in this paper, we focus on efficient testing procedures to assess the effect of an entire gene set through the joint analysis of multiple genomic platforms, such as epigenomic and transcriptomic data.
Joint integrative analyses become substantially challenging when considered on the level of gene sets, where the number of model parameters rapidly increases as the size of the gene set grows. Additionally, correlation structure within a gene on the level of methylation sites, as well as between genes on the transcript expression level, may cause conventional univariate or multivariate tests to perform poorly [10, 20, 21]. We therefore propose a variance component test to assess the total effect of a set of methylation loci and mRNA gene expressions across a gene set on disease outcome. The test statistic for the joint gene set analysis follows a mixture of χ ^{2} distributions, which we may approximate analytically, or empirically using a perturbation procedure, after specifying a disease model for the whole gene set (e.g. epigenetic effect only, or epigenetic effect and gene expression effect, or both epigenetic and gene expression effect as well as their interactions). However, because the true disease models underlying different genes may vary, we also construct two gene set level omnibus tests to accommodate different disease models. A general overview of our approach is presented in Fig. 1.
The biological motivation for this paper lies in the connection between DNA methylation (DNAm) patterns and lung cancer survival. In particular, we are interested in the total joint effect of DNAm and downstream mRNA expression levels for all genes in a related pathway on survival probability in 559 subjects with both epigenomewide DNAm and RNAsequencing data from The Cancer Genome Atlas (TCGA). We demonstrate the utility of our integrative testing procedures by identifying significant gene sets that can be further explored for potential biomarkers of prognosis or even therapeutic targets.
Methods
Our integrative gene set testing approach can be viewed as a variance component test [6, 10] under the generalized linear mixed model framework [22].
Integrated gene model and test of total effects
Huang et al. [6] proposed a method to jointly analyze the effects of a set of genetic markers and a corresponding measure of gene expression within a single candidate gene on disease outcome, which is applicable to the analysis of epigenetic and transcriptomic data. Briefly, let Y _{ i } represent the dichotomous disease outcome of subject i (i=1,…,n) and let X _{ i } represent r covariates of interest for subject i. Further assume that Y _{ i } is associated with the r covariates of interest X _{ i } (with the first covariate set as the intercept), the methylation levels at a set of p CpG loci within the candidate gene (M _{ i }=(M _{1i },…,M _{ pi })^{′}), the corresponding gene expression (G _{ i }), and possibly their interactions. Then, the underlying model for any given candidategene total effect test is:
where $\boldsymbol {\beta }_{X}=\left (\beta _{X_{1}},{\ldots },\beta _{X_{r}}\right)^{\prime }$, $\boldsymbol {\beta }_{\boldsymbol {M}}=\left (\beta _{M_{1}},{\ldots },\beta _{M_{p}}\right)^{\prime }$, β _{ G }, $\boldsymbol {\beta }_{C}=\left (\beta _{C_{1}},{\ldots },\beta _{C_{p}}\right)^{\prime }$ represent the regression coefficients for the covariates, the CpG loci, gene expression, and the interactions between the CpG set and gene expression, respectively. Then, the null hypothesis for a singlegene test of total effect is:
which can be cast into a variance component testing framework by assuming: 1) the elements of β _{ M } are independent and follow an arbitrary distribution with mean 0 and variance τ _{ M } and 2) the elements of β _{ C } are independent and follow an arbitrary distribution with mean 0 and variance τ _{ C }. In other words, the outcome model (1) becomes a logistic mixed model and the null hypothesis may be reexpressed as:
Using the above model specifications, the score statistics may be derived for τ _{ M },β _{ G } and τ _{ C } respectively as:
where $\mathbb {M}=\left (\boldsymbol {M}_{1},{\ldots },\boldsymbol {M}_{n}\right)^{\prime }$, G=(G _{1},…,G _{ n })^{′}, $\mathbb {C}=\left (\boldsymbol {C}_{1},{\ldots },\boldsymbol {C}_{n}\right)^{\prime }$, C _{ i }=G _{ i } M _{ i }, μ̂_{0}=(μ̂_{01},…,μ̂_{0n })^{′}, and $\Hat {\mu }_{0i}=e^{X_{i}^{\prime }\Hat {\boldsymbol {\beta }_{X}}}/\left (1+e^{X_{i}^{\prime }\Hat {\boldsymbol {\beta }_{X}}}\right)$ is the mean Y _{ i } under the null model
where β _{ X } ̂ is the maximum likelihood estimator of β _{ X }. Using a conventional approach to combine the score statistics for each component such that $Q_{conv}=\boldsymbol {U}^{\prime }\mathcal {I}^{1}\boldsymbol {U}$, where $\boldsymbol {U}=(U_{\tau _{M}}, U_{\beta _{G}}, U_{\tau _{C}})$), would involve combining score statistics from different scales and requires the existence of the 8th moment of Y to calculate the efficient information matrix of U, $\mathcal {I}$. Therefore, the component score statistics are instead summed to create a weighted test statistic for the null hypothesis (3), denoted as Q _{∗} statistics:
where Q _{∗}={Q _{ MGC },Q _{ MG },Q _{ M },Q _{ G }} represents the underlying disease models MGC, MG, M, and G which correspond to the model specifications that include 1) CpG, gene expression, and their interactions across the full gene set, 2) the CpG and gene expression effects across the full gene set, 3) only CpG effect, and 4) only gene expression effect respectively, and the weights a _{1}, a _{2}, and a _{3} defined as the inverse square root of the variances for their corresponding score statistics to make $U_{\tau _{M}}$, $U_{\beta _{G}}^{2}$ and $U_{\tau _{C}}$ comparable.
Because $U_{\tau _{M}}$, $U_{\beta _{G}}^{2}$, and $U_{\tau _{C}}$ are all quadratic functions of Y, the null distribution of Q _{∗} may be approximated with a mixture of χ ^{2} distributions, thus we may derive pvalues for Q _{∗} by using the Satterthwaite scaled χ ^{2} approximation [23] or the characteristic function inversion method [24]. Alternatively, one can perform the test by conducting a resamplingbased perturbation procedure [25–27]. The perturbation procedure is used to approximate the null distribution of Q=Q(β _{ X } ̂) by resampling realizations of its asymptotic distribution under H _{0}. Specifically, it can be shown that
where ε is a multivariate normal random variable with mean 0 and covariance $\boldsymbol {D}=\left (\begin {array}{cc} \boldsymbol {D}_{XX} & \boldsymbol {D}_{XV} \\ \boldsymbol {D}_{VX} & \boldsymbol {D}_{VV} \end {array}\right)=n^{1}\boldsymbol {U}^{\prime }\boldsymbol {W}\boldsymbol {U}$, U=(U _{1},…,U _{ n })^{′}, $\boldsymbol {U}_{i}=(\boldsymbol {X}_{i}^{\prime },\boldsymbol {V}_{i}^{\prime })$, $\boldsymbol {V}_{i}=(\sqrt {a_{1}}\boldsymbol {M}_{i}^{\prime }, \sqrt {a_{2}}G_{i},\sqrt {a_{3}}\boldsymbol {C}_{i}^{\prime })^{\prime }$, W=diag{μ _{0i }(1−μ _{0i })}, and A _{ l } is the lth row of $\boldsymbol {A}=\left [\boldsymbol {D}_{XV}^{\prime }\boldsymbol {D}_{XX}^{1},\boldsymbol {\mathrm {I}}_{2p+1}\right ]$ where I is the (2p+1) dimensional identity matrix. In other words, Q _{∗} can be shown to follow a mixture of χ ^{2} distributions. The perturbation procedure then approximates the asymptotic distribution of Q _{∗} by generating realizations of ε, ε̂, repeatedly, where $\hat {\boldsymbol {\epsilon }}=n^{1/2}\sum _{i=1}^{n}\boldsymbol {U}_{i}^{\prime }(Y_{i}\hat {\mu }_{0i})\mathcal {N}_{i}$ and $\mathcal {N}_{i}$ are independent N(0,1). For perturbation b, we generate $\mathcal {N}^{(b)}=\left (\mathcal {N}_{1}^{(b)},{\ldots },\mathcal {N}_{n}^{(b)}\right)$, b=1,…,B (the number of perturbations) to obtain the realization of the distribution of ε, from which we approximate the distribution of Q _{∗}.
Integrated gene set model and test of total effects
We expand our model to extend the singlegene joint test proposed by Huang et al. [6] to a full gene set. Let J×1 vector G _{ i } represent the expression level for j=1,…,J genes for subject i, and $\boldsymbol {M}_{i}=\left (\boldsymbol {M}_{1i}^{\prime },{\ldots },\boldsymbol {M}_{ji}^{\prime },{\ldots },\boldsymbol {M}_{Ji}^{\prime }\right)^{\prime }$, represent the K×1 methylation value vector for the p _{ j } CpG loci of gene j with $\boldsymbol {M}_{ji}=\left (M_{1i},{\ldots },M_{p_{j}i}\right)^{\prime }$, $K=\sum _{j}p_{j}$. Similarly, to allow for interaction effects, let $\boldsymbol {C}_{i}=\left (\boldsymbol {C}_{1i}^{\prime },{\ldots },\boldsymbol {C}_{ji}^{\prime },{\ldots },\boldsymbol {C}_{Ji}^{\prime }\right)^{\prime }$, where $\boldsymbol {C}_{ji}=\left (G_{ji}M_{1i},{\ldots },G_{ji}M_{p_{j}i}\right)^{\prime }$. The model thus underlying a gene set test which includes interactions between the methylation sites and gene expression can be specified as:
where $\boldsymbol {\beta }_{M}=\left (\boldsymbol {\beta }_{M_{1}}^{\prime },{\ldots },\boldsymbol {\beta }_{M_{J}}^{\prime }\right)^{\prime }_{K\times 1}$, $\boldsymbol {\beta }_{G}=\left (\beta _{G_{1}},\beta _{G_{2}},{\ldots },\beta _{G_{J}}\right)^{\prime }_{J\times 1}$, and $\boldsymbol {\beta }_{C}=\left (\boldsymbol {\beta }_{C_{1}}^{\prime },{\ldots },\boldsymbol {\beta }_{C_{J}}^{\prime }\right)^{\prime }_{K\times 1}$ represent the coefficients for all CpG loci, gene expression, and withingene crossproduct interactions across the gene set, and $\boldsymbol {\beta }_{M_{j}}=\left (\beta _{M_{j1}},{\ldots },\beta _{M_{jp_{j}}}\right)^{\prime }_{p_{j}\times 1}$ and $\boldsymbol {\beta }_{C_{j}}=\left (\beta _{C_{j1}},{\ldots },\beta _{C_{jp_{j}}}\right)^{\prime }_{p_{j}\times 1}$. The resulting hypothesis test for the total effect of a gene set is:
As the gene set grows, however, the number of parameters to test becomes intractable under standard likelihoodbased multivariate testing methods. Similar to the above single gene analyses, we resort to an empirical Bayes approach by assuming that the effect parameters β’s share common distributions for each gene j: 1) the elements of $\boldsymbol {\beta }_{M_{j}}$ are independent and follow an arbitrary distribution with mean 0 and variance $\tau _{M_{j}}$ and 2) the elements of $\boldsymbol {\beta }_{C_{j}}$ are independent and follow another arbitrary distribution with mean 0 and variance $\tau _{C_{j}}$. Based on the above assumptions, we construct a test for the following null hypothesis:
We use a modified variance component testing procedure to obtain our test statistic, Q _{Net∗}. For the gene set being tested:
where K _{ j∗} indicates the kernel of the underlying disease model specification for gene j: $K_{j*}=a_{1j}\mathbb {M}_{j}\mathbb {M}_{j}^{\prime }+a_{2j}\boldsymbol {G}_{j}\boldsymbol {G}_{j}^{\prime }+a_{3j}\mathbb {C}_{j}\mathbb {C}_{j}^{\prime }$ for the MGC model, and $K_{j*}=a_{1j}\mathbb {M}_{j}\mathbb {M}_{j}^{\prime }+a_{2j}\boldsymbol {G}_{j}\boldsymbol {G}_{j}^{\prime }$, $K_{j*}=a_{1j}\mathbb {M}_{j}\mathbb {M}_{j}^{\prime }$, and $K_{j*}=a_{2j}\boldsymbol {G}_{j}\boldsymbol {G}_{j}^{\prime }$ for the MG, M, and G only models, respectively; we again chose the weights w _{1},…,w _{ J } to be the inverse of the standard deviation to make each Q _{ j } comparable. In closed form calculations, we assume all genes follow the same model specification: M, G, MG, or MGC such that we obtain as test statistics Q _{NetM }, Q _{NetG }, Q _{NetMG }, or Q _{NetMGC }. We note that the diseasemodel specifying only gene expression effects is in fact equivalent to the singleplatform (i.e. nonintegrative) gene set testing method proposed by Huang and Lin [10] with working independence among the genes. Their approach, called the total effect of a gene set (TEGS), is therefore a special case of the integrative methods presented here.
Under the null, Q _{Net∗} can be shown to follow a mixture of χ ^{2} distributions. Thus, as in the singlegene total effect test, we may calculate pvalues for Q _{Net∗} either by using the characteristic function inversion method (Davies method), the resamplingbased perturbation procedure, or approximate by matching the first two moments of the scaled χ ^{2} distribution (Satterthwaite method). We will refer to this method as the integrated total effect of a gene set (iTEGS) with iTEGSM, iTEGSG, iTEGSMG and iTEGSMGC denoting tests under the M, G, MG, and MGC models, respectively.
Integrated pathwaywide omnibus tests
Omnibus chisquared gene set test
A gene set drawn from a network or pathway is comprised of many genes, and each of these genes may have different underlying disease models wherein causal relationships with disease risk might be best represented by differing models M, G, MG, and MGC. The algorithm to obtain the empirical null distribution of the sum of χ ^{2} statistics of the gene set is as follows:

1.
For each gene j in the gene set:

a.
Calculate the observed Q ̂_{ jM }, then obtain its empirical distribution $\left \{\Hat {Q}_{jM}^{(b)}, b=1,{\ldots },B\right \}$ where B denotes the number of perturbations.

b.
Repeat a.) for Q ̂_{ jG }, Q ̂_{ jMG }, and Q ̂_{ jMGC } respectively.

c.
Obtain pvalues $\text {Pr}\left (\Hat {Q}_{j*}^{(b)}>Q_{j*}\right)$ for Q ̂_{ jM }, Q ̂_{ jG }, Q ̂_{ jMG }, Q ̂_{ jMGC }. Denote these as P ̂_{ jM },P ̂_{ jG },P ̂_{ jMG }, and P ̂_{ jMGC }, respectively, and $\Hat {P}_{j_{\min }}=\min \left (\Hat {P}_{jM}, \Hat {P}_{jG}, \Hat {P}_{jMG}, \Hat {P}_{jMGC}\right)$. Transform $\Hat {P}_{j_{\min }}$ to its corresponding $\chi ^{2}_{1}$ quantile denoted $\Hat {T}_{j_{\min }}$ (the $\chi ^{2}_{1}$ statistic with tail probability $\Hat {P}_{j_{\min }}$).

d.
Obtain the empirical distribution of $\Hat {T}_{j_{\min }}$, $\left \{\Hat {T}_{j_{\min }}^{(b)}\right \}$ where $\Hat {T}_{j_{\min }}^{(b)}$ is the χ ^{2} statistic with tail probability of $\Hat {P}_{j_{\min }}^{(b)}=\min \left (\Hat {P}_{jM}^{(b)},\Hat {P}_{jG}^{(b)},\Hat {P}_{jMG}^{(b)},\Hat {P}_{jMGC}^{(b)}\right)$

a.

2.
Sum the J observed $\Hat {T}_{j_{\min }}$ across the gene set such that $\Hat {T}_{\text {Net}}=\sum _{j=1}^{J}\Hat {T}_{j_{\min }}$. To obtain the empirical null for T ̂_{Net}, calculate $\left \{\Hat {T}_{\text {Net}}^{(b)}=\sum _{j=1}^{J}\Hat {T}_{j_{\min }}^{(b)}\right \}$. Calculate the geneset pvalue by obtaining the proportion of values that are more extreme than the observed T ̂_{Net}.
This approach, which we term the chitransformed integrated network omnibus total effect test (iNOTEchi), should provide a powerful approach for testing gene sets in cases where the true underlying disease models for the genes in a gene set are unknown.
Omnibus uniform network model gene set test
While iNOTEchi provides the flexibility that different genes may follow different disease models (M, G, MG or MGC), its performance may depend on whether the true underlying models for each gene are correctly selected, which introduces another source of uncertainty in model specification. In cases where the disease risk signal is not easily differentiable between the disease risk models, omnibus selection of disease models for each gene may not necessarily improve the power of the method. Therefore, we developed another test that determines a consensus disease model that is most generally applicable across the whole gene set. The complete algorithm is as follows:

1.
For each gene j in the gene set:

a.
Calculate the observed Q ̂_{ jM }, then obtain its empirical distribution $\left \{\Hat {Q}_{jM}^{(b)}, b=1,{\ldots },B\right \}$ where B denotes the number of perturbations.

b.
Repeat a.) for Q ̂_{ jG }, Q ̂_{ jMG }, and Q ̂_{ jMGC } respectively.

a.

2.
Sum the J observed Q ̂_{ j∗} across the gene set under each disease model such that we have three test statistics: Q ̂_{NetM }, Q ̂_{NetG }, Q ̂_{NetMG }, Q ̂_{NetMG C }. Calculate their associated pvalues $\text {Pr}\left (\Hat {Q}_{\text {Net}*}^{(b)}>\Hat {Q}_{\text {Net}*}\right)$, denoted P ̂_{Net∗}, then select as our omnibus network test statistic:
$$\Hat{P}_{\text{Net}_{\min}}=\min\left(\Hat{P}_{\text{Net}M},\Hat{P}_{\text{Net}G}, \Hat{P}_{\text{Net}MG},\Hat{P}_{\text{Net}MGC}\right) $$ 
3.
Obtain the empirical null for $\Hat {P}_{\text {Net}_{\min }}$ by calculating
$\left \{\Hat {P}_{\text {Net}_{\min }}^{(b)} =\min \left (\Hat {P}_{\text {Net}M}^{(b)}, \Hat {P}_{\text {Net}G}^{(b)}, \Hat {P}_{\text {Net}MG}^{(b)},\Hat {P}_{\text {Net}MGC}^{(b)}\right)\right \}$. Calculate the gene set pvalue as above by comparing the observed $\Hat {P}_{\text {Net}_{\min }}$ to $\left \{\Hat {P}_{\text {Net}_{\min }}^{(b)}\right \}$ and obtaining the proportion of values that are more extreme than the observed $\Hat {P}_{\text {Net}_{\min }}$, or by using the Satterthwaite method.
We term this approach the uniform model integrated network omnibus total effect test (iNOTEuni).
Simulation studies
We simulated DNAm based on Infinium HumanMethylation 450K Beadchip data obtained from the lung tissue samples of 681 lung cancer patients in The Cancer Genome Atlas. To realistically simulate disease outcome and gene expression, high correlation CpG blocks were identified across the epigenome to generate CpG sets which were then used to model gene expression. One causal CpG was selected per CpG set and gene expression was simulated for each subject i by the linear regression model: $\phantom {\dot {i}\!}\boldsymbol {G}_{i}=\delta _{0}+\boldsymbol {M}_{j_{\text {causal}}}\delta +\boldsymbol {\epsilon }_{i}$, where $\boldsymbol {\epsilon }_{i}\thicksim \mathcal {MVN}(\boldsymbol {0}, \boldsymbol {\Sigma })$ and Σ is a J×J covariance matrix with diag(1) and betweengene covariance equal to 0.7. Withingene covariance was accounted for by the covariance structure in actual subject data (from which the CpG blocks were drawn). For each simulation, a casecontrol sample of 100 cases and 100 controls were randomly selected from a simulated cohort of 681 subjects.
To evaluate the performance of the proposed omnibus methods, iNOTEchi and iNOTEuni, we conducted power simulations for gene set sizes of 10 and 50 at signal density proportions (i.e. the proportion of genes randomly selected to be causal within the gene sets) of 0.2, 0.5, 0.8, 1.0 across seven different simulation settings. The seven scenarios varied the mixture of underlying disease models for the causal genes in a given gene set as follows: 1) all genes follow Monly models; 2) all genes follow MG models; 3) all genes follow MGC models; 4) 50:50 mixture of Monly and MG models; 5) 50:50 mixture of Monly and MGC models; 6) 50:50 mixture of MG and MGC models; 7) onethird mixture of M, MG, MGC models.
We next compared our proposed methods, iTEGS, iNOTEchi, and iNOTEuni with two existing methods: 1) gene set association analysis (GSAA) [5], an integrative variant of the common gene set enrichment analysis (GSEA) approach to gene set testing, and 2) a more recent estimating equationbased integrative method proposed by Zhao et al. [7] which assumes that any effects of the exposure (e.g., methylation) are fully mediated by a mediator (e.g., gene expression) to produce the outcome which we will simply refer to as the ‘Zhao’ method. The Zhao method requires estimation of parameters and thus struggles to converge if the size of the gene set gets too large (e.g., the number of genes is greater than 5). To accommodate the competing method, we reduced the size of the gene set to three genes, each with 11 corresponding CpG loci, but note that the number of parameters is still quite large (i.e., 36 main effect parameters) relative to our sample size. To compare the power performance of GSAA which tests for a competitive null hypothesis [28], 49 background gene sets of equal size (3 genes per set) and null effect on disease risk were simulated in the same manner as the causal gene set in each simulation.
Application: pathwaywide association scans in TCGA
To illustrate the utility of our method, we obtained an initial sample of preprocessed level 3 genomic data from 681 lung adenocarcinoma (LUAD) and lung squamous cell carcinoma (LUSC) patients in The Cancer Genome Atlas (TCGA) database (http://cancergenome.nih.gov/) with DNAm data assayed on the Illumina Infinium Human Methylation 450K. Among the 681 subjects, 559 also had measured mRNA expression and clinical outcome data. From the 559 patients with both levels of genomic data, we identified a final analytic sample of 249 subjects who had complete information on oneyear survival since cancer diagnosis. Methylation and RNASeq data were adjusted for batch effects using the ComBat method in the Surrogate Variable Analysis (sva) Bioconductor package [29].
To obtain candidate pathways to test, we next scanned the Molecular Signatures Database (MsigDB; version 5.1) [4] for all gene sets that were associated with the keywords “lung” and “(cancer OR carcinomas)” in homo sapiens, and identified 103 gene sets of varying sizes (ranging from as small as 5 to as large as 456 genes in the gene set) for joint testing with integration of epigenomic and transcriptomic data. Among these, four gene sets were excluded due to the absence of methylation probes, mRNA expression data, or both, in all the genes that comprised each gene set, resulting in a final 99 gene sets for our joint analyses. The 99 gene sets were then scanned using iTEGS under the M, MG, and MGC diseaserisk models, as well as with the two iNOTE methods. The iTEGSG test, assuming mRNA gene expression effects only, was calculated to provide a benchmark for assessing the benefits of integrating methylation data, and incorporated in the iNOTE omnibus model selection algorithm. Finally, all gene set tests were adjusted for potential confounding covariates: smoking history (pack years), sex, age at diagnosis, race (white, black, other), pathologic tumor stage at time of initial biopsy, and cell type (adenocarcinoma, squamous cell carcinoma).
Results
Simulation study
Size and power
With the gene set size of 50, type I errors were protected for the variance component test statistics of iTEGS under each of the three gene set models assuming all causal genes within the set follow M, MG, or MGC models (Table 1). The iNOTEuni method was also well protected with a type I error rate close to 0.05. The type I error rate of iNOTEchi was 0.052 under the gene set size of 10 but slightly inflated when the gene set became larger: 0.067 for the gene set size of 25 and 0.08 for the gene set size of 50.
To evaluate the performance of the iNOTE methods with respect to power, we conducted power simulations for a set of 50 genes with signal density of 20% (i.e. 10 genes with one causal CpG locus). Power curves for simulation settings where all causal genes follow 1) M, 2) MG, 3) MGC, and 4) an approximately equal mixture of M, MG, and MGC diseaserisk models are presented in Fig. 2. Other mixtures of disease risk models were also assessed but results were similar to those of the fourth simulation setting (Additional file 1: Figure A.1). Increasing the causal signal density proportion from 20% to 80% resulted in sharp increases in power across all simulation settings, as expected (Additional file 1: Figure A.2).
In the first simulation setting where all 10 causal genes in the gene set follow the M diseaserisk model, iTEGSM demonstrates the greatest power, as expected (Fig. 2 a). The other two model formulations, iTEGSMG and iTEGSMGC, overspecify gene expression and interaction parameters for testing and thus suffer a performance loss in power. Similarly, in the simulation setting under the MG model, iTEGSMG, which correctly specifies the model, has the most optimal power performance, with iTEGSMGC achieving very similar power performance (Fig. 2 b). However, iTEGSM performs considerably worse under settings where both methylation and gene expression effects are present. In the third simulation setting where the methylationbyexpression interaction terms are present (i.e., the MGC model) and the true disease risk model is MGC, iTEGSMGC and iTEGSMG again have similar power performance, but iTEGSM demonstrates a steep drop in power as it tests only for the presence of a portion of the signal (Fig. 2 c). The final simulation setting in which the causal genes are randomly assigned to M, MG, or MGC diseaserisk models in equal proportion, the performance between the different iTEGS statistics is similar to the second simulation setting (Fig. 2 d).
Notably, across all simulation settings, the iNOTEchi and iNOTEuni tests reveal strong power performance that is nearly equivalent to the iTEGS under the correctly specified model, with the exception of the first simulation setting, where they are slightly less powerful. In the first simulation setting, iNOTEuni outperforms iNOTEchi; but in all other simulation settings however, iNOTEchi exhibits a slight power advantage compared to iNOTEuni, particularly in the case of mixtures of different causaldiseaserisk models across different causal genes within a given gene set.
Comparison to existing approaches
We also studied the performance of iTEGS and the two iNOTE tests in comparison to two competing approaches to integrative analysis, GSAA and the Zhao method using the same four simulation settings described in the internal power comparisons (to review power performance for additional mixtures of diseaserisk models, see Additional file 1: Figure B.1) In the 3gene setting, our methods behave as in the 50gene simulations where the correctly specified iTEGS demonstrates optimal power performance. Importantly, both omnibus approaches, iNOTEuni and iNOTEchi, and the correctly specified iTEGS tests consistently outperform GSAA and the Zhao method under various simulation settings (Fig. 3). Our variance componentbased tests especially dominate the Zhao method in the presence of high direct CpG methylation effects and strong correlation between methylation loci and gene expression (Fig. 3 a), which suffers from major power loss due to the presence of only direct methylation effects, rather than mediated effects through gene expression. The power of the Zhao method is somewhat recovered in simulation settings where the gene expression signal exists. The GSAA method, which tests for a competitive null hypothesis, achieved very low power across all of the simulation settings.
Application: lung cancer survival associated gene sets
We next analyzed the TCGA lung cancer data using iTEGS (under each of the Monly, MG and MGC models), iNOTEchi, and iNOTEuni. Among the 99 lung cancer associated MsigDB gene sets that were tested, iTEGS identified 57, 59, and 52 significant gene sets (p<0.05) under the MGC, MG, and M model specifications, and iNOTEchi and iNOTEuni identified 51 and 58 significant gene sets respectively (Table 2). The counts of identified gene sets using our proposed methods all exceeded what we expected under the null, i.e., 5. Gene sets that were identified as significantly associated with oneyear survival after Bonferroni correction at p<5×10^{−4} in at least one of each of the iTEGS and iNOTE tests are reported in Table 3. The pvalues obtained with the Davies method for the iTEGS statistics were generally quite similar to the perturbationbased empirical pvalues when the gene set sizes were small, but tended to vary when the gene sets grew in size (Additional file 1: Table C.1).
A total of 28 gene sets were identified as significant by at least one of the iTEGS tests and by at least one of the omnibus iNOTE tests. There were 23 and 28 gene sets with significant iNOTEchi and iNOTEuni pvalues after Bonferroni correction, respectively. Interestingly, the iTEGSMGC, iTEGSMG, iNOTEchi and iNOTEuni outperformed the iTEGSG in their ability to identify gene sets significantly associated with oneyear survival which were known a priori to be related to lung cancer, despite the fact that many of the gene sets curated by the MsigDB were obtained from gene expression studies. This is supportive of the notion that screening of gene sets using efficiently integrated multiplatform ‘omic data can increase the ability to identify potentially mechanistic disease pathways. Similar patterns supporting the utility of integrative analysis also emerged in additional exploratory gene set screening analyses with different outcomes (e.g. pathological stage of tumor at initial biopsy) and in different pathway databases (e.g. BIOCARTA and KEGG pathways, which include gene sets not specific to lung cancer) can be viewed in Additional file 1: Tables D.1D.3, E.1, and E.2.
The GSAA method only identified 8 significant gene sets, of which only one survived a Bonferroni adjustment. This is a predictable feature of the adapted KolmogorovSmirnov algorithm employed by the GSAA approach, which ignores betweengene correlation among the genes in a gene set and instead uses relative gene rankings among all possible genes under consideration. Thus, the GSAA approach is dependent on not only the size of the gene set being tested, but also the proportion of significantly associated genes belonging to a gene set of interest versus the proportion that does not. Indeed, GSAA may not reliably retrieve diseaseassociated gene sets when the proportion of signal genes in the gene set is small, even if the associations are strong and highly significant.
Among the top gene sets identified by iTEGS and iNOTE in Table 3, we recovered several involving KRAS expression and EGFR signaling, both of which are canonical genes implicated in cancer literature, as well as others related to a microRNA associated with cancer, mirlet7a3. We also retrieved several gene sets previously identified as predictive of lung cancer survival, lending further credibility to both the integrative approach and our findings. For illustrative purposes, we created methylation and mRNA expression heatmaps for one small but interesting gene set which was identified as associated with oneyear survival in our analyses: the Gautschi SRC signaling gene set (pvalues: iTEGSMGC=0.017, iTEGSMG=0.030, iTEGSM=0.653; iTEGSG=0.007; iNOTEchi=0.005, iNOTEuni=0.015; GSAA=0.205) [30], which is comprised of a set of highly downregulated genes in lung cancer cell lines after the application of an SRC inhibitor. Refined characterization of the individual genes viable for testing in the gene set showed that nonsurvivors had generally higher mRNA expression values than survivors (Fig. 4); these findings are biologically consistent with those of Gautschi et al. [30] that SRC inhibition, and therefore reduced expression of genes in the Id family, is associated with decreased cancer cell invasion.
Discussion
Our proposed approach has two advantages: first, it is a variance componentbased score test where the testing procedure is constructed under the null without estimating the large number of effect parameters; second, the omnibus tests approach the optimal performance demonstrated under correct model specification by synthesizing the evidence from three candidate models and are thus robust to model misspecification. In our simulation studies, we found that iTEGS and iNOTE dominated two competing methods, GSAA and the Zhao method. All three tests use information across multiple genomic platforms. However, the GSAA first discards information by using weighted pvalues across individual genes to integrate different genomic data, and then performs an adapted KolmogorovSmirnov test which assesses a competitive null hypothesis [28]. The Zhao method requires strong assumptions that all methylation effects on disease risk are mediated through gene expression, and struggles to converge when the ratio of parameters to the sample size is too large or when there is strong correlation between CpG loci. Although our simulations assumed causal associations between DNAm and gene expression, our testing procedures remain legitimate tests of joint effect even in cases where such associations do not exist. It should also be noted that the original development of GSAA and the Zhao method had slightly different purposes than the multiplatform integrated analyses of a gene set. For example, GSAA focuses on examining signal enrichment within a gene set by testing competitive hypotheses rather than selfcontained null hypotheses; the Zhao method was designed to gain power by exploiting eQTL (expression quantitative trait loci) effects. Their suboptimal performance demonstrates the imperative need for an efficient screening test specifically intended for the joint analysis of gene sets by integration of multiplatform genomic data.
The perturbation procedures used in iNOTEchi and iNOTEuni are the main source of computational burden in the omnibus approaches; however, it is worth noting that perturbation procedures resample from the asymptotic null distributions of the genelevel Q statistics, and thus both 1) preserve the covariance structure within and between genes when conducting gene set tests and 2) are far more efficient than permutation procedures requiring direct reshuffling of the data. It is additionally much easier to adjust for covariates using perturbation procedures than using permutations particularly when there exist associations between the genetic data and the covariates. To run one simulation in the MGonly setting using a 2.60GHz Intel Xeon E52670 CPU to test 50 genes with 11 CpG sites and 1000 perturbations each, the approximate computation time is 29, 31, and 32 seconds for iTEGS (any model specification), iNOTEchi, and iNOTEuni respectively. For the Davies approximation to iTEGS, the computation time is about 22, 22, and 25 seconds for the M, MG, and MGC model specifications respectively.
In our data application to the lung cancer survival data, we were able to recover a sizeable number of significant gene sets. Many of these gene sets tended to be least significant when tested under the iTEGS statistic with only the DNAm diseaserisk model specification, but grew increasingly more significant with the inclusion of mRNA gene expression and interaction specifications. This is biologically plausible in that given a true gene pathway, it is highly unlikely that the CpG sites that biologically map within the bounds of a given gene will behave in a strictly linear manner; the remainder of the significant signal in these gene sets can be deduced to arise from the synergistic or antagonistic interaction effects between DNAm and mRNA expression, which are more properly characterized under MGC models. It is also worth noting that a significant gene set identified by one of our methods could be driven by a small subset of very significant gene members (i.e., signals are sparse), whether the signal arises from the main effects of DNAm or RNA expression, or their interactions. Indeed, this is a distinct advantage of our approach, as sparse signals may nonetheless have high biological significance with respect to disease pathways (for example, in the case of CpG loci in gene promoter regions). In these cases, it is useful to conduct further locusbylocus or genebygene analyses characterizing the gene members in the set, as we did for our TCGA application and the Gautschi SRC signaling pathway.
Conclusion
While the iNOTE approaches make fewer assumptions about the underlying causal disease models in a gene set, the tradeoff is an increase in computational burden. Both iNOTE methods are robust to model misspecification and, importantly, performed with close to optimal power across all simulations settings, particularly those in which the gene set is comprised of mixtures of different disease risk models – a highly likely biological scenario.
We propose two efficient procedures for gene set screening which use selfcontained hypothesis tests, and therefore do not rely on the size or proportion of signals within, compared to without, the tested gene set. Furthermore, iNOTE and iTEGS can easily incorporate the adjustments for potential confounding covariates. Our methods dominated two competing methods with respect to power, and further recovered a much greater number of gene sets known a priori to be associated with lung cancer in our scans for gene sets associated with lung cancer survival. In particular, gene sets related to KRAS, EGFR, mirlet7a3 were found to be significantly associated with lung cancer survival. Finally, our methods are easily extended to include more or different genetic platforms. iTEGS and iNOTE software implemented in R in the present manuscript are available in Additional file 2. For any future updated versions, it may also be downloaded via the inote package at https://cran.rproject.org/.
Abbreviations
 DNAm:

DNA methylation
 iTEGS:

integrated total effect of a gene set
 iNOTE:

Integrated network omnibus total effect test
 GSAA:

Gene set association analysis
 M:

Risk model with only methylation effects
 G:

Risk model with only gene expression effects
 MG:

Risk model with only methylation and gene expression main effects
 MGC:

Risk model for methylation, gene expression, and their interactions
References
 1
Tseng GC, Ghosh D, Zhou XJ. Integrating Omics Data. New York: Cambridge University Press; 2015.
 2
Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, Lin X. Powerful SNPset analysis for casecontrol genomewide association studies. Am J Hum Genet. 2010; 86(6):929–42.
 3
Liu L, Lei J, Sanders SJ, Willsey AJ, Kou Y, Cicek AE, Klei L, Lu C, He X, Li M, Muhle RA, Ma’ayan A, Noonan JP, Sestan N, McFadden KA, State MW, Buxbaum JD, Devlin B, Roeder K. DAWN: a framework to identify autism genes and subnetworks using gene expression and genetics. Mol Autism. 2014; 5(1):22.
 4
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, Paulovich A, Pomeroy SL, Golub TR, Lander ES, Mesirov JP. Gene set enrichment analysis: a knowledgebased approach for interpreting genomewide expression profiles. PNAS. 2005; 102(43):15545–50.
 5
Xiong Q, Ancona N, Hauser ER, Mukherjee S, Furey TS. Integrating genetic and gene expression evidence into genomewide association analysis of gene sets. Genome Res. 2012; 22(2):386–97.
 6
Huang YT, Vanderweele TJ, Lin X. Joint analysis of SNP and gene expression data in genetic association studies of complex diseases. Ann Appl Stat. 2014; 8(1):352–76.
 7
Zhao SD, Cai TT, Li H. More powerful genetic association testing via a new statistical framework for integrative genomics. Biometrics. 2014; 70(4):881–90.
 8
Goeman JJ, van de Geer SA, de Kort F, van Houwelingen HC. A global test for groups of genes: testing association with a clinical outcome. Bioinformatics. 2004; 20(1):93–9.
 9
Liu D, Ghosh D, Lin X. Estimation and testing for the effect of a genetic pathway on a disease outcome using logistic kernel machine regression via logistic mixed models. BMC Bioinformatics. 2008; 9(1):292.
 10
Huang YT, Lin X. Gene set analysis using variance component tests. BMC Bioinformatics. 2013; 14(1):210.
 11
Liu X, Jessen WJ, Sivaganesan S, Aronow BJ, Medvedovic M. Bayesian hierarchical model for transcriptional module discovery by jointly modeling gene expression and ChIPchip data. BMC Bioinformatics. 2007; 8(1):283–0.
 12
Savage RS, Ghahramani Z, Griffin JE, de la Cruz BJ, Wild DL. Discovering transcriptional modules by Bayesian data integration. Bioinformatics. 2010; 26(12):158–67.
 13
Kirk P, Griffin JE, Savage RS, Ghahramani Z, Wild DL. Bayesian correlated clustering to integrate multiple datasets. Bioinformatics. 2012; 28(24):3290–7.
 14
Klein HU, Schäfer M, Porse BT, Hasemann MS, Ickstadt K, Dugas M. Integrative analysis of histone ChIPseq and transcription data using Bayesian mixture models. Bioinformatics. 2014; 30(8):1154–62.
 15
Stingo FC, Chen YA, Vannucci M, Barrier M, Mirkes PE. A Bayesian graphical modeling approach to microRNA regulatory network inference. Ann Appl Stat. 2010; 4(4):2024–48.
 16
Badea L. Extracting gene expression profiles common to colon and pancreatic adenocarcinoma using simultaneous nonnegative matrix factorization.Pac Symp Biocomput. 2008; 290(13):279–90.
 17
Kutalik Z, Beckmann JS, Bergmann S. A modular approach for integrative analysis of largescale geneexpression and drugresponse data. Nat Biotechnol. 2008; 26(5):531–9.
 18
Zhang S, Liu CC, Li W, Shen H, Laird PW, Zhou XJ. Discovery of multidimensional modules by integrative analysis of cancer genomic data. Nucleic Acids Res. 2012; 40(19):9379–91.
 19
Huang YT, Liang L, Moffatt MF, Cookson WOCM, Lin X. iGWAS: integrative genomewide association studies of genetic and genomic data for disease susceptibility using mediation analysis. Genet Epidemiol. 2015; 39(5):347–56.
 20
Kong SW, Pu WT, Park PJ. A multivariate approach for integrating genomewide expression data and biological knowledge. Bioinformatics. 2006; 22(19):2373–80.
 21
Lu Y, Liu PY, Xiao P, Deng HW. Hotelling’s T2 multivariate profiling for detecting differential expression in microarrays. Bioinformatics. 2005; 21(14):3105–13.
 22
Breslow NE, Clayton DG. Approximate Inference in Generalized Linear Mixed Models. JASA. 1993; 88(421):9–25.
 23
Satterthwaite FE. An approximate distribution of estimates of variance components. Biometrics. 1946; 2(6):110–4.
 24
Davies RB. Numerical inversion of a characteristic function. Biometrika. 1973; 60(2):415–7.
 25
Parzen MI, Wei LJ, Ying Z. A resampling method based on pivotal estimating functions. Biometrika. 1994; 81(2):341–50.
 26
Cai T, Lin X, Carroll RJ. Identifying genetic marker sets associated with phenotypes via an efficient adaptive score test. Biostatistics. 2012; 13(4):776–90.
 27
Kline P, Santos A. A score based approach to wild bootstrap inference. J Econ Methods. 2012; 1(1):1–19.
 28
Goeman JJ, Buhlmann P. Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics. 2007; 23(8):980–7.
 29
Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in highthroughput experiments. Bioinformatics. 2012; 28(6):882–3.
 30
Gautschi O, Tepper CG, Purnell PR, Izumiya Y, Evans CP, Green TP, Desprez PY, Lara PN, Gandara DR, Mack PC, Kung HJ. Regulation of Id1 expression by SRC: implications for targeting of the bone morphogenetic protein pathway in cancer. Cancer Res. 2008; 68(7):2250–8.
Acknowledgements
Not applicable.
Funding
This work was supported by the Brown Salomon Research Award and Brown Junior Faculty Research Awards in Genetics and Population Studies to SHC and YTH, the National Cancer Institute [5R03CA18293702 to YTH], and the Ministry of Science and Technology [1052118M001014MY3 to YTH].
Availability of data and material
The datasets generated in the simulation and analyzed in the application portion of this study are in whole or part obtained from The Cancer Genome Atlas data portal: http://cancergenome.nih.gov/. iTEGS and iNOTE implemented in R are freely available through the inote package at https://cran.rproject.org/.
Authors’ contributions
SHC contributed to the experimental design, conducted the simulation studies and real data analyses, wrote the manuscript, and developed the coding software for the networklevel integrative methods. YTH conceived the study methodologies, guided the research, wrote the manuscript and provided critical feedback on the manuscript. All authors read and approved the final manuscript.
Ethics approval and consent to participate
Not applicable.
Consent for publication
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Additional files
Additional file 1
Supplementary information. Figure A Internal power simulation across various diseasemodel settings for moderately sized gene sets Figure B Power simulations comparing variancecomponentbased total effect gene set testing procedures to existing methods under mixture diseasemodel settings Table C : Davies approximation pvalues for gene sets signficantly associated with lung cancer in TCGA subjects after Bonferroni correction Table D Counts of overlapping significant BIOCARTA/ KEGG gene sets associated with oneyear lung cancer survival status by iTEGS, iNOTE, and GSAA Table E Counts of overlapping significant lung cancer gene sets associated with pathological stage of tumor at diagnosis by iTEGS, iNOTE, and GSAA; Table E.2: Variance componentbased total effect test pvalues for lung cancer gene sets significantly associated with pathological stage of tumor after Bonferroni correction. (PDF 2410 kb)
Additional file 2
iNote installable R package. (TAR 319 kb)
Rights and permissions
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver(http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
About this article
Received
Accepted
Published
DOI
Keywords
 Pathway analysis
 Data integration
 Epigenetics
 Gene expression
 Gene set analysis
 Integrative genomics