Frequency-based time-series gene expression recomposition using PRIISM

Abstract

Background

Circadian rhythm pathways influence the expression patterns of as much as 31% of the Arabidopsis genome through complicated interaction pathways, and have been found to be significantly disrupted by biotic and abiotic stress treatments, complicating treatment-response gene discovery methods due to clock pattern mismatches in the fold change-based statistics. The PRIISM (Pattern Recomposition for the Isolation of Independent Signals in Microarray data) algorithm outlined in this paper is designed to separate pattern changes induced by different forces, including treatment-response pathways and circadian clock rhythm disruptions.

Results

Using the Fourier transform, high-resolution time-series microarray data is projected to the frequency domain. By identifying the clock frequency range from the core circadian clock genes, we separate the frequency spectrum to different sections containing treatment-frequency (representing up- or down-regulation by an adaptive treatment response), clock-frequency (representing the circadian clock-disruption response) and noise-frequency components. Then, we project the components’ spectra back to the expression domain to reconstruct isolated, independent gene expression patterns representing the effects of the different influences.

By applying PRIISM on a high-resolution time-series Arabidopsis microarray dataset under a cold treatment, we systematically evaluated our method using maximum fold change and principal component analyses. The results of this study showed that the ranked treatment-frequency fold change results produce fewer false positives than the original methodology, and the 26-hour timepoint in our dataset was the best statistic for distinguishing the most known cold-response genes. In addition, six novel cold-response genes were discovered. PRIISM also provides gene expression data which represents only circadian clock influences, and may be useful for circadian clock studies.

Conclusion

PRIISM is a novel approach for overcoming the problem of circadian disruptions from stress treatments on plants. PRIISM can be integrated with any existing analysis approach on gene expression data to separate circadian-influenced changes in gene expression, and it can be extended to apply to any organism with regular oscillations in gene expression patterns across a large portion of the genome.

Background

Differential gene expression studies typically use the fold change statistic (the ratio of mRNA quantities between two samples) as input, and have been used to discover genes involved in adaptive stress responses which have not been previously characterized (i.e., “novel genes”) [1]. Specifically, to correct for changes in gene expression induced by non-treatment related influences, fold-change values for time-series data are usually calculated using treatment and control data at every timepoint [1]. One of the major factors causing gene oscillations under control conditions is the molecular circadian clock, which influences physiology and metabolism in preparation for predictable changes in light and temperature [2]. However, a wide range of biotic and abiotic stress treatments have been shown to disrupt rhythmic clock patterns through amplitude changes or phase shifts [3–8], resulting in significant fold changes for genes which are clock-influenced but are not involved in direct stress response. Figure1 demonstrates that genes can be differentially regulated due to direct stress responses (I), indirectly differentially regulated through disruption of clock pathways induced by the stress (II) or a combination of both (III). Additional complications in regulation patterns arise from the complexity of transcription factor pathways, in which targets may be regulated by clock components directly or through interactions with their transcription factors (Figure1). For this reason, novel treatment-response gene discovery methods are complicated by the disruption of synchronization of the circadian rhythm pathways, but this complexity is not reflected in existing methods including fold change studies, clustering analysis approaches, and more complex time-serial-based algorithms [1, 2, 5, 6, 9–17].

Figure 1

Biotic and abiotic stresses both directly and indirectly influence target gene expression patterns. Genes found to be differentially expressed may be influenced by (I) only direct treatment influences, (II) only indirect circadian-clock disruption influences, or (III) both direct treatment response and indirect clock influences.

In this paper, we present the PRIISM (Pattern Recomposition for the Isolation of Independent Signals in Microarray data) algorithm to perform novel stress-response gene discovery analyses which correct for differential gene expression patterns induced by the circadian clock. As described previously [6], although core circadian clock gene patterns undergo significant changes in phase and amplitude as a result of stress, they maintain oscillating frequencies which remain similar to each other, and still remain close to the circadian pattern of one cycle per day. It has also been shown that stress results in significantly increased average expression levels for stress-response genes [6], which are reflected in the low-frequency signals (where one oscillation cycle occurs over the course of several days) for these genes. We assume that although circadian clock influences and adaptive stress-response influences can interact with each other (Figure1), they still cycle at very different rates from each other (and therefore maintain separate dominant frequency ranges) under stress conditions. Based on these observations, we have developed PRIISM to project gene expression data to the frequency domain using the Fourier Transform, isolate independent signals, and then project them back to the expression domain to reconstruct independent gene expression patterns representing the effects of different genetic influences. PRIISM is capable of separating one gene expression pattern into three distinct gene expression patterns: (1) The treatment-frequency gene expression pattern, which has much of the complicating circadian influences removed, and consequently can be used to more accurately identify differentially regulated genes which are involved in direct treatment response, (2) the clock-frequency gene expression pattern, representing rhythmic patterns with a period of approximately one cycle per day, and (3) the noise-frequency gene expression pattern (Figure2). By applying PRIISM on a cold-treatment dataset, we demonstrate that it can identify known treatment-response genes with a much lower false-positive rate than the existing methods, and can also identify important regulatory timepoints which are not obvious in the unprocessed data. In addition to improving performance when conducting novel treatment-response gene discovery, PRIISM also provides gene expression data which represent only circadian clock influences, and may be useful for circadian clock analysis studies.

Figure 2

PRIISM separates gene expression data into three independent gene expression datasets. PRIISM separates (A) the original gene expression patterns under control and treatment conditions (used to calculate the fold change pattern) into (B) treatment-frequency, clock-frequency and noise-frequency gene expression patterns. The cold-induced gene COR15A (AT2G42540) is shown as an example.

Biological approaches such as the use of constant light and clock component genetic knockout mutants are applied in order to attempt to remove the influences of the circadian clock on target gene expression. However, constant light is an unnatural condition which reduces the applicability of the results, because natural biotic and abiotic genetic stress-response patterns depend on the time-of-day (the point in the light/dark cycle) at which the treatment is applied [6, 18, 19]. Likewise, the use of genetic knockout mutants of circadian clock genes can reduce disruptions due to circadian input; However, since stress response genes may be regulated by clock components, the results of such a study are also difficult to interpret [7, 19, 20].

Most existing computational approaches for studying differential gene expression in microarray datasets involve clustering algorithms designed to group genes with similar expression profiles, with the goal of identifying potential annotations for unknown genes [10–17]. However, the gene distance measures used by all of these clustering methods are unable to distinguish adaptive-response gene expression patterns from circadian clock disruption gene expression patterns, and so may cluster genes with similar clock influences but very different treatment-response influences. Bar-Joseph et al’s (2003) continuous representation model for finding differentially expressed genes in time series micro array datasets (which has been used to find more cell-cycle response genes in yeast than conventional clustering methods) is also unable to filter clock influences from treatment response influences on gene expression patterns [21].

Several studies have shown that between 6% and 31% of the Arabidopsis genome is influenced by circadian clock genetic components [5, 22, 23]; while another study suggests that there are significant baseline circadian oscillations for nearly 100% of the genome [24]. A number of approaches have been developed for analyzing the circadian rhythms of genes in time-series datasets [5, 25–28]. Fourier analysis (which can be used to identify dominant frequencies in time-series data) has been applied to successfully identify periodic genes by treating time-series microarray datasets as time-domain signals [28–32]. However, these Fourier analysis methods have not been widely used in differential gene expression study methods, because 1) in existing Fourier analysis applications [28–32], a fixed frequency range was used as a priori knowledge to discover genes with similar oscillations, but novel genes may have totally different frequency patterns under different treatment conditions and; 2) to accurately capture oscillating rhythms, high resolution time course gene expression data is essential according to Nyquist sampling theorem[25, 33], but such data have not been available until recently.

As the price of running microarrays and RNA-seq chips continues to fall, high-resolution time-series gene expression datasets that contain enough information to identify and characterize circadian-frequency rhythms for every gene are becoming available [34–36]. Recently, Espinoza et al. (2010) produced one such microarray dataset, which measured 16 timepoints covering a 58-hour time period with a cold treatment in Arabidopsis[7]. Cold-stress genetic responses in Arabidopsis are particularly well-characterized, and have been shown to significantly dampen and phase-shift the oscillations of the core clock genes CCA1 and LHY, which have regulatory influences over some cold-responsive transcription factors, including CBF1 CBF2 and CBF3[20]. Disruption of the expression patterns of other circadian output marker genes due to cold treatment has also been reported, including constant over expression of CAB2 and CCR2, and constant underexpression of CAT3[6, 9]. For these reasons, this is an ideal dataset to test whether the PRIISM algorithm is able to separate the strong circadian-clock influences on cold-response genes from treatment-response influences.

Methods

A wide range of biotic and abiotic stress treatments have been shown to significantly disrupt the cyclic patterns of core circadian clock genes and their downstream target genes [3–8]. When a stress treatment is constantly applied, adaptive stress-response genes are expected to be differentially regulated, while influences from the circadian clock will cause oscillations in target gene expression patterns. In PRIISM, by projecting the gene expression data to the frequency domain using the Fourier transform [37], the resulting amplitude spectra peak at different frequencies, caused by these different influences. The Fourier transform is a mainstream signal processing technique that simplifies period gram analysis by identifying the dominant frequencies in the amplitude spectrum. By distinguishing the clock frequency range from the core circadian clock genes in the frequency domain, we can separate the spectrum to different sections containing treatment-related, clock-related and noise-related influences. Then, we project the amplitude spectra back to the expression domain to reconstruct isolated, independent gene expression patterns representing the effects of different frequency components. This method can be applied to any dataset which has sufficiently high resolution and length to measure frequencies of at least one cycle per day, and which uses a treatment which is applied at a frequency significantly different than the clock frequency.

PRIISM has four steps (Figure3). In the first step, gene expression data are pre-processed to fit the requirements of the Fourier transform, after which the Fourier transform is performed to produce an amplitude spectrum for every gene (Figure3A, 3B). In the second step, a clock vector that defines the frequency range and the amplitudes of the core circadian clock genes is identified based on the spectra of core circadian clock genes (Figure3C). In the third step, the clock vector is used to decompose every gene’s spectrum into three components (treatment, clock and noise; Figure3D). In the final step, the inverse Fourier transform is applied to project each spectrum component back to the expression domain, resulting in three independent expression patterns (Figure3E, 3F).

Figure 3

Workflow of the PRIISM algorithm. The 0 to 26 hour time-frame in the cold for AtgolS3 (AT1G09350) is used as an example.

Pre-processing and Fourier analysis

Time series gene expression data are often unevenly sampled, and the disruption of clock patterns caused by the treatment varies over time. To be able to apply the Fourier transform (which requires steady and evenly sampled input), pre-processing is required. First, the whole time course is divided into overlapped frames. The size of these frames can be changed depending on the experiment; If they are too long, then it may be difficult to capture changes over time, and if they are too short, then it is more difficult to capture the treatment-frequency patterns (particularly for low-resolution data). For this experiment, the first time frame is 26 hours long due to the two-hour light period at the start of the time period, and all the other time frames are 24 hours long, starting and ending at each light/dark transition ( Additional file 1: Figure S1A). Second, within each time frame, the gene expression data is interpolated (making the time points evenly sampled), and the mean of the gene expression data for each gene is shifted to zero (refer to explanation of Eq. 2). The Fourier transform is then applied on each overlapping time frame individually, and the final expression values for each timepoint are calculated using a weighted average for each time frame, where higher weights are used for expression values near the center of each time frame ( Additional file 1: Figure S1B).

Fourier analysis is a signal processing technique [37] for the study of two processes: The Fourier transform (the process of decomposing a signal into a sum of components with different frequencies) and the inverse Fourier transform (the operation of reconstructing the signal from these components). Specifically, the discrete Fourier transform (DFT) and its inverse have been used to transform gene expression signals and to reconstruct the discrete signal, respectively [38]. The Fourier coefficient of DFT (G _n ) measures the contribution of the corresponding frequency component to the original signal and is given in Eq. 1 [37]:

\begin{array}{cc}\hfill G_n=\sum _{k=0}^{K-1}g\left(kT\right)e^{-i2\pi \frac{n}{NT}k}\hfill & \hfill n=0,\dots ,N-1\hfill \end{array}

(1)

where, g(kT) is the sampled signal of K samples with the sampling interval T; i is the imaginary unit. The frequency of the corresponding component n is denoted as f _n (i.e.,\frac{n}{NT}), where N is the number of frequency components. The DFT maps a time course signal into the frequency domain by producing a spectrum. An amplitude spectrum (plotted as the amplitude versus frequency) is a common frequency domain representation of the original signal. Fast Fourier transform (FFT) is an efficient algorithm to compute the DFT and its inverse [37]. Because of its popularity, it has been built into most modern analysis tools including MATLAB and R [39, 40].

The Fourier coefficient of the zero-frequency component (G ₀ ), derived from Eq. (1) where f _n  = 0, is shown in Eq. 2 as given in [37]:

G_0=\sum _{k=0}^{K-1}g\left(kT\right)e^{-i2\pi k*0}=\sum _{k=0}^{K-1}g\left(kT\right)

(2)

Note that there is a dominant peak at zero frequency in the spectrum of the expression value, which may bias the identification of the true dominant peak to frequency zero. To avoid such bias, we shift the mean of the time course gene expression values for each gene to zero (and consequently G ₀  = 0), leading to the removal of the peak at zero frequency. For example, the mean expression value for the gene shown in Figure3A is reduced from 10.6 to 0, and will be added back proportionally to the reconstructed gene expression values during signal recomposition.

Identification of the circadian clock frequency range

The Arabidopsis circadian clock is composed of multiple feedback loops. Three genes, Circadian Clock Associated 1 (CCA1), Late Elongated Hypocotyl (LHY) and Timing of CAB Expression 1 (TOC1) compose the first and most important feedback loop controlling the circadian clock, while Pseudo Response Regulators 7 and 9 (PRR7 and PRR9) form a secondary feedback loop with CCA1 and LHY, and a third feedback loop involving TOC1 is regulated by unknown components [41, 42]. It has been found that through these feedback loops, eight core circadian clock genes (CCA1 LHY PRR7 PRR9 ELF4 GI LUX and TOC1) and their downstream gene targets regulate a wide range of downstream pathways, including germination, leaf development, organelle morphology, photosynthesis, and cell wall development [2, 18, 43–46].

The Fourier transform is performed on these eight core circadian genes (Figure3Ci). The frequency components with relative amplitudes greater than 0.7 (corresponding to half of the maximum value in the spectra) are chosen as dominant frequencies [47]. We define the union of these eight sets of dominant frequencies as C ircadian C lock F requency R ange (CCFR), noted as f _{c_min} f _{c_max} , where f _{c_min} is the lowest frequency, and f _{c_max} is the highest frequency (Figure3Ci). Note that in this example, the dominant clock frequency is significantly lower than one cycle per day, due to the stress-induced disruption of clock patterns. The weight of each frequency component in the CCFR is derived as:

\begin{array}{cc}\hfill w_n=\frac{\sum _{m=1}^8{\left|G_{\mathit{mn}}\right|}^2-min\left(G\right)}{max\left(G\right)-min\left(G\right)}\hfill & \hfill n\in [c\_min,c\_max]\hfill \end{array}

(3)

where |G _mn | is the magnitude of the Fourier coefficient of the n _th frequency component for the m _th core circadian gene, G=\left\{\sum _{m=1}^8{\left|G_{mc\_min}\right|}^2,\sum _{m=1}^8{\left|G_{m(c\_min+1)}\right|}^2,\dots ,\sum _{m=1}^8{\left|G_{mc\_max}\right|}^2\right\} is the set of the summed power of eight core clock genes present at each frequency component within the Circadian Clock Frequency Range (CCFR), and w _n is the weight for the frequency component at frequency f _n . The vector {w _{c_min} w_{c_min+1},…, w _{c_max} } defines the gain-frequency response of a tapering bandpass filter within the CCFR.

Signal decomposition and recomposition

We apply Fourier analysis on each gene, producing the relative amplitude spectrum from which we identify three distinct sections: Treatment-frequency, clock-frequency and noise-frequency components (Figure3Di). For the treatment-frequency decomposition, given a relatively narrow frequency band, we used a low pass filter with a steep cut-off frequency to gain the optimal balance between removing ringing artifact and approximating desired frequency response [48]. This issue is addressed in detail in the “Justification for choosing a steep cut-off frequency for the low-pass filter” section of the ( Additional file 1: Figure S3Figure S4 and Figure S5) [49].

The clock component is derived by bandpass filtering. Fourier coefficients of the clock components of each gene are modulated by the weight of the corresponding frequency components, as given by Eq. 4:

\begin{array}{cc}\hfill {\widehat{G}}_c=w_c{G}_c\hfill & \hfill c\in \left[c\_min,c\_max\right]\hfill \end{array}

(4)

The tapering filtering results in clock-frequency expression patterns that are noise reduced and with less artifacts caused by a discontinuity in the filter function. The reconstructed high frequency expression pattern is considered to be noise, and it is not studied in this paper. Therefore, we simply applied an ideal high pass filter. The reweighted spectra used for the signal reconstruction of the three frequency components sections are shown in Figure3Dii.

The inverse discrete Fourier transform (IDFT) is calculated according to Eq. 5 [37]:

\begin{array}{cc}\hfill g\left(kT\right)=\frac{1}{N}\sum _{n=0}^{N-1}{G}_n{e}^{i2\pi \frac{n}{NT}k}\hfill & \hfill k=0,\dots ,K-1\hfill \end{array}

(5)

The inverse Fourier transform is performed on the full spectrum, including the filtered spectra for each gene. Similar to using the clock vector as a tapering band-pass filter to remove noise, we added a coarse graining process to make sure there is no overlapping between any of the two frequency bands, which may increase the robustness of component selection. The mean of the original gene expression values (which was removed in the pre-processing step), is added back proportionally to each gene expression curve based on the amplitude distribution of each component in the spectra before shifting the mean (Figure3F), according to Eq. 6:

\begin{array}{cc}\hfill g_L^{\text{'}}\left(kT\right)=g_L\left(kT\right)+\frac{\sum _{k=0}^{K-1}g\left(kT\right)\sum {\left|G_L^0\right|}^2}{K\sum {\left|G_n^0\right|}^2}\hfill & \hfill k=0,\dots ,K-1\hfill \end{array}

(6)

where g_L^{\text{'}}\left(kT\right) is the treatment expression level at timepoint kT for a given gene, g_L\left(kT\right) is a result of inverse discrete Fourier transform (Eq. 5) on treatment frequency at timepoint kT, G _n ⁰ is the Fourier coefficient of original gene expression values and G _L ⁰ is the Fourier coefficient of original gene expression values in the treatment-frequency band [0,f _{c_min} -1]. Similarly, we compute g_C^{\text{'}}\left(kT\right) and g_N^{\text{'}}\left(kT\right) i.e., the clock expression level and noise level.

Note that because the entire warm and cold gene expression datasets are mean-shifted based on their relative amplitudes in each component, the reconstructed time-zero fold change values may not necessarily be equal to zero (Figure2B).

Results and discussion

This study analyzes an Arabidopsis Affymetrix ATH1 microarray dataset (containing 22,810 probes) generated by Espinoza et al. (2010), which consists of 16 timepoints collected over the course of 58 hours in both warm (20 °C) and cold (4 °C) conditions under a 16-hour light/8-hour dark cycle starting at ZT14 (14 hours after dawn) [7]. This dataset was chosen for the analysis because it has separate control and treatment arrays, it has sufficiently high resolution (sampled at 2 hours and every 4 hours after that), and cold is a well-studied treatment in Arabidopsis[6, 7, 9, 20, 50, 51].

Gene expression data was RMA normalized using the “affylmgui” program available as part of the Bioconductor software package and annotated using annotation data available from TAIR (version 10, available ftp://ftp.arabidopsis.org/Genes/TAIR10_genome_release/). The gene expression data were interpolated to every 2 hours using B-spline regression, and were segmented into four overlapping gene expression time frames (from both the warm and cold treatments), which were combined using a weighted average ( Additional file 1: Figure S1) [52, 53]. PRIISM was applied on this “original” dataset, resulting in three independent and isolated gene expression datasets (treatment-frequency, clock-frequency and noise-frequency).

Treatment-response gene discovery

In order to show the advantage of PRIISM, we identified known cold-response genes using maximum fold changes and principal component analysis in the treatment-frequency dataset compared to the original dataset. Fold change values were calculated by subtracting the logged gene expression value in the warm from the logged gene expression value in the cold at every timepoint. Lists of Arabidopsis genes upregulated by cold treatment when grown on agar plates or grown in soil were collected from a previous study by Vogel et al [54]. The 302 cold-upregulated genes found in the intersection of these lists were used to define the set of “cold standard” (COS) upregulated genes. Receiver-Operator-Characteristic (ROC) curves (which have been shown to be an effective method for evaluating gene expression data [55]) were generated for these COS-upregulated genes [54] by distinguishing each ranked gene as either a true positive or a false positive (Figure4). A larger area under an ROC curve indicates that more COS-upregulated genes are identified. The line at which the number of true positives is equal to the number of false positives is indicated in Figure4, and only the data above this line are considered biologically relevant. By ranking genes by their maximum fold change values in the treatment-frequency dataset, 52.6% (159/302) of known COS upregulated genes can be identified, compared to only 21.2% (64/302) in the original dataset (Table1) [54]. This difference may be explained by the disruptions contributed by the clock-frequency influences and the noise-frequency influences, which are present in the original dataset. This shows that more COS-upregulated genes can be identified by ranking by the maximum fold change in the treatment-frequency dataset compared to the original dataset.

Figure 4

ROC curves for COS-upregulated genes. ROC Curves for the 26-hour treatment-frequency fold change (dashed black line), the treatment-frequency maximum fold change (solid black line), the original maximum fold change (solid grey line), and original PCA plot distance data (dashed grey line) are shown. The point at which the number of false positives is equal to the number of true positives (dotted grey line) and random gene selection (dotted black line) are also shown.

Table 1 Summary of ROC analysis for genes upregulated by cold treatment

Principal component analysis (PCA) is a linear component composition method that has been applied to summarize different gene expression influences under different conditions, and consequently has been used for differential gene expression studies in microarray datasets [56]. PCA was performed on the original dataset ( Additional file 1: Figure S2A), and the Euclidean distance from the bottom-left of the PCA plot of the first and second component was used to rank the genes, allowing for the construction of an ROC curve based on this data (Figure4). These data show that only 13.9% (42/302) of the cold upregulated genes can be identified in the original PCA plot. The first PCA components of the treatment-frequency data and the clock-frequency data were also plotted ( Additional file 1: Figure S2B) and ranking based on Euclidean distance from the bottom-right was able to identify 46.0% (139/302) of the COS-upregulated genes.

These results showed that, in both maximum fold change and PCA analyses, the ranked treatment-frequency fold change analyses produce fewer false positives than the original methodology by distinguishing more COS-upregulated genes (Table1).

The identification of important gene regulation timepoints using PRIISM

In the previous section, it was shown that gene discovery in the treatment-frequency data produced by PRIISM constantly outperforms the same analyses on the original data. Although these approaches are useful for poorly studied treatment responses, a knowledge-based approach may be used to identify more treatment-response genes with a lower false positive rate.

Cold treatments have been shown to induce the expression of the transcription factors C-repeat/DRE Binding Factor genes CBF1, CBF2 and CBF3[57], which are induced in parallel with the cold transcription factors RAV1 and ZAT12[50]. Some of the important targets of CBF transcription factors include Cold-Responsive (COR) genes COR15A, COR15B COR47, and COR78[20, 50, 58, 59]. All of the cold transcription factors and targets included in these lists have also been shown to be gated by the circadian clock, making them ideal for evaluating PRIISM’s ability to remove clock-frequency influences [6, 20, 23].

In the treatment-frequency data, a peak in the fold change patterns can be observed in the well-studied cold response transcription factors and cold regulated (COR) response genes at the start of the first night (at approximately 26 hours) (Figure5C, 5D). The peaks of the transcription factors can be seen to occur before the peaks of their target genes, as is expected for a TF-target relationship. By contrast, these peaks are not apparent in the original fold change data (Figure5A, 5B). For this reason, an ROC curve was computed using the fold change value at 26 hours in the treatment-frequency fold change data (Figure4, Table1). Table1 shows that 194/302 (64.2%) of the true-positive COS-upregulated genes can be identified with a 50% false positive rate in the treatment-frequency 26-hour fold change data, compared to only 64/302 (21.2%) for the maximum fold change in original data and 42/302 (13.9%) for the PCA plot of the original data .

Figure 5

Fold change patterns of cold transcription factors and target genes before and after PRIISM processing. The original fold change patterns for important cold transcription factors (A) and some of their important target (COR) genes (B) are shown, along with their the treatment-frequency fold change patterns for the same genes (C & D). Thin, dashed vertical lines are used to indicate the 26 hour position on each graph.

This data shows that the fold change value at 26 hours in the treatment-frequency data is the best predictor of whether a gene is involved in adaptive cold response. The top 25 ranked genes based on fold changes at 26 hours in the treatment-frequency dataset are shown in Table2. Included in this table is the “Cold Upregulation Category” for each gene, which indicates whether a gene was upregulated in the cold when plants were grown in soil (“Soil”), on agar plates (“Plate”), on both growth mediums (“COS”), or on neither (“Novel”) in Vogel et al’s study [54]. In this table, 22/25 of the genes belonged to the COS group, 2 belonged to the “Plate” group, and 1 belonged to the “Soil” group, suggesting that the PRIISM method has successfully identified known cold-regulation genes [54]. Table3 shows the top 25 ranked genes which were not part of the COS-upregulated gene list in Vogel et al (2005) [54]. 10/25 of the genes in this list belonged to the “Soil” group, 9 belonged to the “Plate” group, and 6 were novel genes not identified in Vogel et al’s study [54]. All of the novel genes (and all but one of the 25 genes in this list) have been previously identified as being involved in cold response in other studies, suggesting that PRIISM has identified a list of very important cold-response genes (See “Comments” column in Table 8).

Table 2 Genes ranked based on their treatment-frequency fold change values at 26 hours

Table 3 The top 25 ranked non-COS genes based on treatment-frequency fold change values at 26 hours

The results of a case study on ATGolS3 (AT1G09350), the gene with the largest fold change in the treatment-frequency data at 26 hours are shown in Figure6. The logged original gene expression curve under warm conditions has a minimum expression level of approximately 6, which is reflected by a flat treatment-frequency expression curve with a nearly constant value of 6 (Figure6A). The rhythmic pattern of the original data in warm conditions is captured in the clock-frequency gene expression curve, and the sharp peaks and sudden changes in slope are captured in the noise-frequency curve (Figure6A). The original gene expression data under cold conditions peaks quite strongly during the first night but retains some cyclical expression. The PRIISM-processed gene expression data shows that the treatment-frequency gene expression is constantly higher in the cold, with a peak at 26 hours, while the clock-frequency gene expression data is only marginally increased, but is increased more in the first day than in the second day (Figure6B, 6C). The fold change graph shown in Figure6C indicates that most of the increase in gene expression is due to treatment-frequency influences for this gene, but the clock-frequency influences upregulate the gene more strongly early in the cold treatment. The noise-frequency fold change pattern matches many of the sharp peaks and valleys in the original fold change pattern, suggesting that much of the noise has indeed been removed (Figure6C).

Figure 6

A case study examining PRIISM output gene expression and fold change data. The fold change patterns (A), warm gene expression patterns (B) and cold gene expression patterns (C) for the original and PRIISM-processed data for AtgolS3 (AT1G09350), the most highly upregulated gene in response to cold at 26 hours in the treatment-frequency data.

To test the statistical significance of PRIISM’s ability to discover treatment-response genes, P-values were calculated using a Z-test for both the maximum fold change from the original dataset and the fold change values at 26 hours in the treatment-frequency dataset. Figure7 shows the number of genes that were found to be significant (P value < =0.05) in these tests, and how many belonged to the COS-upregulated gene list from Vogel et al [54]. Out of the 161 genes significant in the treatment-frequency data at 26 hours, 98 of them (60.9%) were COS upregulated genes, compared to 154 out of 379 (39.3%) for the original dataset.

Figure 7

Venn diagram showing COS-upregulated genes in original and PRIISM-processed significant gene lists. The number of genes in the overlaps between COS-upregulated genes (A) and the significant genes (P value ≤ 0.05) in both the maximum fold change in the original dataset (B) and the fold change at 26 hours in the treatment-frequency dataset (C) are shown.

Clock-frequency data analysis

The clock vectors calculated by Equation 3 under both warm and cold conditions for each of the time frames are shown in Figure8. The difference between the length and the shape of the warm and cold vectors indicates the circadian rhythm disruption caused by the cold stress. Figure8A shows drastically different frequency profiles for the warm and cold conditions, caused by an abrupt phase shift in the expression data. The clock genes continue to have disrupted frequencies in the second time frame (Figure8B), but appear to return to normal oscillating frequencies, possibly with different phases, in time frames 3 and 4 (Figure8C, 8D).

Figure 8

Clock vectors under warm and cold conditions.

To study whether the clock-frequency data produced by PRIISM successfully isolated cyclic clock influences from treatment-response influences, the clock-frequency gene expression patterns of eight well-studied cold response genes were matched with standard clock patterns according to the pattern-matching algorithm HAYSTACK [26]. This algorithm (the key component of The Diurnal Project) utilizes a model-based pattern matching algorithm to calculate the phase and cyclic pattern type for each gene in a dataset, and also calculates the correlation of each gene to the closest model, which can be used as an indication of how strong the clock influence is on the gene [26]. HAYSTACK provides T-test P-values indicating the probability that an input pattern matches a gene expression model, and provides several types of cyclic clock pattern models to use for comparison [26]. This analysis included the COR genes which have been shown to be under circadian clock control under warm conditions, but gated by cold transcription factors (including the CBF genes) under cold conditions [20]. The results in Table4 indicate that the P values for the clock-frequency gene expression data from PRIISM are substantially lower than the original data (under both warm and cold conditions), often by several orders of magnitude, demonstrating enrichment of clock-frequency gene expression in this data.Note that the remaining portion of the spectrum of the clock-frequency components is simply discarded in PRIISM. In our future work, it will be interesting to further test whether feeding it into the treatment-frequency component will construct more precise results.

Table 4 A comparison of the clock patterns between PRIISM-processed and original gene expression data

Conclusions

Circadian rhythm pathways influence the expression patterns of as much as 31% of the Arabidopsis genome through complicated interaction pathways, and have been found to be significantly disrupted by biotic and abiotic stress treatments, complicating treatment-response gene discovery methods due to clock pattern mismatches in the fold change statistic. The PRIISM algorithm outlined in this paper is designed to separate pattern changes induced by different forces, including treatment pathways and circadian clock rhythm disruptions. By applying PRIISM on a cold-response dataset, we systematically evaluated our method using maximum fold change and PCA analyses. The results of this study showed that the ranked treatment-frequency fold change analyses produce fewer false positives than the original methodology, and the 26 hour timepoint in the PRIISM produced dataset was the best statistic for distinguishing the most known cold-response genes. In addition, PRIISM also provides gene expression data which represent only circadian clock influences, and may be useful for circadian clock studies. In fact, any existing analysis approach on gene expression data can utilize PRIISM to separate circadian-influenced changes in gene expression. In conclusion, PRIISM is a novel approach for overcoming the problem of circadian disruptions from stress treatments on plants. PRIISM can be integrated with any existing analysis approach on gene expression data to separate circadian-influenced changes in gene expression, and it can be extended to apply to any organism with regular oscillations in gene expression patterns across a large portion of the genome. In future work, we will apply the discrete wavelet transforms (DWT) on higher resolution datasets in order to further enhance the ability of PRIISM to distinguish circadian clock disruption influences from treatment-response pathway influences.