We consider a recently obtained transcriptomics data set from Arabidopsis thaliana plants exposed to eight environmental conditions differing in light intensity and/or temperature (Caldana et al., 2011). The response to each condition was followed over 23 time points (including the zero time point) in 20 min intervals up to 360 min starting from the same reference. In addition, samples were taken after 5, 10, 640, and 1280 min. We note that the same set of time points, covering both a linear and a logarithmic scale, was used in all conditions. Arabidopsis thaliana plants, grown in soil at 21°C and 150 μEm⁻²s⁻¹, were either kept at this condition or transferred to seven different environments reflecting a light gradient ranging from darkness to high-light stress and a temperature gradient from 4°C over 21 to 32°C. To identify similarities and differences between these data sets, either one (light or temperature) or both environmental parameters were changed. This resulted in the following seven (combinations of) environmental conditions: (i) 4°C and darkness (abbreviated as 4-D), (ii) 21°C and darkness (21-D), (iii) 32°C and darkness (32-D), (iv) 4°C and 85 μEm⁻²s⁻¹ (light; 4-L), (v) 21°C and 75 μEm⁻²s⁻¹ (low-light; 21-LL), (vi) 21°C and 300 μEm⁻²s⁻¹ (high-light; 21-HL), and (vii) 32°C and 150 μEm⁻²s⁻¹ (light; 32-L). Together with the set of plants kept at the original conditions (21°C and 150 μEm⁻²s⁻¹, abbreviated as 21-L and referred to as control), eight different conditions were considered, as illustrated in Figure 1. Single RNA sample isolated from a pool of the six independent plants per condition and time point was used for transcript profiling with the Affymetrix ATH1 array, resulting in eight data tables including the time-resolved gene expression (Caldana et al., 2011).

As a starting point, we used the 15,089 genes analyzed in the original study after normalization and filtering of the microarray data using standard methods. The data set was subsequently filtered for genes which exhibit fold-change of at least two with respect to the first time point (0 min corresponding to the control condition) in at least two time points for any condition and a coefficient of variation within the considered time domain of at least one. These are reasonable criteria to remove less informative, and possibly non-differentially expressed, genes. This strategy was used as there is only one replicate for each time point, precluding application of more rigorous statistics. The preprocessing step resulted in identification of 2,276 genes which were used in the analysis based on STATIS and dual-STATIS.

Suppose that we are given K data tables X₁, …, X_K corresponding to the investigated conditions, and let they be gathered as blocks of a matrix X = [X₁|…|X_K]. Moreover, let each one of these tables include observations on n genes, corresponding to the rows, over t time points, given by the columns. Furthermore, we will use xijl denote the expression of the i^th gene in the j^th time point of the l^th condition. In addition, each data table is preprocessed, i.e., centered and normalized according to the recommendation of Abdi et al. (2012) and Smilde et al. (2003). Each of the rows is assigned a non-negative mass m, so that the masses over all rows sum up to one. Here we used uniform values for the masses, resulting in a vector of masses m is of unit length.

In the first step of STATIS, the cross-product matrix Si=XiXiT is calculated for every table X_i, 1 ≤ i ≤ K. Each of the resulting cross-product matrices captures the relationship between the rows, also known as within-table structure, as seen by the corresponding table. To investigate the between-table structure, one then determines the matrix C in which the entry c_ij equals the inner product of the cross-product matrices S_i and S_j, 1 ≤ i,j ≤ K. More concretely, cij=∑p=1t∑q=1tSp,qiSp,qj. The optimum weights for combining the data tables are then obtained from the first eigenvector of C, u₁, which corresponds to the eigenvalue of largest weight. The entry ui1 of the first eigenvector captures to the similarity between the i^th table, 1 ≤ i ≤ K, and all other investigated data tables. Therefore, the quantity αi=ui1∕∥u1∥ corresponds to the weights for the i^th table, so that ∑i=1Kαi=1. The compromise cross-product matrix is then given by S=∑i=1KαiSi.

In the second step, generalized eigen value decomposition of the compromise S is performed, so that S = PΛP^T under the constraint P^TMP = I, where I is the identity matrix. It is easy to show that the eigenvalue decomposition of S is equivalent to generalized SVD of X, since S = PAP^T with A = diag(a), where a = [α₁1^T|…|α_K1^T]. The loadings can then be computed as Q=XTMPΔ-1, where Λ = Δ². The compromise factor scores, F, can be computed from S as F=SMPΔ-1.

The loadings of the variables from table i, 1 ≤ i ≤ K, are obtained as Q=XiTMPΔ-1, while the factor scores, F, for table i, 1 ≤ i ≤ K, are given by Fi=SiMPΔ-1.

In addition, similar to PCA, in STATIS one can determine the contribution of rows, columns, and tables to the principal components. The contribution of the i^th row to the component b is given by ctri,b=mifi,b2λb, where λ_b is the l^th largest eigenvalue from the decomposition of the compromise S and f_i,b denotes the factor score of the i^th observation for the b^th dimension. Since λb=∑imifi,b2, all contributions take values between 0 and 1. The larger the contribution, the more the row contributes to the component. In order to base the interpretation of a given component on the rows that have significant contributions, a bootstrap procedure can be applied. Briefly, a bootstrap sample (i.e., a uniform random sample with replacement) is performed on the set of tables. By repeatedly conducting STATIS analyses on random samples of tables bootstrap contributions ctri,b* are obtained. Finally, these bootstrap estimates of contribution are transformed in bootstrap ratios, that can be interpreted in a similar way as a t-statistic (Abdi et al., 2012). Significantly contributing rows whose factor scores have different sign can then be contrasted to interpret the component. Similarly, based on the loading of the j^th variable for the b^th dimension, q_j,b, the contribution of j^th column to the b^th component can be determined as ctrj,b=ajqj,b2.

Likewise, the contribution of the l^th table can then be simply defined as the sum of the contributions of its variables: ctrl,b=∑jctrj,b.

The contribution of all tables of often referred to as the interstructure or inner product map (Lavit et al., 1994; Abdi et al., 2012). The contributions of columns and tables to a component have the same properties as the contributions of rows, described above; hence, they can be used in interpretation in a similar way. It is also important to identify tables that represent well the overall similarity with respect to all considered tables, i.e., to quantify their contribution to the compromise S. To this end, the investigation of table weights α_i as well as the R_V coefficient have been considered (Thioulouse et al., 1997). The R_V coefficient is a measure of similarity of two cross-product matrices S_k and Sk′ both n × n matrices, defined as the inner product (Horn and Johnson, 2006).

RVSk,Sk′=∑it∑jtSi,j,kSi,j,k′, which takes values between 0 and 1 as both matrices are positive semidefinite and normalized, thus, being equal to the cosine of the two matrices (Abdi et al., 2012).

The analysis based on dual-STATIS can be carried our analogously on the transposed data table. In other words, in the first step of dual-STATIS, the cross-product matrix is calculated for every table X_i, 1 ≤ i ≤ K, as Si=XiTXi. Each of the resulting cross-product matrices in the dual-STATIS setting thus captures the relationships between the columns of the original data tables, referred to as within-table structure. The subsequent steps of dual-STATIS are identical to the STATIS analysis, finally providing a compromise, S′, for the variables rather than the observations of the original data tables.

A reference implementation, on which our analysis is based on, is available in the ade4 package for the open source statistical software R (Gentleman and Rober, 1996; Thioulouse et al., 1997).

Following considerations by Caldana et al. (2011) we consider a second normalization strategy for the whole time-course data set. Here, all transcript expression levels at every time point over all experimental conditions are normalized with respect to the corresponding control time point by subtracting the obtained fold-changes. This does not affect the control condition; however, it allows better assessment of the stress similarities and specificities by further eliminating circadian and diurnal responses as measured under control condition. We note that based on this normalization strategy, the previously outlined procedure for identifying informative genes now resulted in 3025 genes. We repeatedly apply both STATIS and dual-STATIS analysis to reveal and emphasize findings pertaining to diurnal changes.

Similar to the classical PCA, analysis based on STATIS allows describing the contributions of tables, variables, and observations. Here, the K = 8 tables correspond to the eight environmental conditions obtained from the Arabidopsis experiment described in Section “Materials and Methods,” where the variables denote the 23 time points measured for each condition and the observations correspond to the 2,276 genes which are obtained using the aforementioned filtering strategy. An analysis of the contributions of the individual data tables to the compromise S then renders it possible to determine the condition which has the strongest influence on the determined principal components and further reflects good overall similarity to the other data tables from the remaining conditions. Figure 2A displays the obtained table weights α_i for all conditions (x-axis), as well as the similarity quantified by the RV2 coefficient of the compromise and the individual tables (y-axis). One can observe that the low-light condition, 21-LL, as well as the 21-D, both at control temperature, have the strongest effect in the derivation of the compromise. This indicates that the transcriptomics changes under normal temperature regimes coupled with darkness/low-light conditions are characteristic for the entire data set, thus, superimposing the changes in the remaining combinations of conditions. Moreover, the light stress (21-HL) and cold stress condition (4-L) have the lowest influence within the complete dataset reflected by the low RV2 coefficient and table weight.

Subsequent investigations of the table interstructure, i.e., the contribution of the tables to the first, second, and third principal components (Figure 2B, left and right) further corroborates this finding: With exception of the aforementioned light stress condition (21-HL), the remaining light regimes at control temperature conditions, namely, 21-LL, 21-D, and 21-L, exhibit strong contributions to the first principal component. Changes in ambient temperature, e.g., cold and heat, contribute reciprocally to the second principal component. Additionally, the second principal component separates the light conditions at control temperatures (21-D and 21-L). Finally, the effect of the diurnal rhythm is captured in the third principal component. Here, the darkness conditions, namely, 4-D, 21-D, and 32-D, contribute positively, while the high-light condition, 21-HL, does so negatively. Interestingly, the data table corresponding to the combination of cold stress and light, i.e., 4-L, exhibits a similar contribution to the third principal component as the other darkness conditions. This observation could be partially attributed to the fact that the light intensity at 4°C is lower (85 μE) than the others “L” controls (150 μE) as an attempt to avoid light stress. Furthermore, this might be related to the attenuating effect of the low temperature as previously reported by Caldana et al. (2011) as well as a general stress response, since both extended dark and cold conditions have been widely reported to affect plant growth.

Furthermore, the first principal component in the interstructure characterizes the dominant transcriptional response within the complete dataset as it corresponds to the largest eigenvalue which is proportional to the variance explained. The strong effect on the covariance structure of transcripts in all datasets under the control temperature condition can be attributed to the experimental design; namely, out of eight environmental conditions, four were conducted at different light regimes and the same temperature of 21°C. Moreover, the second strongest effect, corresponding to the contributions to the second principal component is the contrast of cold and heat stress, while the third strongest influence is ascribed to the diurnal cycle.

Apart from the contributions of individual tables to the compromise, the pairwise cosine similarity, i.e., the R_V coefficient, can be used to further assess the similarities between the tables. Figure 3 illustrates the similarity between conditions as a heatmap (based on the function heatmap.2() in the gqplot R package). It becomes apparent that the type of conditions, in all considered combinations of perturbations, has an effect on the clustering: Primarily, similar temperature conditions (4, 21, and 32°C) cluster together (as derived by the cluster dendrogram) in comparison to conditions in the same light regime (4-D and 4-L and 32-D and 32-L, respectively). Under control, the secondary effect of light intensity can be seen by the co-clustering of the low-light/dark (21-LL, 21-D) and high-light/normal light (21-HL, 21-L) conditions.

In summary, these results demonstrate the advantages of the multi-table analysis based of STATIS as compared to classical PCA. In the classical PCA, distinct dimensions, corresponding to experimental factors, of the investigated data tables cannot be separately investigated. The reason is that individual table structure gets lost, since a combination of environmental conditions and time points, e.g., 21-L 60 min, serves as variables (often referred to as loadings) in a PCA. The analysis by STATIS however, allows the characterization of each of the eight experimental conditions independently of the time-factor. A similar consideration of the same dataset by classical PCA – as conducted in Figure 2A in Caldana et al. (2011) – displays each condition by all of its measured time points resulting in 184 variables in the loadings plot.

The analysis of the contributions of variables reflects the influence of the particular time point with respect to the overall transcriptional adjustments of the Arabidopsis in response to the applied stresses. Figure 4 shows the contributions of each time point per condition on the first (x-axis) and second (y-axis) principal component. Arrows between two successive time points illustrate the sequential progression of the contributions over time with the compromise space defined by the first two principal components. While all the trajectories of time points >0 min show substantial contribution to either the first or second principal component – clearly illustrating the ongoing transcriptional adjustments of Arabidopsis’ transcriptome – three observations are of particular interest:

Characterization of observations is a pivotal step of any descriptive statistical method. We note that in our setting, the observation correspond to the 2,276 genes. To this end, STATIS allows not only the investigation of the contribution of tables (i.e., conditions) and variables (i.e., time points), but also the analysis of the contributions of observations to the principal components. In particular, for the purpose of visualization, contributions of individual observations, rows in matrix X, are obtained by projection onto the principal components of the compromise. Figure 5 displays the contribution of genes to the first two principal components. Based on the projection of observation on the compromise, one can notice that the sign of the contribution of genes to the first principal component accounts either for the repression (red points) or the induction (blue) of the respective transcripts.

In order to determine which key pathways are associated with the genome-wide adjustments of gene expression, a gene set enrichment analysis (GSEA; Subramanian et al., 2005) of ontology terms is typically performed. However, an important intermediate step is to identify the observations (i.e., genes) which contribute significantly to either principal component. Here, bootstrap ratios (cf. Materials and Methods) are used to identify significant contributions at a significance level of 1%. Those genes (1371 in total) which are found to be significantly contributing to either the first or second principal component (illustrated as dark points in the small panel in Figure 5) are subjected to a GSEA by using MapMan terms (Thimm et al., 2004). Table 1 gives an overview of the MapMan terms (Version 1.1, January 2010¹) found to be over-represented at a significance level of 1% [enrichment determined by hypergeometric distribution (Rivals et al., 2007) and a FDR corrected p-value ≤0.01 cutoff (Benjamini and Hochberg, 1995)]. By carefully inspecting Table 1, it is becomes apparent that when plants sense changes in the environment (i.e., light and temperature), several signal transduction pathways are triggered to adjust to the “new” condition as illustrated by the over-enrichment of MapMan bins such as “hormone metabolism” – brassinosteroid, ethylene, and ABA, “receptor kinases,” “DNA synthesis/chromatine structure. histone,” “RNA regulation of transcription” and “micro RNA.”

Following this response, as a result of both diurnal changes and light/dark contrasts, the carbon/nitrogen (C/N) metabolism is clearly affected. This response is reflected by the over-enrichment of MapMan bins associated to photosynthesis-related processes such and “PS light reactions” (darkness; bins 1.1.1; 1.1.30) followed by an over-enrichment of bins involved in sugar metabolism. For instance, starch synthesis (high and normal light) and degradation (darkness) as a result of nutrient limitation, e.g., bin 2.1.2 “major CHO metabolism.synthesis.starch” and 2.2.2 “major CHO metabolism. degradation.starch.” The metabolism of minor sugars, such as members of the raffinose (cold) family and trehalose (darkness) which are considered to act as signal under stress or energy deprivation, are also affected (bins 3.1 and 3.2). High-light (21-HL) results in over-enrichment of bins involved in “PS.photorespiration” and detoxification such as “misc. gluthatione.”

The C/N imbalance at darkness conditions is further supported by the over-enrichment of bins involved in nitrogen metabolism and nutrient recycling such as the bins containing genes related to “Amino acid metabolism” and “Protein degradation.” Next to the C/N balance disruption, under limited nutritional condition, like darkness, remobilization of carbon and nitrogen can occur through the secondary metabolism, as it is reflected by the over-enrichment of the bins involved in “secondary metabolism.phenylpropanoid,” “secondary metabolism.glucosinolates,” “secondary metabolism.flavonoids.” Additionally, bins 14.15 “S-assimilation.AKN” and “misc.nitrilases” are over-enriched as a result of early changes of abundance of the metabolite O-acetylserine (OAS), which is in accordance with existing results analyzing changes of OAS level within this data set as well as providing further evidence that OAS is a signaling molecule [aside from its role in sulfur assimilation (Nikiforova et al., 2005) based on further studies of transgenic plants (Hubberten et al., 2012)]. Altogether, changes in C/N balance affect plant growth. Plant growth is also closely regulated by cell wall expansion which requires constant remodeling of its components. We observed several MapMan bins involved in the cell wall remodeling might play a role in the remobilization of sugar in response to sugar starvation (Lee et al., 2007). Transcriptional changes in genes related to cell wall modification, synthesis, and degradation under starvation or dark treatment has been associated with the putative role of cell wall as alternative energy source and as a means of restricting cell growth and elongation (Baena-Gonzalez and Sheen, 2008).

Finally, it is important to realize that these findings are based on genes of significant contribution to the first and second principal component of the compromise of eight environmental conditions and thus represent the common aspects of Arabidopsis’ transcriptomic response under light and temperature changes.

The previously performed three steps of a STATIS analysis can be complemented by a dual-STATIS analysis. As outlined in Section “Materials and Methods,” the analysis based on dual-STATIS employs the transposed transcriptomic data set (i.e., in each table, a row now corresponds to a time point, while a column stands for a gene). Therefore, it enables to obtain a compromise based on the covariance between time points. Figures 6A,B display the obtained table weights and the R_V coefficients, as well as the compromise of time points. We observe that – aside from minor differences – all eight experimental conditions exhibit comparable table weights (within the range of ∼0.34 to ∼0.37) and R_V coefficients (within the range of 0.83–0.94). This indicates that the covariance of time points is similar under all experimental conditions suggesting similar temporal progression of transcriptional regulation.

Moreover, the contribution of time points, by projecting the corresponding rows onto the principal components of the compromise, displays an almost perfect cycle. With the exception of some intermediate time points (∼120–280 min) as well as the last time point (1280 min), all observations contribute equally to either the first or second principal component as indicated by their equidistant displacement from the origin. Such representation could be explained by a gradually progressing system under tight control of diurnal rhythm. Interestingly, these results contrast the anticipation of stronger contributions of early time points as a result of fast system-wide adaptations to environmental changes.

The obtained findings pointed out that the control condition, 21-L, exhibits a considerable effect with respect to the contribution of conditions as well as time points (cf. Figures 2 and 4). Caldana et al. (2011) suggested to normalize the whole time-course for each stress to the respective time points measure under control condition to better assess the stress similarities and specificities further eliminating circadian and diurnal responses. Therefore, we repeated the analysis based on STATIS with a correspondingly re-normalized data set (cf. Material and Methods section): Figure 7A displays table weights and the coefficient of the compromise and the individual tables. In contrast to Figure 2A), now, the control condition has the least influence and, thus, the smallest similarity to the compromise. Furthermore, heat and control temperature conditions under low-light/darkness exhibit higher similarities to the compromise. A subsequent investigation of the interstructure (Figure 7B) shows that low-light/darkness conditions contribute the most to the first principal component, which captures differences among the light regimes across the eight experimental experiments. Moreover, the second principal component corresponds to temperature differences, since the low temperature conditions, 4-L and 4-D, contribute positively, while the high temperature experiments, 32-L and 32-D, contribute negatively. Finally, the third principal component facilitates the separation between the control, 21-L, and the remaining conditions.

The analysis of the contributions of variables, i.e., time points, additionally emphasizes the diminished effect of circadian regulation within the control condition (Figure 8A): as compared to all other seven conditions, the trajectories of the time point in 21-L display the weakest variations and contributions to either principal component. This observation is further confirmed by the fact that except from the high-light and the control conditions, 21-HL and 21-L, all others do not exhibit the aforementioned closing of the cycle with the last time point (1280 min). Similar to the interstructure, the contribution of the light regime to the first principal component (low-light/darkness positive; high-light negative) becomes more pronounced. Likewise, time points from the experiments in low temperature exhibit negative contributions, whereas those from high temperature experiments show positive contributions to the second principal component.

By performing an analysis based on dual-STATIS on the transposed and normalized data set, in which observations correspond to time points, Arabidopsis’ stress response becomes apparent (Figure 8B). In contrast to the temporal progression observed in the previously used dataset (cf. 4b), early time points, e.g., 5, 10, 20 40, 60 min, are grouped together in a non-sequential manner further separated from successive time points indicating the systems response to perturbation. This is additionally illustrated by considering the length of the trajectories between early time points in Figure 8A, which are generally of greater length than trajectories between two later time points.

From these observations, it becomes apparent that with limited presence of diurnal effects, the influence of observations, i.e., genes, allows a more detailed picture of the pathways affected by stress adjustments than the precious analysis. Table 2 displays the MapMan bins that are additionally enriched by performing a GSEA on highly contributing observations by using the normalized dataset (note that all MapMan bins exhibiting over-enrichment are also found using the normalized dataset). For instance, MapMan (sub) bins 2.1.2.2, 2.1.2.3, and 2.2.1.4 corresponding to “major CHO metabolism.synthesis.starch.starch synthase,” “major CHO metabolism. synthesis.starch.starch branching,” and “major CHO metabolism.degradation. sucrose.hexokinase” characterize in greater detail. Here, the effects of darkness and high-light conditions reciprocally lead to degradation of starch stored transiently in the chloroplasts or increase of synthesis as primary product of photosynthesis in leaves (Zeeman et al., 2007). As a second example, two more effects on transcriptional regulation become visible – MapMan bin 27.3.3 “RNA.regulation of transcription.AP2/EREBP, APETALA2/Ethylene-responsive element binding protein family” and 27.3.32 “RNA.regulation of transcription.WRKY domain transcription factor family.” The WRKY family of transcriptions factors, with up to 100 representatives in Arabidopsis (Eulgem et al., 2000), has been shown to be involved in response to abiotic stresses in plants (Chen et al., 2012). More specifically, one member, WRKY22 (AT4G01250), is exclusively induced over four, three, and fivefold in the three darkness conditions, 21-D, 4-D, and 32-D, respectively, and has been shown to regulate darkness-induced leaf senescence in Arabidopsis (Zhou et al., 2011). As a summary, the analysis of contributions of tables, variables, and observations by application of STATIS provides insights for different dimension of the data tables (i.e., conditions, time points, and genes). Moreover, the STATIS analysis on the re-normalized data allows exposing those biological components which remain conserved across all conditions, by subsequently reducing typically preserved dynamics of the organism, e.g., diurnal behavior, by application of different normalization strategies.