Poisson–Tweedie mixed-effects model: A flexible approach for the analysis of longitudinal RNA-seq data

In the last decade, RNA-sequencing (RNA-seq; Mortazavi et al. (2008)) has progressively replaced microarrays and has become the standard method to measure gene expression (DE). A typical goal of RNA-seq studies is the identification of genes that are differentially expressed across groups of subjects, or are associated with a variable of interest. To this end, a statistical model is postulated to relate the expression levels of a gene to the variable(s) of interests and relevant confounders, and hypothesis testing is performed to assess the significance of the relevant differences or associations. While microarrays yielded continuous measurements that could be analysed using models that relied on normality assumptions (Smyth, 2004), RNA-seq measurents are discrete, and as such they require statistical models suitable for the analysis of count data. While the very first attempts to model RNA-seq data relied on the use of the Poisson distribution (Marioni et al., 2008), it was soon observed that the distribution of RNA-seq counts is typically characterized by phenomena such as overdispersion (the variation of counts is higher than the mean), zero-inflation and heavy tails that the Poisson distribution cannot account for.

To take overdispersion into account, Robinson et al. (2010) and Love et al. (2014) proposed two statistical models based on the negative binomial (NB) distribution, implemented in the R packages edgeR and DESeq2. To deal with small sample sizes, which were a very common feature of the first cross-sectional RNA-seq studies, they proposed to estimate the dispersion parameter of the NB distribution by borrowing information across genes. Esnaola et al. (2013) tested the adequacy of the NB distribution on data from a large-scale RNA-seq experiment. They showed that although the NB fitted sufficiently well the empirical distribution of counts for many genes, such distribution was inadequate for a relevant proportion of genes that exhibit strong zero-inflation or heavy tails. Therefore, they proposed a more flexible generalized linear model (GLM) based on the Poisson–Tweedie (PT) distribution, implemented in the R packages tweeDEseq. The PT is a three-parameter distribution that encompasses the NB as a special case, and differently from the NB it can accommodate different levels of zero-inflation and heavy-tailness for a given amount of overdispersion. A different approach, called limma-voom, was proposed by Law et al. (2014), who combined a linear model for the log-transformed counts with sample-specific precision weights, making it possible to apply the limma pipeline, originally developed for microarray data, to RNA-seq measurements.

At the beginning of the diffusion of RNA-seq, the high sequencing costs discouraged the set-up of longitudinal RNA-seq studies; for this reason, limma-voom, edgeR, DESeq2 and tweeDEseq were developed to deal with cross-sectional RNA-seq datasets. However, interest in the dynamic evolution of gene expression and declining sequencing costs are nowadays motivating more longitudinal studies. One such example is a mice experiment carried out at the Leiden University Medical Center (LUMC) that aimed to identify biomarkers of disease progression for Duchenne muscular dystrophy (DMD). In this experiment, described more in detail in Section 4, blood samples were collected from healthy and dystrophic mice every 6 weeks over a period of 30 weeks, and gene expression in blood was quantified using RNA-seq.

In longitudinal studies, repeated measurements are collected from the same subject at several time points. This introduces correlation between the repeated measurements from the same subject, a correlation that is not present in cross-sectional designs where all samples can be assumed to be independent. Nevertheless, often longitudinal RNA-seq data are analysed using popular approaches (mostly edgeR, DESeq2 and limma-voom) for cross-sectional RNA-seq data, either ignoring the longitudinal nature of the data and assuming samples from the same individual to be independent, or adding for each individual a subject-specific indicator as a covariate (fixed effect). Cui et al. (2016) showed that modelling longitudinal RNA-seq data using such approaches results into increased false positive rates (FPR), or in a loss of power when the FPR is controlled. They argued that longitudinal RNA-seq data should instead be analysed using generalized linear mixed models (GLMMs, McCulloch et al., 2008), and they illustrated their point by comparing those ‘improper’ modelling approaches to a Poisson GLMM. They motivated their choice to use the Poisson distribution instead of a more flexible alternative based on frequent convergence problems that they encountered when attempting to fit a NB GLMM with random intercept. Although their choice of a Poisson GLMM simplifies model estimation, it is clear that such a model cannot account for the different combinations of overdispersion, zero-inflation and heavy tails typical of RNA-seq counts. Therefore, the use of GLMMs based on more flexible distributions for count data, such as the NB and PT distributions, should be preferred even if their estimation is more involved.

Recent proposals tailored to the analysis of longitudinal RNA-seq data include MaSigPro (Nueda et al., 2014), ShrinkBayes (van de Wiel et al,. 2014), DyNB (Äijö et al., 2014) and timeSeq (Sun et al., 2016). MaSigPro fits a negative binomial GLM (NB GLM) that lets the expected value of a gene depend on time polynomials and their interaction with a group indicator; similarly to the methods for cross-sectional RNA-seq, this approach relies on a simple GLM that does not model properly the repeated measurements design. ShrinkBayes can fit Poisson, NB, zero-inflated Poisson (ZIP) and zero-inflated NB (ZINB) GLMMs that can accommodate a random intercept. Owing to the possibility to choose the ZIP and the ZINB as distributions, this approach can be employed to model both zero-inflation and the correlation between repeated measurements. However, ShrinkBayes does not allow to fit GLMMs for heavy-tailed distributions. Moreover, the method requires the user to preselect one specific distributional shape that is fitted to all genes, rather than letting the estimation procedure determine it automatically and separately for each gene as in the case of the PT family of models proposed by Esnaola et al. 2013. Finally, often interest may lie in assessing the joint significance of two or more coefficients, for example that of a main effect and of the corresponding interaction term, or in testing a system of linear combinations of the regression coefficients (an example of this type of hypotheses can be found in Section 4). However, the R package ShrinkBayes allows users to test only hypotheses that involve a single regression coefficient (although with some additional programming one may still be able test more complex hypotheses). DyNB fits a semiparametric model that combines a NB GLM with Gaussian-processes regression. Finally, timeSeq considers a NB GLMM where time is modelled using splines. Both DyNB and timeSeq are semiparametric and restricted to the NB case, thus they do not attempt to model zero-inflation and heavy tails.

In this article we present ptmixed, a PT GLMM that is suitable for the analysis of longitudinal RNA-seq data and that, more in general, can be employed to flexibly model longitudinal counts that are overdispersed, zero-inflated and/or heavy-tailed. Owing to the use of the PT distribution, our approach is capable to capture different combinations of zero-inflation and heavy-tails that typically characterize the distribution of overdispersed counts, and it automatically performs model selection among a range of different discrete distributions separately for each gene, similarly to tweeDEseq (Esnaola et al., 2013). Moreover, our mixed-effects approach allows to properly model within-subject correlations (Cui et al., 2016) and to flexibly test both one-dimensional and multidimensional statistical hypotheses. To the best of our knowledge, our proposal is the first extension of the PT GLM of Esnaola et al. (2013) to the context of random effects models, while a different extension based on generalized estimating equations has been proposed by Bonat et al. (2018). We provide an R implementation of functions for the estimation of the proposed model and related hypothesis testing through the R package ptmixed (Signorelli et al., 2020), which can be freely downloaded from the Comprehensive R Archive Network (CRAN).

The remainder of the article is organized as follows: In Section 2 we introduce our model, describe the use of the adaptive Gauss-Hermite quadrature rule to approximate its likelihood, present our computational approach to maximum likelihood (ML) estimation and discuss the use of Wald and likelihood ratio tests (LRTs) for hypothesis testing. In Section 3 we assess the performance of our method using simulated data, and we compare it to several alternative methods that are commonly used in the analysis of longitudinal RNA-seq studies. In Section 4 we illustrate an application of ptmixed to the data from our motivating experiment. We conclude by discussing advantages and limitations of our approach in Section 5.

In this section, we evaluate the performance of ptmixed through several simulations. We begin by evaluating the accuracy of our approach for ML estimation and hypothesis testing across a range of different distributional shapes (from zero-inflated to heavy-tailed) and sample sizes (Section 3.1). In Section 3.2 we show the added value of approximating the likelihood function with the AGHQ method by comparing it to the LA. Finally, in Section 3.3 we compare the performance of our method to that of several state of the art approaches for the analysis of RNA-seq data. R scripts to reproduce the results presented here are available at https://github.com/m-signo/ptmixed.

3.1 Accuracy of ML estimates and hypothesis testing

To assess the accuracy of the ML estimates and of the Wald and LR tests computed following the procedures outlined in Section 2.3 and 2.4 we simulated data from the PT mixed model in (2.1) considering five different values of the power parameter $a$, so as to cover a wide range of distributional shapes: zero-inflated ($a=-5,-1$), NB ($a=0$), heavy-tailed ($a=0.5$) and Poisson ($a=1$). We fixed the dispersion parameter $D=3$ (except for the case $a=1$, where $D=1$) and the number of repeated measurements per subject $m_i=5 \forall i\in \{1,...,n\}$, and progressively increased the number of subjects $n\in \{10,25,50,100,150\}$. We splitted individuals into two groups, denoted by $d_i\in \{0,1\}$, of equal size. We let the mean ${\mu }_{\mathit{ij}}$ depend on a subject-specific random intercept $v_i$, on group $d_i$ and time $t_i$ in such a way that $\text{log}({\mu }_{\mathit{ij}})={\beta }_0+v_i+{\beta }_1{d}_i+{\beta }_2{t}_{\mathit{ij}}$, with $v_i\sim N(0,{\sigma }^2=0.5)$.

In simulation A, we let $\beta =(2.5,0.3,0)$ and test whether the time effect is null, that is, $H_0:{\beta }_2=0$. Table 2 reports the root mean square error (RMSE) of ${\widehat{\beta }}_2$ for different values of $a$ and $n$. We can observe that the magnitude of the RMSE is comparable across different values of $a$, and that it decreases linearly with $\sqrt{n}$ as expected. Table 3 compares the FPR at $\alpha =0.05$ between the Wald and LR test for each of the distributions considered and increasing sample size; we observe that the two tests exhibit similar performance, producing an FPR close to the nominal type I error ($\alpha =0.05$) already with small sample sizes.

Table 3: False positive rates of the Wald and likelihood ratio tests for different values of $a$ and $n$ in simulation A ($H_0:{\beta }_2=0$, $\alpha =0.05$). Values in each cell are estimated using 5 000 random replicates

Table 3: False positive rates of the Wald and likelihood ratio tests for different values of $a$ and $n$ in simulation A ($H_0:{\beta }_2=0$, $\alpha =0.05$). Values in each cell are estimated using 5 000 random replicates

View larger version

Table 4: RMSE of ${\widehat{\beta }}_1$ in simulation B. Values in each cell are estimated using 5000 random replicates

Table 4: RMSE of ${\widehat{\beta }}_1$ in simulation B. Values in each cell are estimated using 5000 random replicates

View larger version

In simulation B, instead, we let $\beta =(2.5,0,0.2)$ and test whether the group coefficient is null, that is, $H_0:{\beta }_1=0$. Also in this case, we see that the RMSE of the parameter of interest (here ${\widehat{\beta }}_1$) is of comparable magnitude across different values of $a$ and that as expected it decreases linearly with $\sqrt{n}$ (Table 4). Assessment of the FPR at $\alpha =0.05$ (Table 5) points out that for very small sample sizes (namely $n=10$) both the Wald and LR test are anti-conservative. For larger sample sizes, both methods achieve FPRs much closer to $\alpha =0.05$.

3.2 Comparison between the AGHQ and Laplace approximations of the likelihood

To quantify the improvement associated to the use of the AHGQ in the approximation of the likelihood of our model, we illustrate a comparison of its performance to that of the LA, focusing on the case with $a=-1$ and $n=100$ considered in simulation B (for the description of the settings of this simulation, see Section 3.1). Table 6 compares the mean of the ML estimates obtained when using the LA or an AGHQ with 5 quadrature points. The results show that even for $n=100$ the LA can yield rather biased parameter estimates, in particular for $D$, $a$ and ${\sigma }^2$. Use of the AGHQ makes it possible to obtain more accurate parameter estimates for all parameters.

Table 5: False positive rates of the Wald and likelihood ratio tests for different values of $a$ and $n$ in simulation B ($H_0:{\beta }_1=0$, $\alpha =0.05$). Values in each cell are estimated using 5 000 random replicates

Table 5: False positive rates of the Wald and likelihood ratio tests for different values of $a$ and $n$ in simulation B ($H_0:{\beta }_1=0$, $\alpha =0.05$). Values in each cell are estimated using 5 000 random replicates

View larger version

Table 6: Evaluation of the bias of the ML estimates obtained using the LA and the AGHQ. Values in each cell are estimated using 5 000 random replicates

Table 6: Evaluation of the bias of the ML estimates obtained using the LA and the AGHQ. Values in each cell are estimated using 5 000 random replicates

View larger version

Importantly, use of the AGHQ produced considerably higher convergence rates (87.5% versus the 65.2% of the LA) using the default starting values computed by ptmixed (thus reducing the risk that one may need to try different starting points to ensure convergence). Furthermore, the AGHQ reduced the proportion of false positives (5.6% versus 7% for the LA, $\alpha =0.05$) enabling a better type I error control. All these performance improvements are due to the fact that the use of more quadrature points allowed a better approximation of the likelihood function, and they required on average a 12% increase in computing time with respect to the LA (Table 7).

3.3 Comparison with alternative modelling approaches

In this section, we describe the results of two simulation studies (C and D) designed to compare ptmixed to several state of the art approaches that are typically employed in the analysis of longitudinal RNA-seq data. Contrary to Section 3.1, where we assessed the performance of our method in parameter estimation and testing for a fixed $\theta =(\beta ,D,a,{\sigma }^2)$, here we evaluate the overall performance of each method on simulated datasets comprising $G=500$ genes each (i.e., each gene has a different ${\theta }_g$); the reason for this is that some of the methods included in the comparison (DESeq2, edgeR and limma-voom) carry out parameter estimation by borrowing information across different genes, and therefore they require information from all genes to estimate and test the parameters associated to each individual gene.

Table 7: Convergence rate, estimated FPR and mean computing time using the LA and the AGHQ. Values in each cell are estimated using 5000 random replicates

Table 7: Convergence rate, estimated FPR and mean computing time using the LA and the AGHQ. Values in each cell are estimated using 5000 random replicates

View larger version

We simulate datasets comprising gene expression values of $G=500$ genes for a sample of $n$ subjects who are splitted into two groups of equal size. We assume that $m_i=5$ repeated measurements are available for each subject. The data are generated from model (2.1) assuming that

$$\text{log}({\mu }_{\mathit{gij}})={\beta }_{g0}+{\beta }_{g1}{d}_i+{\beta }_{g2}{t}_{\mathit{ij}}+o_{\mathit{ij}}+v_{\mathit{gi}},$$

where $d_i\in \{0,1\}$ is the group indicator, $t_{\mathit{ij}}\in \{0,1,...,4\}$ represents time, $o_{\mathit{ij}}\sim N(0,0.5^2)$ is an offset term whose purpose is to mimic the presence of sample-specific normalization weights, $v_{\mathit{gi}}\sim N(0,{\sigma }_g^2)$ with ${\sigma }_g^2\sim U(0.2,0.8)$ is a random intercept and the dispersion parameter $D$ follows a translated gamma distribution, $D-1\sim \Gamma (\alpha =2,\beta =1)$. Consistently with the results of Esnaola et al. (2013), who showed that the hypothesis that the distribution of RNA-seq counts is NB cannot be rejected for most genes, but that at the same time the distribution is not NB for a relevant proportion of genes, we let the distribution of $Y_{\mathit{gij}}$ be NB ($a_g=0$) for 300 genes, zero-inflated ($a_g\sim U(-10,-1)$) for 100 genes and heavy-tailed ($a_g\sim U(0.3,0.7)$) for the remaining 100 genes.

We have included in this comparison approaches that (a) treat repeated measurements from the same individual as independent, (b) include the subject identifiers as fixed effects into the model, and (c) use random effects to model the correlation induced by the repeated measurements design. An overview of the methods considered is given in Table 8, where we summarize the type of model used and the way in which they handle the repeated measurements design. All methods were fitted specifying as fixed effects the main effect of group and time, with the exception of MaSigPro where also an interaction term was included (in the R package MaSigPro the interaction term is automatically included and cannot be suppressed). As concerns hypothesis testing, for ptmixed we assessed the performance of both the Wald test and the LR test introduced in Section 2.4.

Table 8: List of methods compared in simulations C and D

Table 8: List of methods compared in simulations C and D

View larger version

In simulation $C$, we compare the performance of each method in testing if a gene is up- or down-regulated over time, that is, ${\beta }_{g2}=0$. We let ${\beta }_{g0}\sim N(3,0.5)$, ${\beta }_{g1}\sim N(0,0.3)$, and we assume that 400 genes are not dysregulated over time (i.e., ${\beta }_{g2}=0$), 50 genes are up-regulated with ${\beta }_{g2}\sim U(0.1,0.5)$ and 50 genes down-regulated with ${\beta }_{g2}\sim U(-0.5,-0.1)$. In simulation $D$, instead, we check the performance of each method in testing if a gene is DE between the two groups, that is, ${\beta }_{g1}=0$. We let ${\beta }_{g0}\sim N(3,0.5)$, ${\beta }_{g2}\sim U(-0.1,0.1)$, and we assume that 400 genes are not DE (i.e., ${\beta }_{g1}=0$), 50 genes are up-regulated with ${\beta }_{g1}\sim U(0.5,1)$ and 50 genes down-regulated with ${\beta }_{g1}\sim U(-1,-0.5)$.

First, we assess the precision of the parameter estimates obtained by the different methods by comparing their RMSEs. We report the RMSE for all methods with the exception of MaSigPro, for which the computation was not possible because this method reports parameter estimates only for the significant genes, and not for all genes. Then, we evaluate the extent to which each method controls the type I error by computing the FPR, which is the proportion of truly non-DE genes that are declared as DE by the test at a fixed significance level $\alpha$. Furthermore, for the methods whose FPR appears to be under control we compare the True Positive Rate (TPR), which is the proportion of truly DE genes that are correctly declared as DE by the test at level $\alpha$. Finally, we also present a comparison of the computing times of each method. The results of all methods were corrected for multiple testing. We employed the Benjamini–Hochberg (Benjamini and Hochberg, 1995) correction for all methods list in Table 8, with the exception of ShrinkBayes, for which we employed the Bayesian FDR suggested in van de Wiel et al. (2014). Results for each simulation and sample size are based on 50 random replications (i.e., 50 datasets comprising 500 genes each).

Table 9 compares $\mathit{RMSE}({\widehat{\beta }}_{g2})$ in simulation $C$, where interest lies in the estimation and testing of the time effect ${\beta }_{g2}$, which is a within-subjects effect. In this simulation, the method with the lowest RMSEs is ShrinkBayes, closely followed by tweeDEseq when $n=10$ and $n=20$, and by ptmixed when $n=20$ and $n=40$. The remaining six methods are associated to larger RMSEs across all considered sample sizes. Table 10 compares $\mathit{RMSE}({\widehat{\beta }}_{g1})$ in simulation $D$, where interest lies in the estimation and testing of the group effect ${\beta }_{g1}$, which is a between-subjects effect. Here, ptmixed is by far the method with the lowest RMSEs. The RMSEs of ShrinkBayes, tweeDEseq, edgeR, limma and DEseq2 are between 5 and 20 times higher than those of ptmixed, whereas the three approaches that use subject-specific fixed effects have extremely high RMSEs that do not even decrease with $n$. This last result shows that attempts to model longitudinal data using fixed effects can yield extremely inaccurate parameter estimates of between-subjects effects.

Table 9: RMSE of ${\widehat{\beta }}_2$ in simulation C

Table 9: RMSE of ${\widehat{\beta }}_2$ in simulation C

View larger version

Table 10: RMSE of ${\widehat{\beta }}_1$ in simulation D

Table 10: RMSE of ${\widehat{\beta }}_1$ in simulation D

View larger version

Table 11 reports the estimated FPRs in simulation C. Overall, we can observe that in this case almost all methods control the type I error, achieving an FPR $<5\%$ for all sample sizes. The only exception to this is DESeq2+FE. It can also be noted that the FPRs of the two mixed-effects modelling approaches, ptmixed (both with the Wald test and the LR test) and ShrinkBayes, are much closer to the prespecified $\alpha$ than those of all other methods, which appear to be more conservative.

The estimated TPRs for the methods that control the type I error in simulation C are shown in Table 12. In general, all methods achieve comparable TPRs, with the exception of MaSigPro whose TPR is almost 0 for all $n$. The methods with the highest TPR when $n=10$ are the three fixed-effects approaches, whereas for $n=20$ and $n=40$ edgeR+FE and limma-voom+FE are the most powerful, followed by ptmixed, ShrinkBayes and tweeDEseq.

Table 11: False Positive Rate in simulation C ($H_0:{\beta }_{g2}=0,:\alpha =0.05$)

Table 11: False Positive Rate in simulation C ($H_0:{\beta }_{g2}=0,:\alpha =0.05$)

View larger version

Table 12 True positive rate in simulation C ($H_0:{\beta }_{g2}=0,:\alpha =0.05$) for the methods that control the type I error (FPR $\le \alpha$ in Table 11)

Table 12 True positive rate in simulation C ($H_0:{\beta }_{g2}=0,:\alpha =0.05$) for the methods that control the type I error (FPR $\le \alpha$ in Table 11)

View larger version

Table 13 False positive rate in simulation D ($H_0:{\beta }_{g1}=0,:\alpha =0.05$)

Table 13 False positive rate in simulation D ($H_0:{\beta }_{g1}=0,:\alpha =0.05$)

View larger version

The estimated FPRs in simulation D are presented in Table 13. In this case, ptmixed and MaSigPro control the type I error across all sample sizes, and ShrinkBayes does so for $n=20$ and $n=40$. All other methods fail to control the type I error, producing much higher FPRs: the approaches that assume independence (DESeq2, edgeR, limma-voom and tweeDEseq) yield FPRs that range from 0.184 to 0.342, whereas the approaches that employ fixed effects (DESeq2+FE, edgeR+FE and limma-voom+FE) have FPRs close to 1.

Table 14 True positive rate in simulation D ($H_0:{\beta }_{g1}=0,:\alpha =0.05$) for the methods that control the type I error (FPR $\le \alpha$ in Table 13)

Table 14 True positive rate in simulation D ($H_0:{\beta }_{g1}=0,:\alpha =0.05$) for the methods that control the type I error (FPR $\le \alpha$ in Table 13)

View larger version

The estimated TPRs for the methods that control the type I error in simulation D are shown in Table 14. Once again, MaSigPro has TPRs close to 0, giving an indication that the method has little power to reject $H_0$ both in simulation C and simulation D. As concerns ptmixed, the Wald test appears to be a bit more powerful than the LR test, especially when $n=10$ and $n=20$. Finally, in the two cases in which it controls the type I error ($n=20$ and $n=40$), ShrinkBayes achieves a TPR similar to that of ptmixed's Wald test.

In Tables 1 and 2 of Supplementary File 1, we report the computing time of each method in simulations C and D. Overall, estimation of the GLMM approaches, ptmixed and ShrinkBayes, was considerably slower than that of the other approaches that fit simpler GLMs. For most methods, computing time increased linearly with $n$; this did not hold for the fixed-effects approaches, where computing time increased more than linearly with $n$, and for ShrinkBayes, which took more time to convergence when $n=10$ than when $n=20$ or $n=40$.

In this section, we illustrate how the PT mixed model can be employed to analyse longitudinal RNA-seq data. We consider a longitudinal RNA-seq experiment performed at the LUMC (Leiden, NL) whose aim was to characterize the dynamic evolution of gene expression in mdx mice. The mdx mouse is the most used mouse model for DMD, an X-linked recessive genetic disease caused by protein truncating mutations in the DMD gene that encodes dystrophin. The progression of DMD is characterized by a process of muscle degeneration and regeneration whose consequences are muscle instability and the replacement of muscular tissue with fibrotic and fatty tissues. Eventually, this process leads to loss of muscle function, loss of ambulation and premature death. To date no effective treatment for DMD has been found; one of the reason for drug failure during clinical trials is the use functional tests with high inter- and intra-individual variability, resulting in underpowered studies. The possibility to quantify disease progression non-invasively with objective biomarkers is highly needed. In this experiment, repeated blood measurements were taken in order to obtain estimates of the evolution of blood-based gene expression over time.

Table 15: List of pathways considered in our analysis. The pathways were downloaded form WikiPathways (Slenter et al., 2017)

Table 15: List of pathways considered in our analysis. The pathways were downloaded form WikiPathways (Slenter et al., 2017)

View larger version

The experiment involved five different groups of mice, but for simplicity of illustration here we focus on the comparison of the two most interesting groups: wild-type (wt), comprising 5 healthy mice, and mdx, comprising 3 dystrophic mice that shared the same genetic background of wt mice but carried a nonsense mutation in exon 23 that led to lack of dystrophin. Mice were kept under observation for 30 weeks, during which blood samples were collected from each mouse every 6 weeks starting from 6 weeks of age. Transcriptomic expression in these samples was quantified using RNA-sequencing (GEO accession number GSE132741). We focus our attention on those genes that belong to 5 pathways, listed in Table 15, that are relevant in the DMD pathophysiology. We filter out lowly expressed genes by excluding genes that have less than 1 count per million in more than half of the samples; after filtering, 379 genes are retained in the analysis.

We let the expression value for gene $g$ from mouse $i$ at the $j$-th time point follow a PT GLMM where $Y_{\mathit{gij}}\mid v_{\mathit{gi}}\sim \mathit{PT}({\mu }_{\mathit{gij}},:D_g,:a_g)$ and

where $d_i$ is a dummy variable which is 0 for wt and 1 for mdx mice, $t_{\mathit{ij}}\in \{0,1,2,3,4\}$ denotes time (0 = week 6, 4 = week 30) and $o_{\mathit{ij}}$ is an offset term equal to the logarithm of the effective library sizes computed using the TMM normalization method (Robinson and Oshlack, 2010).

We estimated model 4.1 using an AGHQ with 10 quadrature points and considered up to two different starting values (ptmixed default values and, when necessary, ML estimates of the NB GLMM) for the optimization. ML estimation was successful for 315 genes; for the remaining 64 genes, for which estimation of model 4.1 failed, we attempted to estimate the simpler NB mixed model by fixing $a_g=0$; estimation of such model was successful for 63 of those genes. We do not consider further the remaining gene, for which also estimation of the NB GLMM failed.

Download

Open in new tab

Download in PowerPoint

Figure 1: Histograms showing the distribution of the ML estimates of the power parameter. The left histogram shows the distribution of the ML estimates over the full range of ${\widehat{a}}_g$, whereas the right histogram zooms on the interval $[-5,1]$, where most of the distribution is concentrated

Figure 1 shows the distribution of the ML estimates of the power parameter, ${\widehat{a}}_g$, for the 315 genes for which estimation of the PT GLMM was successful. For 74 genes ${\widehat{a}}_g<-1$, for 62 genes ${\widehat{a}}_g\in [-1,0)$, for 106 genes ${\widehat{a}}_g\in [0,0.5)$ and for the remaining 72 genes ${\widehat{a}}_g\in [0.5,1)$. The spread of ${\widehat{a}}_g$ over a range of values far from 0 is an indication of the fact that for several genes the NB distribution (that is, a PT with $a=0$) does not provide the best fit to the data, in line with the observations made by Esnaola et al. (2013) about the distribution of cross-sectional RNA-seq data.

With the aim of identifying genes whose progression is different between mdx and wt mice, we tested the null hypothesis $H_{01}:{\beta }_{g1}={\beta }_{g3}=0$ of no difference in the average expression level of gene $g$ between wt and mdx at any time point. p-values were adjusted for multiple testing using the Benjamini–Hockberg method (Benjamini and Hochberg, 1995). Detailed results of the test are provided in File 2 of Supplementary material. In short, we found evidence of differential expression (FDR-corrected p-value $<0.05$) for 120 genes; for these genes, we furthermore tested whether there was sufficient evidence that the intercept or the slope are different between groups, $H_{02}:{\beta }_{g1}=0$ and $H_{03}:{\beta }_{g3}=0$ respectively. These tests allowed to further identify 38 genes with a significant intercept (difference between wt and mdx at week 6) and 38 genes where the progression (slope) is different between groups.

Download

Open in new tab

Download in PowerPoint

Figure 2: Plots comparing the trajectories of selected genes in mdx (blue solid line) and wt (red dashed line) mice

In Figure 2, we show the trajectories of four of the significant genes, Igfbp4, Jak1, Stat1 and Tgif1. Igfbp4 is one of the insulin-like growth factor binding proteins, it binds to both IGF-1 and IGF-2 affecting the half-life of IGFs in circulation and their interaction with cell receptors. IGF-1 has been linked to muscle regeneration and multiple articles show a beneficial effect of IGF-1 in DMD (Forcina et al., 2019). IGF binding proteins such as IGFBP3 and IGFBP5 have also been described to be reduced in DMD patients serum compared to healthy controls (Hathout et al., 2019). In this work, we observe a significant elevation of Igfbp4 in mdx mice that is particularly apparent at weeks 6 and 12, when intensive muscle degeneration and regeneration is ongoing. Jak1 and Stat1 are part of the JAK-STAT pathway; inspection of the trajectory plots in Figure 2 shows that the two genes have a similar pattern of expression, with an elevation in mdx mice at week 6 followed by a steep decrease at the later time points. It has been described how the JAK-STAT pathway contributes to muscle repair promoting regeneration in DMD (Moresi et al., 2019). The observed behaviour of Jak1 and Stat1 in mdx mice is particularly interesting for two reasons: the intensive muscle regeneration ongoing at 6 weeks of age in mdx mice, and the previous finding that Stat1 is activated in mdx mice diaphragm at 6 weeks of age compared to healthy mice (Dogra et al., 2008). Finally, for Tgif1 we can observe a downregulation in mdx mice at 6 and 12 weeks; later on, the expression of the gene in mdx mice increases towards the levels observed in wt mice. Tgif1 is part of the TGF$\beta$ pathway, a well-known pathway in the DMD pathophysiology with genetic modifiers mapping to this pathway (Flanigan et al., 2013) and a number of reports involving TGF$\beta$ in DMD (Ceco and McNally, 2013).

In Figure 3, we compare the p-values for $H_{01}:{\beta }_{g1}={\beta }_{g3}=0$ obtained from the PT mixed model with those that would have been obtained if we had used a NB mixed model, focusing on the 311 genes for which both models could be fitted. It can be observed that while for several genes the p-values are quite similar, there are also numerous instances of genes for which the choice between NB and PT matters, as it may entail a different conclusion on the significance of the null hypothesis. In particular, p-values obtained from the PT mixed model are smaller than those from the NB mixed model in 178 cases, and larger 133 times; after application of the Benjamini–Hochberg multiple testing correction, 44 genes are found as significant ($\alpha =0.05$) with both the NB and the PT mixed model; furthermore, 12 genes are identified as significant just by the NB mixed model, whereas 53 genes are found as significant only by the PT mixed model. Overall, these results suggest that while there appears to be a good level of agreement between the two models, identifying for each gene the most suitable distributional shape within the PT family of distributions, rather than arbitrarily fixing $a_g=0$ and choosing the NB distribution, can have important consequences on the identification of DE genes and on downstream analyses.

High-throughput RNA-sequencing technologies have become the standard method to measure gene expression, progressively replacing microarrays. RNA-seq generates count data that exhibit special features such as overdispersion, zero-inflation and heavy tails. The analysis of RNA-seq data arising from longitudinal experiments should be able to deal with such features and, at the same time, to account for the correlation arising from the availability of repeated measurements from the same individual. However, often such data are analysed with methods that either do not account for overdispersion, zero-inflation and heavy tails (Esnaola et al., 2013), or fail to model the dependence between repeated measurements (Cui et al., 2016).

Download

Open in new tab

Download in PowerPoint

Figure 3: Scatter plot comparing the p-values for the hypothesis $H_{01}:{\beta }_{g1}={\beta }_{g3}=0$ computed using the NB ($x$ axis) and PT ($y$ axis) mixed models

In this article, we have introduced ptmixed, a GLMM based on the PT distribution that can be employed to analyse of longitudinal sequencing data. On the one hand, the reliance of ptmixed on the PT distribution allows to flexibly model the typical distributional shapes of RNA-seq counts and, at the same time, to perform automatic model selection thanks to the PT power parameter, which allows to distinguish between several overdispersed distributions that can be obtained as special cases of the PT (Esnaola et al., 2013). On the other hand, ptmixed properly models the correlation between repeated measurements obtained from the same subject through the specification of subject-specific random effects. The applicability of ptmixed extends beyond RNA-seq data, embracing all scenarios where non-independent measurements of a discrete, overdispersed response variable are available; such settings can arise, for example, in ecology (Bolker et al., 2009) and in clinical trials (Albert, 1999).

Our proposal represents the first theoretical development and software implementation of a mixed model based on the PT distribution. An important contribution of our work is to have shown that an accurate evaluation of the integrals in (2.2) makes it possible to model longitudinal overdispersed counts using a GLMM that is more flexible than the commonly used Poisson and NB GLMMs. In our experience, we have observed that the rate of success for the estimation of the PT mixed model can be affected by the sample size and on whether the values of $a$ and ${\sigma }^2$ lie on the boundary of the parameter space; in particular, convergence appears to be easier with larger $n$, and to be more difficult to achieve when the true $a$ is very close to 1 or ${\sigma }^2$ approaches 0.

The simulations presented in Sections 3.1 and 3.2 show that ptmixed allows to accurately perform parameter estimation and hypothesis testing, and demonstrate the advantage connected to the use of the AGHQ instead of the LA. Furthermore, in Section 3.3 we presented a comparison between ptmixed and alternative modelling approaches through two genomewide simulations. Those simulations showed that analysing longitudinal sequencing data with methods that do not properly model the repeated measurements design (assuming independence or introducing subject-specific fixed effects) yield inaccurate parameter estimates (Tables 9 and 10) and fail to control the type I error in tests on between-subjects effects (Table 13). On the other hand, they showed that use of GLMM approaches such as ptmixed and ShrinkBayes make it possible to accurately estimate and test both within- and between-subject effects.

The results of Section 3.3 further allowed to compare ShrinkBayes and ptmixed. As concerns parameter estimation, the RMSEs of ShrinkBayes were found to be 1.5/2 times smaller than those of ptmixed for within-subject effects (Table 9), but 5 times larger than those of ptmixed for between-subjects effects (Table 10). Both methods successfully controlled the FDR at $\alpha =0.05$ in simulation C (Table 11), but ShrinkBayes did not control it when $n=10$ in simulation D. Interestingly, ptmixed was somewhat more powerful than ShrinkBayes in simulation C (Table 12), while the opposite was found in simulation D (Table 14). Finally, the estimation of ptmixed was faster than ShrinkBayes with small sample sizes, but slower with larger sample sizes. Overall, the comparison pointed out that none of the two methods uniformly outperformed the other, suggesting that both methods are valid modelling approaches with their own advantages and weaknesses. Aside from those performance measures, there are a few additional features of ptmixed that make it more flexible than ShrinkBayes: ptmixed is capable of modelling not only zero-inflation but also heavy-tails, it allows to select a different distributional shape for each gene rather than fitting the same distribution to all genes, and it allows to easily test not only simple univariate hypotheses involving a single parameter, but also more general multivariate hypotheses involving linear combinations of parameters and/or the joint significance of two or more parameters.

In Section 4, we illustrated an application of ptmixed to data from a longitudinal RNA-seq experiment involving wt and mdx mice. Through this application, we showed not only that ptmixed is capable to identify as DE many genes known to play a major role in the evolution of DMD, but also that using the more general PT GLMM instead of the simpler NB GLMM can yield substantially different conclusions on which genes are DE. This result emphasizes the importance of properly modelling of zero-inflation and heavy-tails not just from a theoretical point of view but also from an applied perspective, since doing so can substantially change the conclusions of real data analyses.

Supplementary materials for this article are available from http://www.statmod.org/smij/archive.html.

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

The authors gratefully acknowledge funding from the Duchenne Parent Project NL foundation (https://duchenne.nl) and from the Association Francaise contre les Myopathies (grant number 19118).