glmmSeq

The aim of this package is to model gene expression with a general linear mixed model (glmm). 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("KatrionaGoldmann/glmmSeq")

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

# Install CRAN packages
invisible(lapply(c("MASS", "car", "ggplot2", "ggpubr", "lme4", "methods",
                   "parallel", "plotly", "stats", "gghalves"),
                 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 will focus on the synovial data from this cohort.

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 responder", "Moderate responder" )] = "responder"
metadata$EULAR_binary[metadata$EULAR_6m %in% c("Non responder")] = "non_responder"
metadata = metadata[! is.na(metadata$EULAR_binary), ]

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

tpm = tpm[, metadata$SAMID]
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.809307  2.386833 11.790220 10.959627  3.302796  2.386833

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

head(disp)

##     MS4A1    ADAM12 IGHV7-4-1  IGHV3-49  IGHV3-23  ADAM12.1 
## 0.7951463 2.1345465 3.2851952 2.3843452 1.9764054 2.1345465

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

head(sizeFactors)

## SAM20389176 SAM20389177 SAM20389178 SAM20389179 SAM20389185 SAM20389187 
##  4.22168329  0.24675696  0.09137974  3.15511184  0.09391590  0.20458970

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 2 cores:

• 100 genes takes a couple of seconds
• 1000 genes takes about 20 seconds
• All ~20000 genes takes about 9 mins

results <- glmmSeq(~ Timepoint * EULAR_6m + (1 | PATID),
                  id = "PATID",
                  countdata = tpm,
                  metadata = metadata,
                  dispersion = disp,
                  removeDuplicatedMeasures = FALSE,
                  removeSingles=FALSE,
                  progress=TRUE,
                  cores = 1)

## 
## n = 124 samples, 82 individuals
## Time difference of 11.34985 secs

or alternatively using two-factor classification with EULAR_binary:

results2 <- glmmSeq(~ Timepoint * EULAR_binary + (1 | PATID),
                  id = "PATID",
                  countdata = tpm,
                  metadata = metadata,
                  dispersion = disp,
                  removeDuplicatedMeasures = FALSE,
                  removeSingles=FALSE,
                  cores = 1)

## 
## n = 124 samples, 82 individuals
## Time difference of 9.897644 secs

names(attributes(results))

##  [1] "formula"        "stats"          "predict"        "reducedFormula" "countdata"     
##  [6] "metadata"       "modelData"      "optInfo"        "errors"         "variables"     
## [11] "class"

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

stats = data.frame(results@stats)

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

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

Qvalues

The qvalues from each of the pvalue columns can be calculated using glmmQvals. This will output a significance table based on the cut-off (default p=0.05) and add qvalue columns to the @stats slot:

results <- glmmQvals(results, pi0=1)

## 
## q_Timepoint
## -----------
## Not Significant     Significant 
##              26              24 
## 
## q_EULAR_6m
## ----------
## Not Significant     Significant 
##              38              12 
## 
## q_Timepoint:EULAR_6m
## --------------------
## Not Significant     Significant 
##              35              15

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

MS4A1fit <- glmmGene(~ Timepoint * EULAR_6m + (1 | PATID),
                     gene = "MS4A1",
                     id = "PATID",
                     countdata = tpm,
                     metadata = metadata,
                     dispersion = disp['MS4A1'])

MS4A1fit

## Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
##  Family: Negative Binomial(1.2576)  ( log )
## Formula: count ~ Timepoint + EULAR_6m + (1 | PATID) + Timepoint:EULAR_6m
##    Data: data
##  Offset: offset
##       AIC       BIC    logLik  deviance  df.resid 
##  856.8937  879.4559 -420.4468  840.8937       116 
## Random effects:
##  Groups Name        Std.Dev.
##  PATID  (Intercept) 1.087   
## Number of obs: 124, groups:  PATID, 82
## Fixed Effects:
##                          (Intercept)                             Timepoint  
##                              1.81058                               0.03217  
##           EULAR_6mModerate responder                 EULAR_6mNon responder  
##                              0.66051                              -1.06887  
## Timepoint:EULAR_6mModerate responder       Timepoint:EULAR_6mNon responder  
##                             -0.22577                               0.00898

Paired Plots

For variables which are paired according to an ID (the random effect), we can view the model using paired plots. In this case the samples can be paired 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 <- rep(c("dodgerblue3", "goldenrod3", "seagreen4"), each=2)

pairedPlot(glmmResult=results,
           geneName = "IGHV3-23",
           x1Label = "Timepoint",
           x2Label="EULAR_6m",
           xTitle="Time",
           IDColumn = "PATID",
           graphics = "ggplot",
           colours = plotColours,
           modelColour = modColours,
           modelLineColour = modColours,
           fontSize=10,
           x2Offset = 8,
           logTransform=TRUE,
           addViolin = TRUE,
           pairedOnly = FALSE)

Or using base graphics, with or without the model fit overlaid:

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

p1 = pairedPlot(glmmResult=results2,
                geneName = "FGF14",
                x1Label = "Timepoint",
                x2Label="EULAR_binary",
                IDColumn = "PATID",
                graphics="base",
                fontSize=0.65,
                colours=c("coral", "mediumseagreen"),
                modelSize = 1)

p2 = pairedPlot(glmmResult=results,
                geneName = "EMILIN3",
                x1Label = "Timepoint",
                x2Label="EULAR_6m",
                IDColumn = "PATID",
                addModel=TRUE,
                graphics="base",
                fontSize=0.65,
                colours=plotColours)

par(oldpar)

Model plots

Alternatively to plot the model fits alone you can use the modelPlot function:

library(ggpubr)

p1 <- modelPlot(results,
                "ADAM12",
                x1Label="Timepoint",
                x2Label="EULAR_6m",
                xTitle="Time",
                fontSize=8,
                x2Offset=6,
                overlap=FALSE,
                graphics="ggplot",
                colours = plotColours)

p2 <- modelPlot(results,
                "ADAM12",
                x1Label="Timepoint",
                x2Label="EULAR_6m",
                xTitle="Time",
                fontSize=8,
                x2Offset=1,
                addErrorbars = FALSE,
                overlap=TRUE,
                graphics="ggplot",
                colours = plotColours)

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

Fold Change Plots

The comparative fold change (for x1Label variables) between conditions (x2Label and x2Values variables) can be plotted using fcPlots for all genes to highlight significance.

(By setting graphics=“plotly” this can be viewed interactively)

# Genes to label:
labels = c('MS4A1', 'FGF14', 'IL2RG', 'IGHV3-23', 'ADAM12', 'FGFRL1', 'IL36G', 
           'BLK', 'SAA1', 'CILP', 'EMILIN3', 'EMILIN2', 'IGHJ6', 
           'CXCL9', 'CXCL13')

fcPlot(glmmResult=results,
       x1Label="Timepoint",
       x2Label="EULAR_6m",
       x2Values=c("Good responder", "Non responder"),
       pCutoff=0.1,
       labels=labels,
       useAdjusted = FALSE,
       plotCutoff = 1)

Genes on the x-y plane, such as IGHJ6, will have associations in the same direction whereas genes on the x=-y axis have association in opposite directions, such as ADAM12. This allows us to pick out genes of potential interest which we can have another closer look at:

p1 <- pairedPlot(glmmResult=results,
                 geneName = "ADAM12",
                 x1Label = "Timepoint",
                 x2Label="EULAR_6m",
                 IDColumn = "PATID",
                 graphics="ggplot",
                 colours = "grey60",
                 modelColour = rep(plotColours, each=2),
                 modelLineColour =rep(plotColours, each=2),
                 addViolins=FALSE,
                 fontSize=8,
                 logTransform=T) + theme(plot.subtitle=element_text(size=9))

p2 <- pairedPlot(glmmResult=results,
                 geneName = "IGHJ6",
                 x1Label = "Timepoint",
                 x2Label="EULAR_6m",
                 IDColumn = "PATID",
                 graphics="ggplot",
                 addViolins = FALSE,
                 colours = c("blue"),
                 fontSize=8,
                 modelSize=0.1,
                 logTransform=T) + theme(plot.subtitle=element_text(size=9))

ggarrange(p1, p2, ncol=2)

Or flipping the x1 and x2 labels we can look at the fold change between response at different time points. This might be interesting to see if there are differences at certain time points which are not present at others.

fcPlot(glmmResult=results,
       x2Label="Timepoint",
       x1Label="EULAR_6m",
       x1Values=c("Good responder", "Non responder"),
       labels=labels,
       pCutoff=0.1,
       useAdjusted = F,
       plotCutoff = 1,
       graphics="ggplot")

MA plots

An MA plot is an application of a Bland–Altman plot. The plot visualizes the differences between measurements taken in two samples, by transforming the data onto M (log ratio) and A (mean average) scales, then plotting these values.

maPlots <- maPlot(results,
                  x1Label="Timepoint",
                  x2Label="EULAR_6m",
                  x2Values=c("Good responder", "Non responder"),
                  colours=c('grey', 'midnightblue',
                             'mediumseagreen', 'goldenrod'),
                  labels=labels,
                  graphics="ggplot")

maPlots$combined

maPlots <- maPlot(results,
                  x2Label="Timepoint",
                  x1Label="EULAR_6m",
                  x1Values=c("Good responder", "Non responder"),
                  colours=c('grey', 'midnightblue',
                             'mediumseagreen', 'goldenrod'),
                  labels=labels,
                  graphics="ggplot")

maPlots$combined

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:
## 
##   Myles Lewis, Katriona Goldmann and Elisabetta Sciacca (NA). glmmSeq: General Linear
##   Mixed Models for Gene-Level Differential Expression. R package version 0.1.0.
##   https://github.com/KatrionaGoldmann/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.1.0},
##     url = {https://github.com/KatrionaGoldmann/glmmSeq},
##   }

References

Statistical software used in this package:

1. 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.

2. 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/

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

4. 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