Gene expression profiles give a snapshot of mRNA concentration levels as encoded by the genes of an organism under given experimental conditions. Early studies of this data often focused on a single point in time which biologists assumed to be critical along the gene regulation process after the perturbation. However, the static nature of such experiments severely restricts the inferences that can be made about the underlying dynamical system.
With the decreasing cost of gene expression microarrays time series experiments have become commonplace giving a far broader picture of the gene regulation process. Such time series are often irregularly sampled and may involve differing numbers of replicates at each time point [1
]. The experimental conditions under which gene expression measurements are taken cannot be perfectly controlled leading the signals of interest to be corrupted by noise, either of biological origin or arising through the measurement process.
Primary analysis of gene expression profiles is often dominated by methods targeted at static
experiments, i.e. gene expression measured on a single time-point, that treat time as an additional experimental factor [1
]. However, were possible, it would seem sensible to consider methods that can account for the special nature of time course data. Such methods can take advantage of the particular statistical constraints that are imposed on data that is naturally ordered [7
The analysis of gene expression microarray time-series has been a stepping stone to important problems in systems biology such as the genome-wide identification of direct targets of transcription factors [13
] and the full reconstruction of gene regulatory networks [15
]. A more comprehensive review on the motivations and methods of analysis of time-course gene expression data can be found in [17
Testing for Expression
A primary stage of analysis is to characterize the activity of each gene in an experiment. Removing inactive or quiet genes (genes which show negligible changes in mRNA concentration levels in response to treatments/perturbations) allows the focus to dwell on genes that have responded to treatment. We can consider two experimental set ups. Firstly, we may be attempting to measure the absolute level of gene expression (for example using Affymetrix GeneChip microarrays). In this case a quiet gene would be one whose expression level is indistinguishable from noise. Alternatively, we might be may be hybridizing two samples to the same array and quantifying the ratio of the expression levels. Here a quiet gene would be one which is showing a similar response in both hybridized samples. In either case we consider such expression profiles will consist principally of noise. Removing such genes will often have benign effects later in the processing pipeline. However, mistaken removal of profiles can clearly compromise any further downstream analysis. If the temporal nature of the data is ignored, our ability to detect such phenomena can be severely compromised. An example can be seen in Figure , where the temporal information is removed from an experimental profile by randomly reordering its expression samples. Disregarding the temporal correlation between measurements, hinders our ability to assess the profile due to critical inherent traits of the signal being lost such as the speed and scale of variation.
Figure 1 Temporal information removed from the profile of gene Cyp1b1 in the experimental mouse data. (a) The centred profile of the gene Cyp1b1 (probeID 1416612_at in the GSE10562 dataset). The blue crosses represent zero-mean hybridised gene expression in time (more ...)
Failure to capture the signal in a profile, irrespective of the amount of embedded noise, may be partially due to temporal aggregation
effects, meaning that the coarse sampling of gene expression or the sampling rates do not match the natural rates of change in mRNA concentrations [18
]. For these reasons, the classification scheme of differential expression in this paper is focused on reaching a high true positive rate
) and is to serve as a pre-processing tool prior to more involved analysis of time-course microarray data. In this work we distinguish between two-sample
testing and experiments where control
cases are directly-hybridized on the microarray (For brevity, we shall refer to experiments with such setups as one-sample testing
). The two-sample
setup is a common experimental setup in which two groups of sample replicates are used [13
]; one being under the treatment effect of interest and the other being the control group, so to recover the most active genes under a treatment one may be interested in testing for the statistical significance of a treated profile being differentially expressed with respect to its control counterpart. Other studies use data from a one-sample
], in which the control
cases are directly hybridized on a microarray and the measurements are normalized log fold-changes between the two output channels of the microarray [20
], so the analogous goal is to test for the statistical significance of having a non-zero signal.
A recent significant contribution in regards to the estimation and ranking of differential expression of time-series in a one-sample
setup is a hierarchical Bayesian model for the analysis of gene expression time-series (BATS) [11
] which offers fast computations through exact equations of Bayesian inference, but makes a considerable number of prior biological assumptions to accomplish this (cf. Simulated data
Gene Expression Analysis with Gaussian Processes
] offer an easy to implement approach to quantifying the true signal and noise embedded in a gene expression time-series, and thus allow us to rank the differential expression of the gene profile. A Gaussian process is the natural generalisation of a multivariate Gaussian distribution to a Gaussian distribution over a specific family of functions
-- a family defined by a covariance function
, i.e. a metric of similarity between data-points (Roughly speaking, if we also view a function as a vector with an infinite number of components, then that function can be represented as a point in an infinite-dimensional space of a specific family of functions
and a Gaussian process as an infinite-dimensional Gaussian distribution over that space).
In the context of expression trajectory estimation, a Gaussian process coupled with the squared-exponential
covariance function (or radial basis function
, RBF) -- a standard covariance function used in regression tasks -- makes the reasonable assumption that the underlying true signal in a profile is a smooth
], i.e. a function with an infinite degree of differentiability. This property endows the GP with a large degree of flexibility in capturing the underlying signals without imposing strong modeling assumptions (e.g. number of basis functions) but may also erroneously pick up spurious patterns (false positives) should the time-course profiles suffer from temporal aggregation. From a generative viewpoint, the profiles are assumed to have been corrupted by additive white Gaussian noise. This property makes the GP an attractive tool for bootstrapping simulated biological replicates [24
In a different context, Gaussian process priors have been used for modeling transcriptional regulation. For example in [25
], while using the time-course expression of a-priori known direct targets (genes) of a transcription-factor, the authors went one step further and inferred the concentration rates of the transcription-factor protein itself and [26
] extended the same model for the case of regulatory repression. The ever-lingering issue of outliers in time series is still critical, but is not addressed here as there is significant literature on this issue in the context of GP regression, which is complementary to this work.
For example [19
] developed a probabilistic model using Gaussian processes with a robust noise model specialised for two-sample testing to detect intervals
of differential expression, whereas the present work optionally focuses on one-sample
testing, to rank the differential expression and ultimately detect quiet/active
genes. Other examples can also be easily applied; [28
] use a Student-t
distribution as the robust noise model in the regression framework along with variational approximations to make inference tractable, and [29
] employ a Student-t
observation model with Laplace approximations for inference. The standard GP regression framework is straightforward to use here with a limited need for manual tweaking of a few hyper-parameters. We describe the GP framework, as used here for regression, in more detail in the Methods