 Methodology article
 Open Access
 Published:
A permutationbased multiple testing method for timecourse microarray experiments
BMC Bioinformatics volume 10, Article number: 336 (2009)
Abstract
Background
Timecourse microarray experiments are widely used to study the temporal profiles of gene expression. Storey et al. (2005) developed a method for analyzing timecourse microarray studies that can be applied to discovering genes whose expression trajectories change over time within a single biological group, or those that follow different time trajectories among multiple groups. They estimated the expression trajectories of each gene using natural cubic splines under the null (no timecourse) and alternative (timecourse) hypotheses, and used a goodness of fit test statistic to quantify the discrepancy. The null distribution of the statistic was approximated through a bootstrap method. Gene expression levels in microarray data are often complicatedly correlated. An accurate type I error control adjusting for multiple testing requires the joint null distribution of test statistics for a large number of genes. For this purpose, permutation methods have been widely used because of computational ease and their intuitive interpretation.
Results
In this paper, we propose a permutationbased multiple testing procedure based on the test statistic used by Storey et al. (2005). We also propose an efficient computation algorithm. Extensive simulations are conducted to investigate the performance of the permutationbased multiple testing procedure. The application of the proposed method is illustrated using the Caenorhabditis elegans dauer developmental data.
Conclusion
Our method is computationally efficient and applicable for identifying genes whose expression levels are timedependent in a single biological group and for identifying the genes for which the timeprofile depends on the group in a multigroup setting.
Background
Timecourse microarray experiments are widely used to study the temporal profiles of gene expression. In these experiments, the gene expressions are measured across several timepoints, enabling the investigator to study the dynamic behavior of gene expressions over time.
A number of statistical methods have been developed in recent years for identifying differentially expressed genes from timecourse microarray experiments. Park et al. [1] proposed a permutationbased twoway ANOVA to compare temporal profiles from different experimental groups. Luna and Li [2] proposed a statistical framework based on a shapeinvariant model together with a false discovery rate (FDR) procedure for identifying periodically expressed genes based on microarray time course gene expression data and a set of known periodically expressed guide genes. Storey et al. [3] represented gene expression trajectories using natural cubic splines and then compared the goodness of fit of the model under the null hypothesis to that under alternative hypothesis. The null distribution of these statistics was approximated through a bootstrap method. Di Camillo et al. [4] proposed test statistics using the maximum distance between two time trajectories or comparing the areas under two time course curves. Approximating the null distribution of the test statistics using a bootstrap method, they show that their test statistics are more powerful than Storey et al. [3] if the number of measurement time points is small. Hong and Li [5] introduced a functional hierarchical model for detecting temporally differentially expressed genes between two experimental conditions for cross sectional designs, where the gene expression profiles are treated as functional data and modelled by basis function expansions. Angelini et al. [6] modelled timecourse data within a framework of a Bayesian hierarchical model and use Bayes factors for testing purposes.
Permutation resampling methods have been popularly used to derive the null distribution of highdimensional test statistics while preserving the complicated dependence structure among genes in microarray data analysis. In this paper, we present a permutationbased multipletesting method for timecourse microarray experiments when independent subjects contribute gene expression data at different time points. While the method can be generalized to broad class of goodnessoffit test statistics for regression curves, for illustration we use the Ftest type statistic based on natural splines used by Storey et al. [3]. We propose computationally efficient algorithms for identifying the genes whose expression levels are timedependent in a single biological group and for identifying the genes whose timeprofile differs among different groups. For the multiple group setting, we will consider two sets of hypotheses. In the first set, any difference among the curves, including vertically shifted parallel curves, is considered to constitute a discrepancy among the groups. For the second set, only differences in the actual timetrends are considered to be of interest after removing the vertical shift. We shall refer to these as "timecourse" and "timetrend" hypotheses, respectively. Note that if two separated curves can be overlapped by a vertical shift, then they have different timecourses, but the same timetrend. The test on a timetrend hypothesis will remove potential batch effects in microarray experiments.
The rest of the article is organized as follows. We first present a nonparametric test method to identify differential gene expression in a timecourse microarray. We then present simulation results to evaluate the statistical properties of the proposed method. Next, we apply the proposed method to the Caenorhabditis elegans dauer developmental data [7]. Lastly, we give a brief discussion of the methods.
Methods
At first, we briefly review a smoothing method to estimate a gene expression profile over time. Using the smoothing method, we discuss a nonparametric test method for identifying genes whose expression levels are timedependent in a single biological group and for identifying the genes for which the timeprofile depends on the group among multiple groups. We approximate the null joint distribution of the test statistics using a permutation method.
Estimation of the TimeCourse Profile
Suppose that subject i(= 1,..., n) contributes gene expression levels on m genes (y_{i 1},..., y_{ im }) at time t_{ i }. For gene j(= 1,..., m), we consider a time trajectory model E(y_{ ij }t) = μ_{ j }(t), where μ_{ j }(·) is the unknown function that is parameterized by an intercept plus a pdimensional linear basis:
Here [W_{1}(t),..., W_{ p }(t)] is a prespecified pdimensional basis that is common to all m genes, and β_{ j }= [β_{0, j}, β_{1, j},..., β_{p, j}]^{T}is a (p + 1)dimensional vector of unknown parameters for gene j. Similar to Storey et al. [3], we employ a Bspline basis (see chapter IX in de Boor [8]) and place the knots at the 0,1/(p  1), 2/(p  1),...,(p  2)/(p  1), 1 quantiles of the observed time points.
Let
W denotes the design matrix based on the spline model. Then, the least square estimator of β_{ j }is obtained by = (W^{T}W)^{1}W^{T}y_{ j }, where y_{ j }= (y_{1j},..., y_{ nj })^{T}.
One Group Case
In the case of a single biological group (K = 1), we often want to discover genes whose expression levels are timedependent. For gene j(= 1,..., m), we want to test the hypotheses
against
Under H_{ j }, the constant is estimated as . Under , we obtain the estimate , where () is estimated as described in the previous section.
For gene j, the sum of squares of errors (SSE) is expressed as . Let and denote the SSE under H_{ j }and , respectively. Storey et al. [3] employ the Fstatistic
for testing H_{ j }against . It is noted that for the permutationbased multiple testing described below, the (n  p  1)/p factor in the F_{ j }test statistic will have no impact on the results and as such can be omitted from the computations.
In order to generate the null distribution of the vector of test statistics (F_{1},..., F_{ m }) for the m genes, we randomly match the microarray of n subjects {(y_{i 1},..., y_{ im }), i = 1,..., n} with their measurement times {t_{1},..., t_{ n }} at each permutation. Let (,..., ñ) be a permutation of (1,..., n). Then {(, y_{i 1},..., y_{ im }), i = 1,..., n} is a permutation sample of the original data {(t_{ i }, y_{i 1},..., y_{ im }), i = 1,..., n}.
Familywise error rate (FWER) is defined by the probability of rejecting any null hypothesis H_{ j }when all m null hypotheses are true. A singlestep multiple testing procedure controlling the FWER at α level can be described as follows, refer to e.g., Westfall and Young [9] and Jung et al. [10].
Multiple Testing for Time Trend of One Group

1.
Compute the the Ftest statistics (f_{1},..., f_{ m }) from the original data.

2.
From the bth permutation data (b = 1,..., B), compute the Ftest statistics ().

3.
Singlestep procedure to control the FWER

(a)
From the bth permutation data, calculate u_{ b }= max_{1≤j≤m}.

(b)
For gene j, calculate the adjusted pvalue by , where I(·) is an indicator function.

(c)
For a specified FWER level α, discover gene j if <α.
False discovery rate (FDR) is another popular type I error for multiple testing adjustment that is defined by the expected value of the proportion of the number of erroneously rejected null hypotheses among the total number of rejected null hypotheses, refer to Benjamini and Hochberg [11]. A multiple testing procedure to control the FDR at α level can be obtained by replacing Step 3 in above algorithm with Step 3' as described below, refer to Tusher et al. [12] and Storey [13].
3'. Multuple testing controlling the FDR:
(a) For gene j, estimate the marginal pvalue by p_{ j }= B^{1} ≥ f_{ j }).

(b)
For a chosen constant λ ∈ (0, 1), such as 0.95 [13], estimate the qvalue of gene j by

(c)
For a specified FDR level α, discover gene j (or reject H_{ j }) if q_{ j }<α.
The testing algorithm can be considerably simplified during permutations. First, is invariant under permutations, and as such one does not have to recalculate for the permutation samples. Second, suppose that we fix the gene expression data {(y_{i 1},...., y_{ im }), i = 1,..., n} and shuffle the measurement times t_{1},..., t_{ n }in each permutation. Let I denote the n × n identity matrix. Then, noting that = {I  W(W^{T}W)^{1}W^{T}}y_{ j }, permutation replicates of can be obtained by simply permuting the columns of I  W(W^{T}W)^{1}W^{T}. Thus, I  W(W^{T}W)^{1}W^{T}does not have to be recomputed for the permutation samples. Furthermore, given that m is considerably larger than n, permuting the columns of I  W(W^{T}W)^{1}W^{T}, a matrix of dimension n × n, is more efficient than permuting the rows of [y_{1},...,y_{ m }], a matrix of dimension m × n.
K Group Case
In order to compare the timecourse profiles of gene expression measurements among different experimental groups, we assume that a fixed number of measurement times are prespecified commonly among the K groups and at least one subject is assigned to each time point from each group. Let t_{1} < ⋯ <t_{ L }denote the L time points chosen, and n_{ kl }denote the number of patients from group k(= 1,..., K) observed at time t_{ l }(l = 1,..., L). We use the notations n_{k·}= to denote the number of patients from group k and n_{·l}= to denote the number of patients at time point l. So, denotes the total number of subjects in the study. The design and sample size under each condition is summarized in Table 1.
Let (y_{kli 1},..., y_{ klim }) denote the expression measurements for m genes at time t_{ l }(l = 1,..., L) from subject i(= 1,..., n_{ kl }) belonging to group k(= 1,..., K). The expression values are modelled as
where μ_{ kj }(t) = β_{0, kj}+ .
In the Kgroup setting, we want to identify the genes with different time profiles in different groups. The hypotheses for gene j are specified as
against
Under , the estimator is estimated from the group k data, {(t_{ l }, y_{ klij }), 1 ≤ i ≤ n_{ kl }, 1 ≤ l ≤ L}. Let .
Under H_{ j }: μ_{1j}(t) = ⋯ = μ_{ Kj }(t)(= μ_{ j }(t)), the groupfree estimator is estimated using the pooled data, {(t_{ l }, y_{ klij }), 1 ≤ i ≤ n_{ kl }, 1 ≤ k ≤ K, 1 ≤ l ≤ L}. Let denote the estimator of the common time trajectory under H_{ j }.
For gene j, the SSE under H_{ j }is calculated as , where (t) is the estimate of μ_{1j}(t) = ⋯ = μ_{ Kj }(t) from the pooled data. The SSE under is calculated as , where (t) is the estimate of μ_{ kj }(t) from the group k data.
We reject H_{ j }in favor of for a large value of the Fstatistic
The null distribution of the test statistics (F_{1},..., F_{ m }) is approximated using a permutation method. A permutation sample is generated by permuting the gene expression data within each time point: the gene expression data of n_{·l}subjects at time t_{ l }, {(y_{kli 1},..., y_{ klim }), 1 ≤ i ≤ n_{ kl }, 1 ≤ k ≤ K} are randomly partitioned into K groups of size n_{1l},..., n_{ Kl }. The subjects at different time points are not permuted. For each subject, the random vector (y_{kli 1},..., y_{ klim }) is counted as a data point, so that the m genes are not permuted. One permutation sample is obtained by conducting this permutation process for all L time points. Note that there are
different permutations. Table 2 demonstrates a permutation when K = 2. The proposed restricted permutation maintains the time trend in the whole population and allows heteroscedastic error models. Multiple testing to control FDR or FWER is conducted as in one group case, but by using the K group F statistics and permuting the observed expressions within each timepoint. We can save computing time by utilizing the fact that the design matrices of the regression models are invariant to permutations.
Results
Simulations
In this section, we investigate the performance of our method for control of the FWER and power using extensive simulations. We also apply the proposed methods to a real data set.
Simulation Study
The three scenarios considered are based on amplitude variation, phase variation and a homoscedastic versus a heteroscedastic error model. We restrict ourselves to the single and twogroup (i.e., K = 2) cases.
Simulation Settings
We set m = 1, 000. Given a trend μ_{ j }(t) for gene j(= 1,..., m), expression data (y_{1},..., y_{ m }) measured at time t are generated by
Let a_{1},..., a_{1000} and b_{1},..., b_{100} be IID N (0, 1) random variables. Then, heteroscedastic error terms are generated as follows. For l = 1,..., 100 and j = 1,...,10, we generate ϵ_{10(l1)+j}= a_{10(l1)+j} + b_{ l }.
Note that the error terms (ϵ_{1},...,ϵ_{ m }) consist of 100 independent blocks of size 10, and the error terms in block l(= 1,..., 100), (ϵ_{10(l1)+1},...,ϵ_{10l}), have a compound symmetry correlation structure with correlation coefficient ρ, which is set at 0, 0.3 or 0.6. We choose L = 11 measurement times t_{ l }= 0, 1,..., 10, and simulate 4 replications at each time point for each group.
In a singlegroup case, nonprognostic genes (genes under H_{ j }) have model μ_{ j }(t) = 0, and prognostic genes (genes under ) have μ_{ j }(t) = 4 exp(t) in Simulation 1 and μ_{ j }(t) = sin(2πt) in Simulation 2.
In a twogroup case (K = 2), we consider three different simulation models. In Simulation 1 (amplitude variation model), nonprognostic genes have equal time trends for both groups μ_{1j}(t) = μ_{2j}(t) = exp(t), and prognostic genes have μ_{1j}(t) = exp(t) for group 1 and μ_{2j}(t) = 2.5 exp(t) for group 2, see the left panel of Figure 1. In Simulation 2 (phase variation model), nonprognostic genes have equal time trends for both groups μ_{1j}(t) = μ_{2j}(t) = sin(2πt), and prognostic genes have μ_{1j}(t) = sin(2πt) for group 1 and μ_{2j}(t) = sin(2π(t  1/4)) for group 2, see the right panel of Figure 1.
In Simulations 1 and 2, all m = 1, 000 genes are nonprognostic under the global null hypothesis . Under a specific alternative hypothesis , the first m_{1} = 10 genes are prognostic, and the remaining m_{0} = 990 genes are nonprognostic.
In Simulation 3 of a twosample case, we consider heteroscedastic error models. Nonprognostic genes have μ_{1}(t) = μ_{2}(t) = t, and prognostic genes have μ_{1}(t) = t and μ_{2}(t) = 2.5 + t. For both groups (k = 1, 2), the first 100 genes (1 ≤ j ≤ 100) have heteroscedastic error terms t^{1.5} × ϵ_{ kj }, and the remaining 990 genes (101 ≤ j ≤ 1, 000) have homoscedastic error terms ϵ_{ kj }. We generate (ϵ_{1},...,ϵ_{ m }) from the blocked compound symmetric multivariate normal distribution as in a homoscedastic error model. The first 5 genes with heteroscedastic error terms (1 ≤ j ≤ 5) and the first 5 genes with homoscedastic error terms (101 ≤ j ≤ 105) are prognostic, and all the remaining 990 genes are nonprognostic. Figure 2 displays expression levels of a nonprognostic gene under the homoscedastic error model (left panel) and under the heteroscedastic error model (right panel).
Under each setting, N = 1, 000 simulation samples are generated and the singlestep procedure to control the FWER at 5% is applied to each sample. The null distribution of the test statistic is approximated from B = 1, 000 resampling (permutation or bootstrap) replications for each simulation sample. An empirical FWER under H_{0} (or the global power under H_{ a }) is obtained by the proportion of samples that reject any H_{ j }.
The bootstrap method by Storey et al. [3] generates the resampling data under null distribution as follows. We consider one group case here, but cases with K groups are done similarly.

1.
Fit the timecourse model under , and calculate n residuals, e_{ ij }= y_{ ij } (t_{ i }).

2.
Fit the timecourse model under H_{ j }to obtain the fitted population average .

3.
Generate a resampling data set under H_{0}, {(), i = 1,..., n}, by randomly selecting the residual vectors (e_{i 1},..., e_{ im }) among n subjects and adding to the vector of fitted values ().
Simulation Results
Simulation results are reported in Table 3 under H_{0} and in Table 4 under H_{ a }. From Table 3, we observe that both the permutation method (PERM) and the bootstrap method (BOOT) accurately control the FWER under the homoscedastic error models (Simulations 1 and 2). Under the heteroscedastic error model (Simulation 3), however, the bootstrap method is very anticonservative, while the permutation method still control the FWER accurately. From Table 4, we observe that the two methods have almost identical global power in the homoscedastic error models. Power comparison under the heteroscedastic error model is meaningless since the bootstrap method does not control the FWER under H_{0}.
Case Study
In this section, we present the results from applying our method to the analysis of the Caenorhabditis elegans dauer developmental data discussed by Wang and Kim [7] who use cDNA microarrays to profile gene expression differences during the transition from the dauer state to the nondauer state (experimental group) and after feeding of starved L_{1} worms (control group). The cDNA microarray expressions are measured on m = 18,556 genes to examine the transition from dauer into normal development, where dauer animals were harvested at 0, 1.5, 2, 3, 4, 5, 6, 7, 8, 10, and 12 hours after feeding and each time point was repeated three or four times. Wang and Kim [7] perform another cDNA microarray experiment to profile gene expression at 0, 1, 2, 3, 4, 5, 6, 7, 8, 10, and 12 hours after feeding of starved L_{1} worms and each time point was repeated four times. This data set is available for download at http://cmgm.stanford.edu/~kimlab/dauer/. For the purpose of permutation within each measurement time, we need to unify the measurement times between groups. So, we regard the time point t = 1.5 in the experimental group as t = 1.
For this analysis, we will consider both timecourse and timetrend hypotheses. A timecourse hypothesis for a gene is to test any discrepancy in trajectory of its expression level over time as we have considered so far. In contrast, a timetrend hypothesis is to test any discrepancy in time trend of the gene's expression level after removing the difference in overall expression level between groups. For testing a timetrend hypothesis, the testing procedures we have discussed in the methods section can be extended by simply subtracting off groupspecific means at each time point from the observed expressions first. We will contrast the results from our permutation method to those obtained by the bootstrap method suggested by Storey et al. [3]. Each analysis is based on B = 10, 000 resampling replicates and a natural spline basis as the one used in the simulation studies.
The top sixteen genes in terms of the realized value of the F statistic for testing the timecourse and timetrend hypotheses are shown in Figures 3 and 4, respectively. In each case, the estimated time trajectory for each group is superimposed. For the timecourse hypothesis (Figure 3), most of the top genes (e.g., Y59A8C.D, F46F2.3) seemingly fall into the vertical shift category while a few (e.g., B0511.5, K06A4.1) seemingly exhibit differing timetrends. This is perhaps not surprising as the F statistic tends to be largest if the curves are separated by a vertical shift. Timetrend test (Figure 4) identifies genes for which the timetrends differ between the two arms.
Next, we will compare the result from applying our method to those obtained by employing the bootstrap approach. The number of significant genes, at a given FWER level, based on permutation and bootstrap FWER adjusted Pvalues, are shown in Figure 5 for the timecourse (topleft panel) and timetrend (topright panel) hypotheses. The permutation method tends to discover more genes for a FWER level of 0.07 or higher under both timecourse and time trend hypotheses. As illustrated in Figure 5 (bottomleft panel), at the FWER level of 0.05, for the timecourse hypothesis, there are 624 genes selected by the bootstrap method but not by the permutation method for the timecourse hypothesis. For the timetrend hypothesis (bottomright panel), twenty genes are identified by the permutation method but not by the bootstrap method, while 93 genes are selected by the bootstrap method but not by the permutation method. The supplementary material provides the biological properties of 13 genes (out of 20) that are identified only by the permutation method [see Additional file 1].
From each of these three sets of nonempty symmetric differences, the 9 genes with the largest difference in FWER adjusted Pvalues between the permutation and the bootstrap methods are illustrated in Figures 6 to 8. As we have illustrated in the simulation study, the bootstrap method may be severely anticonservative if the errors are heterogeneous over time. This may explain the large set of genes that are significant according to the bootstrap method but not by the permutation method for the timecourse hypothesis. The spline estimator used is not robust estimator of the regression curve in presence of outliers in which case it may give the misleading impression that the time trajectories are time dependent when in fact they are horizontal lines. Another thing to note is that, in some cases, the difference between the two time trajectories is primarily driven by their difference at the baseline, t = 0. It is conceivable that some of these genes would not be prognostic if the observations at baseline were to be omitted from the analysis.
Discussion
We have considered two sets of hypotheses in the multigroup setting. For the timecourse hypothesis, any difference among the groups, including parallel curves shifted vertically, would be considered interesting. For the timetrend, only cases where the timetrend is group dependent would be of interest. This method has several advantages compared to the bootstrap method suggested by Storey et al. [3]. First, as our simulation results have shown, the bootstrap method may not control the FWER if the error variability for each gene is heterogeneous over time. The permutation method, on the other hand, controls the FWER in the heteroscedastic case as it only requires exchangeability within time points under the null hypothesis. The bootstrap method is based on the restrictive assumption that the error model is additive and that the error terms are not only exchangeable within but also across time points.
Second, the bootstrap method requires that, in addition to matrix of observed expressions, the matrix of residuals be stored to avoid recalculating them for each bootstrap replication. Thus, compared to the permutation method, the memory requirement for the bootstrap method is about twice as large. We have illustrated our permutation method by employing the regression goodnessoffit statistic based on natural splines used by Storey et al. [3]. Our method can be extended by using other regression goodnessoffit statistics. More specifically, if one is solely interested in testing for significant genes, but not in estimating the time trajectories, then one could consider using a simple mean trace model where the timetrajectory at each point is estimated by averaging the expressions. This statistic may be more sensible if the number of timepoints is small.
A referee requested that we compare the power between the Fstatistic based on the estimated timetrajectory at each point by averaging the observations with that based on the smoothed timetrajectory proposed in this paper. For simplification, we conducted simulations in a single gene case.
For subject i assigned to group k(= 1, 2) and time t_{ l }(= 0, 1,...,10), the gene expression level was generated by y_{ kli }(t_{ l }) = sin[2π{t_{ l } (k  1)/4}] + ϵ_{ kli }, where ϵ_{ kli }are IID N (0, 1) random variables. Four subjects were assigned to each time point for each group, so that n = 88(= 2 × 11 × 4). We generated 10,000 simulation samples and each sample was permuted B = 1, 000 times. At α = 0.05 level, the empirical power of the statistic based on the smoothed timetrajectories was 0.9572 while that of the standard Fstatistic is 0.8644. We observed similar comparisons under the wide range of simulation settings.
In the methods section, for the onesample case we have proposed an efficient algorithm based on permuting columns of the projection matrix, rather than entire matrix of expressions. To evaluate the gain in efficiency empirically, we compare the two approaches for calculating FWERadjusted Pvalues based on m = 10,000 genes, n = 100 patients and B = 10,000 permutations. For each approach, we replicate this simulation 10 times. The mean processing times on an AMD Opteron 8200 processor are 1,850 seconds based on permuting the matrix of expressions versus 1,764 seconds based on permuting only the columns of the projection matrix. Our approach is not only more elegant, but, as this example illustrates, may provide considerable gain in efficiency for large scale simulations such as those used in empirical power calculations where the number of simulation replicates for each design scenario and the number of markers greatly exceed 10 and 10,000 respectively.
Computer programs in R are available from http://www.duke.edu/~is29/TC/.
Conclusion
In conclusion, our permutationbased multiple testing method for timecourse microarray experiments is computationally efficient and applicable for identifying the genes whose expression levels are timedependent in a single biological group or for identifying the genes whose timeprofiles are different among different groups.
References
 1.
Park T, Yi S, Lee S, Lee SY, Yoo DH, Ahn JI, Lee YS: Statistical tests for identifying differentially expressed genes in timecourse microarray experiments. Bioinformatics 2003, 19: 694–703. 10.1093/bioinformatics/btg068
 2.
Luan Y, Li H: Modebased methods for identifying periodically regulated genes based on time course microarray gene expression data. Bioinformatics 2004, 20: 332–339. 10.1093/bioinformatics/btg413
 3.
Storey JD, Xiao W, LeeK JT, Tompkins RG, Davis RW: Significance analysis of time course microarray experiments. Proc Natl Acad Sci 2005, 102: 12837–12842. 10.1073/pnas.0504609102
 4.
Di Camillo B, Toffolo G, Nair SK, Greenlund LJ, Cobelli C: Significance analysis of microarray transcript levels in time series experiments. BMC Bioinformatics 2007, 8(Suppl 1):S10. 10.1186/147121058S1S10
 5.
Hong F, Li H: Functional Hierarchical Models for Identifying Genes with Different TimeCourse Expression Profiles. Biometrics 2006, 62: 534–544. 10.1111/j.15410420.2005.00505.x
 6.
Angelini C, De CD, Mutarelli M, Pensky M: A Bayesian Approach to Estimation and Testing in Timecourse Microarray Experiments. Statistical Applications in Genetics and Molecular Biology 2007., 6(1):
 7.
Wang J, Kim S: Global analysis of dauer gene expression in Caenorhabditis elegans. Development 2003, 130: 1621–1634. 10.1242/dev.00363
 8.
de Boor C: A Practical Guide to Splines. SpringerVerlag: New York; 2001.
 9.
Westfall PH, Young SS: Resamplingbased Multiple Testing: Examples and Methods for Pvalue Adjustment. Wiley: New York; 1993.
 10.
Jung SH, Bang H, Young S: Sample size calculation for multiple testing in microarray data analysis. Biostatistics 2005, 6: 157–169. 10.1093/biostatistics/kxh026
 11.
Benjamini Y, Hochber Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. JR Statist Soc B 1995, 57: 289–300.
 12.
Tusher V, Tibshirani R, Chu G: Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci 2001, 98: 5116–5121. 10.1073/pnas.091062498
 13.
Storey JD: A direct approach to false discovery rates. Journal of the Royal Statistical Society, Series B 2002, 64(1):479–498.
Author information
Affiliations
Corresponding author
Additional information
Authors' contributions
IS and KO performed statistical analysis and wrote the manuscript. SLG supported the research. SK conducted the biological interpretation of the statistical analysis results. SJ proposed the research project. All authors read and approved the final manuscript.
Authors’ original submitted files for images
Below are the links to the authors’ original submitted files for images.
Rights and permissions
About this article
Cite this article
Sohn, I., Owzar, K., George, S.L. et al. A permutationbased multiple testing method for timecourse microarray experiments. BMC Bioinformatics 10, 336 (2009). https://doi.org/10.1186/1471210510336
Received:
Accepted:
Published:
Keywords
 Bootstrap Method
 Null Distribution
 Prognostic Gene
 Time Trajectory
 Permutation Method