New Research In
Physical Sciences
Social Sciences
Featured Portals
Articles by Topic
Biological Sciences
Featured Portals
Articles by Topic
 Agricultural Sciences
 Anthropology
 Applied Biological Sciences
 Biochemistry
 Biophysics and Computational Biology
 Cell Biology
 Developmental Biology
 Ecology
 Environmental Sciences
 Evolution
 Genetics
 Immunology and Inflammation
 Medical Sciences
 Microbiology
 Neuroscience
 Pharmacology
 Physiology
 Plant Biology
 Population Biology
 Psychological and Cognitive Sciences
 Sustainability Science
 Systems Biology
Significance analysis of time course microarray experiments

Contributed by Ronald W. Davis, June 12, 2005
Abstract
Characterizing the genomewide dynamic regulation of gene expression is important and will be of much interest in the future. However, there is currently no established method for identifying differentially expressed genes in a time course study. Here we propose a significance method for analyzing time course microarray studies that can be applied to the typical types of comparisons and sampling schemes. This method is applied to two studies on humans. In one study, genes are identified that show differential expression over time in response to in vivo endotoxin administration. By using our method, 7,409 genes are called significant at a 1% falsediscovery rate level, whereas several existing approaches fail to identify any genes. In another study, 417 genes are identified at a 10% falsediscovery rate level that show expression changing with age in the kidney cortex. Here it is also shown that as many as 47% of the genes change with age in a manner more complex than simple exponential growth or decay. The methodology proposed here has been implemented in the freely distributed and opensource edge software package.
The identification of genes that show changes in expression between varying biological conditions is a frequent goal in microarray experiments (1). Differential expression can be studied from a static or temporal viewpoint. In a static experiment, the arrays are obtained irrespective of time, capturing only a single moment of gene expression. In a temporal experiment the arrays are collected over a time course, allowing one to study the dynamic behavior of gene expression. A large amount of work has been done on the problem of identifying differentially expressed genes in static experiments (2). Because the regulation of gene expression is a dynamic process, it is also important to identify and characterize changes in gene expression over time. Here we present a general statistical method that identifies genes differentially expressed over time.
Several clustering methods have been applied to time course microarray data, including hierarchical clustering (3, 4), principal componentsbased clustering (5), Bayesian modelbased clustering (6), and Kmeans clustering of curves (7, 8). None of these clustering methods is directly applicable to identifying genes that show statistically significant changes in expression over time. The Kmeans clustering method has been modified to compare expression over time between two groups (9), but this method can only be applied to a few hundred genes at a time because of the computational cost of fitting a single model to all genes simultaneously (7, 8). This approach also requires that the statistical significance be calculated under the assumption that the clustering model estimated for one of the groups is true, which is nonstandard and potentially problematic.
The method that we propose draws on ideas from the extensive statistical literature on time course data analysis (10, 11), particularly splinebased methods (1216). It is applicable to detecting changes in expression over time within a single biological group and to detecting differences in the behavior of expression over time between two or more groups. Individuals may be sampled at multiple time points, or each time point may represent an independently sampled individual. The computational cost is not substantially greater than methods for static experiments, so there is no impeding limit on the number of genes that may be tested. For example, the method is applied here to microarray studies that each measures ≈40,000 genes.
The proposed method is applied to two recent studies carried out on humans. These studies encompass both types of sampling (longitudinal and independent; discussed below) and both types of differential expression over time (between groups and within a single group). In one study, gene expression was monitored over time in control individuals and endotoxintreated individuals. Endotoxin [lipopolysaccharide (LPS)] is widely used to study acute inflammatory and immune response; our goal was to understand the mechanism of endotoxin response by identifying genes with expression that is different over time between the treated and untreated groups (17). In a second study, we examined the effect of age on gene expression in the kidney, where samples were obtained from human subjects ranging in age from 27 to 92 years. The goal here was to identify genes whose expression changes significantly with respect to age in human kidney cortex tissue (18). Genes are identified in both studies that corroborate previous findings and provide insights into these problems.
Existing methods for detecting differential expression from static sampling designs (2) are appropriate for comparing unordered categorical conditions, such as three different cancer tumor types (19) or two different treatments (20). It can be shown that these are not generally applicable to time course experiments (see Supporting Appendix, which is published as supporting information on the PNAS web site). Typically these existing methods fail to properly use the temporal structure present in the data, either leading to a loss in power or incorrect calculation of significance. The method that we propose is specifically designed for time course experiments; expression over time is modeled flexibly, and statistical significance is calculated while accounting for sources of dependence over time.
Materials and Methods
Sample Preparations. Details on the protocols are described in refs. 17 and 18. For the endotoxin study, four subjects were administered endotoxin and four were administered a placebo. Blood samples were collected before infusion and at 2, 4, 6, 9, and 24 h afterward. The leukocytes were isolated from the blood samples by using a modified lysis protocol (17). Total RNA was extracted by using an RNeasy kit (Qiagen, Valencia, CA). Samples from hours 4 and 6 were unavailable for one of the controls. For the kidney aging study, samples were obtained from normal kidney tissues removed at nephrectomy or renal transplant biopsy from 74 patients ranging in age from 27 to 92 years. The samples were dissected into cortex and medulla sections based on histological evaluation. Each frozen tissue section was homogenized, and total cellular RNA was isolated according to TRIzol reagent protocol. Comprehensive evaluations were performed and reported elsewhere on various medical factors of the patients (18), and it is unlikely that these factors have confounded ageregulated changes in gene expression (18, 21). Only cortex samples (72 in total) were used in our analysis.
Microarray Analysis. Messenger RNA was amplified and hybridized onto human U133A and U133B GeneChips according to the protocols recommended by Affymetrix (Santa Clara, CA); 44,924 probe sets on the arrays were analyzed. Normalization was performed by using dchip, and expression levels were calculated by using the perfectmatchonly model (22). Expression values were then transformed by taking log_{2}(Data + 10), where the relatively negligible number 10 was added to stabilize the variance of values close to zero. The human subjects observed in the kidney aging study did not represent a purely random sample. We attempted to alleviate potential confounding from this sampling scheme by selecting 35,068 probe sets whose expression is well explained by the available sex variable, irrespective of any agedependent behavior (Supporting Appendix).
Statistical and Computational Details. Thorough details, including formula derivations, algorithms, and theoretical justifications can be found in Supporting Appendix.
Functional Analysis of Significant Genes. Probe sets on Affymetrix U133 GeneChips were mapped to gene IDs. Among the 4,163 probe sets significant at 0.1% for the endotoxin study, 2,914 unique genes were identified from 3,892 probe sets having mapped gene IDs. Among the 417 probe sets significant at 10% for the kidney study, 300 unique genes were identified from 364 probe sets with mapped gene IDs. In both studies, the Ingenuity Pathways Knowledge Base (Ingenuity Systems, Mountain View, CA) was used for functional analysis of genes. Briefly, the Ingenuity Pathways Knowledge Base consists of >10^{6} individually modeled relationships into an ontology of >550,000 biological concepts. Relationships between genes, proteins, small molecules, complexes, cells, processes, and diseases were manually extracted by Ph.D.level scientists from >200,000 peerreviewed articles (17).
Results
Experimental Objectives and Statistical Formulations. We developed a general statistical method that identifies genes showing temporal differential expression. This method was applied to two human studies encompassing both types of temporal differential expression.
In one study, kidney cortex samples were obtained from 72 human subjects ranging in age from 27 to 92 years. Only one array was obtained per sample, and the age of the subject was recorded (see Materials and Methods). Fig. 1a displays a simulated example of expression measured from a single gene for this type of study. The solid line is the population average time curve for the probe set, which is its true average expression over time with all sources of variation removed. The points are the observed expression values, one per each individual, and these can be thought of as independent random deviations from the solid line due to biological and measurement variation. “Independent sampling” was performed in this study because each sample of cortex tissue represents an independently sampled individual.
To determine whether each gene has expression that changes with age, our method involves performing a hypothesis test on each gene of whether its population average time curve is flat. We call this type of differential expression “withinclass temporal differential expression.” Fig. 1b shows the expression measurements from a highly significant gene, CRABP1, a cellular retinoic acidbinding protein. The gene was tested by first fitting a model under the null hypothesis that there is no differential expression, and then under the alternative hypothesis that there is differential expression. The null model is the dashed flat line that minimizes the sum of squares among all possible flat lines. The alternative model is the solid curve that minimizes the sum of squares among a general class of curves, namely natural cubic splines. A statistic is calculated that compares the goodness of fit of these two models. This statistic is a quantification of evidence for differential expression, and the larger it is the more differentially expressed the gene appears to be. For every gene, a statistic was calculated in this way, and a significance cutoff is applied to them by using a false discovery rate criterion (23). This process involves calculating the null distribution of the statistics when there is no differential expression and is accomplished through a data resampling technique.
In another study, eight human volunteers were randomly divided into endotoxintreated and control groups of equal size (see Materials and Methods). Fig. 2a is a simulated example of expression measurements from a single gene in a group of four individuals. The solid line is the population average time curve for this gene. The dashed lines are the average time curves for the individuals, meaning that these are the true underlying time curves for each individual with the sources of variation removed up to their individual variation. The deviation of an expression value from its corresponding “individual average time curve” can be thought of as an independent random event. The deviation of an individual average time curve from the population average can also be thought of as an independent random event. However, this implies that the deviations of expression values from the population average time curve are correlated within individuals. The sampling performed in this study is called “longitudinal sampling” because each individual is observed at more than one time point.
Here, the goal is to identify genes that show significant differential expression across time between the endotoxintreated and control groups. We call this type of differential expression “betweenclass temporal differential expression.” The approach used here is conceptually similar to that in the kidney study. Fig. 2b shows the expression measurements from a significant gene in the study, IFN regulatory factor 1. Under the null hypothesis of no differential expression, the treated and control groups have the same population average time curve. Therefore, a single curve (a natural cubic spline) is fit to the combined groups, which is represented by the dashed curve in Fig. 2b . The alternative model is formed by fitting a separate curve to each group, as is shown by the solid lines in Fig. 2b . A statistic is computed based on the improvement in goodness of fit in going from a single curve to the separate curves for each group. As before, this statistic is a quantification of evidence for differential expression, and the larger it is the more differentially expressed the gene appears to be. A significance cutoff is applied to these statistics in the same fashion as in the kidney aging study.
In contrast to a static experiment, it is more difficult in the time course setting to form statistics that accurately quantify differential expression. Determining the distribution of the statistics when there is no differential expression is also more challenging. The behavior of expression over time may vary greatly from gene to gene, so a flexible modeling approach must be taken in forming statistics. Under certain study designs (e.g., the endotoxin study) there may be dependence between the expression measurements within a single individual, which complicates the formation of statistics and the simulation of the null distribution. Finally, simple permutation methods cannot be used to simulate null statistics because of the complex structure of time course measurements.
Proposed Statistical Method. Methodology was developed to address these issues in a statistically rigorous fashion. In doing so, a general model for gene expression over time within a single biological group was first formulated. Although a single model can be applied to both studies (Supporting Appendix), a simplified version is possible for the kidney study because of its independent sampling scheme. Let y_{ij} be the relative expression level of gene i in individual j, where there are i = 1, 2,..., 35,068 probe sets and j = 1, 2,..., 72 individuals. Individual j is observed at age t_{j} , which lies somewhere between 27 and 92 years. The expression values are modeled by
where μ _{i} (t_{j} ) is the population average time curve for gene i evaluated at time t_{j} and where ε _{ij} is the random deviation from this curve. In terms of Fig. 1a , μ _{i} (t) is shown by the solid line. The distance between this curve and an observed expression value is ε_{ij}. The ε _{ij} are assumed to be independent random variables with mean zero and genedependent variance .
The following model was developed for the endotoxin study and for longitudinal sampling in general. Let y_{ijk} be the relative expression level of gene i on individual j at the kth time point, where there are i = 1, 2,..., 44,924 probe sets and j = 1, 2,..., 8 different individuals sampled (four within each group). For each individual, there were k = 1, 2,..., 6 time points observed at times t_{jk} , except for the one control, who is missing two time points. The expression values are modeled by
The population average time curve for gene i is again μ _{i} (t). Individual j deviates from μ _{i} (t) by γ _{ij} , implying that μ _{i} (t) + γ _{ij} is the individual average time curve for individual j. The measurement error and remaining sources of random variation are modeled by ε _{ijk} . The solid line in Fig. 2a is represented by μ _{i} (t), and the dashed lines are represented by μ _{i} (t) + γ _{ij} . Each expression value deviates from its corresponding dashed line by ε _{ijk} . The γ _{ij} and ε _{ijk} are assumed to be independent random variables with means equal to zero and genedependent variances and , respectively. The case for which γ _{ij} is modeled as a curve will be dealt with elsewhere; however, the endotoxin study did not contain enough observations to permit this extra level of complexity.
The population average time curve μ _{i} (t) is parameterized in terms of an intercept plus a pdimensional linear basis:
where s(t) = [s _{1}(t), s _{2}(t),..., s _{p}(t)] ^{T} is a prespecified pdimensional basis, α _{i} is the unknown genespecific intercept, and β _{i} = [β _{i} _{1}, β _{i} _{2},..., β _{ip} ] ^{T} is a pdimensional vector of unknown genespecific parameters. A natural choice for the basis is a polynomial of degree p, which would result in modeling μ _{i} (t) = α _{i} + β _{i} _{1} t + β _{i} _{2} t ^{2} + ··· + β _{ip}t^{p} . The polynomial basis was effective in both studies, although a natural cubic spline basis is more flexible and makes fewer assumptions (1216). Furthermore, a natural cubic spline parameterization of μ _{i} (t) resulted in increased power in both studies. For either type of basis, the curve was estimated by minimizing the sum of squares between the curve and the observed expression values, which involves assigning values to only α _{i} and β _{i} . The model fitting procedure for longitudinal sampling is slightly more complicated in that it also takes into account the dependence of measurements within an individual (Supporting Appendix). We also developed a method to automatically choose a single value of the basis dimension p for all genes based on a singular value decomposition of the expression data (Supporting Appendix). We found that allowing p to vary from gene to gene or choosing a single p in a post hoc fashion lead to anticonservative inflation of significance due to overfitting.
Not only does the basis representation of μ _{i} (t) facilitate model fitting, but it also greatly simplifies the hypothesis tests for differential expression. The tests can now be written in terms of the α _{i} and β _{i} for each gene, which implies that the tests do not depend on specific time points, so general sampling schemes may be analyzed as we have claimed. In the kidney aging study, the null hypothesis of no differential expression is equivalent to restricting μ _{i} (t) to be a constant, and the alternative hypothesis of differential expression allows μ _{i} (t) to be a curve. Thus, the null hypothesis model is fit under the constraint that μ _{i} (t) = α _{i} and β _{i} = 0, and the alternative hypothesis model is fit under the general parametrization of μ _{i} (t). For the endotoxin study, the null hypothesis is that the treated and control groups have equal μ _{i} (t) (i.e., equal α _{i} and β _{i} ), and the alternative hypothesis is that they are not equal. The null hypothesis model is obtained by fitting a curve to the two groups combined, and the alternative hypothesis model is obtained by fitting a separate curve to each group. In this particular study, we were not interested in a difference in the intercepts α _{i} because all individuals started out as untreated at time 0. Therefore, the intercept was implicitly assumed to be equal between the two groups under both hypotheses, which comes down to a test for equality of the β _{i} between the two groups. This equality will not always be the case when testing for betweenclass temporal differential expression, so we have developed model fitting methods for both scenarios (Supporting Appendix).
A statistic for each gene was then formed that quantifies differential expression. The statistic was defined to be an analogue of the t and F statistics that are commonly used in static differential expression methods. The statistic compares the goodness of fit of the model under the null hypothesis to that under the alternative hypothesis. First, fitted values from the null and alternative models are calculated for each observed value. The residuals of the model fits are then obtained by subtracting the fitted from the observed values. Calculating to be the sum of squares of the residuals obtained from the null model, and from the alternative model, the statistic for gene i was constructed as
This statistic is proportional to the typical F statistic used in comparing two nested models. The intuition behind the formula is that quantifies the increase in goodness of fit, and dividing this difference by provides exchangeability of the F _{i} across the genes. Justification and exact formulas for all cases can be found in Supporting Appendix.
The null distribution of these statistics is calculated through a data resampling method called the bootstrap (24), where residuals from the alternative model are resampled and added back to the null model to simulate the case where there is no differential expression (Supporting Appendix). We note that our resampling scheme takes into account dependence between time points when the study involves longitudinal sampling. The observed statistics and null statistics were used to estimate a Q value for each gene (Supporting Appendix), which estimates the falsediscovery rate incurred when calling the gene significant (23, 25).
Analysis of Systemic Inflammatory Response Induced by LPS. In the endotoxin study there are 4,163, 7,409, and 12,720 significant probe sets at respective Qvalue cutoffs of 0.1%, 1%, and 5%. These results indicate that an endotoxin injection causes a profound gene expression response in blood leukocytes, which is expected for a number of reasons. First, in vivo administration of LPS invokes an acute systemic inflammation that dramatically perturbs the body's physiology. Second, some differential expression may be due to changing distributions of cell populations in the blood in addition to transcriptional changes (17). Also, measuring differential expression over time is a more sensitive study design than the typical static design; any change over time qualifies as differential expression. A t test (23) and the sam software (20, 23) were used to test for differential expression between 0 and 2 h and between 0 and 4 h. The significance results are displayed in Table 1, where it can be seen that an actual time course analysis offers a sizeable increase in statistical power over a static design analysis.
Under the assumption that the control individuals have flat expression over time, it would have sufficed to simply test the endotoxin individuals for withinclass temporal differential expression. However, we found that, among the controls, at least 6% of the probe sets show withinclass temporal differential expression. We also tested the endotoxintreated individuals for withinclass temporal differential expression and found the significance to be inflated to unrealistic levels. Therefore, it is necessary to include both treated and control individuals in a test for betweenclass temporal differential expression to identify genes differentially expressed because of LPS injection.
To get a broad picture of the behavior of differentially expressed genes, a singular value decomposition was performed on the 4,163 most significant probe sets. Fig. 3 shows the top two eigengenes that explain 66% and 16% of the variance, respectively. The relevant information is the shape of each eigengene, not the magnitude or direction of differential expression. Among the 4,163 most significant probe sets, 27% are upregulated as the eigengene is drawn; 73% are downregulated, which is simply the reflection of the eigengene's shape across the time axis. The second eigengene shows more complex behavior, where changes in expression occur in both directions over the time course. We identified 2,914 unique genes among the 4,136 probe sets; 756 showed upregulation during the time course, and 2,158 showed downregulation. Global functional analysis was performed by using a gene ontology built from experimental evidence compiled in the Ingenuity Pathways Knowledge Base (Fig. 4 and Data Set 1, which is published as supporting information on the PNAS web site).
The most significant functional groups from the 756 upregulated genes take part in immune response (121 genes), inflammatory disease (41 genes), cellular movement (84 genes), tissue morphology (72 genes), and cell death (144 genes). These groups are consistent with an intense response of leukocytes to LPS. The apparent expression of many elements of inflammation are upregulated, including secretive cytokines, chemokines, and associated proteins (CCL20, CCL3, CCL4, CXCL16, CXCL2, IL1A, IL1B, IL8, TNF, TNFSF13B, IL1RAP, IL1RN, TNFAIP3, and TNFAIP6); their membrane receptors (CCR1, IL10RB, IL18R1, IL18RAP, IL1R1, IL1R2, IL4R, IL8RA, IL9R, and TNFRSF6), Tolllike receptors (TLR4, TLR5, TLR6, and TLR8), Fc receptors (FCAR, FCER1G, FCGR1A, FCGR2A, and FCGR2B), IFN receptors (IFNGR1, IFNGR2, and IFNAR1), protein tyrosine phosphatases (PTP4A1, PTPN2, PTPN22, PTPNS1, PTPRC, PTPRJ, PTPRN2, and PTPRO), Janus tyrosine kinases/signal transducers and activators of transcription (JAK2, JAK3, STAT2, STAT3, STAT5A, and STAT5B), and NFκB/IKB proteins (NFκB2, NFκBIA, and NFκBIE). These and other elevated genes participate in the activation of innate immune response and many functions of leukocytes, including cell movement of leukocytes, infiltration, migration, phagocytosis, activation, chemotaxis, proliferation, and recruitment. The most significantly upregulated gene is ORM1, a key acute phase plasma protein.
The most significant functional groups from the genes showing downregulation are protein synthesis (83 genes), including genes of elongation initiation factors, ribosomal proteins, mitochondrial ribosomal protein; and RNA posttranscriptional modification (60 genes), including genes such as heterogeneous nuclear ribonucleoproteins, small nuclear ribonucleoproteins, and splicing factors. Members of RNA polymerase II are also decreased. Correspondingly, genes involved in oxidative phosphorylation, such as mitochondrial complexes I (NDUF), II (SDH), III (UQCR), and IV (COX) are downregulated. This concerted suppression of the cell's protein syntheses, transcription programs, and energy productions may reflect a generalized stress response in blood leukocytes. Interestingly, expression levels of members associated with the MHC class II are also suppressed, consistent with the reduced antigenpresenting capability after LPS shock.
Taking these results together, we observed a comprehensive transcriptional response to LPS stimulation, where leukocyte cells reallocate resources and upregulate transcription of genes involved in innate immune and defense mechanisms.
Analysis of AgeRelated Genes in the Kidney Cortex. We detected substantial differential expression in this study, although not nearly as much as in the previous study. At Qvalue cutoffs of 5% and 10%, there are 187 and 436 significant probe sets, respectively. Among the 436 probe sets significant at Q ≤ 10%, 320 unique genes were identified, with 240 show increasing expression with age and 80 show decreasing expression.
A functional analysis was also performed on these ageregulated genes (Data Set 2, which is published as supporting information on the PNAS web site). The most significant functional group is immune response (69 genes), which includes genes associated with the complement component (C1QA, C1QB, C1QG, C1QR1, C1R, C1S, and C4A), cytokines, chemokines, and receptors (CCL2, CCL20, CCR1, CX3CR1, CXCL9, CXCR4, IL10RB, CSF1R, and CSF2RB), and the MHC II complex (FCER1G, HLADMA, HLADPA1, HLADQB1, HLADRB1, and HLADRB3). All but three of these genes show upregulated expression with age, which suggests an increased abundance of immune response genes in the kidney cortex, perhaps as a result of a low level of inflammation over life or as a compensation for the decreased immune function in older age. Renal infiltration with immunocompetent cells, which is known to correlate with renal diseases and increase over age, can also account for some of these observed changes.
The aging process has previously been shown to involve a widespread and complicated alteration in gene expression (18). In our analysis, significant genes were identified that have a variety of celltocell signaling and interaction (68 genes). For example, 34 genes are involved in the quantity or mobilization of calcium and 14 genes in tyrosine phosphorylation. Fiftynine genes were identified to be involved in apoptosis and cell death. The expression levels of apoptosisenhancer proteins MYC (CMYC), SP1 (Sp1 transcription factor), caspase 1, and ubiquitin D increase with age in the kidney, whereas the antiapoptosis protein HSPA9B (mortalin2) decreases with age. However, the expression of BCL2A1, an inhibitor of the intrinsic apoptosis pathway, is observed to be increasing with age (26).
Gene expression may change with age in response to declining kidney function and increased susceptibility to kidney disease. Interestingly, a group of 12 genes known to be localized in mitochondria (HSPA9B, COX8A, COX7C, AKAP1, BCKDHA, CLPP, FDX1, AMT, DBT, ATP5G3, AK2, and NME4) are identified as significant, and their expression level uniformly declines in older age. A number of ageregulated genes are also known to be involved in renal diseases, including C1QA, CCR1, LYN, PTPRC, and TNFSF13B, which all show increasing expression with age (2731). The expression of TRPM6 (a transient receptor potential cation channel) negatively correlates with age, and mutations of this gene cause hypomagnesemia with secondary hypocalcemia (32).
We used a 4dimensional basis in the significance analysis, although other choices of p are possible. A 1dimensional basis is equivalent to identifying differential expression by fitting a linear regression and testing whether the slope is zero. Applied to logtransformed expression data, a linear regression is only able to identify genes that show simple exponential behavior in their differential expression. By testing the null hypothesis of linear or flat logtransformed expression, we found that at least 25% of the 35,068 probe sets are differentially expressed in a manner more complex than an exponential function of age. Among the 417 most significant probe sets when using the 1dimensional basis, it is estimated that at least 47% of these are more complex than simple exponential behavior. At Q ≤ 5%, just over half of the genes found to be significant using the 4dimensional basis are also significant when using the 1dimensional basis. Therefore, there is potentially a functional class of genes showing simple exponential temporal differential expression and another functional class showing substantially more complex differential expression.
Discussion
We have proposed a significance method to identify differentially expressed genes in time course microarray experiments and applied it to studies involving two types of sampling and both types of temporal differential expression. The method may also be applied to more complicated situations, where three or more groups are compared, for example. We have justified the method in terms of well established statistical concepts (Supporting Appendix).
Temporally differentially expressed genes were identified in two time course microarray studies on humans. The human in vivo model of bacterial endotoxin represents a unique opportunity to examine the onset of systemic inflammation in blood leukocytes. Our analysis suggests that a genomewide transcriptional response takes place where blood leukocytes reallocate resources and upregulate transcription of genes involved in innate immune and defense mechanisms. The inclusion of controls in this study played an important role, in which the controls themselves showed differential expression over time. In the clinical setting, comparing a treatment group to time 0 does not provide a proper control over the entire time course. The significance analysis of the kidney aging study indicates that a large proportion of genes have expression that increases with age, including those involved in signal transduction, cell growth, and genome stability. Genes related to metabolism were found to be decreasing with age. It was also shown that a large proportion of gene expression changing with age does so in a manner more complex than exponential growth. Therefore, in general it is important to use flexible models when identifying genes that are temporally differentially expressed.
We have developed a freely available pointandclick software package called edge that includes this methodology and new methodology for static experiments. edge is available from J.D.S. upon request.
Acknowledgments
We thank Stuart Kim, Steve Calvano, and Stephen Lowry for generously sharing data and James Dai for helping with the calculations used for Table 1. This research was supported by National Institutes of Health Grant R01 HG00291301 (to J.D.S.) and National Institutes of Health LargeScale Collaborative Research Grant U54 GM211903¶ (principle investigator, R.G.T.).
Footnotes

↵ † To whom correspondence may be addressed. Email: jstorey{at}u.washington.edu or dbowe{at}stanford.edu.

Author contributions: J.D.S., W.X., R.G.T., and R.W.D. designed research; J.D.S., W.X., and J.T.L. performed research; J.D.S., W.X., and J.T.L. contributed new reagents/analytic tools; J.D.S., W.X., and J.T.L. analyzed data; and J.D.S. and W.X. wrote the paper.

Abbreviation: LPS, lipopolysaccharide.

↵ ¶ Participating investigators: H. V. Baker, P. E. Bankey, T. R. Billiar, B. H. Brownstein, S. E. Calvano, I. H. Chaudry, J. P. Cobb, C. Cooper, B. Freeman, R. L. Gamelli, N. S. Gibran, B. G. Harbrecht, W. Heagy, D. M. Heimbach, D. N. Herndon, J. W. Horton, J. L. Hunt, J. Johnson, J. A. Lederer, T. Logvinenko, S. F. Lowry, J. A. Mannick, B. A. McKinley, C. MillerGraziano, M. N. Mindrinos, J. P. Minei, L. L. Moldawer, E. E. Moore, F. A. Moore, A. B. Nathens, G. E. O'Keefe, L. G. Rahme, D. G. Remick, Jr., M. B. Shapiro, R. L. Sheridan, G. M. Silver, R. D. Smith, S. Somers, G. Stephanopoulos, M. Toner, H. S. Warren, M. A. West, S. E. Wolf, and V. R. Young.
 Copyright © 2005, The National Academy of Sciences
References

↵
Slonim, D. K. (2002) Nat. Genet. 32 , 502508. pmid:12454645
 ↵

↵
Eisen, M. B., Spellman, P. T., Brown, P. O. & Botstein, D. (1998) Proc. Natl. Acad. Sci. USA 95 , 1486314868. pmid:9843981
 ↵

↵
Alter, O., Brown, P. O. & Botstein, D. (2000) Proc. Natl. Acad. Sci. USA 97 , 1010110106. pmid:10963673

↵
Wakefield, J., Zhou, C. & Self, S. (2003) Bayesian Stat. 7 , 721732.

↵
BarJoseph, Z., Gerber, G., Jaakkola, T. S., Gifford, D. K. & Simon, I. (2003) J. Comp. Bio. 10 , 341356.

↵
Luan, Y. & Li, H. (2003) Bioinformatics 19 , 474482. pmid:12611802

↵
BarJoseph, Z., Gerber, G., Simon, I., Gifford, D. K. & Jaakkola, T. S. (2003) Proc. Natl. Acad. Sci. USA 100 , 1014610151. pmid:12934016

↵
Diggle, P., Heagerty, P., Liang, K. Y. & Zeger, S. (2002) Analysis of Longitudinal Data (Oxford Univ. Press, Oxford), 2nd Ed.

↵
Ramsay, J. O. & Silverman, B. W. (1997) Functional Data Analysis (Springer, Berlin).
 ↵
 ↵

↵
Calvano, S. E., Xiao, W., Richards, D. R., Felciano, R. M., Baker, H. V., Cho, R. J., Chen, R. O., Brownstein, B. H., Cobb, J. P., Tschoeke, S. K., et al. (2005) Nature, in press.
 ↵
 ↵

↵
Tusher, V., Tibshirani, R. & Chu, C. (2001) Proc. Natl. Acad. Sci. USA 98 , 51165121. pmid:11309499

↵
Higgins, J., Wang, L., Kambham, N., Montgomery, K., Maxon, V., Vogelmann, S., Lemley, K., Brown, P., Brooks, J. & van de Rijn, M. (2004) Mol. Biol. Cell 15 , 649656. pmid:14657249

↵
Li, C. & Wong, W. H. (2001) Proc. Natl. Acad. Sci. USA 98 , 3136. pmid:11134512

↵
Storey, J. D. & Tibshirani, R. (2003) Proc. Natl. Acad. Sci. USA 100 , 94409445. pmid:12883005

↵
Efron, B. & Tibshirani, R. J. (1993) An Introduction to the Bootstrap (Chapman & Hall, Boca Raton, FL).
 ↵
 ↵
 ↵

Mitchell, D., Pickering, M., Warren, J., FossatiJimack, L., CortesHernandez, J., Cook, H., Botto, M. & Walport, M. (2002) J. Immunol. 168 , 25382543. pmid:11859149
 ↵
 ↵