glmmSeq

Lifecycle: Maturing License: MIT CRAN status Downloads 2025-02-16 GitHub issues Travis

glmmSeq

glmmSeq logo

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.

Installing from CRAN

install.packages("glmmSeq")

Installing from Github

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

Installing Locally

Or you can source the functions individually:

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

But you will need to load in the additional libraries:

# 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

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()
PATID Timepoint EULAR_6m EULAR_binary
PAT300 6 Good responder
PAT209 6 Good responder
PAT219 6 Moderate responder
PAT211 6 Good responder
PAT216 6 Good responder
PAT212 6 Good responder

Count data

kable(head(tpm)) %>% kable_styling() %>%
  scroll_box(width = "100%")
S9001 S9002 S9003 S9004 S9006 S9007 S9008 S9009 S9010 S9011 S9012 S9013 S9014 S9016 S9017 S9018 S9019 S9020 S9021 S9023 S9025 S9029 S9034 S9035 S9036 S9038 S9039 S9040 S9042 S9043 S9044 S9045 S9047 S9048 S9049 S9052 S9053 S9054 S9056 S9059 S9060 S9063 S9065 S9066 S9067 S9068 S9069 S9070 S9072 S9073 S9074 S9075 S9076 S9077 S9078 S9079 S9080 S9081 S9083 S9084 S9085 S9086 S9087 S9088 S9089 S9090 S9091 S9092 S9093 S9094 S9095 S9096 S9097 S9098 S9099 S9100 S9101 S9102 S9103 S9104 S9105 S9106 S9107 S9108 S9109 S9110 S9111 S9112 S9113 S9114 S9115 S9116 S9117 S9119 S9120 S9121 S9122 S9123 S9124 S9125 S9127 S9128 S9129 S9130 S9131 S9132 S9133 S9134 S9135 S9136 S9137 S9138 S9139 S9140 S9141 S9142 S9143 S9144 S9145 S9146 S9147 S9148 S9149
MS4A1 1.622818 2.024451 0.6995153 15.742876 2.123731 57.280078 2.011160 0.580000 1.006783 0.391420 3.886618 2.746742 31.647637 20.16757 0.925393 3.4271907 3.130648 1.422480 63.14643 0.295537 12.7877279 4.817494 15.874491 0.569317 4.429667 1.060991 2.1850700 5.321541 3.542383 1.676829 5.471954 29.942195 4.5519540 7.252108 1.9540320 12.947362 73.430132 2.364819 5.2372870 12.281007 17.637846 72.557856 2.379427 0.444729 89.042394 2.4190970 226.612930 8.01254 5.328680 0.7414830 4.669566 5.218332 0.0676825 26.0579080 77.100135 23.2918620 50.306168 5.24953 1.342645 0.511810 0.2428008 7.756002 1.020052 0.602854 0.2444728 2.047680 2.110914 0.683498 24.664095 1.9789030 0.1702800 3.910493 6.5513870 0.8342158 31.7618280 1.136854 17.98745 3.369558 8.577294 1.453945 41.031738 42.991580 27.58865 14.876902 21.6405734 42.049748 75.6557000 2.501497 3.560322 32.1091249 7.008288 0.746264 10.0969886 92.3355460 0.5157370 0.8011538 5.075003 1.979931 0.2068106 2.717007 0.984857 0.7771422 0.000000 8.2769093 1.261793 2.361150 2.372850 0.6968002 43.9459380 5.0257346 0.281387 0.4190982 1.716238 16.25508 1.9852750 26.591585 27.994183 1.325606 1.323197 12.47092 2.3049434 1.4901770 0.229867
ADAM12 2.581805 1.110548 0.9573992 5.466742 2.645153 0.972897 0.658545 1.159264 0.462860 3.074193 10.947068 9.678026 1.205472 11.63376 1.772810 0.3506414 3.528288 17.416140 6.36844 4.064450 0.2666618 2.113812 0.473972 1.222569 30.012798 0.181003 2.4554526 4.136035 13.961315 4.195791 3.629029 8.362213 0.6999296 7.894137 0.3268513 2.847538 7.015717 0.717947 0.4396555 1.141057 1.069326 1.895644 0.486478 14.504778 1.897848 0.9943284 4.140489 41.30644 8.812227 0.3317777 8.699836 8.323936 0.0292791 0.9615935 2.512841 0.8960788 3.766372 7.13890 1.917924 8.893462 0.0966324 1.550900 4.760995 1.769894 4.6901815 1.200418 1.392691 0.439077 0.658254 0.3610814 0.2501964 0.610416 0.6842673 2.8749874 0.6717734 16.374299 1.36295 8.725470 0.258386 0.214470 3.573125 2.086027 13.44347 4.007794 0.0338254 1.595552 0.9645215 1.913151 5.780085 0.7098401 5.242844 1.115879 0.6499524 0.5397263 0.2299409 0.7178614 0.535170 0.282374 0.2003997 2.202211 0.310417 2.0991998 0.215329 0.1962535 25.010318 2.970828 1.461817 2.6506093 0.6501276 0.9633052 8.715446 32.4636380 4.281851 0.83041 0.9840191 2.363019 0.314493 1.179129 0.656344 1.57062 0.4597666 0.4707508 4.727155
IGHV7-4-1 0.992885 0.000000 2.2060400 26.381700 0.000000 2158.300000 0.644889 0.000000 0.000000 1.500990 4.945980 38.268000 55.570100 107.91600 0.000000 0.0000000 0.000000 1.080940 27.51650 0.000000 21.4136000 6.121120 15.190300 0.000000 0.597722 0.000000 0.0844527 1228.620000 51.561100 1.149330 1.434460 3728.700000 6.3708500 1.553570 0.0000000 165.928000 30.714000 0.000000 5.5235500 0.770965 1603.290000 53.605000 16.063900 6.209400 1021.480000 0.2965720 85.466200 9.18953 41.191800 5.0514200 47.140500 42.541900 0.4660800 835.2970000 1929.760000 57.2972000 39.990700 0.00000 0.000000 0.411695 0.0000000 26.176100 1.033220 2.273400 0.0000000 11.443300 0.272802 0.000000 81.135100 0.7083360 0.0000000 1.298760 0.5145140 0.1260250 13.9310000 2.345910 77.63540 2.181890 63.096300 0.000000 6.473680 20.211900 347.99200 39.273300 2.3200600 584.705000 2147.2900000 242.387000 14.263100 52.4861000 0.399860 0.000000 0.7025670 11.2577000 0.6126850 4.5281300 155.926000 320.086000 1.3791200 1.587010 0.167711 0.0000000 0.000000 5.6540200 0.000000 0.000000 0.560530 1.2747600 50.6122000 5.3320100 0.000000 0.0000000 0.757454 132.98900 3.4184700 18.130100 3.085550 1.183940 0.277095 135.99100 0.7799880 0.6956290 0.000000
IGHV3-49 0.000000 0.196831 8.7892000 345.493000 2.093480 858.980000 3.201300 0.000000 0.000000 0.000000 114.334000 17.150200 136.605000 1080.73000 0.196674 0.0000000 0.000000 0.931275 140.22000 0.000000 202.4950000 60.114600 145.421000 0.145099 32.808500 1.207290 3.5942600 467.534000 55.301200 0.869369 14.462700 928.660000 64.4575000 0.430536 0.0000000 80.024800 53.191800 0.000000 35.4998000 0.303541 1002.870000 127.201000 634.856000 22.119900 975.855000 0.0000000 623.330000 129.50800 180.819000 29.8933000 183.605000 170.369000 0.2201470 403.4460000 384.597000 24.9283000 426.619000 1.08364 0.591472 4768.210000 0.0000000 93.669200 9.291720 0.607172 0.0000000 9.736200 1.132510 0.000000 175.472000 1.0305600 0.0885169 0.972741 12.5531000 0.8804760 163.1520000 1.764370 37.22090 15.345300 401.861000 0.575801 11.926900 45.198600 317.65400 523.093000 171.4520000 174.163000 5414.9000000 5.045300 386.829000 1190.5200000 4.854510 0.169181 15.7114000 415.2760000 0.1668720 0.1950730 34.623200 26.949400 2.3880600 2.021500 0.000000 0.0000000 0.000000 1.5027300 0.000000 1.783170 49.278100 1.1861700 99.7745000 10.0077000 0.000000 0.1908660 3.835610 22.31070 2.4220200 22.133600 440.657000 0.000000 0.000000 22.10260 5.3278400 0.5799090 0.210089
IGHV3-23 0.715133 3.999940 87.6508000 2349.850000 5.962460 5180.460000 17.278900 0.000000 0.000000 1.437710 985.721000 123.982000 1374.860000 3568.60000 9.857240 0.5625450 2.900430 2.172310 5859.05000 0.000000 502.8400000 366.704000 1793.510000 2.294430 97.123900 11.607200 8.1870600 1080.330000 411.058000 9.718060 114.608000 2306.720000 317.0130000 7.954070 2.7908400 865.088000 3563.990000 0.123914 23.4808000 4.896920 1681.420000 1929.860000 2241.850000 26.191600 5228.450000 4.1899600 2807.770000 1213.59000 387.839000 203.5700000 374.506000 374.388000 0.3392100 1359.9900000 1535.600000 118.7600000 2708.140000 4.59234 0.607162 93.584600 0.3873680 1273.530000 1.666060 5.899820 1.8143100 18.921500 17.654300 0.567110 1593.270000 0.7579240 0.2794190 6.994180 216.3070000 4.8487200 1286.5300000 12.741600 167.65000 115.720000 2142.700000 3.020540 79.259200 6471.240000 2508.41000 4077.530000 802.0710000 1884.330000 5188.7400000 296.555000 1696.260000 3027.5700000 215.969000 0.249956 79.7680000 338.8720000 2.2993700 5.4541000 214.318000 152.085000 9.7046700 11.809000 0.300002 0.0000000 0.000000 15.5774000 0.217214 4.113890 229.482000 14.7361000 1145.9400000 35.9383000 0.613331 5.6297400 4.861760 188.38500 94.7553000 297.408000 3505.310000 19.762300 0.860007 208.00400 54.3620000 52.5728000 0.785236
ADAM12.1 2.581805 1.110548 0.9573992 5.466742 2.645153 0.972897 0.658545 1.159264 0.462860 3.074193 10.947068 9.678026 1.205472 11.63376 1.772810 0.3506414 3.528288 17.416140 6.36844 4.064450 0.2666618 2.113812 0.473972 1.222569 30.012798 0.181003 2.4554526 4.136035 13.961315 4.195791 3.629029 8.362213 0.6999296 7.894137 0.3268513 2.847538 7.015717 0.717947 0.4396555 1.141057 1.069326 1.895644 0.486478 14.504778 1.897848 0.9943284 4.140489 41.30644 8.812227 0.3317777 8.699836 8.323936 0.0292791 0.9615935 2.512841 0.8960788 3.766372 7.13890 1.917924 8.893462 0.0966324 1.550900 4.760995 1.769894 4.6901815 1.200418 1.392691 0.439077 0.658254 0.3610814 0.2501964 0.610416 0.6842673 2.8749874 0.6717734 16.374299 1.36295 8.725470 0.258386 0.214470 3.573125 2.086027 13.44347 4.007794 0.0338254 1.595552 0.9645215 1.913151 5.780085 0.7098401 5.242844 1.115879 0.6499524 0.5397263 0.2299409 0.7178614 0.535170 0.282374 0.2003997 2.202211 0.310417 2.0991998 0.215329 0.1962535 25.010318 2.970828 1.461817 2.6506093 0.6501276 0.9633052 8.715446 32.4636380 4.281851 0.83041 0.9840191 2.363019 0.314493 1.179129 0.656344 1.57062 0.4597666 0.4707508 4.727155

Dispersion

Using negative binomial models requires gene dispersion estimates to be made. This can be achieved in a number of ways. A common way to calculate this for gene i is to use the equation:

Dispersioni = (variancei - meani)/meani2

This can be calculated using:

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

Alternatively, we recommend using edgeR to estimate of the common, trended and tagwise dispersions across all tags:

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

head(disp)

or with DESeq2 using the raw counts:

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

Size Factors

There is also an option to include size factors for each gene. Again this can be estimated using:

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

Or using edgeR these can be calculated from the raw read counts:

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

Similarly, with DESeq2:

sizeFactors <- estimateSizeFactorsForMatrix(counts)

Note the sizeFactors vector needs to be centred around 1, since it used directly as an offset of form log(sizeFactors) in the GLMM model.

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:

gene expression ~ time + response + (1 | patient)

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

Outputs

This creates a GlmmSeq object which contains the following slots:

names(attributes(results))
##  [1] "info"      "formula"   "stats"     "predict"   "reduced"   "countdata" "metadata" 
##  [8] "modelData" "optInfo"   "errors"    "vars"      "class"

The variables used by the model are in the @modeldata:

kable(results@modelData) %>% kable_styling()
Timepoint EULAR_6m
0 Good
6 Good
0 Moderate
6 Moderate
0 Non-response
6 Non-response

The model fit statistics can be viewed in the @stats slot which is a list of items including fitted model coefficients, their standard errors and the results of statistical tests on terms within the model using Wald type 2 Chi square. To see the most significant interactions we can order pvals by Timepoint:EULAR_6m:

stats <- summary(results)

kable(stats[order(stats[, 'P_Timepoint:EULAR_6m']), ]) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
Dispersion AIC logLik meanExp (Intercept) Timepoint EULAR_6mModerate EULAR_6mNon-response Timepoint:EULAR_6mModerate Timepoint:EULAR_6mNon-response se_(Intercept) se_Timepoint se_EULAR_6mModerate se_EULAR_6mNon-response se_Timepoint:EULAR_6mModerate se_Timepoint:EULAR_6mNon-response Chisq_Timepoint Chisq_EULAR_6m Chisq_Timepoint:EULAR_6m Df_Timepoint Df_EULAR_6m Df_Timepoint:EULAR_6m P_Timepoint P_EULAR_6m P_Timepoint:EULAR_6m
IGHV3-23 3.2625571 1587.9550 -785.9775 5.9461793 6.8466618 -0.0213954 0.1369640 -0.8391744 -0.3010409 -0.4047322 0.4780053 0.0906367 0.5132355 0.6319642 0.1277273 0.1623995 7.6196498 14.2612080 8.8446578 1 2 2 0.0057736 0.0008002 0.0120062
CXCL13 4.4416343 953.3999 -468.7000 2.8488723 3.5376960 -0.0405699 0.4921163 -1.4651396 -0.3181046 0.0958175 0.3701901 0.0971787 0.5619716 0.7302592 0.1464648 0.1781329 4.2900182 4.3648750 6.7344828 1 2 2 0.0383368 0.1127663 0.0344846
FGF14 0.2915518 375.0390 -179.5195 1.0784074 -0.0846561 0.1148791 -0.0436742 0.9233238 -0.0035230 -0.1731267 0.2014575 0.0467465 0.3090738 0.3176553 0.0711392 0.0783921 5.5445601 5.0815288 5.6637074 1 2 2 0.0185382 0.0788061 0.0589036
IL24 0.4177215 376.8418 -180.4209 1.0076111 0.3478985 -0.0321855 0.2217543 -0.7615232 -0.1166349 0.1298076 0.1978092 0.0500220 0.2721123 0.4427618 0.0793378 0.1003294 1.9439213 1.0333167 5.5708465 1 2 2 0.1632434 0.5965105 0.0617030
HHIP 1.0353883 537.1812 -260.5906 1.4454128 0.7998605 -0.0399783 -0.1533630 0.5663556 0.1881176 0.0332220 0.2089736 0.0569350 0.3219342 0.3888551 0.0826902 0.0975405 1.0348344 4.7912275 5.5680299 1 2 2 0.3090259 0.0911167 0.0617899
MS4A1 3.7894278 898.6829 -441.3414 2.5917399 2.7363666 0.0490839 0.1933006 -1.7237593 -0.2107974 0.1746154 0.3366877 0.0894116 0.5110370 0.6777978 0.1347790 0.1655775 0.0112465 5.6829474 5.4378669 1 2 2 0.9155431 0.0583396 0.0659451
ADAM12 2.3919119 635.7162 -309.8581 1.6675082 1.6569079 -0.1236820 -0.2683597 -0.5047247 -0.0347932 0.2699943 0.2756081 0.0750725 0.4213671 0.5491945 0.1146527 0.1351884 2.5418566 2.3034269 5.2104279 1 2 2 0.1108644 0.3160947 0.0738873
ADAM12.1 2.3919119 635.7162 -309.8581 1.6675082 1.6569079 -0.1236820 -0.2683597 -0.5047247 -0.0347932 0.2699943 0.2756081 0.0750725 0.4213671 0.5491945 0.1146527 0.1351884 2.5418566 2.3034269 5.2104279 1 2 2 0.1108644 0.3160947 0.0738873
IL2RG 0.8332758 1077.5836 -530.7918 4.3589131 3.4913730 -0.0544556 0.2592244 -0.6099845 -0.1005393 0.0741228 0.1768206 0.0432425 0.2498527 0.3251585 0.0650989 0.0795733 6.8218763 3.0209455 4.9774475 1 2 2 0.0090048 0.2208056 0.0830158
IL16 0.2419416 974.1912 -479.0956 4.5153510 3.2821854 -0.0163063 -0.1351457 -0.3377998 0.0421534 0.0983187 0.0936428 0.0242321 0.1389855 0.1814161 0.0364868 0.0447040 1.0986548 0.2976587 4.9710429 1 2 2 0.2945615 0.8617162 0.0832821
EMILIN3 2.8558506 522.1530 -253.0765 1.2381694 1.0154229 0.1211976 -2.0741731 -0.8821688 0.2905788 0.0936821 0.3076525 0.0803497 0.5615737 0.6368474 0.1317963 0.1513720 15.7455706 9.3862545 4.8709433 1 2 2 0.0000725 0.0091580 0.0875564
BLK 1.0360018 532.3068 -258.1534 1.4610649 0.9405211 0.0152304 0.2889937 -0.6186524 -0.1477599 0.0407170 0.2366606 0.0545010 0.3164346 0.4438536 0.0838779 0.1061424 0.6098285 2.2049332 4.1711544 1 2 2 0.4348524 0.3320510 0.1242354
IGHV1-69 6.1631087 1150.4929 -567.2465 3.6446738 4.5052838 0.0254377 1.0678755 0.3278992 -0.1788186 -0.3840284 0.4260683 0.1132972 0.6469940 0.8341386 0.1701383 0.2070201 2.1904029 4.5564214 3.6033230 1 2 2 0.1388730 0.1024674 0.1650245
PADI4 2.0736630 472.8051 -228.4025 1.0912063 0.2320863 0.1335275 0.3673813 0.1562910 -0.2084779 -0.1075547 0.2903618 0.0735494 0.4303320 0.5600246 0.1128728 0.1359462 0.6800222 0.2066688 3.4338576 1 2 2 0.4095791 0.9018254 0.1796169
LILRA5 0.8432028 787.5680 -385.7840 2.8996321 2.1304753 -0.0337206 0.3077720 -0.0538474 -0.0772630 0.0741083 0.1682072 0.0451020 0.2531459 0.3301617 0.0677630 0.0815717 2.3877871 0.6639472 3.3402216 1 2 2 0.1222866 0.7175063 0.1882262
IGHV5-10-1 5.5424948 1103.9302 -543.9651 3.3881684 4.8972357 -0.0222775 0.3631187 -0.0880881 -0.2792942 -0.0040133 0.6692493 0.1111527 0.6495944 0.8665262 0.1672614 0.2171240 2.8323786 0.3035822 3.0812410 1 2 2 0.0923814 0.8591678 0.2142481
IGHV3OR16-8 4.4604808 728.3564 -356.1782 1.8706896 2.9350513 -0.1123000 -0.2194377 -0.3427260 -0.1967505 0.1245430 0.4521884 0.0999641 0.5628201 0.7216100 0.1486832 0.1890428 5.6767222 2.1922516 3.0625772 1 2 2 0.0171915 0.3341632 0.2162568
PHOSPHO1 1.5924621 697.7595 -340.8798 2.1552363 1.4569268 0.0953195 0.0654148 0.0650941 -0.1432879 -0.0867006 0.2317083 0.0605645 0.3511829 0.4523096 0.0922445 0.1113749 0.4904948 1.0829806 2.4703254 1 2 2 0.4837067 0.5818804 0.2907874
S100B 3.9902943 1068.8305 -526.4153 3.4429454 2.6698463 -0.0027976 0.4409051 -0.0676109 -0.1150445 0.1545318 0.3455394 0.0918901 0.5240098 0.6768203 0.1380287 0.1673427 0.0514062 0.6091472 2.4356538 1 2 2 0.8206342 0.7374378 0.2958724
IGHV2-70D 6.4682516 857.1816 -420.5908 2.3395373 3.7263964 -0.1188104 0.3457869 -0.8728989 -0.2362598 -0.2549742 0.4369532 0.1163295 0.6635574 0.8571650 0.1751086 0.2152746 10.2774075 5.8448827 2.3830046 1 2 2 0.0013467 0.0538022 0.3037646
PPIL4 0.1705260 779.2111 -381.6056 3.2865169 2.2892523 -0.0074728 -0.0643493 -0.0248366 0.0481031 0.0498531 0.0961397 0.0244635 0.1450737 0.1868842 0.0362900 0.0440609 1.5034161 0.5849859 2.2273088 1 2 2 0.2201465 0.7464005 0.3283568
IGHJ3P 4.7990172 1085.2659 -534.6329 3.4861371 4.9831233 -0.0810362 -0.0437041 -1.6508368 -0.2109033 -0.0303276 0.3759482 0.0999974 0.5711212 0.7377153 0.1503448 0.1829902 5.8423310 9.6859013 2.1042886 1 2 2 0.0156451 0.0078838 0.3491882
IGHV7-4-1 11.6864201 1101.3138 -542.6569 3.0586094 5.3345104 -0.0825371 0.3334437 -2.6497950 -0.2951203 -0.3087896 0.5864586 0.1559470 0.8909674 1.1504345 0.2343129 0.2878264 5.5351281 16.9832522 2.0262489 1 2 2 0.0186384 0.0002052 0.3630828
LILRB3 0.4876740 945.3343 -464.6671 3.9474519 2.8953372 -0.0258468 0.2114848 -0.1435780 -0.0203999 0.0663273 0.1263695 0.0337764 0.1909106 0.2488508 0.0505310 0.0614438 0.8158759 1.2121560 1.9058338 1 2 2 0.3663887 0.5454861 0.3856146
IL36RN 32.6380418 500.2509 -242.1254 0.3962182 0.9057894 -0.3373192 -1.8590999 -1.0185520 0.5013853 0.0103691 0.9858262 0.2684549 1.5254572 1.9455415 0.4019638 0.5027282 0.7092212 0.6064456 1.7620524 1 2 2 0.3997025 0.7384345 0.4143575
IGHV3OR16-12 12.1178160 539.6728 -261.8364 0.8508034 1.5129747 -0.0122481 -0.6073245 -0.9450944 -0.3231954 -0.1338251 0.6024031 0.1602493 0.9189713 1.1916242 0.2498543 0.2984656 1.7312008 4.0960551 1.6732438 1 2 2 0.1882576 0.1289891 0.4331713
IGHV5-78 5.5798767 569.8859 -276.9429 1.1452612 1.0630323 0.0397579 -0.0398071 0.6148506 -0.1990246 -0.1458526 0.4174604 0.1105884 0.6345936 0.8092159 0.1691995 0.2013966 0.5850486 1.7295015 1.4883378 1 2 2 0.4443402 0.4211565 0.4751290
IGHV3-49 10.8631558 1271.1690 -627.5845 3.9774695 5.7635601 -0.1677971 -0.4705135 0.4518069 -0.2392901 -0.2599215 0.5651950 0.1503341 0.8586770 1.1064226 0.2259680 0.2742812 9.0581025 2.5620382 1.4850886 1 2 2 0.0026153 0.2777541 0.4759015
IGHV1OR15-4 4.1384009 737.1263 -360.5631 1.9435728 2.4925096 -0.0853148 0.8616510 -0.8866376 -0.1479020 0.0087225 0.3523456 0.0940639 0.5335321 0.6969450 0.1411081 0.1736233 4.7041753 6.3141697 1.3100295 1 2 2 0.0300894 0.0425496 0.5194344
IL2RA 1.1303114 545.7334 -264.8667 1.5164121 1.3234090 -0.1238205 0.0366296 -0.7791813 -0.0168109 0.1045125 0.2320172 0.0578359 0.3155436 0.4380150 0.0864420 0.1086980 8.0196429 2.5247489 1.2479916 1 2 2 0.0046273 0.2829813 0.5357992
EMILIN2 0.2149030 1022.9549 -503.4775 4.8723536 3.4282169 -0.0478693 0.1907900 0.2159081 0.0307185 0.0370946 0.0852932 0.0229787 0.1287306 0.1654152 0.0341139 0.0412203 3.7233248 9.7663101 1.1782536 1 2 2 0.0536574 0.0075731 0.5548115
CXCL11 4.3308929 558.9347 -271.4674 1.2009946 1.4333300 -0.0924072 -0.1022289 -1.4166330 -0.1214451 0.0706083 0.3665961 0.0986191 0.5577861 0.7597657 0.1509614 0.1891935 3.2760657 4.5591537 1.1209216 1 2 2 0.0702973 0.1023275 0.5709459
IGHV1OR21-1 7.6550929 357.7987 -170.8994 0.4791891 -0.7429909 0.1076602 1.1516429 -0.9438997 -0.2080242 -0.1730854 0.5357089 0.1379441 0.7790846 1.1726698 0.2039689 0.2979740 0.0006757 5.1554051 1.1160091 1 2 2 0.9792623 0.0759483 0.5723500
CXCL9 7.7483868 1111.5068 -547.7534 3.2358534 4.1831439 -0.1589647 -0.1505766 -1.9470603 -0.0796054 0.1738691 0.4778537 0.1271857 0.7259892 0.9396534 0.1911696 0.2327512 3.2587813 4.3825236 1.1034645 1 2 2 0.0710419 0.1117756 0.5759513
IL36G 55.6403102 592.4912 -288.2456 0.3500668 0.9983391 -0.3974202 -2.2326598 -1.8715276 0.5421133 0.1857025 1.2833104 0.3479228 1.9795693 2.5460500 0.5217741 0.6460756 0.5065785 0.7179743 1.0879009 1 2 2 0.4766238 0.6983833 0.5804507
SAA1 28.3095252 964.5999 -474.2999 1.7170262 3.3639676 -0.2882037 -1.3730549 -0.6174883 0.0034016 0.4195828 0.9130726 0.2432061 1.3885002 1.7882564 0.3666009 0.4431516 1.5819094 2.4806886 1.0261759 1 2 2 0.2084859 0.2892846 0.5986441
EAF2 0.5728653 724.8370 -354.4185 2.6574912 2.0987539 -0.0752532 0.2348340 -0.3563597 -0.0430360 0.0245402 0.1639909 0.0392425 0.2172068 0.2901849 0.0587959 0.0730532 10.6659078 3.5062721 0.9426974 1 2 2 0.0010913 0.1732298 0.6241599
IGHV7-27 14.5870583 478.5468 -231.2734 0.5865966 0.9619084 -0.4094644 -0.3547719 -0.9275209 0.0653310 0.3182143 0.6635347 0.1904673 1.0110845 1.3177256 0.2856313 0.3366538 6.4310379 0.0762884 0.9125222 1 2 2 0.0112143 0.9625741 0.6336484
CILP 3.4124368 1167.6923 -575.8462 3.9198405 3.5693806 0.1148880 -0.8956630 0.2112706 0.0958713 -0.0338927 0.4055634 0.0914031 0.5413809 0.6810776 0.1336388 0.1646131 5.5763914 2.8852445 0.7728476 1 2 2 0.0182042 0.2363073 0.6794825
IGHJ6 6.8375817 1218.9556 -601.4778 3.8318813 6.0208110 -0.0985737 -0.2121636 -0.0204641 -0.1566502 -0.1585477 0.9116274 0.1323751 0.8044644 1.0778771 0.2004475 0.2554268 2.4568288 0.9107552 0.7333773 1 2 2 0.1170148 0.6342084 0.6930254
FGFR2 4.6484751 484.1123 -234.0562 0.8968766 0.8562658 -0.0725830 -0.9483749 -1.1838893 0.1214805 0.1107619 0.3862805 0.1041300 0.6084193 0.8076024 0.1594249 0.1983454 0.0203936 3.1938923 0.6806885 1 2 2 0.8864432 0.2025140 0.7115253
IGLV4-3 12.6762412 578.4695 -281.2348 0.9284219 1.1862191 0.0748321 -0.4837408 -1.0268217 -0.1506973 -0.2101672 0.6179043 0.1638814 0.9425475 1.2285000 0.2483871 0.3082768 0.0235003 3.2907566 0.6188545 1 2 2 0.8781632 0.1929396 0.7338671
SAA2 59.8447169 871.9089 -427.9545 0.8633720 2.4771414 -0.3747888 -2.0207012 -1.8018214 0.1545014 0.3539453 1.3275257 0.3542867 2.0218950 2.6061812 0.5349216 0.6463870 1.1134766 1.1184541 0.3091783 1 2 2 0.2913275 0.5716507 0.8567671
FGFRL1 0.4980641 705.6527 -344.8264 2.5933500 1.8343637 0.0371596 -0.2538947 -0.0556390 0.0160313 0.0082663 0.1459868 0.0369363 0.2165398 0.2760852 0.0557414 0.0668921 3.1425857 1.7555598 0.0832086 1 2 2 0.0762728 0.4157048 0.9592493
FGF2 0.5313473 506.2867 -245.1434 1.5782689 0.6276046 0.0754127 0.1511292 0.0584818 0.0153345 0.0146872 0.1770026 0.0443837 0.2634762 0.3428861 0.0652600 0.0802985 8.2319662 0.9794955 0.0656300 1 2 2 0.0041159 0.6127809 0.9677176
IGHV7-40 15.9375305 492.4418 -238.2209 0.6102728 1.3965106 -0.2999772 -1.0354265 -1.8729132 0.0202320 0.0876576 0.6899511 0.1877134 1.0563185 1.3922601 0.2887295 0.3614421 4.5933447 2.9065540 0.0589744 1 2 2 0.0320963 0.2338029 0.9709433

Summary statistics on a single gene can be shown specifying the argument gene.

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

Estimated means based on each gene’s fitted model to show fixed effects and their 95% confidence intervals can be seen in the @predict slot:

predict = data.frame(results@predict)
kable(predict) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
y_0_Good y_6_Good y_0_Moderate y_6_Moderate y_0_Non.response y_6_Non.response LCI_0_Good LCI_6_Good LCI_0_Moderate LCI_6_Moderate LCI_0_Non.response LCI_6_Non.response UCI_0_Good UCI_6_Good UCI_0_Moderate UCI_6_Moderate UCI_0_Non.response UCI_6_Non.response RE_var
MS4A1 15.4308165 20.7152501 18.7213987 7.0949707 2.7527690 10.5360538 7.9761734 9.1362857 8.8122864 2.8391654 0.8690196 3.2873828 29.852673 46.968933 39.7729664 17.730073 8.719869 33.768026 0.0000000
ADAM12 5.2430737 2.4963164 4.0090258 1.5491357 3.1650953 7.6145225 3.0548060 1.2417443 2.1464872 0.6925006 1.2474752 2.9782353 8.998877 5.018421 7.4877166 3.465443 8.030483 19.468225 0.0000000
IGHV7-4-1 207.3712020 126.3793288 289.4417585 30.0243900 14.6540308 1.4003960 65.6968833 30.2738485 77.7576114 6.1747931 2.1061275 0.1747709 654.564010 527.575301 1077.4061868 145.990964 101.959933 11.221028 0.0000000
IGHV3-49 318.4801306 116.3703172 198.9486119 17.2968322 500.3796107 38.4383496 105.1910128 29.3400265 56.0316118 3.7578024 77.5451612 5.4682057 964.242010 461.555504 706.3967807 79.615790 3228.824994 270.199551 0.0000000
IGHV3-23 940.7352958 827.4010136 1078.8228597 155.8677956 406.4607522 31.5231393 368.6219183 216.2096233 480.0138327 50.9357383 116.6439723 4.5858142 2400.787509 3166.336572 2424.6358821 476.969030 1416.364171 216.691794 0.0866606
ADAM12.1 5.2430737 2.4963164 4.0090258 1.5491357 3.1650953 7.6145225 3.0548060 1.2417443 2.1464872 0.6925006 1.2474752 2.9782353 8.998877 5.018421 7.4877166 3.465443 8.030483 19.468225 0.0000000
FGFRL1 6.2611492 7.8249818 4.8572335 6.6833182 5.9222990 7.7778511 4.7031366 5.4424407 3.4727930 4.5342782 3.7175543 4.8018105 8.335286 11.250530 6.7935860 9.850905 9.434597 12.598366 0.0058745
IL36G 2.7137709 0.2500280 0.2910324 0.6933905 0.4176179 0.1172444 0.2193802 0.0099191 0.0151667 0.0211411 0.0056096 0.0010308 33.569820 6.302378 5.5845918 22.741997 31.090474 13.335887 0.0000000
BLK 2.5613158 2.8064030 3.4195700 1.5439407 1.3797036 1.9300645 1.6106958 1.6207459 2.1312931 0.8267431 0.6368202 0.8505197 4.072984 4.859428 5.4865560 2.883305 2.989198 4.379850 0.0605014
SAA1 28.9036429 5.1281515 7.3222139 1.3259116 15.5876561 34.2865737 4.8275799 0.5508631 0.9424508 0.1098705 0.7654772 1.4752416 173.051632 47.739514 56.8887158 16.001040 317.416396 796.865491 0.0000000
CILP 35.4946022 70.7185951 14.4937515 51.3298996 43.8445868 71.2805130 16.0302311 29.1752956 5.4851105 19.4940816 13.9477856 21.5621746 78.593176 171.416247 38.2980127 135.156847 137.824587 235.640034 0.2045698
EMILIN3 2.7605305 5.7122239 0.3468891 4.1037785 1.1425403 4.1476108 1.5104748 2.7596693 0.1381293 1.8195270 0.3830255 1.4660145 5.045121 11.823700 0.8711551 9.255701 3.408124 11.734315 0.0000000
EMILIN2 30.8216346 23.1270044 37.3005036 33.6529884 38.2492888 35.8547740 26.0766992 18.7030942 30.8774307 26.7788797 28.9721502 26.7968848 36.429962 28.597318 45.0596936 42.291673 50.497049 47.974413 0.0000000
IGHJ6 411.9125386 228.0053388 333.1681965 72.0459673 403.5688022 86.2818082 68.9941242 24.0959355 61.8136434 6.2157409 64.6730770 4.6567131 2459.223035 2157.477324 1795.7370099 835.076861 2518.324248 1598.670633 0.7559507
CXCL9 65.5716793 25.2633673 56.4055326 13.4791834 9.3566147 10.2318948 25.7015333 7.8711053 19.3240801 3.7026558 1.9161087 1.9545686 167.291386 81.086163 164.6434960 49.069748 45.689600 53.562547 0.0000000
IGHV3OR16-12 4.5402163 4.2185263 2.4735396 0.3305456 1.7645228 0.7345045 1.3941217 0.9711074 0.6347242 0.0547623 0.2352117 0.0838919 14.786058 18.325434 9.6394588 1.995175 13.237184 6.430858 0.0000000
IGHV1OR21-1 0.4756891 0.9075291 1.5047880 0.8240433 0.1850942 0.1250000 0.1664633 0.2635245 0.4965508 0.2081891 0.0239577 0.0120623 1.359339 3.125360 4.5602320 3.261686 1.430015 1.295362 0.0000000
IGHV1-69 90.4940260 105.4155442 263.2645407 104.8860184 125.6101523 14.6090165 39.2593594 37.3270213 101.3620625 33.2848976 30.8036086 3.3414811 208.591502 297.704895 683.7688249 330.512564 512.209806 63.870887 0.0000000
IL2RA 3.7562045 1.7869063 3.8963438 1.6757315 1.7232770 1.5347682 2.3837031 0.9628602 2.4156874 0.8829025 0.8100611 0.6772090 5.918972 3.316197 6.2845447 3.180505 3.666000 3.478266 0.0393846
IGHV2-70D 41.5291831 20.3592259 58.6849465 6.9710229 17.3483512 1.8418858 17.6364566 7.0058075 22.0505583 2.1252119 4.0882652 0.3851986 97.790225 59.164925 156.1830267 22.866031 73.616871 8.807257 0.0000000
IGHV5-10-1 133.9190766 117.1635409 192.5497180 31.5295171 122.6270467 104.7318121 36.0716441 25.2976774 51.2844100 7.4548141 17.3970482 7.9921043 497.186073 542.630657 722.9369300 133.351474 864.364600 1372.448614 0.1424450
S100B 14.4377496 14.1974270 22.4378966 11.0640242 13.4938676 33.5367082 7.3344982 6.1165215 10.3670359 4.3513002 4.3127135 10.2548290 28.420296 32.954504 48.5634668 28.132426 42.220394 109.676212 0.0000000
IL24 1.4160885 1.1674076 1.7676532 0.7237789 0.6612491 1.1878059 0.9609745 0.7059609 1.1821727 0.3847167 0.2982333 0.5877234 2.086743 1.930476 2.6430975 1.361666 1.466135 2.400590 0.0547142
IGHV3OR16-8 18.8224682 9.5950892 15.1138813 2.3662605 13.3608027 14.3792139 7.7582904 2.9292861 6.1523549 0.7459776 3.7615760 2.2120645 45.665384 31.429411 37.1287763 7.505840 47.456452 93.470057 0.0135151
FGFR2 2.3543526 1.5231307 0.9120056 1.2229624 0.7206342 0.9061504 1.1042386 0.5817560 0.3629698 0.4152740 0.1794751 0.2196149 5.019727 3.987801 2.2915247 3.601567 2.893514 3.738856 0.0000000
FGF14 0.9188282 1.8305534 0.8795628 1.7156739 2.3132829 1.6309814 0.6190850 1.2487274 0.5555778 1.1138378 1.4294503 0.9296847 1.363698 2.683473 1.3924795 2.642698 3.743591 2.861293 0.0000000
IGHJ3P 145.9294540 89.7390599 139.6890923 24.2346226 28.0022961 14.3550762 69.8440468 35.8886856 60.1435508 8.7705948 8.0704125 3.8967415 304.899366 224.391023 324.4411452 66.964322 97.160906 52.882187 0.0000000
FGF2 1.8731183 2.9449178 2.1787119 3.7554829 1.9859282 3.4099108 1.3240292 1.9939340 1.4861790 2.4858794 1.1168331 1.9938050 2.649921 4.349462 3.1939527 5.673506 3.531334 5.831810 0.0000000
IL16 26.6339143 24.1515315 23.2670852 27.1702535 18.9989853 31.0767533 22.1679042 19.2600883 18.9118709 21.1415036 13.9463148 22.4649620 31.999660 30.285244 28.6252616 34.918173 25.882210 42.989816 0.0007472
LILRA5 8.4188674 6.8767861 11.4529494 5.8846392 7.9775228 10.1650340 6.0544235 4.5387996 7.9046738 3.6963449 4.5712713 5.7288426 11.706701 10.419096 16.5939864 9.368438 13.921919 18.036438 0.0000000
CXCL11 4.1926376 2.4082175 3.7852093 1.0491572 1.0168372 0.8921720 2.0437813 0.9690135 1.6605085 0.3627796 0.2759238 0.2244880 8.600828 5.984965 8.6285676 3.034159 3.747259 3.545717 0.0000000
IGHV7-40 4.0410745 0.6680766 1.4348842 0.2678342 0.6210134 0.1737184 1.0451981 0.1166984 0.2991986 0.0345152 0.0580424 0.0110776 15.624103 3.824614 6.8813581 2.078364 6.644413 2.724253 0.0000000
PHOSPHO1 4.2927466 7.6052889 4.5829437 3.4367675 4.5814744 4.8246311 2.7258452 4.3937329 2.7321980 1.8230721 2.1396775 2.1835113 6.760352 13.164300 7.6873540 6.478828 9.809847 10.660382 0.0000000
EAF2 8.1560005 5.1926003 10.3148846 5.0725952 5.7110002 4.2127523 5.9140547 3.5052783 7.3784966 3.3609507 3.4581910 2.4307316 11.247841 7.692142 14.4198539 7.655936 9.431383 7.301210 0.0052440
IGHV5-78 2.8951365 3.6750970 2.7821533 1.0699636 5.3542086 2.8329389 1.2773783 1.3373736 1.0902909 0.3290647 1.3759831 0.6715812 6.561733 10.099151 7.0993692 3.479018 20.834231 11.950219 0.0000000
CXCL13 34.3875981 26.9579033 56.2503527 6.5388759 7.9451079 11.0678220 16.6452238 10.9322846 24.5594757 2.3810038 2.3135587 3.0653275 71.041815 66.475451 128.8342720 17.957509 27.284694 39.962022 0.1893666
IGHV7-27 2.6166853 0.2242760 1.8351688 0.2327794 1.0349855 0.5986291 0.7127535 0.0362055 0.4114237 0.0312680 0.1111345 0.0552486 9.606466 1.389283 8.1858315 1.732965 9.638730 6.486256 0.0000000
PPIL4 9.8675569 9.4349003 9.2525855 11.8069120 9.6254985 12.4124171 8.1728568 7.4612810 7.4435364 9.1550634 7.0032562 9.0262468 11.913665 11.930571 11.5012990 15.226893 13.229592 17.068899 0.0341781
IGLV4-3 3.2746765 5.1305435 2.0187495 1.2805468 1.1728039 0.5206835 0.9754334 1.1457207 0.5002922 0.2351748 0.1463546 0.0544524 10.993581 22.974601 8.1459389 6.972687 9.398195 4.978869 0.0000000
IL36RN 2.4738840 0.3268918 0.3854628 1.0315744 0.8933627 0.1256238 0.3582837 0.0269224 0.0393625 0.0708578 0.0333670 0.0029092 17.081719 3.969126 3.7746972 15.018052 23.918772 5.424677 0.0000000
IL2RG 32.8309945 23.6802048 42.5464916 16.7873928 17.8390265 20.0733303 23.2151411 15.4606682 29.0284587 10.4067275 10.2231705 10.8570511 46.429793 36.269590 62.3596301 27.080228 31.128393 37.113079 0.0407158
IGHV1OR15-4 12.0915833 7.2472419 28.6215692 7.0629279 4.9822025 3.1465877 6.0612270 3.0540945 13.0515653 2.7158829 1.5330821 0.9037643 24.121582 17.197410 62.7659755 18.367858 16.191136 10.955305 0.0000000
HHIP 2.2252305 1.7506563 1.9088433 4.6428720 3.9204878 3.7647380 1.4773829 1.0307931 1.1811761 2.7695040 2.0615799 1.9186382 3.351637 2.973242 3.0847922 7.783437 7.455556 7.387142 0.0000000
PADI4 1.2612285 2.8101880 1.8211490 1.1615614 1.4745861 1.7232541 0.7138925 1.4650902 0.9772170 0.5264347 0.5768283 0.6580193 2.228203 5.390219 3.3939070 2.562948 3.769587 4.512945 0.0000000
SAA2 11.9071780 1.2565988 1.5784451 0.4209319 1.9646616 1.7337022 0.8826655 0.0485264 0.0794400 0.0110057 0.0242254 0.0175383 160.628102 32.539848 31.3631500 16.099335 159.332581 171.380238 0.0000000
LILRB3 18.0895994 15.4909569 22.3499024 16.9343137 15.6701735 19.9782091 14.1208572 11.3557724 16.8836563 12.0341978 10.2941966 12.9525109 23.173778 21.131962 29.5858983 23.829672 23.853667 30.814785 0.0000000

Q-values

The q-values from each of the p-value columns can be calculated using glmmQvals. This will print a significance table based on the cut-off (default p=0.05) and add a matrix named qvals to the @stats slot:

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)

Refitting a different model

glmmRefit() allows a different model to be fitted using the original data. The example below shows how to refit the model without the interaction term and then perform a likelihood ratio test using anova. Note for glmmSeq() objects the LHS of the reduced formula must be "count", while for lmmSeq() objects the LHS must be "gene". For glmmTMB analyses, the GLM family can also be changed.

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")

Fold change plots

The comparative fold change (for x1var variables) between conditions (x2var and x2Values variables) can be plotted using an fcPlot for all genes to highlight significance. This type of plot is most suited to look for interaction between time (x1var) and a two-level factor (x2var), looking at change between two timepoints. In the example below from the R4RA study [1], gene expression pre- and post-drug treatment is compared between two drugs (rituximab & tocilizumab), using the design formula gene_counts ~ time * drug. By setting graphics = "plotly" this can be viewed interactively.

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'))

Log2 fold change between the two time points for individuals treated with rituximab on the x axis and individuals treated with tocilizumab on the y axis with each point representing a gene. Genes showing an interaction between time and drug are coloured blue or gold depending on whether their fold change is greater post-rituximab (blue) or post-tocilizumab (gold). Genes without interaction, but changing significantly over time are coloured green and tend to lie along the line of identity. See the Longitudinal tab in https://r4ra.hpc.qmul.ac.uk for an interactive version of the above plot.

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.

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:

  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. glmmTMB: Mollie Brooks, Ben Bolker, Kasper Kristensen, Martin Maechler, Arni Magnusson (2022). glmmTMB: Generalized Linear Mixed Models using Template Model Builder

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

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

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