Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2895927

Formats

Article sections

- Summary
- 1. Introduction
- 2. Notations, Definitions and Problem Formulation
- 3. Multidimensional Directional FDR Controlling Procedures
- 4. A Simulation Study
- 5. An Application to Time-Course Gene Expression Data
- 6. Concluding Remarks
- Supplementary Material
- References

Authors

Related links

Biometrics. Author manuscript; available in PMC 2011 June 1.

Published in final edited form as:

Published online 2009 July 23. doi: 10.1111/j.1541-0420.2009.01292.x

PMCID: PMC2895927

NIHMSID: NIHMS138779

Wenge Guo, Biostatistics Branch, National Institute of Environmental Health Sciences Research Triangle Park, NC 27709, U.S.A. Email: moc.liamg@oug.egnew;

The publisher's final edited version of this article is available at Biometrics

See other articles in PMC that cite the published article.

Microarray gene expression studies over ordered categories are routinely conducted to gain insights into biological functions of genes and the underlying biological processes. Some common experiments are time-course/dose-response experiments where a tissue or cell-line is exposed for different doses and/or durations of time to a chemical. A goal of such studies is to identify gene expression patterns/profiles over the ordered categories. This problem can be formulated as a multiple testing problem where for each gene the null hypothesis of no difference between the successive mean gene expressions are tested and further directional decisions are made if it is rejected. Much of the existing multiple testing procedures are devised for controlling the usual false discovery rate (FDR) rather than the mixed directional FDR, the expected proportion of Type I and directional errors among all rejections. Benjamini and Yekutieli (2005) proved that an augmentation of the usual Benjamini-Hochberg (BH) procedure can control the mixed directional FDR while testing simple null hypotheses against two-sided alternatives in terms of one dimensional parameters. In this article, we consider the problem of controlling the mixed directional FDR involving multidimensional parameters. To deal with this problem, we develop a procedure extending that of Benjamini and Yekutieli based on the Bonferroni test for each gene. A proof is given for its mixed directional FDR control when the underlying test statistics are independent across the genes. The results of a simulation study evaluating its performance under independence as well as under dependence of the underlying test statistics across the genes relative to other relevant procedures are reported. Finally, the proposed methodology is applied to a time-course microarray data obtained by Lobenhofer et al. (2002). We identified several important cell-cycle genes, such as DNA replication/repair gene MCM4 and replication factor subunit C2, which were not identified by the previous analyses of the same data by Lobenhofer et al. (2002) and Peddada et al. (2003). Although some of our findings overlap with previous findings, we identify several other genes that compliment the results of Lobenhofer et al. (2002).

In many applications researchers are interested in identifying trends in mean response over ordered categories in large scale experiments. With the advent of microarray technology such experiments are common in the literature where investigators are routinely conducting experiments to investigate changes in mean gene expressions over time or dose of a chemical or cancer stage, etc. For example, Lobenhofer et al. (2002) studied the effect of 17−*β* estradiol on the gene expression of MCF-7 breast cancer cells as the cells progressed through various phases of cell division cycle. In another experiment, Tamoto et al. (2004) investigated the changes in gene expression with tumor progression in esophagal cancer and identified genes implicated in the early stages of esophagal squamous cell carcinoma. Recently, Bochkina and Richardson (2007) discussed the analysis of a time-course gene expression data where cells from the H2Kb muscle cell line of mouse were treated by insulin (0, 2 or 12 hours of exposure).

In studies such as those described above, identification of statistically significant genes that have similar mean expression profiles over ordered categories is often an important goal to researchers. By identifying such genes, the researchers may potentially discover co-regulated genes belonging to similar pathways and gain insights into biological functions and processes of groups of genes with similar patterns of expressions.

Peddada et al. (2003) introduced an order restricted inference based method for identifying significant genes and group them according to various patterns of inequalities. Implicitly in their methodology, two decisions are being taken for each gene. First, it is decided whether or not a gene is significant using a method exercising a control over gene specific Type I error rate. Then, a suitable inequality pattern is assigned for each selected significant gene based on the values of the underlying test statistics. The directional error that can potentially occur in addition to the usual Type I error due to assigning wrong inequality pattern to a selected significant gene has not been addressed in that paper. In this paper, we take care of both the Type I and directional errors. We do that by taking a multiple testing approach to the main problem where for each gene the mean expressions are successively compared across the ordered categories and a null hypothesis signifying no particular directional pattern is formed to test against the union of all possible directional patterns. We develop a method of simultaneously testing these null hypotheses and determining a directional pattern upon rejection of each of them controlling both the Type I and directional errors in an overall sense.

Although directional error has not been discussed extensively in the literature, it is perhaps a common error that occurs in applications. While testing a null hypothesis *H*_{0} : *θ* = 0 against the two-sided alternative *H*_{1} : *θ* ≠ 0, for some single parameter *θ* of interest, researchers commonly conclude either *θ* > 0 or *θ* < 0 upon rejection of *H*_{0} depending on the sign of the underlying test statistic, keeping the directional error controlled in addition to the Type I error. However, when multiple hypotheses are tested and the number of parameters describing directional patterns is larger, even moderately, than one, as is the case with time course or dose-response microarray data, controlling the directional errors is a problem.

A traditional approach to dealing with directional as well as Type I errors from a multiple testing point of view is to apply a method that controls the so called mixed directional familywise error rate (mdFWER), which is the probability of one or more Type I or directional errors, a variant of the classical familywise error rate (FWER) (Finner, 1994, 1999; Liu, 1996; Sarkar et al., 2004; Shaffer, 1980). However, when the number of null hypotheses is large, as in the context of microarray experiments, the notion of mdFWER, just like the FWER, is too stringent, allowing little chance to make true directional as well as non-directional discoveries. The FDR (False Discovery Rate), due to Benjamini and Hochberg (1995), is a more powerful concept of overall Type I error rate than the FWER in the context of multiple testing and is now most commonly used in large scale scientific investigations, especially in microarray gene expression studies. A variant of it while controlling both Type I and directional errors would be more powerful than the mdFWER. Two such variants have been introduced in the literature (Benjamini et al., 1993), the pure directional FDR which is the expected proportion of directional errors among rejected hypotheses and the mixed directional FDR (mdFDR) which is the expected proportion of Type I and directional errors among rejected hypotheses. In this article, we will focus on procedures controlling the mdFDR.

Benjamini and Yekutieli (2005) gave a method with independent tests that controls the mdFDR when testing multiple simple hypotheses against two-sided alternatives. They proved that the original Benjamini and Hochberg (1995) procedure controlling the FDR at *α* can be augmented to make directional decision upon rejecting a null hypothesis according to the value of the corresponding test statistic without causing the mdFDR to exceed *α*, a result conjectured before by several authors (Benjamini and Hochberg, 2000; Shaffer, 2002; William et al., 1999). Throughout the paper we shall denote Benjamini and Hochberg procedure by BH procedure. Clearly, this method, referred to as the directional BH procedure, can be applied to analyze dose-response microarray data if there are only two ordered categories, but often this is not the case, as such data typically involve more than two ordered categories and the method needs to be suitably extended to accommodate such multiple categories.

We extend the BH directional FDR procedure in this article to develop our proposed multiple testing method that allows us to make a decision on the directional pattern involving multiple parameters once a null hypothesis of no pattern is rejected and maintains a control over the mdFDR. The proposed methodology is then evaluated through a simulation study and applied to the time-course microarray data in Lobenhofer et al. (2002). Our analysis of Lobenhofer's data resulted in the discovery of several cell-cycle genes that were not previously identified by Lobenhoer et al. (2002) and Peddada et al. (2003). Some of our findings complement the previous findings as detailed in Section 5. An important and unique feature of our methodology is that it permits us to specify the time interval of up (or down) regulation of a gene during the 48 hour period of the cell-cycle. One of the usual objectives for conducting cell-cyle time course experiments is to determine the phase of peak expression for a cell-cycle gene and our methodology allows us to make such determinations.

In this section, we present the multiple testing formulation of the problem of identifying expression patterns/trends over ordered categories simultaneously for all the genes, having introduced some notations and definitions related to multiple testing.

Let *μ _{ij}* denote the mean response of the

$${H}_{0j}:{\delta}_{j}=0\phantom{\rule{thinmathspace}{0ex}}\text{against}\phantom{\rule{thinmathspace}{0ex}}{H}_{1j}:{\delta}_{j}\ne 0,$$

(1)

and decide for a rejected *H*_{0j} which component *δ*_{ij}'s are non-zero before declaring their signs to be positive or negative depending on the values of the corresponding test statistics. The declared signs of the *δ*_{ij}'s then determine a possible inequality or directional pattern. For instance, in the case of *q* = 4, suppose for a given gene *j*, the *δ*_{j} = (*δ*_{1j}, … , *δ*_{4j}) is found significantly different from a null vector, with *δ*_{1j} and *δ*_{2j} declared to be positive and negative, respectively, and *δ*_{3j} and *δ*_{4j} zeros. Then, the corresponding directional pattern is *μ*_{1j} < *μ*_{2j} < *μ*_{3j} = *μ*_{4j} = *μ*_{5j}. We can test *H*_{0j} against *H*_{1j} for all the genes applying a suitable multiple testing method. Thus, given *p* ordered categories for each gene, the task of identifying directional patterns of the mean expressions over these categories for all the genes is being formulated as a multiple testing problem where *H*_{0j} is tested against *H*_{1j} simultaneously for all the genes and the signs of the *δ*_{ij}'s are determined subsequent to the rejection of the corresponding *H*_{0j}.

For multiple testing of *H*_{0j} against *H*_{1j}, *j* = 1, … , *m*, we need *p*-values that will provide a valid test for each of these individual testing problems and will allow us to make decisions on the individual *δ*_{ij}'s once a *H*_{0j} is rejected. For that, we consider for each *j* the *p*-value available for testing each component null hypothesis ${H}_{0j}^{i}\phantom{\rule{thinmathspace}{0ex}}:{\delta}_{\mathit{ij}}=0$ against the corresponding component alternative hypothesis ${H}_{1j}^{i}\phantom{\rule{thinmathspace}{0ex}}:{\delta}_{\mathit{ij}}\ne 0$, for *i* = 1, … , *q*, and apply a suitable combination method pooling these *qp*-values by treating *H*_{0j} as an intersection of the subfamily of these *q* component null hypotheses, that is, ${H}_{0j}={\bigcap}_{i=1}^{q}{H}_{0j}^{i}$, and *H*_{1j} as a union of the corresponding *q* alternative hypotheses, that is, ${H}_{1j}={\bigcup}_{i=1}^{q}{H}_{1j}^{i}$. Before we discuss appropriate combination methods to be used, let us explain how to obtain these component *p*-values and state the underlying assumptions.

For every *i* = 1, , *q* and *j* = 1, , *m*, suppose we use the absolute value of a test statistic *T*_{ij} for testing ${H}_{0j}^{i}$ against ${H}_{1j}^{i}$. Let *T*_{ij} ~ *F*_{ij}(*t*, *δ*_{ij}) for some continuous cdf F, which is symmetric about 0 under ${H}_{0j}^{i}$ and gets stochastically larger or smaller as *δ*_{ij} either increases or decreases from 0. In other words, with *F*_{ij}(*t*, *δ*_{ij}) denoting the cdf of *T*_{ij} at t under the parameter *δ*_{ij}, we have *F*_{ij}(*t*, *δ*_{ij}) or *F*_{ij}(*t*, 0) according as *δ*_{ij} > or < 0, and ${F}_{\mathit{ij}}(0,0)=\frac{1}{2}$. Under this setting, a right-tailed test based on the absolute value of *T*_{ij} will be considered for testing ${H}_{0j}^{i}$ against ${H}_{1j}^{i}$, with the corresponding two-sided *p*-value being defined as _{ij} = 2 min {*F*_{ij}(*T*_{ij}, 0), 1 − *F*_{ij}(*T*_{ij}, 0)}. By the assumed distributional property of *T*_{ij}, it is easy to verify that under ${H}_{0j}^{i}$, the two-sided *p*-value _{ij} satisfies

$$\mathit{Pr}\{{\stackrel{~}{P}}_{\mathit{ij}}\le p\}\le p,\phantom{\rule{1em}{0ex}}\text{for any}\phantom{\rule{1em}{0ex}}p\in (0,1).$$

(2)

Given *p*-values for testing ${H}_{0j}^{i}$ against ${H}_{1j}^{i}$, for *i* = 1, … , *q*, a number of combination methods (or methods of pooling the *p*-values) are available in the literature for testing the intersection null hypothesis ${H}_{0j}={\bigcap}_{i=1}^{q}{H}_{0j}^{i}$ against the alternative ${H}_{1j}={\bigcup}_{i=1}^{q}{H}_{1j}^{i}$. Among these, however, the Bonferroni and Simes methods are often used in multiple testing and allow one to make decisions on the individual *δ*_{ij}'s. For a review of these methods, one may see Bernhard et al. (2004). Let _{(1)j} _{(q)j} be the ordered versions of _{ij}, *i* = 1, , *q*, for a fixed *j* = 1, … , *m*. Then, in the Bonferroni test, the pooled (or adjusted) *p*-value is given by *P*_{j} = *q*_{(1)j}; whereas, in the Simes test, it is given by *P*_{j} = min_{1iq} {*q*_{(i)j}/*i*}. While the Bonferroni test does not require any dependence structure in the underlying *p*-values, the Simes test requires a certain type of positive dependence condition that is often satisfied in multiple testing applications (Sarkar and Chang, 1997). Upon rejection of *H*_{0j} using the Bonferroni pooled *p*-value at a level *α*, the *i*th component null hypothesis ${H}_{0j}^{i}$ can be rejected if _{ij}*α*/*q*. For the test based on the Simes pooled *p*-value, ${H}_{0j}^{i}$ corresponding to every _{ij}* _{(Rj)j}* is rejected, where ${R}_{j}=\mathrm{max}\{i:{\stackrel{~}{P}}_{\left(i\right)j}\le \frac{i}{q}\alpha \}$, if the maximum exists; otherwise,

Now, suppose the pooled *p*-value *P _{j}*, based on either Bonferroni or Simes test, is available to us for every

$$\mathit{FDR}=E\left\{\frac{V}{R\vee 1}\right\},$$

(3)

where *R*1 = max(*R*, 1). This method with a control of the FDR at a given level *α* is a stepup test that, given ordered *p*-values *P*_{(1)} *P _{(m)}* with the corresponding null hypotheses

When a *H*_{0j} : * δ_{j}* =

$$\mathit{dFDR}=E\left\{\frac{S}{R\vee 1}\right\},$$

(4)

where *S* denotes the total number of false null hypotheses among *H*_{1}, … , *H _{m}* that are correctly rejected but at least one directional error has been made while deciding upon the signs of the components. In other words,

$$\mathit{mdFDR}=\mathit{FDR}+\mathit{dFDR}=E\left\{\frac{V+S}{R\vee 1}\right\},$$

(5)

the expected proportion of Type I and directional errors among all rejections.

It is important to point out that the goal of this paper is to identify expression patterns of *m* genes over *p* ordered categories. For each gene it is biologically relevant to consider its expression pattern as a whole across *p* ordered categories rather than viewing this to be a problem of testing *qm* separate hypotheses which ignores the intrinsic biological structure present in the problem. Thus, rather than viewing it as a problem of performing *qm* tests, we treat it as a problem of performing a set of *m* tests each involving *q*-dimensional hypothesis. In addition, we want to emphasize that while making directional decisions for the components of * δ_{j}*, no directional errors are being made when

In the next section, we will develop methods to control the mdFDR. This will extend the following directional BH procedure of Benjamini and Yekutieli (2005) from dimension one (i.e., *q* = 1) to a general dimension.

Definition 1 (The level-*α* directional BH Procedure)

- (1)
*Apply the BH method at level α to test H*_{0j}:*δ*_{1j}= 0*against H*_{1j}:*δ*_{1j}≠ = 0*simultaneously for j*= 1, … ,*m, based on the two-sided p-values*_{1j }*j*= 1, …*m*. - (2) Let R denote the total number of null hypotheses rejected.
- (3)
*For every j*= 1, ,*m, with*${\stackrel{~}{P}}_{1j}\le \frac{R}{m}\alpha $, declare*δ*_{1j}>*or*< 0*according as T*_{1j}> 0 or < 0.

It controls the mdFDR at *α* under independence of the underlying test statistics.

We introduce in this section our proposed method of controlling the mdFDR while testing *H*_{0j} : * δ_{j}* =

- (1)
*Apply the BH method at level α to test H*_{0j}against*H*_{1j }*simultaneously for j*= 1, … ,*m, based on the Bonferroni pooled p-values P*= 1, ,_{j}, j*m*. - (2) Let R denote the total number of null hypotheses rejected.
- (3)
*For every i*= 1, ,*q and j*= 1, ,*m with*${\stackrel{~}{P}}_{\mathit{ij}}\le \frac{R}{\mathit{qm}}\alpha $,*if T*> 0,_{ij}*declare δ*>_{ij}*or*< 0*according as T*> 0_{ij}*or*< 0.

*With independent q-dimensional test statistics* **T**_{j} = (*T*_{1j}, , *T _{qj}*),

We offer a proof of Theorem 1 in the Appendix. Benjamini and Yekutieli (2005) gave an indirect proof of this theorem in the special case when *q* = 1 using an approach that relates to the FDR-adjusted confidence intervals for selected parameters they developed in the same paper. However it is not apparent how one could adapt their proof to the present case involving multiple parameters. So, we provide a direct proof in a more general setting.

In Theorem 1, we assume that *q*-dimensional test statistics **T**_{j}'s are independent. However, within each **T**_{j}, we do not impose any restriction on *T _{ij}*'s.

It would be tempting to develop an alternative method based on the Simes pooled *p*-values as follows:

- (1)
*Apply the BH method at level α to test H*_{0j }*simultaneously for j*= 1, … ,*m, based on the Simes pooled p-values P*= 1, … ,_{j}, j*m*. - (2) Let R denote the total number of null hypotheses rejected.
- (3)
*For every j*= 1, ,*m, let*_{(1)j}= 1, ,_{(q)j}be the ordered values of_{ij}, i*q*.

*Let* ${R}_{j}=\mathrm{max}\{i:{\stackrel{~}{P}}_{\left(i\right)j}\le \frac{i}{q}\cdot \frac{R}{m}\alpha \}$, *if the maximum exists; otherwise R _{j}* = 0.

As Simes test is known to be more powerful than the Bonferroni test (Simes, 1986), Procedure 2 would be more powerful than Procedure 1. Unfortunately, however, it would not control the mdFDR, as the Associate Editor pointed out. Consider, for instance, *m* = 1. The augmented test in this procedure in this case reduces to the step-up test with Simes critical values for the q hypotheses. Assume that *q* = 10 and that for half of the hypotheses *δ*_{i1} = 0 and for the remaining *δ*_{i1} is very large. Then the familywise error rate (FWER) of the step-up test with Simes critical values (for the test of the *q* hypotheses) is not controlled; see also Hommel (1988). However, in this scenario, mdFDR FWER. Therefore, Procedure 2 loses the control of the mdFDR in this situation. So, we do not formally propose it in this article as a multidimensional directional FDR controlling procedure, though we will consider it along with Procedure 1 in our simulation studies in the next section.

A simulation study was performed to evaluate the performance of our proposed method, Procedure 1. Specifically, we wanted to investigate the following three questions:

- (i) How does it perform in terms of its control of the FDR, dFDR and mdFDR and also power under independence as well as under types of dependence across the genes?
- (ii) How dose it perform in terms of the same operating characteristics under the independence across the genes when we benchmark it against Procedure 2 (based on the Simes pooled
*p*-values) and the procedure that makes no adjustment to the gene specific*p*-values, that is, simply uses_{(1)j}as the pooled*p*-value? - (iii) How does the performance of Procedure 1 under the independence across the genes change as the the dimension
*q*increases?

We generated *q* + 1 independently distributed *m*-dimensional random normal vectors **Z**_{1}, … , **Z**_{q+1}, where the components *Z _{ij}*,

Figure 1 presents the simulated FDR, dFDR and mdFDR, Figure 2 presents standard deviation of the simulated mdFDR, and Web Figure 1 presents the simulated average power (the proportion of false null hypotheses that are correctly rejected with correct assigned signs) of Procedure 1 plotted against the number of false null hypotheses for *m* = 1000, *q* = 5, *α* = 0.05 and *ρ* = 0 (independence), 0.2, 0.5 and 0.8.

Performance of Procedure 1 under dependence across genes in terms of its control of the FDR, dFDR and mdFDR for *m* = 1000, *q* = 5, *α* = 0.05, and *ρ* = 0, 0.2, 05 and 0.8.

Standard deviation of the mdFDR of Procedure 1 under dependence across genes for *m* = 1000, *q* = 5, *α* = 0.05, and *ρ* = 0, 0.2, 05 and 0.8.

Some interesting observations can be made from Figure 1. With increasing number of truly significant genes, the FDR steadily decreases to zero as long as the dependence across the genes is low or moderately high, while the dFDR slowly increases from zero to a value slightly less than 0.01, no matter what the dependence across the genes is, as long as it is non-negative. Consequently, when the genes are not too highly dependent, with increasing number of truly significant genes although the mdFDR decreases, implying that Procedure 1 as an mdFDR controlling procedure becomes more conservative, it does not however reach zero (see Figure 1(a)-(c)). When the genes are highly dependent, as we see from Figure 1(d), Procedure 1 becomes less conservative as the number of truly significant genes begins to increase from zero, but eventually it becomes more conservative as this number becomes larger.

Also, as seen from Figure 2, the standard deviation of the estimated mdFDR is very small. From Web Figure 1 we see that as the dependence across genes increases, the change in power is small. Overall, the effect of dependence across genes on the performance of the proposed procedure is relatively small. As suggested by the associate editor, under the above simulation settings, we also evaluated the performance of Procedure 1 for the * δ_{j}* = (100, 0, … , 0), in which one component is extremely large and the rest are zero. For such non-null

Web Figure 4 presents an answer to question (ii). As we see from this figure, Procedures 1 and 2 behave quite similarly, at least when the dependence across the genes is not of concern, in terms of controlling the FDR, dFDR and mdFDR and the power, though Procedure 2 is slightly more liberal as expected. Also as expected, if no adjustment is made to gene-specific p-values, we lose the control of the FDR and mdFDR, with the maximum reaching 0.2. It seems surprising that, even without any adjustment to gene specific *p*-values, the dFDR always remains low, though it becomes larger compared to that for Procedures 1 and 2 as the number of false nulls increases.

Web Figures 5-6 provide an answer to question (iii). It is interesting to note that the performance of Procedure 1 in terms of controlling the FDR, dFDR and mdFDR is unaffected by the dimension q when the dependence across the genes is not present. The power, of course, increases with increasing dimension.

Lobenhofer et al. (2002) investigated the effect of estrogen on the expression of cell-cycle genes as MCF-7 breast cancer cells go through the cell division cycle. A normal cell division cycle consists of four major phases, namely, the G1 (or Gap 1), S (Synthesis), G2 (or Gap 2) and M (Mitosis) phase. Genes involved in the cell cycle (known as *cell-cycle genes*) are expected to attain peak gene expression during the phase in which they have a specific biological function in the cell cycle.

According to Lobenhofer et al. (2002), most estradiol treated MCF-7 cells are expected to go through S, G2 and M phases in 12 to 36 hours after treatment and complete the cycle in 48 hours. Genes involved in cell growth and related activities are expected to have maximum expression (or minimum expression if they are anti-growth) during 1 or 4 hours and then monotonically decrease (or increase) in expression as cells go through the remaining phases. On the other hand, genes involved in DNA synthesis, repair and mitosis would have maximum (or minimum) expression during 12 to 36 hours. Thus, such genes may have an *Umbrella* (or *Inverted umbrella*) shaped pattern with a peak or trough during 12 to 36 hours time period. However, according to Lobenhofer et al. (2002), the cells may be asynchronous as they complete the cell division cycle at 48 hours after exposure. For this reason, the expression of some of the cell-cycle genes may not return to their baseline values at 48 hours but may attain a plateau.

Before exposing the MCF-7 breast cancer cells to estrogen, Lobenhofer et al. (2002) first synchronized all the cells to G1 phase by depriving the cells of serum for 24 hours. Synchronization of cells to the same phase at the beginning of the experiment is important for obtaining reliable gene expression data. They then harvested estradiol treated cells after 1, 4, 12, 24, 36 or 48 hours of treatment. Gene expressions using cDNA microarray chips were obtained at each time point. Each cDNA microarray chip consisted of 1900 gene probes. With 8 replicates at each time point, there were a total of 48 microarray chips across the 6 time points.

Motivated by the above observations, in this section we apply the proposed methodology to identify some cell-cycle genes by considering four ordered categories of time points, namely, 1 hour (*T*1), 4 hours (*T*2), mid group (i.e., the union of 12, 24, 36 hours) (*T*3) and 48 hours (*T*4) after treatment. Thus the sample sizes in the four groups are 8, 8, 24 and 8 respectively. Since the major cell division related activity takes place during the 12 to 36 hours time interval, we combined those time periods together to contrast that period from initial cell growth period (1, 4 hours) and the end of mitosis (48 hours).

Suppose *θ*_{T1,j}, *θ*_{T2,j}, *θ*_{T3,j}, and *θ*_{T4,j} denote the mean gene expression of gene *j*, *j* = 1, 2, … , 1900, during time periods *T*1, *T*2, *T*3 and *T*4, respectively. Using notations from the previous section, we let *δ*_{1j} = *θ*_{T1,j} − *θ*_{T2,j}, *δ*_{2j} = *θ*_{T2,j} − *θ*_{T3,j} and *δ*_{3j} = *θ*_{T3,j} − *θ*_{T4,j}.

Let *T _{ij}* denote the test statistic associated with the parameter

By applying our proposed method, Procedure 1 to the list of the pooled *p*-values *P _{j}*'s, we identified 86 differentially expressed genes at level

Comparing our results with those of Lobenhofer et al. (2002) and Peddada et al. (2003), we found that of the 86 genes we identified, 39 were also identified in at least one of the two previous papers. This included 8 of 13 DNA replication/repair genes identified by Lobenhofer et al. (2002). Among the 5 that were not identified by our procedure, we note that except for MCM7, which may be significant at *α* = 0.10, all others had large *p*-values that were not significant even at *α* = 0.20. Interestingly, in addition to MCM3 that was identified by both Lobenhofer et al. (2002) and Peddada et al. (2003), we identified a well known cell-cycle gene MCM4 (http://www.cyclebase.org).

An important step in DNA synthesis during the S phase is the binding of complex proteins to DNA for recruiting other proteins necessary for DNA synthesis. One such complex protein is the replication factor C. Lobenhofer et al. (2002) identified one subunit of this protein, known as replication factor C3. Later the order-restricted inference based methodology of Peddada et al. (2003) identified two additional subunits of this protein, namely, replication factors C4 and C5. Interestingly, the newly introduced methodology identified subunits C2, C3 and C5 as significant genes, thus reinforcing the earlier findings and adding one more subunit to the previous list of replication factor C. Furthermore, based on the proposed methodology it is possible to conclude that the subunits C2, C3 and C5 have peak expression during the 12, 24 or 36 hours time period where the DNA synthesis and replication takes place.

Furthermore, similar to the order-restricted inference procedure of Peddada et al. (2003), the proposed methodology identified the cyclin-dependent kinase inhibitor 1 A (p21 and Cip 1) as repressed during 12 to 36 hours. This gene was not identified by Lobenhofer et al. (2002).

A complete list of all genes identified by this procedure is provided in the Web Supplementary Materials accompanying this article.

In microarray gene expression studies, researchers are often not only interested in identifying differentially expressed genes under different biological conditions, but are also interested in detecting trends in mean response over ordered categories. For instance, in the simple case of two categories (normal versus tumor tissue), researchers are not only interested in identifying significant genes across these two categories, but they are also interested in further identifying the down and up-regulated genes. As the number of ordered categories increases, the trends or directional patterns become complex and the number of directional patterns increases. Except for the usual Type I errors, this also potentially results in a relatively high frequency of directional errors. Hence, it is important to develop statistical methods of identifying trends in mean response over ordered categories while maintaining a control over both the Type I and directional errors.

The approach proposed in this article provides such a methodology. Differently from existing statistical methods (Peddada et al., 2003; Lin et al., 2007), we have formulated the problem of identifying trends in mean response over ordered categories as a multiple testing problem involving successive comparisons and further directional decisions on the multidimensional parameter of each gene. To deal with this problem, we have first suggested a general multidimensional BH-type directional procedure using the Bonferroni test for controlling the mixed directional FDR (mdFDR), an overall measure of both Type I and directional errors within the framework of the FDR, and theoretically proved that the proposed procedure controls the mdFDR at a pre-specified level when the underlying test statistics are independent across the genes. We evaluated the performance of the introduced procedure in the case of dependence through a simulation study. Finally, the whole proposed methodology has been applied to analyze a time-course microarray data and some interesting results have been obtained.

Although our focus was on identifying individual gene expression profile or trend over the ordered categories in some common microarray experiments such as time-course or dose-response experiments, the proposed methodology can be also applied in National Toxicology Program (NTP) studies, where researchers are interested in determining whether for a given tumor type, there is a significant dose effect and then identifying its dose-response profile.

The methodology proposed in this article provides an interesting starting point towards addressing the complex yet important problem of controlling both Type I and directional errors in multiple testing involving multidimensional parameters. The mdFDR controlling property of the proposed directional BH procedure has been established under the assumption that the underlying test statistics are independent across the genes. When gene expressions are obtained by drawing samples from same subjects over time, such an assumption need not be valid. In such cases, not only do we have dependence among gene expressions at a given time point but there may be temporal dependence among gene expressions at different time points. It will be interesting to theoretically investigate the performance of the proposed directional BH procedure under such complex dependence structures. In addition, it will also be interesting to develop more powerful adaptive BH directional FDR procedure by exploiting knowledge of the proportion of true null hypotheses.

The research of the first and third authors is supported by the Intramural Research Program of the NIH, National Institute of Environmental Health Sciences(Z01 ES101744-04) and the research of the second author is supported by a grant from U.S. National Science Foundation (DMS-0603868). We thank Grace Kissling and David Dunson for several useful comments that led to improved presentation of the manuscript. We also thank the editor, the associate editor (AE) and the referees for their comments which have led to substantial improvement in the presentation of the manuscript. We are particularly thankful to the AE for spotting an error in the proof of Theorem 1, for providing an interesting example showing that Procedure 2 cannot control the mdFDR (described in Remark 3), and for making several other important comments that led to improved simulation studies and better exposition of the paper.

7. Supplementary Materials

Web Appendices and Figures referenced in Sections 3 and 4 are available under the Paper Information link at the *Biometrics* website http://www.biometrics.tib.org. A complete list of all genes identified in Section 5 are also available at the Biometrics website.

Wenge Guo, Biostatistics Branch, National Institute of Environmental Health Sciences Research Triangle Park, NC 27709, U.S.A. Email: moc.liamg@oug.egnew.

Sanat K. Sarkar, Department of Statistics, Temple University, Philadelphia, PA 19122, U.S.A. Email: ude.elpmet@tanas.

Shyamal D. Peddada, Biostatistics Branch, National Institute of Environmental Health Sciences Research Triangle Park, NC 27709, U.S.A. Email: vog.hin.shein@adaddep..

- Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J. Roy. Statist. Soc. Ser. B. 1995;57:289–300.
- Benjamini Y, Hochberg Y. On the adaptive control of the false discovery rate in multiple testing with independent statistics. Journal of Educational and Behavioral Statistics. 2000;25:60–83.
- Benjamini Y, Hochberg Y, Kling Y. Working paper 93-2. Tel Aviv University, Dept. of Statistics and Operation Research; 1993. False discovery rate control in pairwise comparisons.
- Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann. Statist. 2001;29:1165–1188.
- Benjamini Y, Yekutieli D. False discovery rate-adjusted multiple confidence intervals for selected parameters. J Amer. Statist. Assoc. 2005;100:71–93.
- Bernhard G, Klein M, Hommel G. Global and multiple test procedures using ordered p-values – A review. Statist. Pap. 2004;45:1–14.
- Bochkina N, Richardson S. Tail Posterior Probability for Inference in Pairwise and Multiclass Gene Expression Data. Biometrics. 2007;63:1117–1125. [PubMed]
- Finner H. Testing multiple hypotheses: general theory, specific problems, and relationships to other multiple decision procedures. Habilitationsschrift, Fachbereich IV Mathematik, Univ. Trier; 1994.
- Finner H. Stepwise multiple test procedures and control of directional errors. Ann. Statist. 1999;27:274–289.
- Hommel G. A stagewise rejective multiple test procedure based on modified Bonferroni test. Biometrika. 1988;75:383–386.
- Lin D, Shkedy Z, Yekutieli D, Burzykowski T, Göhlmann H, Bondt A, Perera T, Geerts T, Bijnens L. Testing for trends in dose-response microarray experiments: a comparison of several testing procedures, multiplicity and resampling-based inference. Statist. App. Gen. Mol. Bio. 2007;6(1) Article 26. [PubMed]
- Liu W. Control of directional errors with step-up multiple tests. Statist. Probab. Lett. 1996;31:239–242.
- Lobenhofer E, Bennett L, Cable P, Li L, Bushel P, Afshari C. Regulation of DNA replication fork genes by 17 beta-estradiol. Molecular Endocrinology. 2002;16:1215–1229. [PubMed]
- Peddada S, Lobenhofer E, Li L, Afshari C, Weinberg C, Umbach D. Gene selection and clustering for time-course and dose response microarray experiments using order-restricted inference. Bioinformatics. 2003;19:834–841. [PubMed]
- Sarkar SK. Some results on false discovery rate in stepwise multiple testing procedures. Ann. Statist. 2002;30:239–257.
- Sarkar SK, Chang C-K. The Simes method for multiple hypothesis testing with positively dependent test statistics. J. Amer. Statist. Assoc. 1997;92:1601–1608.
- Sarkar SK, Sen PK, Finner H. On two results in multiple testing. In: Benjamini Y, Bretz F, Sarkar S, editors. Recent Developments in Multiple Comparisons. Institute of Mathematical Statistics; Beachwood: 2004. pp. 89–99. (IMS Lectures Notes-Monograph Series, 47).
- Shaffer JP. Control of directional errors with stagewise multiple test procedures. Ann. Statist. 1980;8:1342–1347.
- Shaffer JP. Multiplicity, directional (type III) errors, and the null hypothesis. Psychological Methods. 2002;7:356–369. [PubMed]
- Simes RJ. An improved Bonferroni procedure for multiple tests of significance. Biometrika. 1986;73:751–754.
- Tamoto E, Tada M, Murakawa K, Takada M, Shindo G, Teramoto K, Matsunaga A, Komuro K, Kanai M, Kawakami A, Fujiwara Y, Kobayashi N, Shirata K, Nishimura N, Okushiba S, Kondo S, Hamada J, Yoshiki T, Moriuchi T, Katoh H. Gene-expression profile changes correlated with tumor progression and lymph node metastasis in esophageal cancer. Clinical Cancer Research. 2004;10:3629–3638. [PubMed]
- Williams VS, Jones V, Tukey JW. Control errors in multiple comparisons, with examples from state-to-state differences in educational achievement. Journal of Educational and Behavioral Statistics. 1999;24:42–69.

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |