- Methodology article
- Open Access
A permutation-based multiple testing method for time-course microarray experiments
- Insuk Sohn^{1},
- Kouros Owzar^{1, 2},
- Stephen L George^{1, 2},
- Sujong Kim^{3, 4} and
- Sin-Ho Jung^{1, 2}Email author
https://doi.org/10.1186/1471-2105-10-336
© Sohn et al; licensee BioMed Central Ltd. 2009
- Received: 18 March 2009
- Accepted: 15 October 2009
- Published: 15 October 2009
Abstract
Background
Time-course microarray experiments are widely used to study the temporal profiles of gene expression. Storey et al. (2005) developed a method for analyzing time-course 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 time-course) and alternative (time-course) 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 permutation-based 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 permutation-based 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 time-dependent in a single biological group and for identifying the genes for which the time-profile depends on the group in a multi-group setting.
Keywords
- Bootstrap Method
- Null Distribution
- Prognostic Gene
- Time Trajectory
- Permutation Method
Background
Time-course microarray experiments are widely used to study the temporal profiles of gene expression. In these experiments, the gene expressions are measured across several time-points, 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 time-course microarray experiments. Park et al. [1] proposed a permutation-based two-way ANOVA to compare temporal profiles from different experimental groups. Luna and Li [2] proposed a statistical framework based on a shape-invariant 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 time-course 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 high-dimensional test statistics while preserving the complicated dependence structure among genes in microarray data analysis. In this paper, we present a permutation-based multiple-testing method for time-course microarray experiments when independent subjects contribute gene expression data at different time points. While the method can be generalized to broad class of goodness-of-fit test statistics for regression curves, for illustration we use the F-test 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 time-dependent in a single biological group and for identifying the genes whose time-profile 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 time-trends are considered to be of interest after removing the vertical shift. We shall refer to these as "time-course" and "time-trend" hypotheses, respectively. Note that if two separated curves can be overlapped by a vertical shift, then they have different time-courses, but the same time-trend. The test on a time-trend hypothesis will remove potential batch effects in microarray experiments.
The rest of the article is organized as follows. We first present a non-parametric test method to identify differential gene expression in a time-course 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 non-parametric test method for identifying genes whose expression levels are time-dependent in a single biological group and for identifying the genes for which the time-profile depends on the group among multiple groups. We approximate the null joint distribution of the test statistics using a permutation method.
Estimation of the Time-Course Profile
Here [W_{1}(t),..., W_{ p }(t)] is a pre-specified p-dimensional 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 B-spline 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.
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
Under H_{ j }, the constant is estimated as . Under , we obtain the estimate , where ( ) is estimated as described in the previous section.
for testing H_{ j }against . It is noted that for the permutation-based 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}.
Family-wise error rate (FWER) is defined by the probability of rejecting any null hypothesis H_{ j }when all m null hypotheses are true. A single-step 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 F-test statistics (f_{1},..., f_{ m }) from the original data.
- 2.
From the b-th permutation data (b = 1,..., B), compute the F-test statistics ( ).
- 3.
Single-step procedure to control the FWER
- (a)
From the b-th permutation data, calculate u_{ b }= max_{1≤j≤m} .
- (b)
For gene j, calculate the adjusted p-value 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:
- (b)
- (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 re-calculate 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 re-computed 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
Design and sample sizes for a K group case.
Time | ||||
---|---|---|---|---|
Group | t _{ 1 } | ⋯ | t _{ L } | Total |
1 | n _{11} | ⋯ | n _{1L} | n _{1·} |
⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
K | n _{K 1} | ⋯ | n _{ KL } | n _{K·} |
Total | n _{·1} | ⋯ | n _{·L} | n |
where μ_{ kj }(t) = β_{0, kj}+ .
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 group-free 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.
Illustration of a permutation for a K = 2 group case
(a) Original Data | ||||||
---|---|---|---|---|---|---|
Time | ||||||
t _{ 1 } | ⋯ | t _{ l } | ⋯ | t _{ L } | ||
Group | 1 | y _{11} | y _{l 1} | y _{L 1} | ||
y _{12} | y _{l 2} | y _{L 2} | ||||
y _{13} | ||||||
2 | y _{14} | y _{l 3} | y _{L 3} | |||
y _{15} | y _{l 4} | y _{L 4} | ||||
y _{l 5} | ||||||
(b) A permuted Sample | ||||||
Time | ||||||
t _{ 1 } | ⋯ | t _{ l } | ⋯ | t _{ L } | ||
Group | 1 | y _{14} | y _{l 2} | y _{L 3} | ||
y _{12} | y _{l 3} | y _{L 1} | ||||
y _{15} | ||||||
2 | y _{13} | y _{l 1} | y _{L 4} | |||
y _{11} | y _{l 4} | y _{L 2} | ||||
y _{l 5} |
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 two-group (i.e., K = 2) cases.
Simulation Settings
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(l-1)+j}= a_{10(l-1)+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(l-1)+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 single-group case, non-prognostic 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 Simulations 1 and 2, all m = 1, 000 genes are non-prognostic 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 non-prognostic.
Under each setting, N = 1, 000 simulation samples are generated and the single-step 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 }.
- 1.
Fit the time-course model under , and calculate n residuals, e_{ ij }= y_{ ij }- (t_{ i }).
- 2.
Fit the time-course 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
Empirical FWER level for a nominal two-sided FWER of 0.05
Two-group case | ||||||||
---|---|---|---|---|---|---|---|---|
One-group case | Simulation 1 | Simulation 2 | Simulation 3 | |||||
ρ | PERM | BOOT | PERM | BOOT | PERM | BOOT | PERM | BOOT |
0 | 0.060 | 0.046 | 0.057 | 0.067 | 0.042 | 0.056 | 0.062 | 0.458 |
0.3 | 0.048 | 0.041 | 0.049 | 0.050 | 0.057 | 0.065 | 0.050 | 0.438 |
0.6 | 0.047 | 0.037 | 0.051 | 0.047 | 0.051 | 0.047 | 0.045 | 0.474 |
Empirical global power at a two-sided FWER level of 0.05
One-group case | Two-group case | |||||||||
---|---|---|---|---|---|---|---|---|---|---|
Simulation 1 | Simulation 2 | Simulation 1 | Simulation 2 | Simulation 3 | ||||||
ρ | PERM | BOOT | PERM | BOOT | PERM | BOOT | PERM | BOOT | PERM | BOOT |
0 | 0.822 | 0.810 | 0.814 | 0.802 | 0.978 | 0.974 | 0.962 | 0.962 | 0.996 | 1.000 |
0.3 | 0.742 | 0.736 | 0.708 | 0.714 | 0.868 | 0.880 | 0.892 | 0.892 | 0.976 | 1.000 |
0.6 | 0.610 | 0.606 | 0.608 | 0.602 | 0.718 | 0.724 | 0.790 | 0.804 | 0.956 | 1.000 |
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 non-dauer 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 time-course and time-trend hypotheses. A time-course 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 time-trend 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 time-trend hypothesis, the testing procedures we have discussed in the methods section can be extended by simply subtracting off group-specific 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.
Discussion
We have considered two sets of hypotheses in the multi-group setting. For the time-course hypothesis, any difference among the groups, including parallel curves shifted vertically, would be considered interesting. For the time-trend, only cases where the time-trend 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 goodness-of-fit statistic based on natural splines used by Storey et al. [3]. Our method can be extended by using other regression goodness-of-fit 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 time-trajectory at each point is estimated by averaging the expressions. This statistic may be more sensible if the number of time-points is small.
A referee requested that we compare the power between the F-statistic based on the estimated time-trajectory at each point by averaging the observations with that based on the smoothed time-trajectory 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 time-trajectories was 0.9572 while that of the standard F-statistic is 0.8644. We observed similar comparisons under the wide range of simulation settings.
In the methods section, for the one-sample 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 FWER-adjusted P-values 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 permutation-based multiple testing method for time-course microarray experiments is computationally efficient and applicable for identifying the genes whose expression levels are time-dependent in a single biological group or for identifying the genes whose time-profiles are different among different groups.
Declarations
Authors’ Affiliations
References
- Park T, Yi S, Lee S, Lee SY, Yoo D-H, Ahn J-I, Lee Y-S: Statistical tests for identifying differentially expressed genes in time-course microarray experiments. Bioinformatics 2003, 19: 694–703. 10.1093/bioinformatics/btg068View ArticlePubMedGoogle Scholar
- Luan Y, Li H: Mode-based methods for identifying periodically regulated genes based on time course microarray gene expression data. Bioinformatics 2004, 20: 332–339. 10.1093/bioinformatics/btg413View ArticlePubMedGoogle Scholar
- 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.0504609102PubMed CentralView ArticlePubMedGoogle Scholar
- 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/1471-2105-8-S1-S10PubMed CentralView ArticlePubMedGoogle Scholar
- Hong F, Li H: Functional Hierarchical Models for Identifying Genes with Different Time-Course Expression Profiles. Biometrics 2006, 62: 534–544. 10.1111/j.1541-0420.2005.00505.xView ArticlePubMedGoogle Scholar
- Angelini C, De CD, Mutarelli M, Pensky M: A Bayesian Approach to Estimation and Testing in Time-course Microarray Experiments. Statistical Applications in Genetics and Molecular Biology 2007., 6(1):Google Scholar
- Wang J, Kim S: Global analysis of dauer gene expression in Caenorhabditis elegans. Development 2003, 130: 1621–1634. 10.1242/dev.00363View ArticlePubMedGoogle Scholar
- de Boor C: A Practical Guide to Splines. Springer-Verlag: New York; 2001.Google Scholar
- Westfall PH, Young SS: Resampling-based Multiple Testing: Examples and Methods for P-value Adjustment. Wiley: New York; 1993.Google Scholar
- Jung SH, Bang H, Young S: Sample size calculation for multiple testing in microarray data analysis. Biostatistics 2005, 6: 157–169. 10.1093/biostatistics/kxh026View ArticlePubMedGoogle Scholar
- 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.Google Scholar
- 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.091062498PubMed CentralView ArticlePubMedGoogle Scholar
- Storey JD: A direct approach to false discovery rates. Journal of the Royal Statistical Society, Series B 2002, 64(1):479–498.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.