MCMSeq: Bayesian hierarchical modeling of clustered and repeated measures RNA sequencing experiments

Abstract

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.

NEBULA is a fast negative binomial mixed model for differential or co-expression analysis of large-scale multi-subject single-cell data

“…The glmmTMB package utilizes the Template Model Builder automatic differentiation engine, and the INLA package uses an efficient Bayesian framework based on integrated nested Laplace approximations (LAs). We did not include methods based on adaptive Gaussian quadrature (AGQ) ( glmer.nb with nAGQ>1) 24 , 25 or Bayesian methods using Markov chain Monte Carlo (MCMC) 26 because of their computational intensity. The Gaussian variational approximation method 14 can be promising but has not been implemented for the NBMM.…”

“…Although the time complexity is under the assumption of the independent random effects, the two-layer procedure often requires hundreds of iterations to converge, which becomes the major computational bottleneck. On the other hand, Bayesian methods such as MCMC 26 offer better accuracy, but are even slower than the above methods.…”

Compositional Data Analysis using Kernels in Mass Cytometry Data

“…The newer methods such as edgeR or voom do not provide theoretical guarantee of type-I error control either. These methods have been shown to have inflated type-I error rate when applied to other types of data (Hawinkel et al, 2019;Rocke et al, 2015;Vestal et al, 2020;Datta and Nettleton, 2014). In this paper, we propose a multivariate statistical framework, 'Compositional Data Analysis using Kernels' (CODAK), based on kernel distance covariance (KDC) (Hua and Ghosh, 2015) to quantify and test the association of predictors such as grouping or application of drugs with the composition profile of cell types.…”

A comparison of methods for multiple degree of freedom testing in repeated measures RNA-sequencing experiments

“…baySeq [ 13 , 14 ], PairedFB [ 15 ]), or can only perform single DF tests (e.g. MCMSeq [ 16 ], ShrinkBayes [ 17 ]).…”

“…The package rmRNAseq [ 18 ] utilizes the voom normalization method on log-transformed counts and then models the transformed data using a linear model with a continuous auto-regressive structure to account for the correlation in the data. Vestal et al [ 16 ] tested a similar method by using a variance stabilizing transformation (VST) on raw RNA-seq counts and then fitting linear mixed models (LMMs) to the transformed data. They found that this method performed similarly to their hierarchical Bayesian MCMSeq method in terms of T1E and FDR control, but many models failed to converge in small sample size situations.…”