An empirically adjusted weighted ordered p -values meta-analysis method for large-scale simultaneous significance testing in genomic experiments

Abstract

Background

Meta-analysis is a popular approach for combining results from multiple studies investigating the same questions. Meta-analysis has gained wide popularity in genomic analysis due to the availability of large volumes of genomic study results from public databases. In genomic meta-analysis, researchers, often, tend to combine p-values related to significance testing of a gene from multiple studies where thousands of genes are tested simultaneously. The traditional p-value combination approaches aim to find genes which are differentially expressed in at least one of studies. An alternative form of meta-analysis has, recently, gained popularity where the aim is to find genes that are consistently differentially expressed in a large number, possibly a majority, of studies. An approach based on weighted ordered p-values (WOP) has been developed, in the recent past, to perform the latter type of meta-analysis.

Methods

In this article, we discuss the limitations of the WOP meta-analysis method due to its adherence to the standard null distributional assumptions of classical meta-analysis that can lead to incorrect significance testing results. Moreover, we propose a robust meta-analysis method for simultaneous significance testing of multitude of genes that improves the WOP approach using an empirical modification.

Results

Through simulation studies, we demonstrate the superiority of our proposed method over the existing WOP meta-analysis by substantially reducing false discoveries of significant genes and controlling type-I error rates especially in the presence of unobserved confounding variables. We illustrate the utility of our proposed method through a variety of meta-analysis of genomic studies in different diseases.

Keywords

Meta-analysis Fisher's test Stouffer's Z test weighted ordered p-values empirical estimation multiple testing

Introduction

Meta-analysis of multiple genomic studies has become a common practice in recent years due to the fast developments in high-throughput technology and the availability of huge amount of data in public databases.^1–3 The individual genomic studies, usually, have much smaller sample sizes compared to the number of genes which might result in loss of power of statistical detection after adjusting for multiplicities. Meta-analyzing the results from these individual studies can greatly increase the power of the tests and, therefore, has been recognized as an appropriate method for statistical detection and widely adopted in genomic studies.⁴ Although several meta-analysis methods have been developed over the past few years, there are two main approaches for classical meta-analysis.⁵ The first approach is to directly combine p-values of statistical testing from multiple studies. This approach includes popular methods such as the Fisher’s combined probability test,⁶ the Stouffer’s Z test,⁷ and the weighted variations of these methods.⁸ The second approach combines the model-based effect sizes from different studies which includes method such as GeneMeta.⁹ Even though both the approaches have their advantages and disadvantages, the first approach combining p-values is more flexible as it requires less information and assumption from individual component studies. In this article, we focus on the meta-analysis methods that directly combine p-values from multiple studies.

The aim of the classical meta-analysis methods combining p-values from multiple studies, such as the Fisher’s⁶ or the Stouffer’s method,⁷ is to test the alternative hypothesis that at least one of the studies is non-null. For example, in the genomic significance testing context, these methods test the alternative hypothesis that a gene is differentially expressed in at least one study against the null hypothesis that it is not differentially expressed in any of the studies. However, in meta-analysis of genomic studies involving thousands of genes, often it is more reasonable to identify a small subset of genes that are differentially expressed consistently in a majority of studies as opposed to only one or a few studies.⁵ In recent years, a few methods have been developed to address this issue.^5,10,11 These methods order the p-values for the same gene obtained from the different studies and use either a specific $r^{t h}$ ordered p-value¹⁰ or a small subset of ordered p-values leaving out the extreme ones.¹¹ However, such approaches require pre-specification of the choice of $r$ or the proportion of ordered p-values to be combined which can be a subjective choice.^5,10 Moreover, leaving out a number of p-values, during a p-value combination, can lead to considerable loss of information for the purpose of meta-analysis. An alternative approach, proposed by Li and Ghosh,⁵ considers ordered p-values from the multiple studies and combines all of them by weighting them based on their order. Instead of considering only a single ordered p-value or a small subset of p-values, this weighted ordered p-values (WOP) method combines all the ordered p-values giving the maximum weight to the median p-values and down-weighting the extreme ones.⁵ The rationale behind considering this unimodal shape of the distribution of weights is that it is believed that the behavior of the majority of studies is best explained by the p-values that are closer to the center of the distribution.⁵ The WOP method is known to be more robust than the previous approaches since it uses all the available p-value information and one does not need to pre-specify which of the ordered p-values is or are to be used.⁵

Under the WOP framework, although many classical p-value combination approaches can be expressed, Li and Ghosh mainly focused on the two popular tests - Fisher’s combined probability test and the Stouffer’s Z test. One of the key assumptions of the classical p-value combination approaches, including the Fisher’s and the Stouffer’s methods, is that the p-values obtained from the individual studies are uniformly distributed under the null hypothesis of no differential expression. However, in large-scale multiple testing problems, this distributional assumption of the p-values may be violated¹² which can raise questions on the p-value combination approaches of Fisher, Stouffer, and hence the WOP method as well. One of the main differences that exists between a single hypothesis testing framework and a large-scale testing of multiple hypotheses, e.g., 10,000 hypotheses, is in their goals. In a single hypothesis testing framework, one is interested to reject a null hypothesis in favor of an alternative hypothesis with high power, say 90%, whereas in large-scale multiple hypothesis testing of thousands of genes, one often expects to identify a small number of significant genes that can be carried forward for further investigation. That is, it is not desirable to reject 90% of the 10,000 null hypotheses in multiple hypotheses testing of 10,000 genes unlike that in a single hypothesis testing. Rejecting a large proportion of null hypotheses of no significance in simultaneous testing of multiple genes, often, indicates false discoveries which can get multiplied when combining results through meta-analysis. This implies that, although a meta-analysis is usually performed to gain power and decrease type-II error rate, caution should be practiced in large-scale meta-analysis involving thousands of simultaneous hypotheses testing so that the type-I error rates are not inflated, and false discoveries are controlled. This issue of inflated type-I error rate and false discoveries is widespread in large observational genomic studies due to the presence of unobserved variable and hidden confounder effects, e.g., unmeasurable technical artifacts during sequencing, and biological effects such as cell cycle status.¹³ This is because, presence of hidden confounders or unobserved variables dilate the null (hypothesis) distribution, hence violating the underlying theoretical distributional assumptions of the test statistics and the p-values. The consequence can be more serious in the meta-analysis when the theoretical distributional assumptions are violated in many of the individual studies being meta-analyzed.^14,15 Therefore, when p-values are the only sources of information available from the individual studies for meta-analysis, the distributional assumption of the p-values from individual studies needs to be valid or corrected for accurate statistical inference using the WOP meta-analysis method.

One way to solve this problem is to take advantage of the large-scale multiple hypothesis testing framework that allows us to estimate the null distribution empirically and, therefore, eliminating the need to rely on a theoretical asymptotic null distribution.¹² Empirical null distribution can take into account the variation and moderate bias caused by unobserved variables and hidden confounders in large observational studies. In this article we propose a meta-analysis method for large-scale genomic experiments that implements the WOP method while estimating the empirical null distribution parameters through an empirical Bayes approach.¹² Like the WOP method, our proposed method has advantage over the classical meta-analysis methods as it aims to identify genes that are differentially expressed consistently in a majority of studies. Additionally, our proposed method, being able to obtain the null distribution empirically, has much robust performance than the WOP method. Through a variety of simulation studies, we show that our proposed method successfully controls the type I error rate and greatly reduces the false discovery rate (FDR) in comparison to the WOP method, especially in presence of hidden confounding effects in the individual studies. We also show the utility of our proposed method through meta-analysis on three sets of genomic studies on lung cancer, brain cancer, and diabetes.

The rest of the article is organized as follows. In Methods section, we describe our hypothesis setting, the proposed empirically adjusted weighted ordered p-values meta-analysis method, and method for empirical estimation of null distribution. In the Results section, we present results from various simulation scenarios comparing the performances of our proposed method and the WOP method. We also illustrate the application of our proposed meta-analysis method on three different sets of genomic data. We end the article with a Discussion section.

Methods

Hypothesis setting and weighted ordered p-values statistic

Suppose there are $K$ independent studies where each study consists of $G$ genes. Let $θ_{i j}$ denotes the underlying true effect size for the $i^{t h}$ gene in the $j^{t h}$ study, $i = 1,2, \dots, G; j = 1,2, \dots, K$ . That is, $θ_{i j} = 0$ indicates that the $i^{t h}$ gene is not differentially expressed, while a non-zero $θ_{i j}$ indicates that it is differentially expressed, in the $j^{t h}$ study. The goal of our method is to detect genes that are differentially expressed in a majority of studies against the null hypothesis that their true effect sizes are zero in all studies. As a general rule, we target those genes that are differentially expressed in at least half of the studies. That is, for the $i^{t h}$ gene, the hypothesis setting for our meta-analysis method is

{H S}_{m} : {H_{0} : \sum_{j = 1}^{K} I (θ_{i j} \neq 0) = 0 v e r s u s H_{1}^{m} : \sum_{j = 1}^{K} I (θ_{i j} \neq 0) \geq m}

where $m = ⌈ K / 2 ⌉$ , i.e., $m$ is the smallest integer that is not lower than $K / 2$ .

Note that, the hypothesis setting under our meta-analysis method can be generalized for any choice of $m$ ranging from $⌈ K / 2 ⌉ + 1$ to $K$ . Since the WOP meta-analysis method is mostly focused on testing ${H S}_{m}$ for $m = ⌈ K / 2 ⌉$ for simplicity,⁵ we also consider the same choice for $m$ in this article.

Next, we briefly describe the WOP meta-analysis statistic for testing ${H S}_{m}$ for $m = ⌈ K / 2 ⌉$ . Suppose $p_{i j}$ denotes the p-value for testing the null hypothesis $θ_{i j} = 0$ against the alternative that $θ_{i j} \neq 0$ for the $i^{t h}$ gene in the $j^{t h}$ study. For a gene $i$ , the WOP method computes the following summary statistic:

T_{i} = \sum_{j = 1}^{K} w_{j} H (p_{i (j)})

where

p_{i (j)} s

denote the list of p-values, corresponding to the

i^{t h}

gene, ordered over the

K

studies. Here,

w_{j}

denotes the weight corresponding to the

j^{t h}

ordered p-value of the

i^{t h}

gene. Two different weighting schemes based on the binomial distribution are considered – binomial weighting and half-binomial weighting.⁵ The binomial weighting scheme is defined as

w_{j}^{b} = f (j - 1; K - 1,0.5), j = 1,2, \dots, K

, where

f (x; n, π)

denotes the probability mass function of the binomial distribution

B i n (n, π) f o r x = 0,1, \dots, n

. This weighting scheme gives the maximum weight on the median p-value and down-weight the largest and smallest p-values. To further reduce the influence of the smallest p-values on the WOP summary statistic, we considered an alternative weighting scheme called the half binomial weighting scheme. This scheme considers the same weights as in binomial weighting scheme for the median and larger p-values but gives zero weights to the smallest p-values, defined as

w_{j}^{h b} = w_{j}^{b}

for

m \leq j \leq K

and

0

for

j < m

. Further details on the two weighting schemes can be found in Li and Ghosh.⁵ The function

H (.)

in the WOP statistic depends on the choice of the p-value combination method. In particular, two popular p-value combination methods are considered – Fisher’s method⁶ where

H (p_{i (j)}) = - 2 \log (p_{i (j)})

, and Stouffer’s method⁷ where

H (p_{i (j)}) = Φ^{- 1} (1 - p_{i (j)})

Proposed empirically adjusted weighted ordered p-values method

In this section, we describe our proposed meta-analysis method for testing ${H S}_{m}$ for $m = ⌈ K / 2 ⌉$ which empirically modifies the raw p-values from multiple studies and computes multiple testing corrected p-values after appropriately combining them across the studies. This empirical modification of the raw p-values will ensure that the p-values from all studies are uniformly distributed under the null hypotheses, so that the key assumption of the p-value combination methods such as the Fisher’s⁶ and the Stouffer’s⁷ methods is satisfied. Next, we provide the detailed steps of our proposed meta-analysis method:

Assuming that there are $K$ independent studies and $G$ genes in each study,

Step 1

We obtain the p-value $p_{i j}$ for testing the null hypothesis $θ_{i j} = 0$ against the alternative that $θ_{i j} \neq 0$ for gene $i$ in study $j$ , $i = 1,2, \dots, G; j = 1,2, \dots, K$ .

Step 2

We consider the inverse z-transformation to get the corresponding z-scores as follows:

z_{i j} = Φ^{- 1} (p_{i j})

Step 3

The z-scores in step 2 may not follow a standard normal distribution under the null hypotheses. Assuming that the null distribution of the z-scores is normal with mean $δ_{0}$ (not necessarily 0) and standard deviation $σ_{0}$ (not necessarily 1), we estimate the parameters empirically using an empirical Bayes method,¹² details of which can be found in the next section. Suppose ${\hat{δ}}_{0}$ and ${\hat{σ}}_{0}$ are the estimated mean and standard deviation of the null distribution. We modify the z-scores, obtained in step 2, using the estimated parameters as:

z_{i j}^{'} = \frac{z_{i j} - {\hat{δ}}_{0}}{{\hat{σ}}_{0}}

These modified z-scores $z_{i j}^{'}$ are expected to follow a standard normal distribution under the null hypotheses.

Step 4

We convert the empirically adjusted z-scores into corresponding p-values as:

p_{i j}^{'} = Φ (z_{i j}^{'})

Step 5

For a gene $i$ , we order the p-values over the $K$ independent studies. Let $p_{i (j)}^{'}$ denote the $j^{t h}$ ordered p-value for gene $i$ , $i = 1,2, \dots, G; j = 1,2, \dots, K$ . We calculate the WOP summary statistic, as described in the previous section, as follows:

T_{i} = \sum_{j = 1}^{K} w_{j} H (p_{i (j)}^{'})

where

w_{j}

represents the weight corresponding to the

j^{t h}

ordered p-value and

H (.)

denotes the p-value combination method. The choices for the weights

w_{j}

and

H (.)

follow from the previous section.

Step 6

For gene $i$ , since the exact null distribution of the summary statistic, $T_{i}$ , is not readily available, we obtain the p-value, $p^{i},$ by comparing the statistic $T_{i}$ to the numerical distribution by simulating $U (0,1)$ random variables as described below:

(i) We randomly generate p-values from $U (0,1)$ distribution for all the $G$ genes in the $K$ studies. We repeat this data generation process $B$ times. Let $p_{i j}^{(b)}$ denotes the p-value for the $i^{t h}$ gene in the $j^{t h}$ study in the $b^{t h}$ dataset, $i = 1,2, . . ., G; j = 1, 2, \dots, K; b = 1,2, \dots, B$ .

(ii) We calculate the summary statistic $T_{i}^{(b)} = \sum_{j = 1}^{K} w_{j} H (p_{i (j)}^{(b)})$ , for gene $i$ , where $p_{i (j)}^{(b)}$ denotes the $j^{t h}$ ordered p-value for gene $i$ in the $b^{t h}$ dataset.

(iii) For gene $i$ , the p-value corresponding to the summary statistic $T_{i}$ is computed as

p^{i} = \frac{\sum_{b = 1}^{B} I {T_{i}^{(b)} \geq T_{i}}}{B}

Finally, we apply the Benjamini-Hochberg (BH) multiplicity correction method on the set of $G$ p-values, to account for multiple testing^16, and obtain our empirically null-adjusted weighted ordered p-values for the $G$ genes. From this point onwards, this proposed empirically null-adjusted weighted ordered p-values method is referred to as ENWOP method.

Note that, there is an alternative way for obtaining the p-value of the WOP statistic based on permutation analysis.⁵ But this approach requires the original data for each study which, in most situations, are not readily available. Therefore, we focus on the more practical solution which involves computation of the p-values based on the numerical distribution of the WOP statistic as described above.

Empirical estimation of null distribution

Suppose the p-values corresponding to $G$ genes in a study be denoted as $p_{1}, p_{2}, \dots, p_{G}$ . These p-values can be converted into z-scores as $z_{i} = Φ^{- 1} (p_{i})$ , $i = 1,2, \dots, G$ . The null distribution of the z-scores is supposed to be $N (0,1)$ theoretically. However, in large-scale testing situations empirical and theoretical null might differ. The large-scale multiple testing situation enables us to estimate the null distribution of the z-scores. In this section, we will discuss an empirical Bayes method, proposed by Efron,¹² for estimating the null distribution empirically.

The z-scores, corresponding to the $G$ genes, can be categorized into two groups – the “uninteresting” group if the $z_{i}$ is obtained from the null distribution, and the “interesting” group if the $z_{i}$ is obtained from the non-null distribution. Let $π_{0}$ denotes the prior probability of the z-scores belonging to the “uninteresting” group and $π_{1} = 1 - π_{0}$ denotes the prior probability of the z-scores belonging to the “interesting” group. Suppose $f_{0} (z)$ and $f_{1} (z)$ be the densities of the z-scores in the “uninteresting” and the “interesting” groups respectively. The mixture density of the z-scores is defined as $f (z) = π_{0} f_{0} (z) + π_{1} f_{1} (z) .$ Following Bayes theorem, the a posteriori probability of belonging to the “uninteresting” group given $z$ is $\Pr [“ u n i n t e r e s t i n g ” | z] = π_{0} f_{0} (z) / f (z)$ . The null density, $f_{0}$ is estimated from the central peak of the histogram of the z-scores. Assuming that $f_{0}$ is a normal distribution with mean $δ_{0}$ and standard deviation $σ_{0}$ , for z-scores close to zero, we can write $\log (f (z)) = - \frac{1}{2} {(\frac{z - δ_{0}}{σ_{0}})}^{2} + c o n s t a n t$ . The parameters of $f_{0}$ are estimated as: $δ_{0} = \arg \max {f (z)}$ and $σ_{0} = {[- \frac{d^{2}}{d z^{2}} \log f (z)]}_{δ_{0}}^{- 1 / 2}$ .

However, the above estimate of $σ_{0}$ can be unstable.¹² Therefore, a smoothing step is applied where a quadratic curve $a_{0} + a_{1} x_{k} + a_{2} x_{k}^{2}$ is fitted by ordinary least squares to the estimated $\log (f (x_{k}))$ values, for $x_{k}$ within 1.5 units of the maximum $δ_{0}$ , yielding the final estimate of $σ_{0}$ as ${[- 2 a_{2}]}^{- 1 / 2}$ . This method of estimation of the parameters of null distribution is called the method of “central-matching”. More details about this method can be found in Efron¹² and Efron.¹⁷

Results

Simulation studies

We conducted simulation studies to evaluate the performance of our proposed method ENWOP for accurate identification of significant genes in a majority of studies. We simulated continuous gene expression datasets for multiple independent studies. Details of the data generation process are given below.

We considered 10 independent studies each involving continuous gene expression levels for 3000 genes, i.e., $K = 10$ and $G = 3000$ . We considered two groups of subjects in each study where each group consists of 20 subjects, i.e., $n_{1} = n_{2} = 20 .$ We considered 50 genes as differentially expressed between the two subject groups in $1,2, \dots, 10$ studies respectively. That is, in total, 500 genes are differentially expressed between the subject groups in at least one study. Since our alternative hypothesis for a gene is that it is differentially expressed in at least five studies, we aim to identify only the 300 genes $(10 %)$ that are differentially expressed in at least five of the studies.

We generated the (log) expression level for the $i^{t h}$ gene, $l^{t h}$ subject in the $k^{t h}$ group for each study separately using the following model:

y_{i k l} = μ + G_{i} + V_{k} + {G V}_{i k} + W_{i k l} + e_{i k l}

Here

μ

denotes the overall mean effect,

G_{i}

denotes the effect due to the

i^{t h}

gene,

V_{k}

denotes the effect due to the

k^{t h}

subject group, and

{G V}_{i k}

denotes the interaction effect between the

i^{t h}

gene and the

k^{t h}

subject group,

W_{i k l}

denotes the effect of a hidden variable and

e_{i k l}

denotes the error component corresponding to the

i^{t h}

gene,

l^{t h}

subject in the

k^{t h}

group,

i = 1,2, \dots, G; k = 1,2; l = 1,2, \dots, n_{k}

For our simulations, we considered $μ$ , $G_{i}$ , and $V_{k}$ as zero for all $i, k$ and $l$ , for simplicity. Note that, we considered 50 genes to be differentially expressed between the two subject groups in $1,2, \dots, 10$ studies respectively. The differences in magnitudes of (log) expression values of these genes between the two groups are considered as 8, which are obtained through the generation of the interaction terms between the genes and the groups, ${(G V)}_{i k} s$ , as follows:

For study $j$ ,

{(G V)}_{i 1} = - 4, {(G V)}_{i 2} = 4 f o r i = 1, \dots, 25 j

{(G V)}_{i 1} = 4, {(G V)}_{i 2} = - 4 f o r i = 25 j + 1, \dots, 50 j

{(G V)}_{i 1} = {(G V)}_{i 2} = 0 f o r i = 50 j + 1, \dots, G

where

j = 1, \dots, K .

In our simulations, we assumed the presence of a hidden variable which acts as a confounder. The effect of the hidden confounder for the $i^{t h}$ gene, $l^{t h}$ subject in the $k^{t h}$ group was generated such that it varied over the two subject groups, different groups of genes as well as over different studies. We considered $W_{i k l} = u_{i k l} I (s_{i k l} = 1)$ , where $s_{i k l} \sim B e r n o u l l i (0.4)$ and $u_{i k l}$ are generated depending on the gene, subject group and the study ID $j$ as given below:

u_{i k l} = {\begin{array}{c} N (- 1 + j + δ . I, {0.01}^{2}) f o r i = 1, \dots, 25 j; l = 1, \dots, n_{k} \\ N (2 + j + δ . I, {0.01}^{2}) f o r i = 25 j + 1, \dots, 50 j; l = 1, \dots, n_{k} \\ N (5 + j + δ . I, {0.01}^{2}) f o r i = 50 j + 1, \dots, G; l = 1, \dots, n_{k} \end{array}

where $i = 1,2, \dots, G$ ; $k = 1,2; j = 1,2, \dots, K$ . Here, $I = 1$ for group $1$ (i.e., for $k = 1$ ) and $I = 0$ for group $2$ (i.e., for $k = 2$ ). The magnitude of the difference between the means of the distributions of $u_{i k l}$ between the two subject groups is given by $δ$ . In our simulations, we considered $δ = 4$ .

We introduced correlations among some of the genes through the generation of the error terms $e_{i k l} s$ as described below:

We considered four groups of correlated genes given by $C_{1} = \{1, 2, \dots, 30\}, C_{2} = \{121, 122, \dots, 180\}, C_{3} = \{1501,1502, \dots, 1560\}$ and $C_{4} = {2671, 2672, \dots, 2730}$ .

The error term $e_{i k l}$ for the $i^{t h}$ gene, $l^{t h}$ subject in the $k^{t h}$ group is generated as:

e_{i k l} = {\begin{array}{c} \frac{1}{\sqrt{2}} e_{i k l}^{1} + \frac{1}{\sqrt{2}} e_{i k l}^{2} & i f i \in \{C_{1}, C_{2}, C_{3}, C_{4}\} \\ e_{i k l}^{2} & o . w \end{array}, i = 1,2, \dots, G; k = 1,2; l = 1,2, \dots, n_{k}

where $e^{1}$ are generated independently from $N (0, 1)$ in such a way that the values of $e^{1}$ are same for all the genes belonging to the same group, and $e^{2}$ are generated independently from $N (0, 2^{2})$ .

For each study, after generating the (log) gene expression values for all the subjects, we tested whether the genes are differentially expressed between the two subject groups using “limma” in Bioconductor¹⁸ and stored the raw p-values. We then applied the ENWOP method and obtained the list of differentially expressed genes with a BH adjusted p-value cutoff of 0.05. We obtained the type I error rate for the ENWOP method as well as evaluated its performance using sensitivity, specificity, and FDR based on 500 Monte-Carlo iterations. Because of the hypothesis setup, as described in the Methods section, false positives and type I error are not the same. Type I error is rejecting the null hypothesis for a gene that is not differentially expressed in any of the studies while a false positive is rejecting the null hypothesis for a gene that is differentially expressed in less than $m$ studies.⁵ We, additionally, compared the performance of the ENWOP method with the WOP method (without any empirical adjustment).⁵

Table 1 summarizes the simulation results for the ENWOP method and the corresponding WOP method with the two choices for the p-value combination approach (Fisher⁶ and Stouffer⁷) and the two weighting schemes (binomial and half-binomial), as discussed in the Methods section. The type I error rates for the ENWOP method are controlled at 0.05 but the WOP method has extremely high type I error rates in all settings. The ENWOP method also has significantly lower FDR values compared to the WOP method. Although the proposed ENWOP method has slightly lower sensitivity values, the specificity values are much higher compared to the WOP method. In general, the half binomial weighting scheme have slightly lower sensitivity and slightly higher specificity values compared to the binomial weighting scheme. The FDR values for both methods and the type I error rates for the WOP method are also slightly lower for half binomial weighting scheme. The results did not vary significantly between the choices of the p-value combination approach.

Table 1.

Performances of the proposed method (ENWOP) and the WOP method in presence of hidden confounder. Type I error rate, sensitivity, specificity, and FDR values are obtained based on 500 Monte-Carlo iterations. The proportion of differentially expressed genes between the two subject groups is 10%.

Weighting scheme	p-value combination approach	Method	Type I error	Sensitivity	Specificity	FDR
Binomial	Fisher	ENWOP	0.038	0.955	0.983	0.138
	Fisher	WOP	0.530	0.999	0.767	0.677
	Stouffer	ENWOP	0.040	0.945	0.984	0.130
	Stouffer	WOP	0.518	0.997	0.784	0.661
Half-binomial	Fisher	ENWOP	0.041	0.926	0.986	0.120
	Fisher	WOP	0.489	0.992	0.815	0.627
	Stouffer	ENWOP	0.042	0.917	0.987	0.116
	Stouffer	WOP	0.483	0.990	0.822	0.618

We, additionally, varied the proportion of differentially expressed genes between the two subject groups, ranging from 5% to 20%. Figure 1 shows the performances of the ENWOP method as well as the WOP method with varying proportion of differentially expressed genes in presence of hidden confounder in the studies. The ENWOP method has type I error rates controlled at 0.05 consistently in all settings, while the type I error rates of the WOP method are unacceptably high in the same simulated settings (see Figure 1). The FDR values of the ENWOP method are also much lower than those of the WOP method in all settings. Although the FDR values for the WOP method tend to decrease with the increase in the proportion of differentially expressed genes, they are substantially worse (higher) than our proposed ENWOP in every situation. The ENWOP method has slightly lower sensitivity values compared to the WOP method but the values increased when the proportion of differentially expressed genes is increased. The specificity values of the ENWOP method are close to one in all settings but the WOP method has much lower specificity values which further decreased as the proportion of differentially expressed genes is increased. Overall, our proposed ENWOP method significantly outperforms the WOP approach in all situations of hidden confounder presence by accurately identifying truly differential genes.

Figure 1.

Performances of the ENWOP method and the WOP method in presence of hidden confounder with varying proportion of differentially expressed genes between two subject groups. Type I error rate, sensitivity, specificity, and FDR values are obtained based on 500 Monte-Carlo iterations. WOP-BF: WOP method with binomial weighting scheme and Fisher’s p-value combination approach; WOP-BS: WOP method with binomial weighting scheme and Stouffer’s p-value combination approach; WOP-HBF: WOP method with half-binomial weighting scheme and Fisher’s p-value combination approach; WOP-HBS: WOP method with half-binomial weighting scheme and Stouffer’s p-value combination approach; ENWOP-BF: ENWOP method with binomial weighting scheme and Fisher’s p-value combination approach; ENWOP-BS: ENWOP method with binomial weighting scheme and Stouffer’s p-value combination approach; ENWOP-HBF: ENWOP method with half-binomial weighting scheme and Fisher’s p-value combination approach; ENWOP-HBS: ENWOP method with half-binomial weighting scheme and Stouffer’s p-value combination approach.

In addition to the simulated sceanrios with hidden confounders, we also considered some other variations in our simulation models. In particular, we considered a simulation scenario where we assumed the presence of a hidden variable that does not act as a confounder as well as a simulation scenario where there did not exist any effect of a hidden variable/confounder in the studies. Both scenarios are decribed below.

Presence of a hidden variable that does not act as a confounder

In this simulation scenario, we assumed the presence of a hidden variable which affects the outcome but does not vary between the two subject groups. We generated the distribution of the hidden variable for the $i^{t h}$ gene, $l^{t h}$ subject in the $k^{t h}$ group, as $W_{i k l} = u_{i k l} I (s_{i k l} = 1)$ , where $s_{i k l} \sim B e r n o u l l i (0.4)$ and $u_{i k l}$ are generated as given below:

u_{i k l} = N (- 4 + j, {0.1}^{2}) f o r i = 1, \dots, G; k = 1, 2; l = 1, \dots, n_{k}; j = 1, \dots, K

We considered 10% of the genes as differentially expressed between the two subject groups in at least five studies as considered before. The differences in magnitudes of (log) expression values of these differentially expressed genes are considered as two. All the other terms in the model for simulation are generated in the same way as described previously.

Supplementary Table 1 shows the results for the ENWOP method and the WOP method. In this simulation scenario, both the methods have controlled type I error rates under all settings. The FDR values of the ENWOP method are slightly lower compared to the WOP method. Both methods have very similar sensitivity and specificity values. The methods with half binomial weighting scheme have lower sensitivity as well as FDR values compared to those with binomial weighting scheme, similar to what we observed in presence of hidden confounder.

The performances of the two methods with varying proportion of differentially expressed genes, in the presence of hidden variable that does not act as confounder, are shown in Supplementary Figure 1. The type I error rates remained controlled at 0.05 for both methods consistently in all settings. Both the methods have similar FDR values for smaller proportion of differentially expressed genes, but the FDR values of the WOP method slightly increased with increase in the proportion of differentially expressed genes. The sensitivity values are very similar for both methods with half binomial weighting scheme having slightly lower values than binomial weighting scheme. Both methods have very similar specificity values.

No effect of any hidden variable

In this simulation scenario, we assumed that there is no effect of any hidden variable/confounder in the studies. Therefore, we set $W_{i k l} = 0$ , for all $i, k$ and $l$ . Here also, we considered 10% of the genes as differentially expressed between the two subject groups in at least five studies and the differences in magnitudes of (log) expression values of these genes are considered as two. The random error term ( $e_{i k l}$ ) are generated as before with $e^{1}$ drawn independently from $N (0, {0.5}^{2})$ and $e^{2}$ generated independently from $N (0, {5.5}^{2})$ . All the other terms in the model for simulation are generated in the same way as described previously.

Supplementary Table 2 shows the results for the ENWOP method as well as the WOP method. The type I error rates are controlled for both methods under all settings. The sensitivity and FDR values are slightly lower for the ENWOP method compared to the WOP method. Both methods have very similar specificity values. The half binomial weighting scheme have slightly lower sensitivity as well as FDR values compared to the binomial weighting scheme, consistent with what observed in the previous simulation scenarios.

Supplementary Figure 2 shows the simulation results for the two methods with varying proportion of differentially expressed genes when there is no effect of hidden variable/confounder in the studies. The performances of the methods are very similar to what we observed in the previous scenario in the presence of hidden variable that does not act as confounder.

An application to lung cancer studies

We conducted meta-analysis using the ENWOP method on five lung cancer gene expression datasets.^14,19 These five studies will be referred to as Bhattacharjee,²⁰ GSE11969,²¹ GSE29016,²² GSE30219,²³ and GSE43580,²⁴ respectively. Each of the datasets contain normalized expression levels for 7200 genes, and subjects with different types of lung cancer. We aimed to identify the genes that are differentially expressed between two lung cancer types - adenocarcinoma (AD) and squamous cell carcinoma (SQ) in at least three out of the five studies. The sample size of each study is given in Supplementary Table 3(a). To obtain the original p-values for the genes for each dataset, we tested for differential expression between AD and SQ subjects using “limma”.¹⁸ We applied the ENWOP method following the steps described in the Methods section. The empirically estimated mean and the standard deviation of the original z-scores are –0.83 and 1.99, respectively. That is, the empirical null distribution of the original z-scores is much different from the theoretical null distribution of $N (0,1)$ . For comparison, we also applied the WOP method to identify the differentially expressed genes between the two lung cancer types in at least three studies.

We considered a gene significant if the BH adjusted p-value is less than 0.05. Table 2 summarizes the number of differentially expressed genes, identified by the ENWOP method as well as the WOP method, for the two choices of p-value combination approach (Fisher⁶ and Stouffer⁷) and two weighting schemes (binomial and half binomial). The WOP method identified much higher number of significant genes compared to the proposed ENWOP method at each combination of p-value combination approach and weighting scheme. For the WOP method, the maximum number of significant genes (68.3%) is identified by the Fisher’s p-value combination approach with binomial weighting scheme. The ENWOP method with Stouffer’s p-value combination approach and binomial weighting scheme identified the maximum number of significant genes, which is 1474 (20.4%). In all settings, the WOP method identified more than 58% of the genes as significant which is unlikely to be accurate and indicates a possibility of large number of false discoveries. Consistent with our simulation results, the methods with binomial weighting scheme identified higher number of significant genes compared to the methods with half binomial weighting scheme.

Table 2.

The number of significant genes (percentage) identified by our proposed method ENWOP and the WOP method with two choices of p-value combination approaches and two weighting schemes for the lung cancer meta-analysis.

Weighting scheme	p-value combination approach	Method	Number of significant genes (percentage)
Binomial	Fisher	ENWOP	1406 (19.5%)
	Fisher	WOP	4921 (68.3%)
	Stouffer	ENWOP	1474 (20.4%)
	Stouffer	WOP	4672 (64.9%)
Half-binomial	Fisher	ENWOP	1317 (18.3%)
	Fisher	WOP	4286 (59.5%)
	Stouffer	ENWOP	1371 (19.0%)
	Stouffer	WOP	4208 (58.4%)

Figure 2 shows the overlap between the number of significant genes identified by the ENWOP method and the WOP method with both weighting schemes using the (a) Fisher’s, and (b) Stouffer’s p-value combination approaches for the lung cancer meta-analysis. For the Fisher’s p-value combination approach, there are 1234 genes which are identified by both methods with both weighting schemes (see Figure 2(a)). Similarly, for the Stouffer’s p-value combination approach, 1312 genes are identified by both methods with both weighting schemes (see Figure 2(b)). Using the same pair of p-value combination approach and weighting scheme, all the genes identified by the ENWOP method are also identified by the WOP method. In order to identify biological pathways associated with the gene lists identified by the ENWOP method with both p-value combination approaches and weighting schemes, we performed functional annotation clustering using Database for Annotation, Visualization and Integrated Discovery (DAVID) software.²⁵ The ENWOP method identified several biologically relevant pathways including cell cycle, DNA replication, and p53 signaling pathway based on the BH adjusted p-value cutoff of 0.05. We further investigated some genes which are identified as significant by the WOP method but not by our ENWOP method. For example, the gene with entrez ID 11131 was identified by the WOP method with both weighting schemes and p-value combination approaches, but not by the ENWOP method. Supplementary Figure 3 shows the box plots of the expression levels of the gene for the two cancer types (AD and SQ) in each of the five studies. From the boxplots, we observe that there is no significant difference in the expression levels of the gene between the two cancer types. Furthermore, based on the p-values from the differential expression analysis in the individual studies, the gene was not significant at p-value cutoff of 0.05 in any of the studies. This suggests that this gene is unlikely to be an important factor for the difference between AD and SQ patients, and therefore, it is reasonable that the proposed ENWOP method did not identify the gene. Since the genes identified by the WOP method is a superset of those identified by the ENWOP method, these findings are consistent with our simulation results that showed that the sensitivity of the WOP method can be similar or slightly higher than that of the ENWOP method but at the cost of very high false discovery rates and type-1 error rates.

Figure 2.

Venn diagram showing the overlaps between the number of significant genes identified by the ENWOP method and the WOP method with both weighting schemes using the (a) Fisher’s, and (b) Stouffer’s p-value combination approaches for the lung cancer meta-analysis. WOP-BF: WOP method with binomial weighting scheme and Fisher’s p-value combination approach; WOP-BS: WOP method with binomial weighting scheme and Stouffer’s p-value combination approach; WOP-HBF: WOP method with half-binomial weighting scheme and Fisher’s p-value combination approach; WOP-HBS: WOP method with half-binomial weighting scheme and Stouffer’s p-value combination approach; ENWOP-BF: ENWOP method with binomial weighting scheme and Fisher’s p-value combination approach; ENWOP-BS: ENWOP method with binomial weighting scheme and Stouffer’s p-value combination approach; ENWOP-HBF: ENWOP method with half-binomial weighting scheme and Fisher’s p-value combination approach; ENWOP-HBS: ENWOP method with half-binomial weighting scheme and Stouffer’s p-value combination approach.

Further applications

To further compare the performance of the ENWOP method with the WOP method, we considered two other micro-array meta-analyses. The first meta-analysis consists of seven studies on brain cancer and the second consists of 16 studies on diabetes. These datasets were previously analyzed by Li and Ghosh⁵ and Song and Tseng¹⁰ and we obtained the p-value for differential expression for each gene in each study from the supplementary material of Song and Tseng.¹⁰ The seven studies in the brain cancer meta-analysis tested for differential expression for 5836 genes between two subtypes of brain tumors – anaplastic astrocytoma (AA) and glioblastoma multiforme (GBM). The studies in the diabetes meta-analysis tested for differential expression for 6645 genes between different tissues in mice and human subjects. We excluded the two studies from our analysis which included human subjects and meta-analyzed the remaining 14 studies on mice for consistency. Details about the sample sizes of the studies can be found in Supplementary Table 3(b)-(c).

For the brain cancer meta-analysis, we aimed to identify the differentially expressed genes between AA and GBM in at least four out of the seven studies. The empirically estimated mean and the standard deviation of the null distribution of the z-scores are −0.36 and 1.42, respectively, which indicates that the empirical null is much different from the theoretical null distribution of

N (0,1)

. Table 3 shows the numbers of differentially expressed genes, identified by the ENWOP method and the WOP method, for the two different choices of p-value combination approach and the two weighting schemes. In all settings, the WOP method identified much higher number of significant genes compared to the ENWOP method, similar to our previous findings with lung cancer meta-analysis. Supplementary Figure 4 shows the overlap between the number of significant genes identified by the ENWOP method and the WOP method with both weighting schemes using the (a) Fisher’s, and (b) Stouffer’s p-value combination approaches. For the Fisher’s p-value combination approach, there are 498 genes which are identified by both methods with both weighting schemes (see Figure 2(a)), and for the Stouffer’s p-value combination approach, 482 genes are identified by both methods with both weighting schemes (see Figure 2(b)). Functional annotation clustering with genes identified by the ENWOP method identified important biologically relevent pathways at BH adjusted p-value cutoff of 0.05, such as PI3K-Akt signaling pathway, p53 signaling pathway, cell cycle, TNF signaling pathway, and Wnt signaling pathway.

Table 3.

The number of significant genes (percentage) identified by the ENWOP method and the WOP method with two choices of p-value combination approaches and two weighting schemes for the brain cancer meta-analysis.

Weighting scheme	p-value combination approach	Method	Number of significant genes (percentage)
Binomial	Fisher	ENWOP	635 (10.9%)
	Fisher	WOP	2462 (42.2%)
	Stouffer	ENWOP	613 (10.5%)
	Stouffer	WOP	2260 (38.7%)
Half-binomial	Fisher	ENWOP	518 (8.9%)
	Fisher	WOP	1885 (32.3%)
	Stouffer	ENWOP	503 (8.6%)
	Stouffer	WOP	1796 (30.8%)

For the diabetes meta-analysis, we aimed to identify the genes that are differentially expressed in at least seven out of 14 studies. The estimated mean and the standard deviation of the empirical null distribution of the z-scores are −0.21 and 1.12, respectively, which are much deviated from the theoretical null distribution. Table 4 includes the numbers of differentially expressed genes, identified by the ENWOP method and the WOP method, for all the settings. Consistent with the previous analyses, the ENWOP method identified much lesser number of differentially expressed genes compared to the WOP method. Supplementary Figure 5 shows the overlap between the number of significant genes identified by the ENWOP method and the WOP method with both weighting schemes using the (a) Fisher’s, and (b) Stouffer’s p-value combination approaches. Functional annotation clustering with genes identified by the ENWOP method identified several biological processes related to diabetes, at BH adjusted p-value cutoff of 0.05, including response to lipopolysaccharide, response to estradiol and response to glucocorticoid.

Table 4.

The number of significant genes (percentage) identified by ENWOP method and the WOP method with two choices of p-value combination approaches and two weighting schemes for the diabetes data.

Weighting scheme	p-value combination approach	Method	Number of significant genes (percentage)
Binomial	Fisher	ENWOP	152 (2.3%)
	Fisher	WOP	1341 (20.2%)
	Stouffer	ENWOP	157 (2.4%)
	Stouffer	WOP	1258 (18.9%)
Half-binomial	Fisher	ENWOP	140 (2.1%)
	Fisher	WOP	1085 (16.3%)
	Stouffer	ENWOP	143 (2.2%)
	Stouffer	WOP	1084 (16.3%)

Discussion

Meta-analysis is a popular statistical technique for combining results from multiple studies. Such meta-analyses are frequently performed in various domains of scientific studies including medicine, epidemiology, and psychology among others. In genomic significance testing studies related to a disease, meta-analysis is often considered an important avenue to overcome the aspect of low power of detection from individual studies due to the smaller sample sizes compared to the number of genes or hypotheses to be tested. In these scenarios of large-scale hypotheses testing, combination of p-values obtained from individual studies has been one of the most popular approaches for meta-analysis. The classical p-value combination methods have been developed with the goal of declaring a gene to be overall significant if it is found to be statistically significant in at least one of the combining studies. However, such p-value combination methods have been recently criticized for not being able to capture the consistent pattern of signals from genes across multiple studies^5,10 and an alternative form of meta-analysis has been proposed with the aim of finding genes that are consistently identified significant in multiple studies. The weighted ordered p-values (WOP) method is one of the important approaches that serve the latter form of meta-analysis. In this article, we have shown that the existing WOP method of meta-analysis can suffer from inflated type-I error rates and substantial false discoveries, even after multiplicity corrections, in large-scale simultaneous hypotheses testing of thousands of genes especially in the presence of unobserved variables and confounder effects. As a remedy to this problem we have proposed a meta-analysis approach, called ENWOP method, that empirically estimates the null distribution of the test statistic using the large-scale nature of the simultaneous hypotheses testing, as opposed to assuming a theoretical null distribution which can be violated in the presence of hidden confounder effects. The ENWOP method caters to the same meta-analysis aim as the WOP method but has superior performances than the WOP method especially in the presence of potentially hidden variable effects as seen through numerous simulated scenarios. Using a variety of real genomic data meta-analysis, we have demonstrated the usefulness of the proposed ENWOP method in identifying significant and biologically releveant genes for different diseases.

Our proposed ENWOP method is aimed at those situations where p-value is the only source of information that is consistently available for all the genes in all of the studies considered for meta-analysis. There exists another type of meta-analysis where model based approaches are employed for pooling results from multiple studies.²⁶ Since the latter type of meta-analysis requires substantially more information, e.g. effect sizes, than just the p-values in consistent pattern from all the individual studies, we have not focused on that type of meta-analysis in this article. Note that, in the rare event where the raw expression values are readily available for all the genes from every study, researchers have the option of performing the existing WOP meta-analysis through a permutation based hypotheses testing⁵ as an alternative to our proposed ENWOP method. However, such situations of availability of complete expression levels from all the studies are not commonly encountered in a genomic meta-analysis which makes the permutation based WOP approach infeasible in most cases. This, again, highlights the importance of our proposed method in producing an accurate meta-analysis in the absence of abundant information from the component studies.

Supplemental Material

Supplemental Material - An empirically adjusted weighted ordered p-values meta-analysis method for large-scale simultaneous significance testing in genomic experiments

Supplemental Material for An empirically adjusted weighted ordered p-values meta-analysis method for large-scale simultaneous significance testing in genomic experiments by Wimarsha T Jayanetti, N Rao Chaganty, and Sinjini Sikdar in Research Methods in Medicine & Health Sciences.

Footnotes

Declaration of conflicting interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) received no financial support for the research, authorship, and/or publication of this article.

Data availability statement

The that support the findings of this study are publicly available.^10,14,¹⁹ An R code for the proposed method is available as a supplementary material.

Supplemental Material

Supplemental material for this article is available online.

References

Kröger

Mapiye

Entfellner

JBD

, et al. A meta-analysis of public microarray data identifies gene regulatory pathways deregulated in peripheral blood mononuclear cells from individuals with systemic lupus erythematosus compared to those without. BMC Med Genom 2016; 9(1): 66.

Sikdar

Joehanes

Joubert

, et al. Comparison of smoking-related DNA methylation between newborns from prenatal exposure and adults from personal smoking. Epigenomics 2019; 11(13): 1487–1500.

Karim

Bradburn

Roberts

ACCEPTS Study , et al. First-trimester ultrasound detection of fetal heart anomalies: systematic review and meta-analysis. Ultrasound Obstet Gynecol 2022; 59(1): 11–25.

Panagiotou

Willer

Hirschhorn

, et al. The power of meta-analysis in genome wide association studies. Annu Rev Genom Hum Genet 2013; 14: 441–465.

Ghosh

. Meta-analysis based on weighted ordered p-values for genomic data with heterogeneity. BMC Bioinf 2014; 15: 226.

Fisher

. Statistical methods for research workers. London, UK: Oliver & Boyd, 1932.

Stouffer

Suchman

DeVinney

, et al. The American soldier: adjustment during army life. Princeton, NJ: Princeton University Press, 1949.

Willer

Abecasis

. METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 2010; 26(17): 2190–2191.

Choi

Kim

, et al. Combining multiple microarray studies and modeling interstudy variation. Bioinformatics 2003; 19(1): i84–90.

10.

Song

Tseng

. Hypothesis setting and order statistic for robust genomic meta-analysis. Ann Appl Stat 2014; 8(2): 777–800.

11.

Olkin

Saner

. Approximations for trimmed Fisher procedures in research synthesis. Stat Methods Med Res 2001; 10(4): 267–276.

12.

Efron

. Large-scale simultaneous hypothesis testing: the choice of a null hypothesis. J Am Stat Assoc 2004; 99(465): 96–104.

13.

Chen

Zhou

. Controlling for confounding effects in single cell RNA sequencing studies using both control and target genes. Sci Rep 2017; 7: 13587.

14.

Sikdar

Datta

. EAMA: empirically adjusted meta-analysis for large-scale simultaneous hypothesis testing in genomic experiments. PLoS One 2017; 12(10): e0187287.

15.

Sikdar

. Robust meta-analysis for large-scale genomic experiments based on an empirical approach. BMC Med Res Methodol 2022; 22(1): 43.

16.

Benjamini

Hochberg

. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol 1995; 57(1): 289–300.

17.

Efron

. Size, power and false discovery rates. Ann Stat 2007; 35(4): 1351–1377.

18.

Ritchie

Phipson

, et al. Limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res 2015; 43(7): e47.

19.

Hughey

Butte

. Robust meta-analysis of gene expression using the elastic net. Nucleic Acids Res 2015; 43(12): e79.

20.

Bhattacharjee

Richards

Staunton

, et al. Classification of human lung carcinomas by mRNA expression profiling reveals distinct adenocarcinoma subclasses. Proc Natl Acad Sci USA 2001; 98(24): 13790–13795.

21.

Takeuchi

Tomida

Yatabe

, et al. Expression profile–defined classification of lung adenocarcinoma shows close relationship with underlying major genetic changes and clinicopathologic behaviors. J Clin Oncol 2006; 24(11): 1679–1688.

22.

Staaf

Jönsson

, et al. Relation between smoking history and gene expression profiles in lung adenocarcinomas. BMC Med Genom 2012; 5: 22.

23.

Rousseaux

Debernardi

Jacquiau

, et al. Ectopic activation of germline and placental genes identifies aggressive metastasis-prone lung cancers. Sci Transl Med 2013; 5(186): 186ra66.

24.

Tarca

Lauria

Unger

IMPROVER DSC Collaborators , et al. Strengths and limitations of microarray-based phenotype prediction: lessons learned from the IMPROVER diagnostic signature challenge. Bioinformatics 2013; 29(22): 2892–2899.

25.

Dennis

Jr Sherman

Hosack

, et al. DAVID: database for annotation, visualization, and integrated discovery. Genome Biol 2003; 4(5): P3.

26.

Han

Eskin

. Random-effects model aimed at discovering associations in meta-analysis of genome-wide association studies. Am J Hum Genet 2011; 88(5): 586–598.

Supplementary Material

Please find the following supplemental material available below.

For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.

For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.

0.00 MB

23.95 MB