rmRNAseq: differential expression analysis for repeated-measures RNA-seq data

Abstract

Motivation: With the reduction in price of next-generation sequencing technologies, gene expression profiling using RNA-seq has increased the scope of sequencing experiments to include more complex designs, such as designs involving repeated measures. In such designs, RNA samples are extracted from each experimental unit at multiple time points. The read counts that result from RNA sequencing of the samples extracted from the same experimental unit tend to be temporally correlated. Although there are many methods for RNA-seq differential expression analysis, existing methods do not properly account for within-unit correlations that arise in repeated-measures designs.

Results: We address this shortcoming by using normalized log-transformed counts and associated precision weights in a general linear model pipeline with continuous autoregressive structure to account for the correlation among observations within each experimental unit. We then utilize parametric bootstrap to conduct differential expression inference. Simulation studies show the advantages of our method over alternatives that do not account for the correlation among observations within experimental units.

Availability and implementation: We provide an R package rmRNAseq implementing our proposed method (function TC_CAR1) at https://cran.r-project.org/web/packages/rmRNAseq/index.html. Reproducible R codes for data analysis and simulation are available at https://github.com/ntyet/rmRNAseq/tree/master/simulation.

Fig. 1.

Estimated correlations across all 11 911 genes for each pair of time points using the procedure in Section 2.3.2 applied to the log-transformed LPS RNA-seq data

Fig. 2.

Bar diagrams showing numbers of DE genes (FDR is nominally controlled at 0.05) with respect to Line and Time main effects when analyzing the LPS RNA-seq dataset using our method (rmRNAseq), voom-limma, edgeR, DESeq2, splineTimeR and ImpulseDE2

Fig. 3.

Boxplots of the incurred FDR when FDR is nominally controlled at 0.05 for all methods and all contrasts in two simulation scenarios. Each boxplot has 50 data points representing 50 simulated datasets

Fig. 4.

Boxplots of the PAUC when false positive rate is less than or equal to 0.05 for all methods and all contrasts in two simulation scenarios. Each boxplot has 50 data points representing 50 simulated datasets

Fig. 5.

Boxplots of the NTP when false positive rate is less than or equal to 0.05 for all methods and all contrasts in two simulation scenarios. Each boxplot has 50 data points representing 50 simulated datasets