In high-dimensional gene expression experiments such as microarray and RNA-seq experiments, the number of measured variables is huge while the number of replicates is small. As a consequence, hypothesis testing is challenging because the power of tests can be very low after controlling multiple testing error. Optimal testing procedures with high average power while controlling false discovery rate are preferred. Many methods were constructed to achieve high power through borrowing information across genes. Some of these methods can be shown to achieve the optimal average power across genes, but only under a normal assumption of alternative means. However, the assumption of a normal distribution is likely violated in practice. In this paper, we propose a novel semiparametric optimal testing (SPOT) procedure for high-dimensional data with small sample size. Our procedure is more robust because it does not depend on any parametric assumption for the alternative means. We show that the proposed test achieves the maximum average power asymptotically as the number of tests goes to infinity. Both simulation study and the analysis of a real microarray data with spike-in probes show that the proposed SPOT procedure performs better when compared to other popularly applied procedures.

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 š¹SS 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 š¹SS 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 š¹SS test. Both the moderated š‘”-statistic and the š¹SS 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 š¹SS 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, š‘ 2š‘”, for each gene. For each š‘”=1,ā€¦,šŗ, š‘‹š‘” and š‘ 2š‘” are related to true parameters, šœ‡š‘” and šœŽ2š‘”, by š‘‹š‘”āˆ£šœ‡š‘”,šœŽ2š‘”āˆ¼š‘(šœ‡š‘”,šœˆš‘”šœŽ2š‘”) and š‘ 2š‘”āˆ£šœŽ2š‘”āˆ¼(šœŽ2š‘”/š‘‘š‘”)šœ’2š‘‘š‘”, where šœ‡š‘” is the contrast mean for gene š‘”, šœŽ2š‘” 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 š‘”. š‘ 2š‘” is the sample variance for gene š‘”. So šœˆš‘”=1/š‘› and š‘‘š‘”=š‘›āˆ’1.

Example 2.2. Affymetrix microarray experiment with two independent samples. Assume sample sizes are š‘›1 and š‘›2 for treatment A and treatment B, respectively. The statistic š‘‹š‘” is the difference in sample means of normalized expression measurements between two groups for gene š‘”. š‘ 2š‘” is the pooled sample variance. Then šœˆš‘”=1/š‘›1+1/š‘›2 and š‘‘š‘”=š‘›1+š‘›2āˆ’2.
Given the data š‘‹š‘” and š‘ 2š‘”, an ordinary t-test with statistic š‘”š‘”=š‘‹š‘”/āˆššœˆš‘”š‘ š‘” may be used to test the null hypothesis š»0š‘”: šœ‡š‘”=0. 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 šœŽ2š‘” follow some distributions. The residual variances of genes, šœŽ2š‘”, 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 šœ‹1(ā‹…). Based on the Neyman-Pearson fundamental lemma, for a randomly selected gene š‘”, the most powerful test statistic for testing š»0š‘”āˆ¶šœ‡š‘”=0 versus š»1š‘”āˆ¶šœ‡š‘”āˆ¼šœ‹1(šœ‡š‘”) is given by š‘‡NPš‘”=āˆ¬š‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”,šœŽ2š‘”ī€øšœ‹1ī€·šœ‡š‘”ī€øšœ‹ī€·šœŽ2š‘”ī€øš‘‘šœ‡š‘”š‘‘šœŽ2š‘”āˆ«š‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”=0,šœŽ2š‘”ī€øšœ‹ī€·šœŽ2š‘”ī€øš‘‘šœŽ2š‘”,(3.1) where šœ‹(ā‹…) denote the prior distributions of šœŽ2š‘”. And the test rejects the null hypothesis š»0š‘” when š‘‡NPš‘” is large. The simultaneous testing procedure where all genes are tested using the most powerful statistics š‘‡NPš‘”,š‘”=1,2,ā€¦,šŗ, 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 šœŽ2š‘” with the prior distribution: 1šœŽ2š‘”āˆ¼1š‘‘0š‘ 20šœ’2š‘‘0,(3.2) where šœ’2š‘‘0 denotes a chi-square distribution with degrees of freedom š‘‘0 and š‘ 20 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 Ģƒš‘”š‘”=š‘„š‘”āˆššœˆš‘”Ģƒš‘ š‘”,(3.3) where Ģƒš‘ 2š‘”=š‘ 2š‘”š‘‘š‘”+š‘ 20š‘‘0š‘‘š‘”+š‘‘0(3.4) is a shrinkage estimator of šœŽ2š‘” by shrinking š‘ 2š‘” toward š‘ 20.

In practice, the unknown hyperparameters š‘‘0 and š‘ 20 for the distribution of the variance šœŽ2š‘” can be estimated consistently by the method of moments, that is, equating the empirical and expected first two moments of logš‘ 2š‘” [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 š‘(0,šœˆ0šœŽ2š‘”), where šœˆ0 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 šœŽ2š‘”. 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: šœ‹1ī€·šœ‡š‘”ī€øī€·āˆ¼š‘šœƒ,šœ2š‘”ī€ø,(3.5) 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 (š¹SS test) that is computationally faster. The š¹SS statistic shrinks both the estimate of mean šœ‡š‘” and the estimate of variance šœŽ2š‘”.

Both the moderated š‘”-test and the š¹SS 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 šœ‹1(šœ‡š‘”) 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 šœ‹1(šœ‡š‘”) 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 šœŽ2š‘”, 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 (š‘‹š‘”,š‘ 2š‘”), under the alternative and null hypothesis, respectively. By denoting the marginal distributions by š‘š1ī€·š‘‹š‘”,š‘ 2š‘”ī€ø=ī€š‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”,šœŽ2š‘”ī€øšœ‹1ī€·šœ‡š‘”ī€øšœ‹ī€·šœŽ2š‘”ī€øš‘‘šœ‡š‘”š‘‘šœŽ2š‘”,š‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø=ī€œš‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”=0,šœŽ2š‘”ī€øšœ‹ī€·šœŽ2š‘”ī€øš‘‘šœŽ2š‘”,(4.1) statistic (3.1) becomes š‘‡NPš‘”=š‘š1ī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø.(4.2) The null marginal distribution š‘š0(š‘‹š‘”,š‘ 2š‘”) only involves integration with respect to variance šœŽ2š‘”. With consistent estimators of hyperparameters as proposed in Smyth [2], we can estimate š‘š0(š‘‹š‘”,š‘ 2š‘”) consistently. For the alternative marginal distribution š‘š1(š‘‹š‘”,š‘ 2š‘”), 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 š‘š1(š‘‹š‘”,š‘ 2š‘”) nonparametrically with observed values of (š‘‹š‘”,š‘ 2š‘”) 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 š‘šš‘š(š‘‹š‘”,š‘ 2š‘”), can be estimated consistently by nonparametric density estimators with observed (š‘‹š‘”,š‘ 2š‘”) for all genes š‘”=1,ā€¦,šŗ. Can this consistent estimator of š‘šš‘š(š‘‹š‘”,š‘ 2š‘”) help us construct a most powerful test statistic, together with a consistent estimator of š‘š0(š‘‹š‘”,š‘ 2š‘”)?

Suppose that š‘0 and š‘1 are proportions of EE and DE genes, respectively, with 0ā‰¤š‘0,š‘1ā‰¤1 and š‘0+š‘1=1, then the mixture marginal density is š‘šš‘šī€·š‘‹š‘”,š‘ 2š‘”ī€ø=š‘0š‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø+š‘1š‘š1ī€·š‘‹š‘”,š‘ 2š‘”ī€ø.(4.3)

The ratio of mixture marginal density š‘šš‘š(š‘‹š‘”,š‘ 2š‘”) and the null marginal density š‘š0(š‘‹š‘”,š‘ 2š‘”) is a monotonic function of the statistic š‘‡NPš‘” expressed in formula (4.2) because š‘šš‘šī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø=š‘0š‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø+š‘1š‘š1ī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø,=š‘0+š‘1š‘š1ī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø.(4.4) Thus the test that rejects the null hypothesis when š‘šš‘š(š‘‹š‘”,š‘ 2š‘”)/š‘š0(š‘‹š‘”,š‘ 2š‘”) is large is also a most powerful test. Note that to calculate this statistic, we only need to estimate š‘šš‘š(š‘‹š‘”,š‘ 2š‘”) and š‘š0(š‘‹š‘”,š‘ 2š‘”) but do not have to estimate the proportions š‘0 and š‘1.

Let īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”) denote any consistent density estimator of š‘šš‘š(š‘‹š‘”,š‘ 2š‘”), and let īš‘š0(š‘‹š‘”,š‘ 2š‘”) denote any consistent estimator of š‘š0(š‘‹š‘”,š‘ 2š‘”), such that īš‘šš‘šī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘ƒāŸ¶š‘šš‘šī€·š‘‹š‘”,š‘ 2š‘”ī€øasšŗā†—āˆž,īš‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€øš‘ƒāŸ¶š‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€øasšŗā†—āˆž,(4.5) where š‘ƒāˆ’ā†’ denotes convergence in probability. Then the statistic īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”)/īš‘š0(š‘‹š‘”,š‘ 2š‘”) 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 š‘š0(š‘‹š‘”,š‘ 2š‘”) through estimating the hyperparameters š‘‘0 and š‘ 20 of šœŽ2š‘” in Section 3. For š‘šš‘š(š‘‹š‘”,š‘ 2š‘”), any theoretically consistent density estimator īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”) of joint data (š‘‹š‘”,š‘ 2š‘”) 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 īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”)/īš‘š0(š‘‹š‘”,š‘ 2š‘”) 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 š‘š0(š‘‹š‘”,š‘ 2š‘”)

The null marginal density š‘š0ī€·š‘‹š‘”,š‘ 2š‘”ī€ø=ī€œš‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”=0,šœŽ2š‘”ī€øšœ‹ī€·šœŽ2š‘”ī€øš‘‘šœŽ2š‘”=ī€œš‘’āˆ’š‘„2š‘”/(2š‘£š‘”šœŽ2š‘”)ī€·2šœ‹š‘£š‘”šœŽ2š‘”ī€ø1/2īƒ©š‘‘š‘”2šœŽ2š‘”īƒŖš‘‘š‘”/2š‘ 2(š‘‘š‘”/2āˆ’1)š‘’āˆ’š‘‘š‘”š‘ 2š‘”/(2šœŽ2š‘”)Ī“ī€·š‘‘š‘”ī€øīƒ©š‘‘/20š‘ 202īƒŖš‘‘0/2šœŽāˆ’2(š‘‘0š‘”/2+1)š‘’āˆ’š‘‘0š‘ 20/2šœŽ2š‘”Ī“ī€·š‘‘0ī€ø/2š‘‘šœŽ2š‘”=š¶2ā‹…š‘ š‘”2(š‘‘/2āˆ’1)īƒ©š‘„2š‘”/š‘£š‘”+š‘‘0š‘ 20+š‘‘š‘”š‘ 2š‘”2īƒŖāˆ’(1+š‘‘0+š‘‘š‘”)/2,(5.1) where š¶2 is a constant. As in Smyth [2] and Hwang and Liu [4], we assume that the distribution of šœŽ2š‘” does not depend on whether a gene is DE or EE. Then, all genes are used to estimate the parameters š‘‘0 and š‘ 20. We apply the method of moments proposed in Smyth [2] to get estimates of š‘‘0 and š‘ 20. Replacing unknown parameters š‘‘0 and š‘ 20 in š‘š0(š‘‹š‘”,š‘ 2š‘”) by their consistent method of moments estimates leads to a consistent estimator īš‘š0(š‘‹š‘”,š‘ 2š‘”) of š‘š0(š‘‹š‘”,š‘ 2š‘”).

5.2. A Hybrid Method for Estimation of š‘šš‘š(š‘‹š‘”,š‘ 2š‘”)

Although any consistent estimator īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”) can be used to construct a SPOT statistic of the form īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”)/īš‘š0(š‘‹š‘”,š‘ 2š‘”), 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 īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”), 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 šœŽ2š‘” to help improving the accuracy.

In constructing this estimator, we first estimate the marginal density of š‘‹š‘” by the typical kernel density estimate: īš‘“ī€·š‘„š‘”ī€ø=1šŗšŗī“š‘–=11ā„Žš¾ī‚µš‘„š‘”āˆ’š‘„š‘–ā„Žī‚¶,(5.2) where ā„Ž is a positive value known as bandwidth. We estimate the conditional density of š‘“(š‘ 2š‘”āˆ£š‘„š‘”) by using š‘“ī€·š‘ 2š‘”āˆ£š‘„š‘”ī€ø=ī€œš‘“ī€·š‘ 2š‘”āˆ£šœŽ2š‘”,š‘„š‘”ī€øš‘“ī€·šœŽ2š‘”āˆ£š‘„š‘”ī€øš‘‘šœŽ2š‘”=ī€œš‘“ī€·š‘ 2š‘”āˆ£šœŽ2š‘”ī€øš‘“ī€·šœŽ2š‘”āˆ£š‘„š‘”ī€øš‘‘šœŽ2š‘”,(5.3) where the second equality is a result of the independence between š‘ 2š‘” and š‘„š‘” given the parameter šœŽ2š‘”. The distribution of š‘ 2š‘”āˆ£šœŽ2š‘” is (šœŽ2š‘”/š‘‘š‘”)šœ’2š‘‘š‘” for normal-distributed observations. Now we need to estimate š‘“(šœŽ2š‘”āˆ£š‘„š‘”). Denote the set of genes that lie within bandwidth distance to gene š‘” as {š“š‘”āˆ¶š‘–āˆˆš“š‘” if and only if |š‘„š‘–āˆ’š‘„š‘”|<ā„Ž}. We estimate š‘“(šœŽ2š‘”āˆ£š‘„š‘”) by the following approximation that is based on the neighborhood of š‘„š‘”, š“š‘”: 1#ī€½š“š‘”ī€¾ī“š‘–āˆˆš“š‘”š‘“ī€·š‘ 2š‘–āˆ£šœŽ2ī€øšœ‹ī€·šœŽ2ī€øāˆ«š‘“ī€·š‘ 2š‘–āˆ£šœŽ2ī€øšœ‹ī€·šœŽ2ī€øš‘‘šœŽ2.(5.4)

The #{š“š‘”} in formula denotes the number of genes in set š“š‘”. Substituting exact parametric form of š‘“(š‘ 2š‘”āˆ£šœŽ2) and šœ‹(šœŽ2) into above formulas leads to the explicit form īš‘“ī€·š‘ 2š‘”āˆ£š‘„š‘”ī€ø=š¶3ā‹…1#ī€½š“š‘”ī€¾ī“š‘–āˆˆš“š‘”š‘ š‘”2(š‘‘/2āˆ’1)īƒ©š‘‘š‘”š‘ 2š‘–+š‘‘0š‘ 20+š‘‘š‘”š‘ 2š‘”2īƒŖāˆ’(š‘‘0+š‘‘š‘”)/2,(5.5) where š¶3 is a constant. The product between the kernel estimate īš‘“(š‘„š‘”) and the conditional estimate īš‘“(š‘ 2š‘”āˆ£š‘„š‘”) provides us a joint density estimator of the mixture īš‘šš‘šī€·š‘‹š‘”,š‘ 2š‘”ī€ø=īš‘“ī€·š‘„š‘”ī€øā‹…īš‘“ī€·š‘ 2š‘”āˆ£š‘„š‘”ī€ø.(5.6) With this approximation, we cannot theoretically show that the resulting estimator, īš‘šš‘š(š‘‹š‘”,š‘ 2š‘”), is consistent but it works better in practice than the consistent kernel density estimator of the joint density š‘šš‘š(š‘‹š‘”,š‘ 2š‘”).

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ā€‰ā€‰(šœ‡š‘”š‘–,šœŽ2š‘”) for observations of gene š‘” in treatment group š‘–. The way to sample šœ‡š‘”š‘– and šœŽ2š‘” 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 š‘”, š¹SS, 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, šŗ1=2,500 are DE genes whereas the other šŗ0=7,500 are EE. In the second setting, only šŗ1=1,800 are DE while the other šŗ0=8,200 are EE. Gene expression means šœ‡š‘”š‘– and variances šœŽ2š‘” were simulated as follows: šœ‡š‘”1šœ‡=0āˆ€š‘”;š‘”2ī€·āˆ¼Normal0.5,0.32ī€øforš‘”=1to0.3šŗ1;šœ‡š‘”2ī€·āˆ¼Normal1,0.32ī€øī€·forš‘”=0.3šŗ1ī€ø+1to0.9šŗ1;šœ‡š‘”2āˆ¼š‘”1ī€·(0.5)forš‘”=0.9šŗ1ī€ø+1tošŗ1;šœ‡š‘”2ī€·šŗ=0forš‘”=1ī€øšœŽ+1to10000;2š‘”āˆ¼Gamma(2,4)āˆ€š‘”.(6.1)

For each simulated data, SPOT, moderated š‘”, š¹SS, 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 š‘”, š¹SS, 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 š¹SS 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.

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 š¹SS 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.

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 šœŽ2š‘” under which the š¹SS 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, šœŽ2š‘” were sampled from a log-normal distribution with parameters āˆ’0.96 and 0.8. Figure 3 shows that the SPOT procedure and the š¹SS 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 š¹SS 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.

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 š¹SS 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 š‘”, š¹SS, and SPOT procedures over a range of FDR āˆˆ[0,0.15] 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 š¹SS test only detects 442 true positives (Table 1).

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 š¹SS 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 š‘šš‘š(š‘‹š‘”,š‘ 2š‘”) and š‘š0(š‘‹š‘”,š‘ 2š‘”).


Proof of the Claim That the Moderated t-Test Achieves the Optimal Average Power Asymptotically under the Assumptions That šœ‡š‘”āˆ¼š‘(0,šœˆ0šœŽ2š‘”) and 1/šœŽ2š‘”āˆ¼1/š‘‘0š‘ 20šœ’2š‘‘0

Under Smyth's [2] model assumptions, the most power test statistic formula (3.1) derived under the Neyman-Pearson lemma becomes š‘‡NPš‘”=āˆ¬š‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”,šœŽ2š‘”ī€øšœ‹1ī€·šœ‡š‘”ī€øšœ‹1ī€·šœŽ2š‘”ī€øš‘‘šœ‡š‘”š‘‘šœŽ2š‘”āˆ«š‘“ī€·š‘‹š‘”,š‘ 2š‘”āˆ£šœ‡š‘”=0,šœŽ2š‘”ī€øšœ‹0ī€·šœŽ2š‘”ī€øš‘‘šœŽ2š‘”=āˆ¬ī‚€š‘’āˆ’(š‘„š‘”āˆ’šœ‡š‘”)2/2š‘£š‘”šœŽ2š‘”/ī€·2šœ‹š‘£š‘”šœŽ2š‘”ī€ø1/2ī‚ī€·š‘‘š‘”/2šœŽ2š‘”ī€øš‘‘š‘”/2ī‚€š‘ š‘‘š‘”āˆ’2š‘’āˆ’š‘‘š‘”š‘ 2š‘”/2šœŽ2š‘”ī€·š‘‘/Ī“š‘”ī€øī‚š’œ/2āˆ«ī‚€š‘’āˆ’š‘„2š‘”/2š‘£š‘”šœŽ2š‘”/ī€·2šœ‹š‘£š‘”šœŽ2š‘”ī€ø1/2ī‚ī€·š‘‘š‘”/2šœŽ2š‘”ī€øš‘‘š‘”/2ī€·š‘ 2(š‘‘š‘”/2āˆ’1)ī€·š‘‘/Ī“š‘”š‘’/2ī€øī€øāˆ’š‘‘š‘”š‘ 2š‘”/2šœŽ2š‘”ā‹…ā„¬=š¶1ā‹…īƒ©š‘„2š‘”/ī€·š‘£0+š‘£š‘”ī€ø+š‘‘0š‘ 20+š‘‘š‘”š‘ 2š‘”š‘„2š‘”/š‘£š‘”+š‘‘0š‘ 20+š‘‘š‘”š‘ 2š‘”īƒŖāˆ’(1+š‘‘0+š‘‘š‘”)/2,(A.1)whereā€‰š’œā€‰denotesā€‰(š‘’āˆ’šœ‡2š‘”/(2š‘£0šœŽ2š‘”)/(2šœ‹š‘£0šœŽ2š‘”)1/2)(š‘‘0š‘ 20/2)š‘‘0/2(šœŽāˆ’š‘‘0š‘”+2/Ī“(š‘‘0/2))š‘’āˆ’š‘‘0š‘ 20/(2šœŽ2š‘”)š‘‘šœ‡š‘”š‘‘šœŽ2š‘”,ā€‰ andā€‰ā€‰ā„¬ā€‰denotes (š‘‘0š‘ 20/2)š‘‘0/2(šœŽāˆ’2(š‘‘0š‘”/2āˆ’1)/Ī“(š‘‘0/2))š‘’āˆ’š‘‘0š‘ 20/(2šœŽ2š‘”)š‘‘šœŽ2š‘”,ā€‰ which is a monotonic function of Smyth's [2] moderated t-statistic, with š¶1 being some constant. Thus the claim follows with existence of consistent estimates of š‘‘0 and š‘ 20, which has been shown in Smyth [2].