MOTIVATION: As single-cell RNA-sequencing (scRNA-seq) becomes more widely used in transcriptomic research, complex experimental designs, such as paired or longitudinal studies, become increasingly feasible. Paired/longitudinal scRNA-seq enables the study of transcriptomic changes over time within specific cell types, yet guidance on analytical approaches and resources for study planning, such as power analysis, remains limited. Data simulation is a valuable tool for evaluating analysis method performance and informing study design decisions, including sample size selection. Currently, most scRNA-seq simulation methods simulate cells for a single sample, thus ignoring the between-sample and between-subject variability inherent to paired/longitudinal scRNA-seq data. RESULTS: Here, we introduce rescueSim (REpeated measures Single Cell RNA-seqUEncing data SIMulation), a novel method that simulates paired/longitudinal scRNA-seq data using a gamma-Poisson framework and incorporates additional variability between samples and subjects. We demonstrate our method's ability to reproduce important data properties and demonstrate its application in study planning. AVAILABILITY AND IMPLEMENTATION: rescueSim is implemented as an R package and is available at https://github.com/ewynn610/rescueSim.

BACKGROUND: Studies that utilize RNA Sequencing (RNA-Seq) in conjunction with designs that introduce dependence between observations (e.g. longitudinal sampling) require specialized analysis tools to accommodate this additional complexity. This R package contains a set of utilities to fit linear mixed effects models to transformed RNA-Seq counts that properly account for this dependence when performing statistical analyses. RESULTS: In a simulation study comparing lmerSeq and two existing methodologies that also work with transformed RNA-Seq counts, we found that lmerSeq was comprehensively better in terms of nominal error rate control and statistical power. CONCLUSIONS: Existing R packages for analyzing transformed RNA-Seq data with linear mixed models are limited in the variance structures they allow and/or the transformation methods they support. The lmerSeq package offers more flexibility in both of these areas and gave substantially better results in our simulations.

BACKGROUND: As the barriers to incorporating RNA sequencing (RNA-Seq) into biomedical studies continue to decrease, the complexity and size of RNA-Seq experiments are rapidly growing. Paired, longitudinal, and other correlated designs are becoming commonplace, and these studies offer immense potential for understanding how transcriptional changes within an individual over time differ depending on treatment or environmental conditions. While several methods have been proposed for dealing with repeated measures within RNA-Seq analyses, they are either restricted to handling only paired measurements, can only test for differences between two groups, and/or have issues with maintaining nominal false positive and false discovery rates. In this work, we propose a Bayesian hierarchical negative binomial generalized linear mixed model framework that can flexibly model RNA-Seq counts from studies with arbitrarily many repeated observations, can include covariates, and also maintains nominal false positive and false discovery rates in its posterior inference. RESULTS: In simulation studies, we showed that our proposed method (MCMSeq) best combines high statistical power (i.e. sensitivity or recall) with maintenance of nominal false positive and false discovery rates compared the other available strategies, especially at the smaller sample sizes investigated. This behavior was then replicated in an application to real RNA-Seq data where MCMSeq was able to find previously reported genes associated with tuberculosis infection in a cohort with longitudinal measurements. CONCLUSIONS: Failing to account for repeated measurements when analyzing RNA-Seq experiments can result in significantly inflated false positive and false discovery rates. Of the methods we investigated, whether they model RNA-Seq counts directly or worked on transformed values, the Bayesian hierarchical model implemented in the mcmseq R package (available at https://github.com/stop-pre16/mcmseq ) best combined sensitivity and nominal error rate control.