glmmSeq

The aim of this package is to model gene expression with a general linear mixed model (GLMM) as described in the R4RA study [1]. The most widely used mainstream differential gene expression analysis tools (e.g Limma, DESeq2, edgeR) are all unable to fit mixed effects linear models. This package however fits negative binomial mixed effects models at individual gene level using the negative.binomial function from MASS and the glmer function in lme4 which enables random effect, as as well as mixed effects, to be modelled.

install.packages("glmmSeq")

devtools::install_github("myles-lewis/glmmSeq")

functions = list.files("./R", full.names = TRUE)
invisible(lapply(functions, source))

# Install CRAN packages
invisible(lapply(c("MASS", "car", "ggplot2", "ggpubr", "lme4","lmerTest",
                   "methods", "parallel", "plotly", "pbapply", "pbmcapply"),
                 function(p){
                   if(! p %in% rownames(installed.packages())) {
                     install.packages(p)
                   }
                   library(p, character.only=TRUE)
                 }))

# Install BioConductor packages
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
invisible(lapply(c("qvalue"), function(p){
  if(! p %in% rownames(installed.packages())) BiocManager::install(p)
  library(p, character.only=TRUE)
}))

Overview

To get started, first we load in the package:

library(glmmSeq)
set.seed(1234)

This vignette will demonstrate the power of this package using a minimal example from the PEAC data set. Here we focus on the synovial biopsy RNA-Seq data from this cohort of patients with early rheumatoid arthritis.

data(PEAC_minimal_load)

This data contains:

• metadata: which describes each sample. Including patient ID, sample time-point, and six-month EULAR response. Where EULAR is a rheumatoid arthritis response metric based on composite DAS28 scores.
• tpm: the transcript per million RNA-seq count data

These are outlined in the following subsections.

metadata$EULAR_binary  = NA
metadata$EULAR_binary[metadata$EULAR_6m %in%
                        c("Good", "Moderate" )] = "responder"
metadata$EULAR_binary[metadata$EULAR_6m %in% c("Non-response")] = "non_responder"

kable(head(metadata), row.names = F) %>% kable_styling()

kable(head(tpm)) %>% kable_styling() %>%
  scroll_box(width = "100%")

disp <- apply(tpm, 1, function(x){
  (var(x, na.rm=TRUE)-mean(x, na.rm=TRUE))/(mean(x, na.rm=TRUE)**2)
  })

head(disp)

##     MS4A1    ADAM12 IGHV7-4-1  IGHV3-49  IGHV3-23  ADAM12.1 
##  3.789428  2.391912 11.686420 10.863156  3.262557  2.391912

disp  <- setNames(edgeR::estimateDisp(tpm)$tagwise.dispersion, rownames(tpm))

head(disp)

dds <- DESeqDataSetFromTximport(txi = txi, colData = metadata, design = ~ 1)
dds <- DESeq(dds)
dispersions <- setNames(dispersions(dds), rownames(txi$counts))

sizeFactors <- colSums(tpm)  
sizeFactors <- sizeFactors / mean(sizeFactors)  # normalise to mean = 1

head(sizeFactors)

##      S9001      S9002      S9003      S9004      S9006      S9007 
## 0.20325747 0.09330435 0.09737100 3.13456671 0.24515015 4.19419297

sizeFactors <- calcNormFactors(counts, method="TMM")

sizeFactors <- estimateSizeFactorsForMatrix(counts)

Fitting Models

To fit a model for one gene over time we use a formula such as:

gene expression ~ fixed effects + random effects

In R the formula is defined by both the fixed-effects and random-effects part of the model, with the response on the left of a ~ operator and the terms, separated by + operators, on the right. Random-effects terms are distinguished by vertical bars (“|”) separating expressions for design matrices from grouping factors. For more information see the ?lme4::glmer.

In this case study we want to use time and response as fixed effects and the patients as random effects:

To fit this model for all genes we can use the glmmSeq function. Note that this analysis can take some time, with 4 cores:

• 1000 genes takes about 10 seconds
• 20000 genes takes about 4 mins

results <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                   countdata = tpm,
                   metadata = metadata,
                   dispersion = disp,
                   progress = TRUE)

## 
## n = 123 samples, 82 individuals
## Time difference of 3.412316 secs
## Errors in 4 gene(s): IL12A, FGF12, VIL1, IL26

or alternatively using two-factor classification with EULAR_binary:

results2 <- glmmSeq(~ Timepoint * EULAR_binary + (1 | PATID),
                    countdata = tpm,
                    metadata = metadata,
                    dispersion = disp)

## 
## n = 123 samples, 82 individuals
## Time difference of 3.078482 secs
## Errors in 4 gene(s): IL12A, FGF12, VIL1, IL26

names(attributes(results))

##  [1] "info"      "formula"   "stats"     "predict"   "reduced"   "countdata" "metadata" 
##  [8] "modelData" "optInfo"   "errors"    "vars"      "class"

kable(results@modelData) %>% kable_styling()

stats <- summary(results)

kable(stats[order(stats[, 'P_Timepoint:EULAR_6m']), ]) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

summary(results, gene = "MS4A1")

## Generalised linear mixed model
## Method: lme4
## Family: Negative Binomial
## Formula: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##  Dispersion         AIC      logLik     meanExp 
##    3.789428  898.682864 -441.341432    2.591740 
## 
## Fixed effects:
##                                Estimate Std. Error
## (Intercept)                     2.73637    0.33669
## Timepoint                       0.04908    0.08941
## EULAR_6mModerate                0.19330    0.51104
## EULAR_6mNon-response           -1.72376    0.67780
## Timepoint:EULAR_6mModerate     -0.21080    0.13478
## Timepoint:EULAR_6mNon-response  0.17462    0.16558
## 
## Analysis of Deviance Table (Type II Wald chisquare tests)
##                      Chisq Df Pr(>Chisq)
## Timepoint          0.01125  1    0.91554
## EULAR_6m           5.68295  2    0.05834
## Timepoint:EULAR_6m 5.43787  2    0.06595

predict = data.frame(results@predict)
kable(predict) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")

results <- glmmQvals(results)

## 
## Timepoint
## ---------
## Not Significant     Significant 
##              29              17 
## 
## EULAR_6m
## --------
## Not Significant     Significant 
##              41               5 
## 
## Timepoint:EULAR_6m
## ------------------
## Not Significant 
##              46

glmmTMB models

The glmmTMB package is an alternative to lme4 for fitting negative binomial GLMM models. Note, that this package optimises the dispersion parameter for each GLMM model. This has the advantage that a predefined list of dispersions for each gene is not required. But it also means that model fitting is significantly slower than a negative binomial GLMM fit by lme4::glmer with family function MASS::negative.binomial for which the dispersion parameter must be pre-specified. Each glmmTMB model takes about 0.8 seconds, so that even with 8 cores, a typical transcriptome of 20,000 genes would take about 30 minutes, although this will vary depending on the complexity of the design formula. In comparison, the standard glmer pipeline with known dispersions supplied takes about 2 minutes on 8 cores.

Although glmmTMB is slower than lme::glmer with known dispersion, glmmTMB is still faster than the equivalent lme4 function glmer.nb when the dispersion is unknown.

Setting method = "glmmTMB" when calling glmmSeq() switches from using lme4::glmer (the default) to the glmmTMB package. The dispersion argument is then no longer required. The family argument can be used to specify different family functions and glmmTMB provides nbinom1 and nbinom2 (glmmSeq uses the latter by default). Other family functions e.g. poisson can be used.

Gaussian mixed effects models

An alternative to fitting negative binomial GLMM is to use transformed data and fit gaussian LMM. The lmmSeq function will fit many LMM across genes (or other data) where random effects/mixed effects models are useful. It is almost 3x faster than glmmSeq. For example, gaussian normalised DNA methylation data can be analysed using this method. In the example below the RNA-Seq gene expression data is log transformed so that it is approximately Gaussian. Alternatively variance stabilising transformation (VST) can be applied to count data through DESeq2 or the voom transformation in limma voom.

logtpm <- log2(tpm + 1)
lmmres <- lmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                   maindata = logtpm,
                   metadata = metadata,
                   progress = TRUE)

## 
## n = 123 samples, 82 individuals
## Time difference of 0.6729643 secs

summary(lmmres, "MS4A1")

## Linear mixed model
## Formula: gene ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##         AIC      logLik     meanExp 
##  503.415984 -243.707992    2.590358 
## 
## Fixed effects:
##                                 Estimate Std. Error
## (Intercept)                     2.614919    0.30012
## Timepoint                       0.024789    0.06841
## EULAR_6mModerate                0.860136    0.45592
## EULAR_6mNon-response           -0.983157    0.58504
## Timepoint:EULAR_6mModerate     -0.256884    0.10233
## Timepoint:EULAR_6mNon-response -0.009612    0.12407
## 
## Analysis of Deviance Table (Type II Wald chisquare tests)
##                    Chisq Df Pr(>Chisq)
## Timepoint          2.321  1    0.12767
## EULAR_6m           6.098  2    0.04740
## Timepoint:EULAR_6m 7.134  2    0.02823

Hypothesis testing

By default, p-values for each model term are computed using Wald type 2 Chi-squared test as per car::Anova(). The underlying code for this has been optimised for speed.

For LMM via lmmSeq(), an alternative to the Wald type 2 Chi-squared test is an F-test with Satterthwaite denominator degrees of freedom, which has been implemented using the lmerTest package and is enabled using test.stat = "F". This is not available for GLMM.

However, if a reduced model formula is specified by setting reduced, then a likelihood ratio test is performed instead using anova. This will double computation time since two (G)LMM have to be fitted for each gene. For further information on inference testing on GLMM, we recommend Ben Bolker’s GLMM FAQ page.

glmmLRT <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                   reduced = ~ Timepoint + EULAR_6m + (1 | PATID),
                   countdata = tpm,
                   metadata = metadata,
                   dispersion = disp, verbose = FALSE)

summary(glmmLRT, "MS4A1")

## Generalised linear mixed model
## Method: lme4
## Family: Negative Binomial
## Formula: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##  Dispersion         AIC      logLik     meanExp 
##    3.789428  898.682864 -441.341432    2.591740 
## 
## Fixed effects:
##                                Estimate Std. Error
## (Intercept)                     2.73637    0.33669
## Timepoint                       0.04908    0.08941
## EULAR_6mModerate                0.19330    0.51104
## EULAR_6mNon-response           -1.72376    0.67780
## Timepoint:EULAR_6mModerate     -0.21080    0.13478
## Timepoint:EULAR_6mNon-response  0.17462    0.16558
## 
## Likelihood ratio test
## Reduced formula: count ~ Timepoint + EULAR_6m + (1 | PATID)
##   Chisq Df Pr(>Chisq)
##    4.85  2    0.08846

Extracting individual Genes

Similarly you can run the script for an individual gene (make sure you use drop = FALSE to maintain countdata as a matrix).

MS4A1glmm <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                     countdata = tpm["MS4A1", , drop = FALSE],
                     metadata = metadata,
                     dispersion = disp,
                     verbose = FALSE)

The (g)lmer or glmmTMB fitted object for a single gene can be obtained using the function glmmRefit():

fit <- glmmRefit(results, gene = "MS4A1")

## boundary (singular) fit: see help('isSingular')

fit

## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) [
## glmerMod]
##  Family: Negative Binomial(0.2639)  ( log )
## Formula: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##    Data: data
##  Offset: offset
##       AIC       BIC    logLik  deviance  df.resid 
##  898.6829  921.1803 -441.3414  882.6829       115 
## Random effects:
##  Groups Name        Std.Dev.
##  PATID  (Intercept) 2e-07   
## Number of obs: 123, groups:  PATID, 82
## Fixed Effects:
##                    (Intercept)                       Timepoint  
##                        2.73637                         0.04908  
##               EULAR_6mModerate            EULAR_6mNon-response  
##                        0.19330                        -1.72376  
##     Timepoint:EULAR_6mModerate  Timepoint:EULAR_6mNon-response  
##                       -0.21080                         0.17462  
## optimizer (bobyqa) convergence code: 0 (OK) ; 0 optimizer warnings; 1 lme4 warnings

These (g)lmer objects can be passed to other packages notably emmeans for visualising estimated marginal means. Note these are based on the model, not directly on the data.

library(emmeans)

## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

emmeans(fit, ~ Timepoint | EULAR_6m)

## EULAR_6m = Good:
##  Timepoint emmean    SE  df asymp.LCL asymp.UCL
##          0   2.74 0.337 Inf      2.08      3.40
##          6   3.03 0.418 Inf      2.21      3.85
## 
## EULAR_6m = Moderate:
##  Timepoint emmean    SE  df asymp.LCL asymp.UCL
##          0   2.93 0.384 Inf      2.18      3.68
##          6   1.96 0.467 Inf      1.04      2.88
## 
## EULAR_6m = Non-response:
##  Timepoint emmean    SE  df asymp.LCL asymp.UCL
##          0   1.01 0.588 Inf     -0.14      2.17
##          6   2.35 0.594 Inf      1.19      3.52
## 
## Results are given on the log (not the response) scale. 
## Confidence level used: 0.95

emmip(fit, ~ Timepoint | EULAR_6m)

fit2 <- glmmRefit(results, gene = "MS4A1",
                  formula = count ~ Timepoint + EULAR_6m + (1 | PATID))

anova(fit, fit2)

## Data: data
## Models:
## fit2: count ~ Timepoint + EULAR_6m + (1 | PATID)
## fit: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##      npar    AIC    BIC  logLik deviance  Chisq Df Pr(>Chisq)  
## fit2    6 899.53 916.41 -443.77   887.53                       
## fit     8 898.68 921.18 -441.34   882.68 4.8505  2    0.08846 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Model Plots

For variables such as time, which are matched according to an ID (the random effect), we can examine the fitted model using plots which show estimated means and confidence intervals based on coefficients for the fitted regression model, overlaid upon the underlying data. In this case the samples are matched longitudinally over time.

Plots can be viewed using either ggplot or base graphics. We can start looking at the gene with the most significant interaction IGHV3-23:

plotColours <- c("skyblue", "goldenrod1", "mediumseagreen")
modColours <- c("dodgerblue3", "goldenrod3", "seagreen4")
shapes <- c(17, 19, 18)

ggmodelPlot(results,
            geneName = "IGHV3-23",
            x1var = "Timepoint",
            x2var="EULAR_6m",
            xlab="Time",
            colours = plotColours,
            shapes = shapes,
            lineColours = plotColours, 
            modelColours = modColours,
            modelSize = 10)

Alternatively plots can be generated using base graphics, here with or without the model fit overlaid. By default p-value labels are taken from the column names of the pvals object in the stats slot of the S4 result object. These can be relabelled using the plab argument.

oldpar <- par(mfrow=c(1, 2))

modelPlot(results2,
          geneName = "FGF14",
          x1var = "Timepoint",
          x2var="EULAR_binary",
          fontSize=0.65,
          colours=c("coral", "mediumseagreen"),
          modelColours = c("coral", "mediumseagreen"),
          modelLineColours = "black",
          modelSize = 2)

modelPlot(results,
          geneName = "EMILIN3",
          x1var = "Timepoint",
          x2var = "EULAR_6m",
          colours = plotColours,
          plab = c("time", "response", "time:response"),
          addModel = FALSE)

par(oldpar)

To plot the model fits alone set addPoints = FALSE.

library(ggpubr)

p1 <- ggmodelPlot(results,
                  "ADAM12",
                  x1var="Timepoint",
                  x2var="EULAR_6m",
                  xlab="Time",
                  addPoints = FALSE,
                  colours = plotColours)

p2 <- ggmodelPlot(results,
                  "EMILIN3",
                  x1var="Timepoint",
                  x2var="EULAR_6m",
                  xlab="Time",
                  fontSize=8,
                  x2Offset=1,
                  addPoints = FALSE,
                  colours = plotColours)

ggarrange(p1, p2, ncol=2, common.legend = T, legend="bottom")

r4ra_glmm <- glmmSeq(~ time * drug + (1 | Patient_ID), 
                       countdata = tpmdata, metadata,
                       dispersion = dispersions, cores = 8, removeSingles = T)
r4ra_glmm <- glmmQvals(r4ra_glmm)
labels = c(..)  # Genes to label
fcPlot(r4ra_glmm, x1var = "time", x2var = "drug", graphics = "plotly",
       pCutoff = 0.05, useAdjusted = TRUE,
       labels = labels,
       colours = c('grey', 'green3', 'gold3', 'blue'))

labels = c('MS4A1', 'FGF14', 'IL2RG', 'IGHV3-23', 'ADAM12', 'IL36G', 
           'BLK', 'SAA1', 'CILP', 'EMILIN3', 'EMILIN2', 'IGHJ6', 
           'CXCL9', 'CXCL13')
maPlots <- maPlot(results,
                  x1var="Timepoint",
                  x2var="EULAR_6m",
                  x2Values=c("Good", "Non-response"),
                  colours=c('grey', 'midnightblue',
                             'mediumseagreen', 'goldenrod'),
                  labels=labels,
                  graphics="ggplot")

maPlots$combined

Troubleshooting

Mixed effects models are tricky to fit and lme4::glmer sometimes returns errors. In fact, a significant amount of code in glmmSeq is devoted to catching and handling errors to allow parallelisation to continue. Errors in genes are stored in the errors slot.

results@errors[1]   # first gene error

##                                                             IL12A 
## "Error in eval(expr, envir) : PIRLS loop resulted in NaN value\n"

Sometimes glmmSeq returns errors in all genes. This usually means a problem with not enough samples in each timepoint or a mistake in the formula. Since version 0.5.5 glmmSeq now returns a vector of the error messages for all genes, which can be useful for debugging. In the example below, only timepoint 0 is specified.

results3 <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                    countdata = tpm[, metadata$Timepoint == 0],
                    metadata = metadata[metadata$Timepoint == 0, ],
                    dispersion = disp)

## 
## n = 72 samples, 72 individuals

## All genes returned an error. Check call. Check sufficient data in each group

##                                                                              MS4A1 
## "fixed-effect model matrix is rank deficient so dropping 3 columns / coefficients"

If there are errors which are not caught by the error checking core mechanism, this can lead to problems with grouping results after the core models have been fit. Setting returnList=TRUE when calling glmmSeq returns the list output direct from mclapply (or parLapply on windows). This can be helpful for debugging unforeseeen problems in the core loop.

Citing glmmSeq

glmmSeq was developed by the bioinformatics team at the Experimental Medicine & Rheumatology department and Centre for Translational Bioinformatics at Queen Mary University London.

If you use this package please cite as:

citation("glmmSeq")

## To cite package 'glmmSeq' in publications use:
## 
##   Lewis M, Goldmann K, Sciacca E (????). _glmmSeq: General Linear Mixed Models
##   for Gene-Level Differential Expression_. R package version 0.5.6,
##   https://github.com/myles-lewis/glmmSeq,
##   <https://myles-lewis.github.io/glmmSeq/>.
## 
## A BibTeX entry for LaTeX users is
## 
##   @Manual{,
##     title = {glmmSeq: General Linear Mixed Models for Gene-Level Differential
## Expression},
##     author = {Myles Lewis and Katriona Goldmann and Elisabetta Sciacca},
##     note = {R package version 0.5.6, 
## https://github.com/myles-lewis/glmmSeq},
##     url = {https://myles-lewis.github.io/glmmSeq/},
##   }

References

1. Felice Rivellese, Anna Surace, Katriona Goldmann, Elisabetta Sciacca, Cankut Cubuk, Giovanni Giorli, … Michael Barnes, Myles J. Lewis, Costantino Pitzalis, R4RA collaborative group. Rituximab versus tocilizumab in rheumatoid arthritis: synovial biopsy-based biomarker analysis of the phase 4 R4RA randomized trial. Nature medicine 2022; 28(6): 1256-68. doi:10.1038/s41591-022-01789-0

Statistical software used in this package:

2. lme4: Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1-48. doi: 10.18637/jss.v067.i01.

3. car: John Fox and Sanford Weisberg (2019). An {R} Companion to Applied Regression, Third Edition. Thousand Oaks CA: Sage. URL: https://socialsciences.mcmaster.ca/jfox/Books/Companion/

4. MASS: Venables, W. N. & Ripley, B. D. (2002) Modern Applied Statistics with S. Fourth Edition. Springer, New York. ISBN 0-387-95457-0

5. glmmTMB: Mollie Brooks, Ben Bolker, Kasper Kristensen, Martin Maechler, Arni Magnusson (2022). glmmTMB: Generalized Linear Mixed Models using Template Model Builder

6. qvalue: John D. Storey, Andrew J. Bass, Alan Dabney and David Robinson (2020). qvalue: Q-value estimation for false discovery rate control. R package version 2.22.0. https://github.com/StoreyLab/qvalue

7. lmerTest: Alexandra Kuznetsova, Per Brockhoff, Rune Christensen, Sofie Jensen. lmerTest: Tests in Linear Mixed Effects Models

8. emmeans: Russell V. Lenth, Paul Buerkner, Maxime Herve, Jonathon Love, Fernando Miguez, Hannes Riebl, Henrik Singmann. emmeans: Estimated Marginal Means, aka Least-Squares Means