Denis Agniel

Boris P Hejblum
email: bhejblum@hsph.harvard.edu

Variance component score test for time-course gene set analysis of longitudinal RNA-seq data

Keywords: Gene Set Analysis, Longitudinal data, RNA-seq data, Variance component testing, Heteroscedasticity

As gene expression measurement technology is shifting from microarrays to sequencing, the statistical tools available for their analysis must be adapted since RNA-seq data are measured as counts. It has been proposed to model RNA-seq counts as continuous variables using nonparametric regression to account for their inherent heteroscedasticity. In this vein, we propose tcgsaseq, a principled, model-free, and efficient method for detecting longitudinal changes in RNA-seq gene sets defined a priori. The method identifies those gene sets whose expression varies over time, based on an original variance component score test accounting for both covariates and heteroscedasticity without assuming any specific parametric distribution for the (transformed) counts. We demonstrate that despite the presence of a nonparametric component, our test statistic has a simple form and limiting distribution, and both may be computed quickly. A permutation version of the test is additionally proposed for very small sample sizes. Applied to both simulated data and two real datasets, tcgsaseq is shown to exhibit very good statistical properties, with an increase in stability and power when compared to state-of-the-art methods ROAST, edgeR, and DESeq2, which can fail to control the type I error under certain realistic settings. We have made the method available for the community in the R package tcgsaseq.

Introduction

Gene expression is a dynamic biological process of living organisms, whose dysfunction and variation can be related to numerous diseases. During the past two decades, gene expression measurements have developed rapidly, thanks to wide dissemination of microarray technology. In recent years, gene expression measurement technology has been shifting from microarrays to sequencing (RNAseq) technology. The higher resolution of RNA-seq technology provides a number of advantages over microarrays, among which are the ability to make de novo discoveries and an increased sensitivity to low-abundance variants [START_REF] Marioni | RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays[END_REF][START_REF] Malone | Microarrays, deep sequencing and the true measure of the transcriptome[END_REF]. Indeed, RNA-seq is not restricted to a predefined set of probes like microarrays but can instead measure the genome in its entirety.

A large body of statistical methods have been developed to analyze microarray data. But as technology for measuring gene expression is transitioning to RNA-seq, new methodological challenges arise. RNA-seq produces count data, while microarray analysis techniques generally assume continuity. Due to their underlying count nature, RNA-seq data are intrinsically heteroscedastic. Various approaches have been proposed to deal with these issues, mostly relying on modeling the underlying count nature of the data through the use of Poisson or negative binomial distributions [START_REF] Marioni | RNA-seq: An assessment of technical reproducibility and comparison with gene expression arrays[END_REF][START_REF] Anders | Differential expression analysis for sequence count data[END_REF]Robinson and others, 2010). Recently, [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] have instead proposed to use normal-based methods to analyze RNA-seq data by explicitly modeling the heteroscedasticity and accounting for it by weighting. All of these methods currently make potentially restrictive assumptions about the distribution of the data.

While most methods for gene expression data focus on univariate differential gene expression analysis, it has been shown that gene set analysis (GSA) can be a more powerful and interpretable alternative [START_REF] Subramanian | Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles[END_REF]Hejblum and others, 2015). GSA uses a priori defined gene sets annotated with biological functions and investigates their potential association with biological conditions of interest. There are many different approaches to GSA, and a GSA method is typically defined by the type of hypothesis tested as well as how information across genes is aggregated. There are two primary approaches to specifying tested hypotheses: self-contained hypotheses and competitive hypotheses [START_REF] Goeman | Analyzing gene expression data in terms of gene sets: Methodological issues[END_REF]. Rahmatallah and others (2016) recently showed that self-contained GSA tests tend to be more powerful and more robust than competitive ones [START_REF] Goeman | Analyzing gene expression data in terms of gene sets: Methodological issues[END_REF]. Furthermore, some GSA tests rely on univariate gene-level statistics as a first step, aggregating them afterwards in a bottom-up enrichment approach. But when signal strength is weak, single-step, top-down GSA methods relying on direct multivariate modeling are better than those enrichment based ones at leveraging the additional power of GSA over differential gene expression analysis (Hejblum and others, 2015).

As costs keep decreasing for RNA-seq experiments, more complex study designs, such as time-course experiments, have become more common [START_REF] Dorr | Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients[END_REF][START_REF] Leng | EBSeq-HMM: A Bayesian approach for identifying gene-expression changes in ordered RNA-seq experiments[END_REF][START_REF] Baduel | Habitat-Associated Life History and Stress-Tolerance Variation in Arabidopsis arenosa[END_REF]. However, very few GSA approaches can properly accommodate and test hypotheses more complex than simple differential expression, such as change over time. The ROAST method [START_REF] Wu | ROAST: Rotation gene set tests for complex microarray experiments[END_REF] is a linear-model-based gene set testing procedure. [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] have proposed to use it in conjunction with their heteroscedasticity weighting method voom on RNA-seq data, and this combination of voom and ROAST has been identified in the recent review by Rahmatallah and others (2016) as one of the top-performing GSA methods for RNA-seq data. Among the eleven methods compared in the review, only three are derived from top-down models incorporating all the genes in the set at once, and ROAST is the only one that is not limited to comparisons of two groups. Additionally, DESeq2 (Love and others, 2014) and edgeR (Robinson and others, 2010;McCarthy and others, 2012), are currently the most prominent approaches used for gene-level differential analysis of RNA-seq data. They both rely on the assumption that gene counts from RNA-seq measurements follow a negative binomial distribution. edgeR can use the ROAST method to propose a built-in framework for self-contained gene set testing, while DESeq2 can only perform gene-wise tests. Nueda and others (2014) have also developed a genewise method for dealing with time-course course RNA-seq data, considering Poisson and negative binomial distributions through generalized linear models.

Another concern, particularly in the longitudinal setting, is adequately accounting for heterogeneity of effects (Cui and others, 2016). GSA approaches often overlook the potential longitudinal heterogeneity within a gene set while such heterogeneity is not infrequent [START_REF] Ackermann | A general modular framework for gene set enrichment analysis[END_REF] and can be of biological interest [START_REF] Prieto | Algorithm to find gene expression profiles of deregulation and identify families of disease-altered genes[END_REF][START_REF] Hu | Detection of deregulated modules using deregulatory linked path[END_REF], especially if considered gene sets correspond to biological pathways or networks. Hejblum and others (2015) showed the potential statistical power gain when accounting for this heterogeneity in longitudinal microarray studies.

In this article, we propose tcgsaseq, a method to analyze RNA-seq data at the gene-set level, with a particular focus on longitudinal studies. We derive a variance component score test, similar to those that have been proposed in other testing situations [START_REF] Lin | Variance component testing in generalised linear models with random effects[END_REF][START_REF] Wu | Rare-variant association testing for sequencing data with the sequence kernel association test[END_REF][START_REF] Huang | Gene set analysis using variance component tests[END_REF], that facilitates testing both homogeneous and heterogeneous gene sets simultaneously. Variance component tests offer the speed and simplicity of standard score tests, but potentially gain statistical power by using many fewer degrees of freedom, and have been shown to have locally optimal power in some situations (Goeman and others, 2006). Inspired by the voom approach from Law and others (2014), we propose to estimate the mean-variance relationship in a more principled way using local linear regression, to account for the inherent heteroscedasticity of the data. Despite this nonparametric step, we demonstrate that the test statistic has a simple limiting distribution that is valid regardless of any model specification. We also propose a permutation version of the test to deal with small sample sizes. Our method is implemented in the R package tcgsaseq, available on the Comprehensive R Archive Network at https://cran.r-project.org/web/packages/tcgsaseq.

Our general approach to GSA in longitudinal RNA-seq studies has three pri-mary advantages over existing approaches. First, unlike ROAST, our variance component approach remains valid under model misspecification and does not rely on ad hoc aggregation of information across genes. Second, unlike previous approaches to variance component testing in microarray data [START_REF] Huang | Gene set analysis using variance component tests[END_REF], our approach can accommodate the intrinsic mean-variance relationship in RNA-seq data, while remaining fast to compute. Third, our test remains powerful even when patients display heterogeneous trajectories over time.

The remainder of the paper is organized as follows. Section 2 describes our variance component score test, while its asymptotic properties are derived in Section 3, and an estimation strategy and practical recommendations are detailed in Section 4. In Section 5, we present numerical studies assessing the potential impact of ignoring the mean-variance relationship in RNA-seq data, and comparing our approach to state-of-the-art methods of gene set analysis. In Section 6, we apply tcgsaseq to the analysis of real data. Final remarks and comments are discussed in Section 7.

2 Variance component score test for longitudinal gene set analysis

Problem setup

Consider r i = (r T i1 , ..., r T iP ) T , r ij = (r ij1 , ..., r ijni ) T a vector of sequence read counts for subject i, that have been mapped to each of P genes measured at times t i = (t i1 , ..., t ini ) T . Meanwhile x i = (x i1 , ..., x iq ) T is a vector of baseline covariates describing experiment design conditions, all measured on individual i. So the full data considered for analysis are n independent realizations from random vectors

V = {V i = (r T i , t T i , x T i ) T } n i=1
. Typically, the counts r i are normalized in some way in a preprocessing step (Hansen and others, 2012). We take y i to be a normalized version of r i . For example, one standard normalization procedure accounts for the library size L it = P j=1 r ijt by computing the log-counts per million as

y ijt = log 2 10 6 × 0.5 + r ijt 1 + L it . (1)
Based on these data, our interest is in identifying gene sets that have longitudinally changing expression patterns. The rest of this section will develop the proposed test statistic and its limiting distribution.

The test statistic

We are interested in testing for longitudinal changes in a pre-specified gene set, which for the purposes of illustration we take to be the first p genes, y i = (y i1 , ..., y ip ). To develop a variance component score test statistic, we start from the following working model, which is a linear mixed effect model [START_REF] Laird | Random-effects models for longitudinal data[END_REF]Fitzmaurice and others, 2012)

y ijt = α 0j + a 0ij + x T i α j + φ T it β j + φ T it ξ ij + ijt (2)
and which can be rewritten as

y i = α 0 + a 0i + X i α + Φ i β + Φ i ξ i + ε i (3)
where ε i ∼ N (0, Σ εi ), a 0i ∼ N (0, Σ a ), α 0 is a pn i × 1 vector of gene-specific intercepts α 0j , a 0i is a pn i × 1 vector of random intercepts a 0ij , X i is a pn i × pq block diagonal matrix with each block consisting of n i rows equal to the fixed effect baseline covariates x i , α is a pq × 1 vector of gene-specific fixed effects α j , Φ i = Φ(t i ) is a pn i × pK block-diagonal matrix encoding the effect of time with each block having t th row φ it = {φ 1 (t it ), ..., φ K (t it )} for some set of K basis functions {φ k (•)}, β is a pK × 1 vector of gene-specific fixed effects of time β j , ξ i is a pK × 1 vector of gene-specific individual-level effects of time, and Σ εi is a pn i × pn i covariance matrix of measurement errors. Note that Σ εi is indexed by i and may depend on the mean of y i . We assume that a 0i ⊥ ε i . It is important to note that, in practice, the above model is unlikely to hold. Fortunately the testing procedure we propose is entirely robust to its misspecification. We are interested in testing the null hypothesis of no longitudinal change in the normalized gene expression for any genes in the gene set H 0 : The mean of y i does not depend on t i (4)

or, that is, E(y i |t i ) = E(y i ) and E(y i |t i , x i ) = E(y i |x i ). Under the model (3), the null hypothesis (4) corresponds to the following null: H 0 : β = ξ i = 0, which tests for both homogeneous and heterogeneous dynamics within the gene sets. In section 3, we demonstrate that the corresponding variance component score test can be written as

Q = q T q, q T = n -1 2 n i=1 y T µi Σ -1 i Φ i Σ 1/2 ν ( 5
)
where Σ i = Σ a + Σ εi , y µi = y i -α 0 -X i α, Σ ν is a working covariance matrix defined in (10), and Σ

1/2 ν is the symmetric half matrix such that Σ

1/2 ν Σ 1/2 ν = Σ ν .
We also show in the Supplementary Material that, given the dimension of q (pK) is small relative to the number of individuals (n), the asymptotic distribution of the test statistic is a mixture of χ 2 1 random variables, Q -→ pK l=1 a l χ 2 1 , where the mixing coefficients {a l } depend on the covariance of q. This asymptotic distribution holds regardless of the distribution of y.

In the end, p-values may be computed by comparing the observed test statistic Q to the distribution of pK l=1 a l χ 2 1 where a l is an estimate of a l . Details are developed further in Section 3 and the Supplementary Material. When the entries of q are correlated -i.e. there is correlation between genes in the gene set -then the degrees of freedom for the test based on Q ( pK l=1 a l ) may be much lower than the degrees of freedom for a similar Wald or score test (pK), yielding more power to detect departures from the null hypothesis.

In practice, the parameters α 0 , α, Σ i , and Σ ν are unknown and must be estimated. In particular, accounting for the mean-variance relationship encoded in Σ i is vital in RNA-seq data. Hence we estimate the test statistic as

Q = q T q, q T = n -1 2 n i=1 (y i -α 0 -X i α) T Σ -1 i Φ i Σ 1/2 ν (6)
and we further demonstrate in section 3 that plugging in standard estimates of Σ ν , α 0 , and α and a nonparametric estimator of Σ i for the estimated test statistic still yields a similar asymptotic distribution:

Q -→ pK l=1 a l χ 2 1 (7)
for mixing coefficients { a l } given in the Supplementary Material. The strength of our approach is that this simple limiting distribution holds even when the model (3) may be misspecified and despite the presence of the nonparametric estimator of Σ i in (6). In this way, we may account very flexibly for the mean-variance relationship in y while maintaining a simple, powerful test statistic that does not require any particular model to hold.

The convergence in (7) relies on the central limit theorem. There may of course be situations where the central limit theorem fails to kick in, either when the number of genes in the gene set p is quite large, or else when the number of individuals n is small. We now discuss two modifications for these scenarios.

Testing for homogeneous gene sets

When the number of genes p is large relative to n, the convergence of q to a limiting normal distribution may be in doubt. In these cases, it might be better to begin from a working model that assumes that all genes share a common trajectory over time as a useful approximation. To wit, we may reduce the number of parameters in the model by taking X i in model (3) to be a pn i × q counterpart of X i , and similarly α to be q × 1, Φ i to be pn i × K, and ξ i to be K × 1 counterparts, respectively, of the original quantities.

The derivation of the test statistic follows precisely the same lines with updated dimensions. This test will be powerful to detect alternatives where all genes in the set develop in a similar longitudinally. However, this test should be expected to have low power to detect situations where some genes, for example, decline over time while other genes increase over time, as in heterogeneous gene sets.

Using permutation

When n is very small (or alternatively when p is relatively large but the homogeneous strategy described above seems unwise), relying on the limiting distribution (7) may not be accurate. In these cases, permutation may be used to estimate the empirical distribution of Q under the null (4). To perform the permutations, we simply shuffle the time labels within each individual to get permuted observations for the ith individual y * ijt = y ijIt where {I 1 , ..., I ni } is a permutation of {1, ..., n i }. Indeed, under the null, observations of a given gene j for a given individual i are exchangeable, regardless of sampling time. A large number B of permutation-based test statistics Q = {Q * 1 , ..., Q * B } can thus be generated where each Q * b is computed using the bth set of permuted data. P-values may then be computed as:

p = 1 B B b=1 I{ Q ≤ Q * b } (8)
3 Properties of the test statistic

In this section, we derive the test statistic and demonstrate its asymptotic distribution assuming all parameters known. We then take up the distribution of the test statistic when estimating all relevant parameters.

Test statistic derivation

Under the working model (3),

y i |a 0i , ξ i , x i , t i ∼ N (α 0 + a 0i + X i α + Φ i β + Φ i ξ i , Σ εi ).
Integrating over the random intercepts a 0i , we can rewrite the model as

y i |ξ i , x i , t i ∼ N (µ i + Φ i θ i , Σ i ) (9)
where µ i = α 0 +X i α denotes time-independent fixed effects,

θ i = (β T 1 + ξ T i1 , . . . , β T p + ξ T ip
T denotes combined effects of time. The test for no longitudinal changes in expression corresponds to the modelbased hypothesis H 0 : θ i = 0. We write θ i as θ i = ην i and we consider the working assumption that the {ν i } are independently distributed such that E(ν i ) = 0 (under the null) and

Var(ν i ) = Σ ν . (10)
Under this working assumption, H 0 :

θ i = 0 is equivalent to H 0 : η = 0.
To obtain the variance component test statistic, rewrite the model as:

y µi = Φ i θ i + ε i for centered outcome y µi = y i -µ i . Then y µi | ν i , ∼ N (Φ i ην i , Σ i )
and the log-likelihood for y µi can be written as

log L(η) = - 1 2 n i=1 log |Σ i | + (y µi -ηΦ i ν i ) T Σ -1 i (y µi -ηΦ i ν i ) . ( 11
)
Because the target of inference is η, we marginalize over the nuisance parameter ν conditional on the observed data to obtain L * (η) = E{L(η)|V} where the expectation is taken over the distribution of ν. We follow the argument in [START_REF] Commenges | Score Test of Homogeneity for Survival Data[END_REF] and note that the score at the null value is 0:

lim η→0 ∂ log L * (η)/∂η = E n i=1 y T µi Σ -1 i ην T i φ i | V = 0. So we instead consider the score with respect to η 2 , lim η→0 ∂ log L * (η)/∂(η 2 ), which is equal to E ∂ log L(η) ∂η η=0 | V 2 + constant + o p (1) (12) = n i=1 y T µi Σ -1 i Φ i Σ ν n i=1 y T µi Σ -1 i Φ i T + constant + o p (1) (13)
Thus, after normalizing by n -1 , the variance component score test statistic can be written Q = q T q, as described in ( 5). If the dimension of q is small relative to the sample size, q is asymptotically normal by the central limit theorem, which means that the limiting distribution of Q is a mixture of χ 2 1 s (see Supplementary Material for more details).

Estimated test statistic limiting distribution

In practice, estimates need to be supplied for many of the quantities in the test statistic to be estimated in (6). Following the argument for the limiting distribution of Q, Q has a limiting distribution of a mixture of χ 2 1 s so long as q has a limiting normal distribution. To establish this, we define

q = n -1 2 n i=1 (y i -µ i ) T T i to be a version of q where Σ -1 i Φ i Σ 1/2 ν is replaced with its limit, T i = lim n→∞ Σ -1 i Φ i Σ 1/2
ν . The quantities q and q are asymptotically equivalent, and the central limit theorem ensures that q has a limiting normal distribution. Therefore the asymptotic distribution of Q is a mixture of χ 2 1 s (see Supplementary Material for details), and the corresponding testing procedure is valid under any data-generating function and robust to misspecification of the working model 3.

Estimation

The test statistic (6) contains the estimated quantities α 0 , α, Σ ν , and Σ i . In this section, we take up practical issues involved in the estimation of these quantities.

Estimating model parameters

In this section, we will discuss estimation of ( α 0 , α, Σ ν ), assuming that we have a consistent estimator of Σ i in hand (we leave discussion of estimating Σ i to the following section). A natural way to estimate ( α 0 , α, Σ ν ), taking into account the heteroscedasticity in y, is to fit a weighted mixed effects model corresponding to

y ijt = α 0j + x T i α j + φ T it θ ij + ijt , (14)
where the random intercepts a 0ij have been integrated out and θ ij is taken to be N (0, Σ ν ). The weights are taken to be w i = diag( Σ i ) -1 . This model corresponds to the model ( 9). In practice, fitting the full mixed effects model for each of many gene sets may be computationally demanding. In these cases, the simpler fixed effects model

y ijt = α 0j + x T i α j + φ T it β j + ijt (15)
may be fit to estimate α 0 and α. Σ ν may not be estimated directly from (15).

In this case, it may be specified using a working estimate or taken to be the identity matrix. In general, correct specification of Σ ν may be difficult and may not have much impact on the test statistic. Furthermore, its misspecification does not affect the validity of the test.

Estimating the mean-variance relationship

One key feature of our approach is to estimate Cov(y i |x i , t i , ξ i ) = Σ i in such a way as to account for the heteroscedasticity in y. It is important to note that accurate estimation of Σ i will likely increase power but will not affect the validity or type I error of the testing procedure.

To approximate the mean-variance relationship in y, we use information from all P genes. We assume that the diagonal elements of Σ i , which we will denote v = {v ijt } = {Var(y ijt |x i , t i , ξ i )}, may be modeled as a function of their means m = {m ijt } = {E(y ijt |x i , t i , ξ i )}. Specifically, v ijt = ω{m ijt } + e ijt for some unknown function ω(•) and errors which follow the moment conditions E(e ijt ) = 0, Var(e ijt ) = τ 2 , τ > 0.

We follow Law and others (2014) in using local linear regression in estimating ω(•). We will first write out the form of the estimator if all parameters were known:

ω n (x) = ijt ˜ ijt (x)v ijt , ˜ ijt (x) = bijt(x) ijt bijt(x) , b ijt (x) = K{(m ijt -x)/h} S n2 (x) -(m ijt -x) S n1 (x) where S nd (x) = ijt K{(m ijt -
x)/h}(m ijt -x) d for some kernel function K(•) and bandwidth h > 0.

The means m and variances v could in principle be estimated using some parametric model like ( 14) for all P genes. Let ( m ijt , v ijt ) denote the estimated mean and variance for y ijt . Then the mean-variance relationship may be estimated as

ω n (x) = ω n (x)| mijt= mijt,vijt= vijt (16)
As with any smoothing procedure, choice of the bandwidth h is paramount in producing a reliable estimator. Standard cross-validation techniques may be used to select h in practice. Then the diagonal entries of Σ i can be estimated as ω n ( m i ).

Practical considerations Under the model ( 14), we could take

m ijt = α 0j + x T it α j +φ T it β j +φ T it ξ ij and v ijt = (y ijt -α 0j -x T it α j -φ T it β j -φ T it ξ ij ) 2
. However, because P n, fitting ( 14) for such a large number of genes is not practical. Instead, the simpler model ( 15) could be used. Then m ijt could be taken to be α 0j + x T it α j + φ T it β j and v ijt to be (y ijt -α 0j -x T it α j -φ T it β j ) 2 . A similar approach was used to model heteroscedasticity in the context of linear regression in [START_REF] Carroll | Adapting for heteroscedasticity in linear models[END_REF].

On the other hand, fitting a local linear regression on all of m and v -a total of P n i=1 n i observations -may be computationally difficult. In order to reduce the number of points used in the nonparametric fit ( 16), one could follow Law and others (2014) and model the mean-variance relationship at the gene level

v j = ω(m j ) + e j . ( 17
)
The gene-level mean may be estimated as

m j = n -1 n i=1 n -1 i ni t=1 α 0j + x T
it α j + φ T it β j and the gene-level variance as

v j = n -1 n i=1 n -1 i ni t=1 (y ijt - α 0j -x T it α j -φ T it β j ) 2 .

Numerical study

In this numerical study section, we illustrate both the importance of accounting for RNA-seq data's heteroscedasticity and the superiority of tcgsaseq in terms of statistical power and robustness. We have performed simulations under two different settings. The first one demonstrates the good behavior of the asymptotic test under various scenarios in synthetic data. The second one focuses on a very realistic situation where sample size is small and gene counts are generated from a distribution of real observed RNA-seq data. Throughout, for the purposes of testing, we take Σ a and Σ ν to be the identity for simplicity. In practice, more precise estimation of Σ a and Σ ν would serve to increase power. Finally, throughout we use Davies's approximation method [START_REF] Davies | Algorithm AS 155: The Distribution of a Linear Combination of chi-2 Random Variables[END_REF] to compute p-values for the mixture of χ 2 s, implemented in the CompQuadForm R package [START_REF] Duchesne | Computing the distribution of quadratic forms: Further comparisons between the Liu-Tang-Zhang approximation and exact methods[END_REF].

Synthetic data

In this subsection, we use synthetic data to illustrate the behavior of our estimator when the distribution of the data is known. We first look at the importance of accounting for the mean-variance relationship in Σ i in tcgsaseq, and then we compare tcgsaseq to other competing methods, including ROAST, edgeR-ROAST, and DESeq2. In order to adapt DESeq2 to self-contained gene set testing, we use the minimum p-value test to adequately combine univariate p-values for testing a whole gene set while taking into account gene correlation [START_REF] Moskvina | On multipletesting correction in genome-wide association studies[END_REF][START_REF] Lin | Kernel machine SNP-set analysis for censored survival outcomes in genome-wide association studies[END_REF] and we refer to it as DESeq2-min test. We demonstrate that the tcgsaseq testing procedure is robust to even heavy misspecification of the mean model and the mean-variance relationship, while ROAST, edgeR-ROAST and DESeq2-min test may suffer from extreme lack of power compared to tcgsaseq or inflate type I errors under misspecification.

Mean-variance relationship

To illustrate the importance of estimating the mean-variance relationship, we generated data under the model y ijt = log (µ ijt + 0.5)10

6 i µ ijt + 1
where

µ ijt = η ijt i,t η ijt /P + b ij t it + βt it , η ijt = a ij + α j + 3 k=1 x itk + ijt , (18)
with ijt ∼ N (0, α 2 j ), α j ∼ exp(100), a ij ∼ N (0, α 2 j /100), x itk ∼ N (100, 2500), t it ∼ U (0, 1), b ij ∼ N (0, 1). We set n = 200, p = 10, P = 1000, n i = 3, and β was allowed to be set to 6 different values within a range from 0 to 2: (0, 0.5, 0.75, 1, 1.25, 1.5, 2). We considered three methods of accounting for the mean-variance relationship: (i) gene-level estimates using equation ( 17); (ii) voom-type estimation Law and others (2014); (iii) specifying Σ i to be the identity. Throughout this section we test at the 0.05 level. The model ( 15) was used to generate the test statistic and adjust for covariates, and the gene-level equation ( 17) was used to estimate the mean-variance relationship. Note that the mean model ( 15) is misspecified.

When applying our variance component test, type I error was at the nominal level using either strategy (i) or (iii), but was deflated to 0.035 when using strategy (ii). On the other hand, ROAST, which similarly relies on estimating the mean-variance relationship, inflated type I error to 0.56 when using (ii), and had adequate type I errors using the other weighting strategies. Such results indicate that ROAST may drastically inflate type I error rates in some cases, and that using (ii), the voom-type estimator, to account for heteroscedasticity led to worse performance in general.

Figure 1 demonstrates the effect of accounting for heteroscedasticity under the alternative hypothesis. The more accurate modeling of the mean-variance relationship with method (i) yields large power gains over method (ii). Conversely, by drastically mis-modeling the mean-variance relationship, the naive method (iii) yields almost no power for detecting longitudinal changes.

Comparison to competing methods

In this subsection, we compare tcgsaseq to ROAST, edgeR-ROAST, and DESeq2min test using a highly misspecified model and using a negative binomial model. For the misspecified model, we again generate data under model ( 18) with one modification to ensure that the data are positive integers: y ijt = max log β took values between -2 and 2. Because DESeq2 is computationally intractable at higher sample sizes with large P , we set the total number of genes P to be 100. We used model ( 17) to compute the mean-variance relationship for both tcgsaseq and ROAST. Gene-based estimates were used for dispersions in DE-Seq2, and likelihood ratio tests were used to produce test statistics for both edgeR and DESeq2.

The results are depicted in figure 2. We see that tcgsaseq has the highest power at all sample sizes and at all values of β. Though the model is highly misspecified for all methods, the negative-binomial-based methods edgeR-ROAST and DESeq2min test suffer greatly in terms of power, in particular DESeq2min test. In Figure 2, we show ROAST using the mean-variance relationship estimated from model ( 17), but it's important to note that if the voom-type strategy were used, the type I error for ROAST rises to more than 0.1 at all sample sizes. Asymmetrical power curves around the null are expected for both edgeR and DESeq2 as they both assume the data have a negative binomial distribution (which is asymmetric with a heavier right tail). Secondly, we generated data under the negative binomial model, a distributional assumption that edgeR and DESeq2 make. The mean of the negative binomial was specified as

µ ijt = max{1001 + a 0i + x i + (b i + β + β j )t it x i , 0} (19)
where a 0i ∼ N (0, 1), x i ∼ N (µ xi , 1), µ xi ∼ exponential(1/10), b i ∼ N (0, 1), β j ∼ N (0, 1), and the negative binomial dispersion parameter was set to 1. We again let β vary between -2 and 2 and considered sample sizes of n = 50, 100, 150. In this setting, the default procedure to estimate dispersions was used for DESeq2.

The results depicted in Figure 3 show that tcgsaseq and ROAST outperform edgeR-ROAST and DESeq2-min test despite the fact that the distributional assumptions of edgeR and DESeq2 are true. The performance of tcgsaseq and ROAST is in general comparable. Furthermore, it is important to highlight that DESeq2-min test does not control the type I error, and the type I error gets worse as sample size increases. To reinforce this point, we generated another dataset in which the mean model is correctly specified (which might be more amenable to ROAST), but with an error variance that depends on time. Specifically,

y ijt = a ij + 5x i + β j t it + ijt τ j t it , (20)
where a ij ∼ N (0, 1),

x i ∼ N (µ i , 1), µ i ∼ exp(10), β j = I{j = 1}, t ij ∼ U (0, 1),
ijt ∼ logistic, and τ j ∼ exp(10). In this case, applying the ROAST method using strategy (ii) again produced a wildly inflated type I error of 0.95, while using method (iii) also produced inflated type I errors of 0.10. On the other hand, the type I error rate was maintained at the nominal level using method (iii) for tcgsaseq but was highly deflated (0.003) using method (i).

Realistic small samples simulations

In this subsection, we simulated another dataset to illustrate the good behavior of the tcgsaseq permutation test in the realistic setting of a small sample size. We generated data for 6 individuals each measured at 3 time points according to the scenario described in Law and others (2014), using the script provided in their Supplementary Material:

m ij = 0.2 + 1/ µ ij + βt it and y ijt = log (a ij + 0.5)10 6 j a ij + 1 (21)
where

t it ∈ {1, 2, 3}, a ij ∼ Pois(λ ij ), λ ij ∼ Gamma( κij 40mij , 40mij κij µij ), κ ij ∼ χ 2 (
40), and µ ij follows an empirical baseline distribution derived from real RNA-seq counts data, provided in supplementary material of [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF]. This simulation scheme allows us to generate data that realistically resemble real RNA-seq count data. We set n = 6 with n i = 3, and gene sets were constructed such that p ∈ (30, . . . , 400) and for every gene pair (j, j ) in the set cor(µ j , µ j ) > 0.8. β was allowed to be set to 8 different values within a range from 0 to 0.5: (0, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5).

As shown in Figure 4, under this setting, edgeR-ROAST, ROAST and the proposed permutation test all control the type I error at nominal rate. However DESeq2-min test fails to control the type I error, while it is highly deflated for the asymptotic version of tcgsaseq (0.001) which is not surprising given that with such a small sample size it is unlikely the central limit theorem would have come into effect yet. Nonetheless, we observe a steady and consistent increase of power for our method using the permutation test over the compared state-of-the-art approaches. The deviation from the negative binomial distribution hypothesis (which both edgeR and DESeq2 rely upon), as well as the univariate genewise step first needed, explain those poor performances compared to tcgsaseq. Similar results are obtained when adding a random effect of time simulating heterogeneous gene sets (see Supplementary Material). In addition, as gene sets were constructed with correlated genes, these results also indicates that tcgsaseq is robust to and can take advantage of inter-gene correlation.

Application to two real datasets

In this section we present analyses of two real datasets. In both analyses, Σ a and Σ ν were again taken as the identity while Σ i was estimated at the gene level using equation ( 17) for tcgsaseq. The ROAST method was applied in combination with precision weights estimated through the voom approach. Given the sample sizes available in both studies, we used the tcgsaseq permutation test. In both cases, only transformed data were available and not original counts data, preventing us from applying either edgeR-ROAST or DESeq2-min test.

6.1 Longitudinal RNA-seq measurements in successful kidney transplant patients

We analyzed a RNA-seq dataset from [START_REF] Dorr | Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients[END_REF], in which gene expression was measured in the peripheral blood mononuclear cells of 32 kidney transplant patients. Gene expression was measured at 4 time points: before transplantation, 1 week after transplantation, 3 months after transplantation and 6 months after transplantation. The patients had no graft rejection at the time of each sample. We investigated custom gene sets targeted specifically towards kidney transplant. Those kidney-oriented gene sets were derived by the Alberta Transplant Applied Genomics Center from specific pathogenesis-based transcripts [START_REF] Halloran | The molecular phenotype of kidney transplants: Personal viewpoint[END_REF][START_REF] Sellarés | Molecular diagnosis of antibodymediated rejection in human kidney transplants[END_REF][START_REF] Broin | A pathogenesis-based transcript signature in donor-specific antibody-positive kidney transplant patients with normal biopsies[END_REF], and their definition is available at http://atagc.med.ualberta.ca/Research/ GeneLists. We tested for a linear change in gene expression over time while adjusting for patient's age and gender.

Figure 5 shows the p-values for those 9 gene sets for both tcgsaseq and ROAST. At a 5% threshold, our approach tcgsaseq identifies three significant gene sets while the combination of voom and ROAST identifies none. Among those significant gene sets, two relate to T-cell gene expression, corroborating the original results from [START_REF] Dorr | Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients[END_REF]. The gene set annotated as "T-cells" notably includes the gene CD3D, previously highlighted in [START_REF] Dorr | Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients[END_REF]. On the other hand, gene sets related to transplant rejection such as "Donor-specific antibody" or "Gamma-IFN and Rejection" have much higher p-values here, which is what one would expect as these data include only successful transplant patients. Besides, tcgsaseq also detects a significant change in expression for the gene set related to Mast cells, which have recently been highlighted as playing an ambiguous role in kidney transplant [START_REF] Papadimitriou | Mast cell quantitation in renal transplant biopsy specimens as a potential marker for the cumulative burden of tissue injury[END_REF]. Mast cells have been linked both to peripheral tolerance (de Vries and others, 2009), as well as to late graft loss [START_REF] Jevnikar | Late kidney allograft loss: What we know about it, and what we can do about it[END_REF]. These results both reinforce and broaden the original findings from [START_REF] Dorr | Differentially expressed gene transcripts using RNA sequencing from the blood of immunosuppressed kidney allograft recipients[END_REF]. (2016) were especially interested in the difference between two populations of Arabidopsis arenosa, respectively denoted KA and TBG, that have adopted different flowering strategies. In addition, half of the plants were exposed to cold and short day photoperiods (vernalization) between week 4 and week 10 in order to study the corresponding effect on flowering in both populations. Different plant siblings were sampled at each time point in order to avoid the important stress effect of leaf removal on the plants.

Dealing with this complex experiment design, we used tcgsaseq to address two separate biological questions: i) which gene sets have a different activation between the two populations, adjusted for the plant age and the cold exposure, ii) which gene sets have a different activation due to cold exposure, adjusted for population differences and plant age. Using two different modeling strategies, we found that one data-driven gene set constructed by [START_REF] Baduel | Habitat-Associated Life History and Stress-Tolerance Variation in Arabidopsis arenosa[END_REF] from the top 1% differentially expressed genes between the two populations was significant at a 5% threshold for i) and for ii). On the contrary two other datadriven gene sets, again constructed by [START_REF] Baduel | Habitat-Associated Life History and Stress-Tolerance Variation in Arabidopsis arenosa[END_REF] and identified as the respective population-specific responses to cold exposure, were both significant at this 5% threshold for ii) but not for i). In addition, two gene sets from Gene Ontology associated with salt and cold response pathways, respectively, were also investigated and found significant at the 5% level for ii) but not for i). In comparison, ROAST gave similar results for i), but lacked power for ii) identifying only 1 out of 4 significant gene sets. This analysis corroborates the results obtained by [START_REF] Baduel | Habitat-Associated Life History and Stress-Tolerance Variation in Arabidopsis arenosa[END_REF], and it further illustrates the good behavior of tcgsaseq in complex time-course RNA-seq studies.

Discussion

The proposed method detailed in this article constitutes an innovative and flexible approach for performing gene set analysis of longitudinal RNA-seq gene expression measurements. The approach relies on a principled variance component score test that accounts for the intrinsic heteroscedasticity of RNA-seq data, and for which we derive a simple limiting distribution without requiring any particular model to hold. As illustrated in the previous sections, the good performance of the method when applied to various datasets constitutes a major strength of the method. Our numerical study shows the importance of accurately accounting for heteroscedasticity when analyzing RNA-seq data. We also demonstrate the robustness of our testing procedure to model misspecification. When comparing our proposed approach to ROAST, edgeR, and DESeq2, state-of-the-art methods, we illustrate the competing methods' sensitivity to model misspecification, failure to control Type I error, and reduced statistical power compared to tcgsaseq. Goeman and others (2006) derived a multi-dimensional score-test that has locally optimal power given that i) all of the tested parameters are independent, and ii) there are no nuisance parameters. Their second assumption especially is unrealistic in the context of longitudinal RNA-seq studies targeted by this paper. However, our test being related to theirs, it suggests good performance in statistical power can be expected -just as we demonstrated.

Of particular biological interest in longitudinal studies, the proposed approach can test for both homogeneous and heterogeneous gene sets simultaneously. This is especially relevant for gene sets constructed from biological pathways encompassing both regulators and targets, where different expression dynamics are expected across genes. We also provide a homogeneous version of the testing procedure when heterogeneity of gene dynamics within a gene set is not of interest. It also results in a power gain compared to GSA methods sensitive only to homogeneous sets. In addition, the proposed solution for accounting for heteroscedasticity can also deal with non-count data and could in principle be applied widely.

As demonstrated by the simulation studies, the proposed approach is also robust to inter-gene correlation within tested gene sets. While in our simulations we do not estimate inter-gene correlations, one could account for them more formally by following [START_REF] Wang | A unified mixed effects model for gene set analysis of time course microarray experiments[END_REF] and estimating a working correlation matrix with the residuals from an initial gene-wise modeling of y. The resulting estimates could be incorporated into the structure of Σ i to increase power.

This work features 4 important novelties: i) we propose a statistically sound nonparametric estimation of the inherent mean-variance relationship in RNAseq data; ii) we present an original variance component score test statistic for which we derive an asymptotic distribution, and we show that it possesses the conventional root-n convergence rate despite the presence of a nonparametric estimator and is robust to model misspecification; iii) we conduct a numerical study comparing the proposed method and illustrating its good statistical performances under various settings; iv) we provide an implementation of the proposed method as an R package available to the community.

These innovations contrast with previous works from either [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] or Hejblum and others (2015). While the nonparametric estimation of heteroscedasticity from Law and others (2014) is more of a heuristic procedure, our statistically rigorous approach does not include unmotivated transformations or extrapolation outside the range of the estimating function. Furthermore, we show that our test statistic has the conventional root-n rate despite the presence of the nonparametric estimator. As for Hejblum and others (2015), they focus exclusively on microarray data and therefore completely ignore the heteroscedasticity issue. In addition, while the modeling proposed in this article is similar to theirs, our test statistic is completely different. While they rely on a likelihood ratio test, we propose an original variance component score test that remains valid even if the model was misspecified and is very fast to compute.

Finally, in this paper we focus primarily on longitudinal measurements of RNA-seq. However, our approach directly applies to virtually all RNA-seq study designs, including traditional case-control and more complex studies. Our approach allows researchers to incorporate the natural heteroscedasticity in the data into a powerful test statistic that makes no modeling assumptions. Of note, the inclusion of time-varying covariates would require further assumptions concerning the model to be made to ensure the limiting distribution of the test statistic. Evaluating tcgsaseq's performance in a broader array of studies is an area for future research.

Software

Software in the form of R code is available on the Comprehensive R Archive Network as an R package tcgsaseq.

where t it ∈ {1, 2, 3}, ijt ∼ N (0, 0.05), a ij ∼ P ois(λ ij ), λ ij ∼ Gamma( κij 40mij , 40mij κij µij ), κ ij ∼ χ 2 (40), m ij = 0.2+1/ √ µ ij , γ j ∼ N (0, σ γ ), and µ ij follows an empirical baseline distribution derived from real data, provided in supplementary material of [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF]. We set n = 18 with n i = 3, and gene sets were constructed such that p ∈ (30, . . . , 400) and for every gene pair (j, j ) in the set cor(a j , a j ) > 0.8. σ γ ) varied among (0, 0.02, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3).

The results are very similar to those for homogeneous gene sets, see Figure 6.

Figure 1 :

 1 Figure1: Power evaluation in synthetic data according to how heteroscedasticity is accounted for, based on 1,000 simulations.

Figure 2 :

 2 Figure 2: Power evaluation in synthetic data comparing tcgsaseq, ROAST, edgeR-ROAST, and DESeq2-min test, based on 1,000 simulations.

Figure 3 :

 3 Figure 3: Power evaluation in negative binomial data comparing tcgsaseq, ROAST, edgeR-ROAST, and DESeq2-min test, based on 1,000 simulations.

Figure 4 :

 4 Figure 4: Power evaluation in realistically simulated data with a small sample size, based on 500 simulations.

Figure 5

 5 Figure 5: P-values from testing the 9 kidney oriented gene sets investigated

Figure 6 :

 6 Figure 6: Power evaluation for heterogeneous gene sets in realistically simulated data with a small sample size, based on 500 simulations.

Acknowledgments

This work was supported by the National Institutes of Health (NIH) [U54 HG007963 to B.P.H.]. The authors express their deepest gratitude to Professor Tianxi Cai, Harvard University, for her help and support in this work. They also thank Pierre Baduel for his help in analyzing the plant physiology data.

Conflict of Interest: None declared.

Supplementary Material

S.1 Properties of the test statistic

In this section, we establish the limiting distribution of the quantity Q and its estimated counterpart Q.

S.1.1 Limiting distribution of Q

Recall that Q = q T q. Let Γ = Cov(q). Then

where U is an orthonormal matrix of the eigenvectors of Γ, A is a diagonal matrix of the eigenvalues of Γ, and u = Γ -1/2 q is asymptotically standard multivariate normal by the central limit theorem. Noting that, because U is orthonormal,

k is an element of the asymptotically standard normal U T u and a k an eigenvalue of Γ. It finally follows that Q ∼ pK l=1 a l χ 2 1 .

S.1.2 Limiting distribution of Q

Recall that Q = q T q. We first show that q and q are asymptotically equivalent.

for the processes T (

Here we use the fact that Σ i only depends on x i .

We can thus write

Using that fact and T ( x) → p T ( x) it can be shown that

This means that q and q are asymptotically equivalent, and we can examine the limiting behavior of q to understand the limiting behavior of q. Now, ignoring the intercept α 0 for notational simplicity,

where U i (α) corresponds to the influence function for α from a null parametric model and

Thus, clearly q has a limiting normal distribution by the central limit theorem, which further suggests that the asymptotic distribution of Q is a mixture of χ 2 1 s where the mixing probabilities are given by the eigenvalues of Cov(Q i ).

S.2 The voom procedure

Let r ijt be the read count from individual i for gene j at time t. L i = G j=1 r ijt is then the library size, that can vary from one sample to another. [START_REF] Law | voom: Precision weights unlock linear model analysis tools for RNA-seq read counts[END_REF] propose the following transformation to obtain log-counts per million value from sample i for gene g:

Counts are offset away from zero by 0.5 to avoid taking the log of zero (and this also reduces the variability of log-cpm for low expression genes) while the library size is offset by 1 to ensure that any log-cpm is below 1 million. Let X be the design matrix of factors that are suspected to influence gene expression.

The voom procedure to estimate mean-variance weights proposed by Law and others (2014) is the following:

1. α is the OLS estimate from the linear model:

-log 2 (10 6 ) be the average log-count value.

4. f (•) is the predictor obtained from the LOWESS regression [START_REF] Cleveland | Robust Locally Weighted Regression and Smoothing Scatterplots[END_REF] of s j over rj . The final precision weights are then defined as follows:

S.3 Additional simulations

Here we provide additional simulation results in the realistic setting of a small sample size. We generated data for 6 individuals, each measured at 3 time points according to the scenario described in Law and others (2014), using the script provided in their Supplementary Material:

y ijt = log (a ij + 0.5)10 6 j a ij + 1 + γ j t it + ijt (36)