Order restricted inference in chronobiology
Abstract
Producción Científica
Full text
Received: Added at production Revised: Added at production Accepted: Added at production DOI: xxx/xxxx ARTICLE TYPE Order Restricted Inference in Chronobiology Yolanda Larriba 1 | Cristina Rueda 1 | Miguel A. Fernández* 1 | Shyamal D. Peddada 2 1 Departamento de Estadística e Investigación Operativa, Universidad de Valladolid, Valladolid, Spain 2 Department of Biostatistics, Public School of Health, University of Pittsburg, Pittsburgh, USA Correspondence *Miguel A. Fernández. Email: [email protected]a.es Summary This paper is motivated by applications in oscillatory systems where researchers are typically interested in discovering components of those systems that display rhythmic temporal patterns. The contributions of the paper are twofold. First, a methodology is developed based on a circular signal plus error model that is dened using order restrictions. This mathematical formulation of rhythmicity is simple, easily interpretable and very exible, with the latter property derived from the non-parametric formulation of the signal. Second, we address various commonly encountered problems in the analysis of oscillatory systems data. Specically, we propose a methodology for (a) detecting rhythmic signals in an oscillatory system, (b) estimating the unknown sampling time which occurs when tissues are obtained from subjects whose time of death is unknown. The proposed methodology is computationally ecient, outperforms the existing methods and is broadly applicable to address a wide range of questions related to oscillatory systems. KEYWORDS: Constrained Inference, Circular Data, Rhythmicity Detection, Timing Estimation, Oscillatory Systems 1 INTRODUCTION Locomotor activity, blood pressure or body temperature are just a few of the physiological and biological phenomena exhibiting rhythmic or oscillatory processes in nature. Such oscillatory systems contain one or more components that display periodic or rhythmic patterns over each observed period. For example, genes (i.e. the components) participating in a circadian clock often display a rhythmic pattern of expression as shown in Figure 1. The study of such components with temporal rhythmic patterns, and how these patterns change under dierent conditions, is called chronobiology. 1,2,3,4 Chronobiology has been an active area of research during the past two decades, with major impact on treating cardiovascular disorders like hypertension, 5 detecting genes associated with neurodegenerative disorders 6 or depression, 7 and improving the eectiveness of cancer treatments. 8 For instance, Haus 9 demonstrated that the timing of radiation according to host and/or tumour rhythms improves the toxic/therapeutic ratio of the treatment. These and other ndings in biomedical sciences have increased interest in chronobiological experiments, particularly in identifying and/or characterizing rhythmic processes. Although this article is motivated by gene expression studies associated with circadian clock, the proposed methodology is very general and is broadly applicable to other oscillatory systems such as the cell-cycle, 10,11 the endocrinology, 12 vascular processes, 13 etc. In fact, there is considerable interest in pharmacology, psychiatry and other areas of medical sciences to discover genes belonging to oscillatory systems that have rhythmic temporal expression.
2 Larriba et al Time Gene expression 10.0 10.2 10.4 10.6 10.8 Eif4b FIGURE 1 Observed gene expression (blue dots) along 24 hours for gene Eif4b. Green line represents the underlying estimated rhythmic signal From a statistical point of view, the modelling of chronobiological rhythms in biomedical sciences is a challenge because, unlike the time course nancial data in the stock market or heart rate data in intensive care units, the density of time points is generally low 14,15 and the number of periods of data is usually very small. 16,17 For these reasons, standard time series or Fourier models are not convenient. 18,19,20 Another challenge with these chronobiological data is that not all components display the same pattern of expression over time. 21,22,23 Despite this heterogeneity, rhythmic patterns of components of an oscillatory system usually display up-down-up patterns. Yet, it is important to note that in many cases the shape is not exactly sinusoidal or even symmetric. Figure 1 shows an example of these up-down-up asymmetric rhythmic patterns. Models based on parametric functions of time, such as Cosinor, 24,3 , have been proposed in chronobiology to model these patterns. 25 However, these parametric functions are too rigid, as other patterns (e.g. asymmetric ones) frequently appear in biological systems. The rst statistical problem to solve in this context is to determine if the observed pattern is rhythmic or not. There are a wide variety of procedures in the literature to detect rhythmicity including, among others, those based on autocorrelation 26 ; cosine curve-tting 10,27,3 or Fourier analysis. 28 Some non-parametric methods such as JTK_Cycle 29 (JTK) and RAIN, 21 that use Jonckheere- Terpstra test and the Kendall's tau correlation, have also been proposed in the literature. However, these approaches do not detect asymmetric rhythmic patterns properly. Recently, Larriba et al 22 designed an algorithm that successfully identify and classify circadian clock gene expression patterns. This latter approach can be considered as a precursor to the methods we present here. A fundamental assumption made in the above discussion is that the time corresponding to each biological sample is known or can be ascertained. However, in many instances, such as when dealing with samples obtained from human cadavers 6 or human organ biopsies, 30,31 the timing, i.e. the exact time corresponding to each biological sample may be unknown. In such cases, one needs to rst estimate or determine the time associated with each sample before investigating rhythmicity. There are, to our knowledge, two main procedures in the literature to cope with this problem of timing estimation, namely Oscope 32 and CYCLOPS. 33 Both procedures present limitations and do not address the problem in a general context. Oscope, is specically designed to recover cell cycle dynamic in single cell RNA-Seq experiments and works on the single transcript level, so that it is computationally intensive and highly sensitive to the inter-subject variability, inherent in human experiments. CYCLOPS is based on developing neural networks from the rst eigenvectors (called eigengenes) from a singular value decomposition (SVD) analysis. Notice that eigengenes suggest the fundamental (gene) expression patterns across the samples, which in turn represent a biological theme if the data are well organized, 34,35 see Figure S1 in the Supporting Information. CYCLOPS overcomes the drawbacks in Oscope but suers form its own weaknesses. First one is that the data should cover the entire periods, a downside in human studies, since target population must be substantially increased to ll in under-represented times of the day. Second, it uses additional information, such as gene rhythmicity evolutionary information, that it not always available. Third,
Larriba et al 3 the optimization problem solved is quite far, mathematically speaking, from a close-tting formulation. And fourth, in a neural network framework (which is like a black box), it is dicult to assess the inuence of outliers, due to biological or sample noise, or of other artifacts of the data. In this paper, we develop a general methodology, based on Order Restricte Inference (ORI), that solves the two main questions of rhythmic pattern detection and timing estimation commented above besides addressing other scientic questions related to oscillatory systems. In Section 2 we introduce a (circular) signal plus error model which is the basis of the methodology, we develop estimation and testing procedures that can be used to solve the problem of rhythmic pattern detection and the interesting question arising in real practice of the estimation of peaks. Also in this Section, a solution to the problem of timing estimation under the framework of circular signal models is proposed reformulating the statistical problem as that of deriving the optimal circular order . Section 3 is devoted to simulation experiments that illustrates the good performance of the new solutions by comparing them to those provided with alternative approaches, while Section 4 shows the results obtained using several real data sets. Finally, concluding remarks are provided in Section 5. 2 METHODS The key of our methodology is the denition of what we call circular signals . A circular signal can be graphically mapped as a function displaying a temporal up-down-up pattern. Such patterns are commonly seen in biological rhythmic processes as in cell-cycle and the circadian clock (see Figure 1). This up-down-up pattern over a discrete number of values can be described using mathematical inequalities that establish order restrictions among those values. We refer to these signals as circular , since periodic events in the Euclidean space can be mapped as circular processes in the Circular space. 36 Moreover, circular signals can be equivalently formulated both in the Euclidean as well as in the Circular space (see Section 2 in the Supporting Information). For simplicity of exposition, throughout this paper we shall use the term gene to describe the response variable of interest and the term gene expression for the outcome. 2.1 Circular signal model For each gene, suppose its expression is obtained at time points ti, i = 1, . . . ,n in each of the p periods of data, with T being the length of each period. In our set-up both p and T are known. Let Xij be the observed data collected at time point ti, i = 1, . . . , n, within the j th period, Xj= (X1j, . . . , Xnj)0 denote the vector of data at the n time points in the j th period j= 1, . . . , p and Y= (X1., . . . , Xn.)0, where Xi. denotes the average of data collected at time point ti across the p periods, for i= 1, . . . , n . We assume that, for each given time point ti , the data collected across p periods has same expected value, that the covariance matrix of Xj is a diagonal matrix, and that the period samples are independent from one another. Independence, is a reasonable assumption in many applications if p is not too large. Thus, we assume that, for each j= 1, . . . , p , the data satisfy the following signal plus error model: Xj=µ+j, (1) where the signal term µ has an up-down-up shape which is dened more precisely below. For convenience, we shall refer to the signal as circular signal and it should not be confused that the components of µ are angular parameters. These values are points in the Euclidean space that have a rhythmic pattern as shown in Figure 2. It is important to notice that no distributional assumption is needed for the results obtained below unless otherwise stated in the text. Denote L= arg min i=1,...,n µi and U= arg max i=1,...,n µi , the indices on the time point vector for which the signal reaches its minimum and maximum respectively. Common signal shapes in chronobiology describe up-down-up patterns, that is, the signal monotonically increases up to µU, and then decrease up to µL before increasing again so that these indices are usually unique. A typical pattern is provided in the left hand panel of Figure 2. However, as illustrated in this work, the real data exhibit more irregular patterns, such as those in Figure 1 and Figure S2 in the Supporting Information. Notice that L>U or L<U is not usually known to the analyst. Without loss of generality we assume that L>U and we maintain in the rest of the paper, while in practise both options are tested.
4 Larriba et al Now we provide the Euclidean space representation of the up-down-up signals (µ) and their equivalent circular ordered representation (φ) in the Circular space. Proposition 1 (Supporting Information) checks that these two representations are equivalent using a transformation TLU between the Euclidean and Circular spaces. One may refer to Section 2 in the Supporting Information for full details. Given this equivalence, for convenience, we refer to both representations as circular signals . A signal µ in the Euclidean space is said to be up-down-up i µ∈C=SLU CLU , where L, U ∈ {1, . . . , n} , CLU ={µ∈Rn:µ1≤ · · · ≤ µU≥ · · · ≥ µL≤ · · · ≤ µn≤µ1} . A signal φ in the Circular space is said to be circular ordered i φ∈Co={φ∈[0,2π)n:φ1 · · · φnφ1} where can be read as is followed by". φ it is said to follow the circular order o. FIGURE 2 Equivalent formulation of circular signal with L > U . Left: Euclidean space. Right: Circular space. The equivalence between circular signal in the Euclidean space and circular signal on the unit circle is illustrated in Figure 2. The utility of the equivalence of the two formulations to solve the problem of temporal order estimation is discussed in Section 2.4. Moreover, other methodological strategies will arise from the equivalence. Some of them are described in the Discussion Section. 2.2 Circular signal estimation We propose solution to the following mean squares optimization problem as estimator of the circular signal: Y?= arg min Z∈C n X i=1 (Yi−Zi)2. (2) The vector Y? is called the Isotonic Regression (IR) of Y with respect to C with equal weights. If variances are unequal then we use weighted least squares where the weights are inverse of variances at each time point. The IR estimator is a step function, with sets of consecutive components for which the estimator takes the same value, called level sets. See Robertson et al 37 for applications and algorithms to solve IR problems. Note that the order dened by C is not a closed convex cone but a union of closed convex cones (recall that C=SLU CLU ). Consequently the derivation of Y? is non-trivial for C , unlike when the cone of interest is a closed convex cone, see Robertson et al. 37 In order to derive the IR estimator dened in (2), we have designed a computationally ecient algorithm based on theoretical results on the IR estimator. Both the algorithm and the theoretical results are given in Section 3 of the Supporting Information.
Larriba et al 5 Implicit in the estimation of µ∈C , is the problem of estimating tU and tL , the times at which the peak and the trough occur. In the case of genes with periodic expression, these peak time points correspond to the biological/functional activity of the genes. Point estimators for tU and tL can be immediately derived from the point estimation of the indices U and L obtained in Algorithm 1, given in Subsection 3.2 in the Supporting Information. 2.3 Inference for circular signal normal models In chronobiology it is reasonable to assume that the random errors j , j= 1,2, . . . , p , are independently and normally distributed and that the variability in the random errors are not dependent on the value of the signal. Thus, we assume that the variances are homoscedastic. It is well-known that, under this assumption, IR yields the maximum likelihood estimator (MLE) of the corresponding parameter. 37,38 Therefore, in this case Y? as dened in (2), is the MLE of the circular signal µ . In particular, under the assumption of normality, it is straightforward to obtain condence intervals for tU and tL using standard parametric bootstrap. 39 Subsections 3.1 and 4.1 below compare the performance of this approach with Cosinor, an extended methodology to model rhythms in chronobiology, in simulations and real data, respectively. Testing circular signals We formulate the problem of identifying circadian genes using the following testing problem: H0:µ1=· · · =µn (at signal) (3) H1:µ∈C (circular signal) . Detection of rhythmic signals has been considered in the literature by several authors 2,29,21 as a testing problem. In Larriba et al 22 , the authors even classied signals into one of four dierent patterns. As remarked in Section 2.2, it is important to recognize that (3) is not a standard testing problem. Again, the standard ORI theory is not applicable directly because C=SLU CLU in (3) is not a convex cone. 37,38 In fact, the non-convexity issue arises even in the case of simpler alternative hypotheses such as the umbrella order where the location of the peak (or trough) is unknown, i.e. µ1≤ · · · ≤ µr≥ · · · ≥ µn with r unknown. 40,41,42 For this reason, as done in the case of umbrella alternative, we consider a two-step approach. First, we estimate L and U by L? and U? using the IR algorithm (see Subsection 3.2 in the Supporting Information) and we dene the testing problem assuming L and U are known. This way, we do not propagate the uncertainty estimates associated with L? and U? when dealing with Type 1 errors associated with the proposed methodology. Thus, the hypothesis testing problem of interest is formulated as follows: H0:µ1=· · · =µn (at signal) (4) H1:µ∈CLU (circular signal. L, U known) . We test the above hypotheses using the conditional test based on the likelihood ratio (LR) statistic. The use of conditional tests is not new, see Bartholomew, 43 Menéndez and Salvador 44 and Fernández et al 45 among others. It is basically a conditional version of the classical LR test for the above hypotheses where the critical value depends upon the number of level sets m used to derive Y? . The conditional test is computationally simple and often more powerful for interesting alternatives than the LR test (see the above references). The conditional α− level test is given by: Reject H0 if R≥c(m), where the LR statistic is given by R= n P i=1 ωi(Y? i−Y)2 , Y= ( n P i=1 ωiYi)/ n P i=1 ωi and c(m) is chosen so that Pr(χ2 m−1≥ c(m)) = α . Recall that in this work we are assuming that ωi=p σ2 and that σ2 is assumed initially known. In the case σ2 unknown, the R statistic changes and includes an estimator of σ2 . The new statistic has a beta distribution with parameters ((m−1)/2,(2n−m)/2) so that percentiles coming from that distribution are used, instead of the ones from χ2 m−1 , in the determination of the critical value. See Robertson et al 37 pages 69-70.
6 Larriba et al The performance of the above conditional test is compared with JTK, a well known method used in chronobiology to detect rhytmic patterns. 46,47 In Section 3.1 both methods are compared using synthetic data and in Section 4.1 using real data. 2.4 Temporal order estimation As noted earlier, in some applications, time points at which samples were obtained may be unknown. For example, when dealing with autopsy data, the time of death of a set of patients is usually unknown. Yet, using data from those patients, chronobiologists are interested in determining the temporal order among samples (patients' time of death) regarding gene expression data. This section provides a mathematical solution to the problem of temporal order estimation by deriving the optimal circular order among the time point indices of an observed data set. Note that, in a second stage, this method can be combined with the one developed in Subsection 2.3 to identify rhythmic genes. The proposed methodology provides a biologically interpretable solution and overcomes some of the shortcomings of recently introduced CYCLOPS methodology. We shall evaluate the performance of the two methods using measures based on mean squared error (MSE). Suppose Yk= (Yk 1, . . . , Y k n)0 represents the (gene expression) data corresponding to the kth gene, k= 1, . . . , K and D={Yk}K k=1 denotes set of all data vectors from the K genes. Let Π be the set of all possible circular orderings of all indices S={1, . . . , n} around a unit circle. Notice that for each circular order o∈Π , there is a circular signal model so that µ∈Co. For a given circular order o∈Π , we dene a measure of the distance between o and D as follows: d(o, D) = K X k=1 n X i=1 νkYk i−Y?k o,i 2, (5) where Y?k o= (Y?k o,i )n i=1 denotes the IR of Yk o= (Yk o,i)n i=1 under the circular signal model that generates o, and νk denotes the weight associated with the kth element in the data set. For instance, when the experiments are subject to dierent variability, then νk=1 σ2 k . The problem of determining the temporal order among the specimens using the gene expression data is to solve the following optimization problem: arg min o∈Π d(o, D). (6) The above problem (6) is a NP-hard problem, 48 since the unknown order is one among #Π = (n−1)! possible orders. An optimization problem that resembles (6) is formulated in Barragán et al, 49 although in this latter case the statistical problem is one dened in a Circular space. We obtain an approximate solution to the optimization problem (6) by formulating it as a traveling salesman problem (TSP) as follows. The data on gene k are represented by a weighted directed graph where the nodes represent the items (or points) to be ordered. Each pair of nodes (i, j) is connected by an edge of length Lk ij that represents the intensity of the relationship or the distance between i to j in gene k . The information is aggregated resulting in a matrix L of aggregated edge lengths, Lij =PK k=1 δkLk ij. Some details on the choice of the weights are given below. In the problem at hand, if Lij measures the temporal distance between the data at time points ti and tj, then we propose to use the L1 distance for the expression values, Lij = K P k=1 δk|Yk i−Yk j| for i, j ∈S . Let the binary matrix Γ , satisfying restrictions (i) and (ii) below, represent a tour that goes exactly once through all nodes in the graph, starting and ending at the same node with Γij = 1 i the edge (i, j) is active in the tour. There is an obvious one to one relationship between Γ and orders o among the set of indices S . Therefore, the problem of nding a circular order using the representation of an aggregated directed graph, dened by L , is reduced to nding a tour that goes exactly once through all nodes in the graph, starting and ending at the same node. The tour that minimizes the total length is the solution of the well-known TSP which, in our case, is mathematically formulated as follows, b Γ = arg min ΓX ij ΓijLij (7) restricted to
Larriba et al 7 (i) Γij, is a doubly stochastic matrix, (ii)X i,j∈V Γij ≤ |V| − 1∀V⊂S, |V|>1. We conducted several heuristic procedures to provide a set of approximate solutions for the tour b Γ dened in (7). Among those tours, we choose the one ( ˆo ) that minimizes (5). One may refer to Subsection 4.1 in the Supporting Information for a detailed description of the temporal order estimation methodology including a owchart. The above optimization algorithm is exible as dierent weights can be chosen. For example, in the case of gene expression data, information regarding rhythmicity of a set of genes is often available, for instance because they are known to be rhythmic in other organs or species so that the weights can be assigned accordingly. Moreover, the above procedure can be combined with a previous SVD on the initial data matrix D , as in CYCLOPS, using only the rst eigengenes proposed in Ana et al. 33 In Subsections 3.2 and 4.2, we compare the performances of ORI and CYCLOPS approaches when determining the temporal orders among time points for simulations and real data, respectively. Validation measures We dene several measures of agreement between a circular order and a data set based on IR. Relative eciency rates ( RRE ) are similar to measures used in linear regression and also resemble those proposed in Ana et al. 33 The measure compares the total sum of squared errors of a given order ( o ), that may correspond to a linear or a circular (up-down-up) relationship, in reconstructing characteristic expression patterns, relative to the total sum of errors as follows: RRET(o) = PK k=1 SREt(o, Yk) PK k=1 SREt(·,Yk), (8) where, SREt(·,Yk) = n P i=1 Yk i−Yk Yk i2 ;Yk=1 nPn i=1 Yik , SREt(o, Yk) = n P i=1 Yk i−Y?k o,i Yk i2 and Y?k o is the IR of Yk under the circular signal model that generates the order o for k= 1, . . . , K . RRE is a positive measure of the percentage of variability not explained by order o . Smaller values of the RRE measure indicate that the order generates estimators that are closer to the observed values suggesting a more reliable order reconstruction. On the other hand, for experiments where the real timing is known, a measure of concordance ( CRE ) between the real order and the circular order o is dened as follows: CRET(o, REAL) = Pn k=1 SREt(o, Y?k REAL) Pn k=1 SREt(·,Y?k REAL). (9) Again, notice that smaller values of the CRE measure indicates a higher concordance among the real timing and the circular orders considered. 3 SIMULATIONS We generate a data set combining four signal shapes that represent real gene expression patterns, called Cosine, Cosine Two, Asymmetric and Flat (see Section 5 in the Supporting Information for full pattern denitions). The rst three ones represent patterns from rhythmic genes and the Flat pattern represents non-rhythmic genes. A data set with 15,000 genes is generated, 20 % corresponding to rhythmic patterns ( 1,000 genes from each of the three rhythmic signals) and the rest to non-rhythmic patterns, imitating real scenarios. Corresponding to each pattern, we simulate data Xj= (X1j, . . . , Xnj)0 for n= 24 time points and j= 1,2 periods using the simulated data set equation Xj∼N24(µ, σ2I) where σ2 is xed to be 1, so that Y∼N24(µ,σ2 2I) . The values of µ are chosen so as to represent the four dierent signal shapes considered. Moreover, following Wu et al 50 and Larriba et al 22 , in this simulation study we sampled every hour for 2 full days (denoted by 1h/2 days), phase shift ( tU ) is chosen from an uniform distribution in [0,24), and a median amplitude level of 2.5 is also xed. We performed exhaustive simulations using a wide range of values of error variance, phase angles and amplitudes. Since the conclusions were
8 Larriba et al similar to the ones summarized here, we do not present those results in the paper. More details on the simulation design are given in Section 5 in the Supporting Information. 3.1 Results for circular signal inference First, we consider the estimation of µ , tU , and tL . In this case, we compare ORI with Cosinor, in terms of the MSE of the estimators. We drive condence intervals (CI) using parametric bootstrap. 39 For these simulations we have computed 500 percentile condence intervals. Each of these intervals is based on 200 bootstrap replications. Since data periodicity is 24 hours, we computed MSE taking into account the equivalence 0≡24 hours. Table 1 illustrates the average MSE for ORI and Cosinor estimators for the signal ( µ ) and for peak and troughs times ( tU , tL ) for each of the three simulated rhythmic patterns. For each of these three patterns, the average MSE is computed over the corresponding genes in the 15,000 genes data set. Apart from the case when the data are generated according to Cosine function, which is the underlying assumption in the Cosinor model, in all other cases, ORI outperforms Cosinor by having smaller MSE (see Table 1). Notice that, although Cosinor gives very good ts for the Cosine function, this methodology is unable to give reasonable ts for data coming from models, as the Cosine Two pattern, that deviate slightly from that function. 50 This fact is quite important in practice where the Cosine model will not hold perfectly in almost any case. In Table 2, we also compared the procedures in terms of estimated coverage probabilities of 95% condence intervals. As expected, neither procedure does a good job of achieving the nominal 95% level for all patterns. However, between the two procedures, ORI performs substantially better by getting coverage probability closer to the true level of 95%. On the other hand Cosinor performs disastrously when the the underlying model is not Cosine shaped function. In fact, in some cases, the coverage probability of Cosinor can be as low as zero. TABLE 1 Mean MSE for µ , tU and tL Cosine Cosine Two Asymmetric ORI Cosinor ORI Cosinor ORI Cosinor µ 0.32 0.06 0.32 0.84 0.23 1.77 tU 2.10 0.08 1.30 4.15 0.07 1.80 tL 2.25 0.08 1.26 4.11 12.12 25.07 TABLE 2 95 % CI bootstrap coverage percentages (average lengths) Cosine Cosine Two Asymmetric ORI Cosinor ORI Cosinor ORI Cosinor tU 91 (3.16) 100 (0.68) 92 (2.46) 0 (0.75) 100 (0.44) 71 (2.01) tL 93 (3.20) 100 (0.69) 94 (2.36) 0 (0.76) 88 (4.19) 0 (2.00) We compared the performance of ORI with the commonly used JTK procedure for testing hypotheses regarding rhythmicity of a gene. Since hypotheses regarding a large number of genes is being performed, to control for false discovery rate (FDR) we applied the Benjamini-Hochberg (BH) procedure. In Table 3 we provide the FDR as well as the false negative rate (FNR), i.e. a gene with a rhytmic pattern is declared to be non-rhytmic. We performed simulations at nominal FDR α= 0.01 . ORI controlled both the false discovery rate (FDR) and the false negative rate (FNR) for dierent patterns of true signal (Table 3). JTK algorithm fails to detect Asymmetric signal patterns with FNR =0.956, while ORI works well in that, although it has a slightly higher FDR value for Flat pattern than expected (0.025 instead 0.01).
Larriba et al 9 TABLE 3 FNR and FDR comparisons at nominal level of α= 0.01 . False Negative Rate False Discovery Rate Cosine Cosine Two Asymmetric Flat ORI JTK ORI JTK ORI JTK ORI JTK 0.000 0.000 0.000 0.000 0.000 0.956 0.025 0.000 3.2 Results for temporal order estimation In this section we consider the problem of temporal order estimation. We compare the performance of 4 dierent methods for estimating the true temporal order. In addition to ORI methodology described earlier, we consider the recently developed neural network based methodology called CYCLOPS. 33 This methodology performs a reduction of dimensionality of the full gene expression data using those eigengenes 34 that contribute to 85% of the total variability. The third methodology, that we denote as SVD 85% +ORI, is a variation to ORI methodology where we apply ORI on the same 85% CYCLOPS eigengenes. The fourth method we consider is LINEAR which temporally orders the data using the rst eigengene. The above four methods are compared in terms of RRE and CRE described in Section 2.4. As for the data sets considered the true temporal order among samples is known, as a measure of comparison, we computed these two metrics for the true real order (REAL). Notice that for these cases, CRE is obviously zero. In (8) and (9) we dened these agreement ( RRE ) and concordance measures ( CRE ) for the orders derived from the full data set. Since CYCLOPS and SVD 85% +ORI are based on the rst eigengenes gathering for 85% of the variability of the full data set, we also considered, for comparison purposes, similar measures based on these eigengenes, denoted as RRE85% and CRE85% . These measures are fully detailed in Subsection 4.2 of the Supporting Information. Figure 3 shows that the rst eigengene from the simulated data set captures well the rhythm, i.e. the fundamental expression pattern across samples. This fact is revealed by the high RRE values for LINEAR ordering in comparison with circular counterparts (see Table 4) which suggests the existence of a temporal order in this simulated data set. The reconstruction of the temporal order is much better using ORI or SVD 85% +ORI than using CYCLOPS as the latter gives the worst RRE and CRE values in all sets considered, as it is shown in Table 4. In particular, the CRE value is around three times higher for CYCLOPS than for the ORI based solutions. The good performance of ORI and SVD 85% +ORI approaches is also illustrated graphically in Figure S3 in the Supporting Information, where the pattern of dierent simulated rhythmic genes are plotted under four dierent orders. TABLE 4 RRE and CRE values for the entire simulated data set and for the set of the rst eigengenes accounting for 85% of data variability Measure REAL LINEAR ORI SVD 85% +ORI CYCLOPS RRET 0.1971 0.6653 0.1965 0.1986 0.2520 RRE85% 0.1104 0.5197 0.1048 0.0970 0.1626 CRET 0.0000 0.0453 0.0491 0.1471 CRE85% 0.0000 0.0242 0.0406 0.1309 4 REAL DATA APPLICATION We apply the ORI methodology to four well-known data sets in chronobiology, see Hughes et al, 29 Thaben and Westermark 21 and Larriba et al, 22 which are available online at NCBI GEO, (http://www.ncbi.nlm.nih.gov/geo/). The mouse liver and pituitary gland as well as the NIH3T3 cell lines data consisted of 45,101 genes each, whereas the U2OS human cell lines data consisted of 32,321 genes. Each data had 48 time points representing two periods of data, i.e. a sampling frequency of 1 h / 2 days. As in simulation study, results for circular signal estimation and detection, and for temporal order estimation are compared with Cosinor, JTK and CYCLOPS, respectively.
16 Larriba et al 31. Bossé Y, Postma D, Sin D, et al . Molecular signature of smoking in human lung tissues. Cancer Res 2012; 72(15): 3753-3763. 32. Leng N, Chu L, Barry C, et al . Oscope identies oscillatory genes in unsynchronized single-cell RNA-seq experiments. Nat Methods 2015; 12(10): 947-950. 33. Ana R, Francey L, Hogenesch J, Kim J. CYCLOPS reveals human transcriptional rhythms in health and disease. Proc Natl Acad Sci USA 2017; 114(20): 5312-5317. 34. Alter O, Brown P, Botstein D. Singular value decomposition for genome-Wide expression data processing and modeling. Proc Natl Acad Sci USA 2000; 97(18): 10101-10106. 35. Zhang W, Edwards A, Fan W, Zhu D, Zhang K. SvdPPCS: An eective singular value decomposition-based method for conserved and divergent co-expression gene module identication. BMC Bioinformatics 2010; 11. doi: 10.1186/1471-2105-11-338 36. Winfree A. The Geometry of Biological Time . Springer Science Business Media . 2001. 37. Robertson T, Wright F, Dykstra R. Order Restricted Statistical Inference . John Wiley & Sons . 1988. 38. Silvapulle M, Sen P. Constrained Statistical Inference: Inequality, Order and Shape Restrictions . Wiley Series in Probability and StatisticsJohn Wiley & Sons . 2005. 39. Efron B, Tibshirani R. An introduction to the bootstrap . New York : Chapman & Hall . 1993. 40. Peddada S, Harris S, Zajd J, Harvey E. ORIOGEN: Order restricted inference for ordered gene expression data. Bioinformatics 2005; 21(20): 3933-3934. 41. Rueda C, Ugarte M, Militino A. Checking unimodality using isotonic regression: an application to breast cancer mortality rates. Stoch Env Res Risk A 2016; 30(4): 1277-1288. 42. Wylupek G. An Automatic Test for the Umbrella Alternatives. Scand J Stat 2016; 43(4): 1103-1123. 43. Bartholomew D. A test of homogeneity of means under restricted alternatives. J Roy Stat Soc B Met 1961; 23(2): 239-281. 44. Menéndez J, Salvador B. Anomalies of the likelihood ratio tests for testing restricted hypothesis. Ann Statist 1991; 19(2): 889-898. 45. Fernández M, Rueda C, Peddada S. Identication of a core set of signature cell cycle genes whose relative order of time to peak expression is conserved across species. Nucleic Acids Res 2012; 40(7): 2823-2832. 46. Bekker dC, Will I, Hughes D, Brachmann A, Merrow M. Daily rhythms and enrichment patterns in the transcriptome of the behavior-manipulating parasite Ophiocordyceps kimemingiae. PLoS ONE 2017; 12(11): 1-20. doi: 10.1371/journal.pone.0187170 47. Ferrari C, Proost S, Janowski M, et al . Kingdom-wide comparison reveals the evolution of diurnal gene expression in Archaeplastida. Nat Commun 2019; 10(1): 387316. doi: 10.1038/s41467-019-08703-2 48. Bartholdi III J, Tovey C, Trick M. The computational diculty of manipulating an election. Soc Choice Welfare 1989; 6(3): 227-241. 49. Barragán S, Rueda C, Fernández M, Peddada S. Determination of temporal order among the components of an oscillatory system. PLoS ONE 2015; 10(7): 1-14. doi: 10.1371/journal.pone.0124842 50. Wu G, Zhu J, Yu J, Zhou L, Huang J, Zhang Z. Evaluation of ve methods for genome-wide circadian gene identication. J Biol Rhythms 2014; 29(4): 231-242.
Larriba et al 17 51. Larriba Y, Rueda C, Fernández M, Peddada S. A bootstrap based measure robust to the choice of normalization methods for detecting rhythmic features in high dimensional data. Front Genet 2018; 9: 24. doi: 10.3389/fgene.2018.00024 52. Zhang R, Lahens N, Ballance H, Hughes M, Hogenesch J. A circadian gene expression atlas in mammals: Implications for biology and medicine. P Natl Acad Sci USA 2014; 111(45): 16219-16224. 53. Deng Z, Wang L, Hou H, Zhou J, Li X. Epigenetic regulation of IQGAP2 promotes ovarian cancer progression via activating Wnt/β-catenin signaling. Int J Oncol 2016; 48(1): 153-160. 54. Follmann D, Proschan M. A Simple Permutation-Type Method for Testing Circular Uniformity with Correlated Angular Measurements. Biometrics 1999; 55(3): 782â791. How to cite this article: Larriba Y, Rueda C, Fernández MA, and Peddada SD (), Order Restricted Inference in Chronobiology, Stat Med , XX .