1. Introduction

The problem of statistically testing mean difference for each of thousands of variables is commonly encountered in genomic studies. For example, the popularly applied microarray technology allows the gene expression study of tens of thousands of genes simultaneously. The recent advance of next-generation sequencing technology allows the measurement of gene expression in an even higher dimension. These high-throughput technologies have revolutionized the way genomic studies progress and provided rich data to explore. However, these experiments are expensive, and as a consequence, such experiments typically involve only a few samples for each treatment group. This results in the “large , small ” problem for hypothesis testing, and the power of the statistical tests can be very low after controlling the multiple testing error, such as the false discovery rate (FDR).

The normalized signal intensities from microarray experiments are generally assumed to follow normal distributions [1–4]. The recently emerging next-generation sequencing data may also be modeled approximately using normal distributions, when the number of reads are large or under certain transformation [5]. Thus multiple testing problem for normal means has wide applications in genetic and genomic studies, and it is also a general statistical question of interest.

Several testing procedures have been proposed in the context of microarray study, including the SAM test [6], Efron's -test [7], the regularized -test [8], the -statistic [1] and its multivariate counterpart, the -statistic [9], the test of Wright and Simon [10], the moderated -test [2], the test [3] and the test of [11] which is similar to the test, the test [4], and the LEMMA test [12]. Although numerous procedures have been proposed, very few can be justified to achieve the optimal power. Among these procedures, Hwang and Liu [4] proposed a framework and showed that an optimal testing procedure can be derived within such a framework. They also proposed a test with maximum average power (the MAP test) and an approximated version, the test. Here the optimality was defined in terms of maximizing the power averaged across all tests for which the null hypotheses are false while controlling FDR. This method provides theoretical guide for developing optimal multiple testing procedures. The popularly applied moderated -statistic developed by Smyth [2] can also be shown to achieve optimal power asymptotically under different distributional assumptions from the test. Both the moderated -statistic and the test assume that the mean expression levels (or the mean of interesting contrasts) of all genes follow a normal distribution although the parameters for this distribution vary between the two tests. Yet in practice such distribution depends on the population of genes selected in a particular study and often does not follow the prespecified parametric distribution. This raises concerns about the robustness of the moderated and the tests.

The objective of this paper is to develop an optimal and robust multiple testing procedure without any distributional assumptions on the mean. As in Hwang and Liu [4], the optimality is defined in terms of maximizing the power averaged across all tests for which the null hypotheses are false while controlling FDR. We develop a semiparametric optimal testing procedure which we abbreviate as the SPOT procedure. The distribution of the mean expression across genes is not assumed to follow a parametric model which makes our method robust to violations to normal assumptions. We find that the SPOT procedure works very well in simulation studies and in an analysis of real microarray data with spike-in probes.

The remaining of this paper is organized as follows. We first introduce necessary notations in Section 2. Then, in Section 3, we describe the general concepts of optimal testing procedures. We propose our semiparametric optimal testing (SPOT) procedure in Section 4 and describe its implementation in Section 5. Section 6 presents simulation studies. Section 7 shows the analysis result of a real microarray dataset. Section 8 provides a summary of this paper.

2. Notations

An appropriate linear model is typically fitted for each gene based on the design of a microarray experiment. Section 2 of Smyth [2] provides a nice description of this topic. Given the linear model, suppose that we have an interesting contrast to test for each gene. This contrast may be the difference between the means of two treatment groups or linear combination of means from several treatment groups. For the simplicity of description, we call the genes whose contrast means are not zero as the differentially expressed (DE) genes and the genes whose contrast means equal to zero as equivalently expressed (EE) genes. After fitting the linear model for each gene, we obtain an estimate for the contrast for each gene, . In addition, we get the estimate of the sample residual variance, , for each gene. For each , and are related to true parameters, and , by and , where is the contrast mean for gene , is the true residual variance for gene , and the coefficients and are determined by the design of the experiment. Two examples are given as follows.

Example 2.1. Two-channel microarray experiment to compare two treatments. Assume that each sample from treatment A is paired randomly with a sample from treatment B and each pair of samples is cohybridized onto one slide. After normalization and appropriate transformation, the difference of normalized expression measurements between the two samples on each slide is analyzed for each gene. Hence, this is a paired sample case and the number of data points for each gene is , the number of slides. We are interested in identifying DE genes. In this case, is the mean difference of the paired samples for gene . is the sample variance for gene . So and .

Example 2.2. Affymetrix microarray experiment with two independent samples. Assume sample sizes are and for treatment A and treatment B, respectively. The statistic is the difference in sample means of normalized expression measurements between two groups for gene . is the pooled sample variance. Then and .
Given the data and , an ordinary t-test with statistic may be used to test the null hypothesis : . However, the power of such tests is low after controlling multiple testing error. So statistical methods with higher power are in demand for such high-dimensional testing problem as encountered in gene expression studies.

3. Optimal Testing Procedures

In the analysis of high-dimensional gene expression data such as microarray data, we are more interested in the average behavior of the tests across all genes rather than the performance of an individual test. Because the dimension of tests is huge, multiple testing errors should be controlled to avoid too many type I errors. Controlling FDR is an important method for controlling multiple testing errors and is widely used for genomic studies. Although many testing procedures have been developed as reviewed in Section 1, the paper by Hwang and Liu [4] provides some theoretical guide on how to derive optimal testing procedures within an empirical Bayes framework. The optimal tests are defined to be the ones that maximize the power averaged across all genes for which the null hypotheses are false while controlling FDR. Such optimal tests have been called MAP tests, where MAP stands for maximum average power [4].

In a Bayesian framework, we assume model parameters like and follow some distributions. The residual variances of genes, , have been modeled by prior distribution like inverse gamma [2, 10] or log-normal [4] distribution independent of whether the null hypothesis is true or false. For EE genes, the mean of contrast , , is equal to 0. For DE genes, the mean is not 0 almost surely. Denote the alternative distribution of by . Based on the Neyman-Pearson fundamental lemma, for a randomly selected gene , the most powerful test statistic for testing versus is given by where denote the prior distributions of . And the test rejects the null hypothesis when is large. The simultaneous testing procedure where all genes are tested using the most powerful statistics , achieves the highest average power while controlling FDR, as proved in Hwang and Liu [4].

One popular multiple-testing method for microarray data is the moderated  -test proposed by Smyth [2]. Smyth proposed to model the residual variance with the prior distribution: where denotes a chi-square distribution with degrees of freedom and is another hyperparameter. This prior distribution is equivalent to an inverse-gamma distribution and has been shown to fit real data well. Compared to a standard t-test statistic , Smyth's moderated t-statistic takes the form of where is a shrinkage estimator of by shrinking toward .

In practice, the unknown hyperparameters and for the distribution of the variance can be estimated consistently by the method of moments, that is, equating the empirical and expected first two moments of [2]. Smyth [2] showed that the moderated -test is equivalent to the statistic proposed in Lönnstedt and Speed [1] which was derived as the posterior odds under the assumption that the distribution of under the alternative hypothesis follows , where is a constant. In fact, we can prove the claim that the moderated -test achieves the optimal average power asymptotically under their assumptions for and . The proof is in the appendix.

However, note that the assumption that of DE genes follows a normal distribution with mean zero is restrictive. It is likely that, for example, there are more upregulated genes than downregulated genes for some studies which suggests that the mean of should be positive. Hwang and Liu [4] have proposed a more general normal prior distribution of for DE genes: where the mean of this distribution is not necessarily zero but to be estimated based on data. In addition, the variance for this distribution does not depend on the residual variance. Under this model, they have derived an optimal test and an approximated version of the test statistic ( test) that is computationally faster. The statistic shrinks both the estimate of mean and the estimate of variance .

Both the moderated -test and the test have been shown to achieve optimal power asymptotically under the assumption of normal distribution for the alternative means. Simulation studies also confirm that the power of the tests is superior under the model assumptions. However, a single normal distribution assumption on for DE genes may not be appropriate for all cases and the distribution of may consist of a mixture of different subgroup distributions, for example, a mixture of two normal distributions with one having a negative mean and the other having a positive mean. If the parametric distributional assumptions of are violated, the power of an optimal test built under those assumptions will suffer.

4. Semiparametric Optimal Testing (SPOT) Procedure

To obtain a more robust procedure, we propose to model the distribution of the mean nonparametrically while still deriving the optimal procedure. For the variance , the inverse gamma distributional assumption is reasonable and works well in practice, so we still keep this assumption. Hence, we will derive a semiparametric optimal testing procedure that we call the SPOT procedure.

Note that the numerator and denominator of the most powerful test statistic (3.1) are the joint marginal distributions of , under the alternative and null hypothesis, respectively. By denoting the marginal distributions by statistic (3.1) becomes The null marginal distribution only involves integration with respect to variance . With consistent estimators of hyperparameters as proposed in Smyth [2], we can estimate consistently. For the alternative marginal distribution , it is hard to find a consistent estimator without any distributional assumption on . If we were to know which genes are DE, then we could estimate nonparametrically with observed values of () from the DE gene population. Many nonparametric density estimators are consistent, for example, the histogram estimators and the kernel density estimators with proper choices of bandwidths [13]. However, the knowledge of differential expression is the research question of the study and of course is not available for all genes. Considering all genes without separating those that are differentially expressed from those that are not, we have a mixture distribution of differentially expressed and nondifferentially expressed genes. The mixture density of the marginal distributions, denoted by , can be estimated consistently by nonparametric density estimators with observed () for all genes . Can this consistent estimator of help us construct a most powerful test statistic, together with a consistent estimator of ?

Suppose that and are proportions of EE and DE genes, respectively, with and , then the mixture marginal density is

The ratio of mixture marginal density and the null marginal density is a monotonic function of the statistic expressed in formula (4.2) because Thus the test that rejects the null hypothesis when is large is also a most powerful test. Note that to calculate this statistic, we only need to estimate and but do not have to estimate the proportions and .

Let denote any consistent density estimator of , and let denote any consistent estimator of , such that where denotes convergence in probability. Then the statistic has the optimal testing power asymptotically. Notice the convergence with respect to , which is usually huge in the microarray and RNA-seq studies.

We have already discussed the availability of a parametric consistent estimator of through estimating the hyperparameters and of in Section 3. For , any theoretically consistent density estimator of joint data () can be used to construct the test statistic (4.4) with asymptotically optimal average power. For example, nonparametric estimators such as histograms, kernel density estimates, and local polynomial estimators can all be utilized. As our test statistic involves both parametric and nonparametric parts, we name it the semiparametric optimal test (SPOT).

5. Implementation of SPOT

In this section, we discuss details in implementation of the proposed SPOT procedure.

5.1. Estimation of

The null marginal density where is a constant. As in Smyth [2] and Hwang and Liu [4], we assume that the distribution of does not depend on whether a gene is DE or EE. Then, all genes are used to estimate the parameters and . We apply the method of moments proposed in Smyth [2] to get estimates of and . Replacing unknown parameters and in by their consistent method of moments estimates leads to a consistent estimator of .

5.2. A Hybrid Method for Estimation of

Although any consistent estimator can be used to construct a SPOT statistic of the form , in practice, a density estimator that converges fast would be always preferred. It is known that the accuracy of the density estimators goes down quickly as the dimension increases [13]. We have tried a few two-dimensional density estimators for , including the kernel estimators. Due to the curse of dimensionality, the direct two-dimensional density estimators do not perform as satisfactory as a hybrid estimator that we develop and would suggest to use. This hybrid estimator has a component that is similar to kernel estimators, whereas it also utilizes the prior information on variances to help improving the accuracy.

In constructing this estimator, we first estimate the marginal density of by the typical kernel density estimate: where is a positive value known as bandwidth. We estimate the conditional density of by using where the second equality is a result of the independence between and given the parameter . The distribution of is for normal-distributed observations. Now we need to estimate . Denote the set of genes that lie within bandwidth distance to gene as if and only if . We estimate by the following approximation that is based on the neighborhood of , :

The in formula denotes the number of genes in set . Substituting exact parametric form of and into above formulas leads to the explicit form where is a constant. The product between the kernel estimate and the conditional estimate provides us a joint density estimator of the mixture With this approximation, we cannot theoretically show that the resulting estimator, , is consistent but it works better in practice than the consistent kernel density estimator of the joint density .

6. Simulation Study

In order to evaluate the performance of our proposed SPOT procedure, we performed three simulation studies. The gene expression data were simulated from Normal  (,) for observations of gene in treatment group . The way to sample and differs across simulation studies. We assume that there are two treatment groups and 3 replicates per treatment group. For each simulation setting, one hundred sets of gene expression data were independently simulated, and each dataset included 10,000 genes. The performances of the SPOT, moderated , , and ordinary -test statistics were evaluated for by comparing their average behavior averaged across the 100 datasets.

6.1. Simulation Study I

In the first simulation study, we have two settings that differ in the number of DE genes. For the first setting, are DE genes whereas the other are EE. In the second setting, only are DE while the other are EE. Gene expression means and variances were simulated as follows:

For each simulated data, SPOT, moderated , , and ordinary -test statistics were calculated and evaluated using the number of selected true positives at various FDR levels. The plots of number of true positives versus FDR for SPOT, moderated , , and ordinary -test statistics are shown in Figure 1. Simulation settings 1 and 2 generated similar results. The ordinary -test is the poorest method under comparison. The moderated -test is considerably better than the ordinary -test although it is worse than test. Our proposed SPOT test is superior to all other three methods, with the largest number of true positive findings than the other three statistics at the same FDR levels.

Figure 1: Simulation study I: plots of number of true positives (# TP) versus false discovery rate (FDR) from analyses using SPOT, moderated t, and methods.

6.2. Simulation Study II

To check how the variance distribution affects the relative ranking of the SPOT procedure, we did another simulation study the same as the setting 1 of simulation study I except that the variances were simulated from a log-normal distribution, which is the assumption under which the test was derived. As Figure 2 shows, the results are similar to those from simulation I. The SPOT procedure still performs much better than all the other three methods.

Figure 2: Simulation study II: plots of number of true positives (# TP) versus false discovery rate (FDR) from analyses using SPOT, moderated t, and methods.

6.3. Simulation Study III

Typically, the parametric test achieves higher power than the nonparametric test if the parametric assumption is appropriate. To check the robustness of the SPOT procedure, we simulated data under the parametric assumption for both and under which the test was derived. Specifically, for the 2,500 differentially expressed genes, were drawn from a normal distribution with mean 1.2 and standard deviation 0.3, were sampled from a log-normal distribution with parameters −0.96 and 0.8. Figure 3 shows that the SPOT procedure and the test are comparable to each other when FDR is small (less than 0.05) and they are both much better than the moderated -test and the ordinary -test. When FDR is between 0.05 and 0.15, the test is the best while the SPOT procedure is the next best performing procedure, which is still much better than the moderated -test and the ordinary -test.

Figure 3: Simulation study III: plots of number of true positives (# TP) versus false discovery rate (FDR) from analyses using SPOT, moderated t, and methods.

7. Evaluation Using the Golden Spike Microarray Data

In this section, we compare the performances of different methods using a real microarray dataset from experiments conducted using Affymetrix GeneChip in the Golden Spike Project. The Golden Spike Project generated microarray datasets comparing two replicated groups in which the relative concentrations of a large number of genes are known. The two groups are the spike-in group and the control group, each with three chips. Data and information related to this project are available through the website http://www2.ccr.buffalo.edu/halfon/spike/. More specifically, the Golden Spike dataset included 1309 individual cRNAs “spiked in” at known relative concentrations between the two groups. The fold-changes between the spike-in and control group were assigned at different levels for different cRNAs, and the levels ranged from 1.2 to 4. Hence, these cRNAs were truly “differentially expressed” between groups and we consider them as DE genes. In addition, a background sample of 2551 RNA species was present at identical concentrations in both samples. So these 2551 RNA species were not differentially expressed between the two groups. With the knowledge of the true differential expression status, this real microarray dataset provides an ideal case to evaluate the performances of different methods without imposing any distributional assumption for variances and means as usually is done in simulation studies.

With the summary dataset downloaded from the Golden Spike Project website, we calculated the SPOT, the moderated , the ordinary , and the statistics and evaluated their performances using the true statuses of RNA based on the design. Figure 4 shows the plots of number of true positives versus FDR for the ordinary , moderated , , and SPOT procedures over a range of FDR which is of most practical interest. It can be observed that the performance of the SPOT procedure improves over the performances of the other three methods throughout the whole range of FDR in these plots. In addition, the improvement is substantial at lower FDR levels. For example, the SPOT procedure detects 754 true positives at the FDR level of 1% while the moderated -test only detects 466 and the test only detects 442 true positives (Table 1).

Table 1: Golden Spike data: number of true positives selected by three testing procedures at critical FDR levels.

Figure 4: Golden spike data: plots of number of true positive (TP) genes versus false discovery rate (FDR) from analysis using SPOT, moderated ,  and tests.

8. Summary

In this paper, we have derived a semiparametric optimal testing (SPOT) procedure for high-dimensional gene expression data analysis. Although the method is illustrated for analyzing microarray data, it can be applied to any high-dimensional testing problem with normal model. Our test statistic is justified to be asymptotically most powerful, without any assumption on the mean parameter of differential expression. The asymptotic property is derived when the number of genes is large, which is reasonable for high-dimensional gene expression studies. We also provided an approximate version to implement the SPOT procedure in practice and evaluated the performance of our proposed test statistic using both simulation studies and real microarray data analysis. The proposed SPOT method is shown to outperform the popularly applied moderated   and the statistics, which are optimal only under certain normality conditions of the mean. There is still potential in improving the performance of SPOT procedure if better density estimates can be found for the marginal distributions and .

Appendix

Proof of the Claim That the Moderated t-Test Achieves the Optimal Average Power Asymptotically under the Assumptions That and

Under Smyth's [2] model assumptions, the most power test statistic formula (3.1) derived under the Neyman-Pearson lemma becomes where  denotes ,  and   denotes ,  which is a monotonic function of Smyth's [2] moderated t-statistic, with being some constant. Thus the claim follows with existence of consistent estimates of and , which has been shown in Smyth [2].