glmmSeq

Myles Lewis, Katriona Goldmann, Elisabetta Sciacca, Cankut Cubuk, Anna Surace

2021-03-30

Lifecycle: Maturing License: MIT CRAN status Downloads 2021-03-30 GitHub issues Travis

glmmSeq

glmmSeq logo

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.

Installing from CRAN

install.packages("glmmSeq")

Installing from Github

devtools::install_github("KatrionaGoldmann/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", "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:

These are outlined in the following subsections.

Metadata

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()
SAMID PATID Timepoint EULAR_6m EULAR_binary
SAM20389187 PAT300 6 Good responder responder
SAM20389185 PAT209 6 Good responder responder
SAM20389181 PAT219 6 Moderate responder responder
SAM20389179 PAT211 6 Good responder responder
SAM20389177 PAT216 6 Good responder responder
SAM20389176 PAT212 6 Good responder responder

Count data

tpm = tpm[, metadata$SAMID]
kable(head(tpm)) %>% kable_styling() %>%
  scroll_box(width = "100%")
SAM20389176 SAM20389177 SAM20389178 SAM20389179 SAM20389185 SAM20389187 SAM24297965 SAM24297969 SAM24297977 SAM24297980 SAM24297982 SAM24297983 SAM24297984 SAM24297986 SAM24297987 SAM24297991 SAM24297994 SAM24297995 SAM24297996 SAM24297997 SAM24298029 SAM24298035 SAM24298052 SAM24298065 SAM24298074 SAM24298089 SAM24298092 SAM24298098 SAM24298100 SAM24298111 SAM24298112 SAM24298118 SAM24298120 SAM24298123 SAM24298125 SAM24298127 SAM24298133 SAM24298135 SAM24298136 SAM24298143 SAM24298148 SAM24298152 SAM24298157 SAM24298159 SAM24298160 SAM24298162 SAM24298163 SAM24298165 SAM9103802 SAM9103804 SAM9103805 SAM9103806 SAM9103807 SAM9103808 SAM9103809 SAM9103810 SAM9103811 SAM9103812 SAM9103813 SAM9103815 SAM9103816 SAM9103817 SAM9103818 SAM9103819 SAM9103820 SAM9103821 SAM9103822 SAM9103823 SAM9103824 SAM9103825 SAM9103826 SAM9103827 SAM9103828 SAM9103829 SAM9103830 SAM9103832 SAM9103833 SAM9103834 SAM9103835 SAM9103836 SAM9103837 SAM9103838 SAM9103839 SAM9103840 SAM9103843 SAM9103844 SAM9103845 SAM9103847 SAM9103849 SAM9103850 SAM9103851 SAM9103852 SAM9103853 SAM9103854 SAM9103856 SAM9103859 SAM9103860 SAM9103861 SAM9103862 SAM9103863 SAM9103865 SAM9103867 SAM9103868 SAM9103869 SAM9103870 SAM9103875 SAM9103876 SAM9103878 SAM9103880 SAM9103881 SAM9103885 SAM9103887 SAM9103888 SAM9103889 SAM9103890 SAM9103891 SAM9103892 SAM9103896 SAM9103897 SAM9103898 SAM9103901 SAM9103902 SAM9103904 SAM9103905
MS4A1 57.280078 2.123731 1.200281 15.742876 2.024451 1.622818 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 2.957004 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 0.972897 2.645153 1.166610 5.466742 1.110548 2.581805 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 2.232818 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 2158.300000 0.000000 1.030500 26.381700 0.000000 0.992885 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.205779 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 858.980000 2.093480 8.066680 345.493000 0.196831 0.000000 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 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 5180.460000 5.962460 23.445900 2349.850000 3.999940 0.715133 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 1.102910 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 0.972897 2.645153 1.166610 5.466742 1.110548 2.581805 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 2.232818 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.809307  2.386833 11.790220 10.959627  3.302796  2.386833

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)
##     MS4A1    ADAM12 IGHV7-4-1  IGHV3-49  IGHV3-23  ADAM12.1 
## 0.7951463 2.1345465 3.2851952 2.3843452 1.9764054 2.1345465

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

head(sizeFactors)
## SAM20389176 SAM20389177 SAM20389178 SAM20389179 SAM20389185 SAM20389187 
##  4.22168329  0.24675696  0.09137974  3.15511184  0.09391590  0.20458970

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

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

Similarly, with DESeq2:

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:

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

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

Outputs

This creates a GlmmSeq object which contains the following slots:

names(attributes(results))
##  [1] "formula"        "stats"          "predict"        "reducedFormula" "countdata"     
##  [6] "metadata"       "modelData"      "optInfo"        "errors"         "variables"     
## [11] "class"

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

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

The model fit statistics can be viewed in the @stats slot. To see the most significant interactions we can order by P_Timepoint.EULAR_6m:

stats = data.frame(results@stats)

kable(stats[order(stats$P_Timepoint.EULAR_6m), ]) %>%
  kable_styling() %>%
  scroll_box(width = "100%", height = "400px")
Dispersion AIC logLik X.Intercept. Timepoint EULAR_6mModerate.responder EULAR_6mNon.responder Timepoint.EULAR_6mModerate.responder Timepoint.EULAR_6mNon.responder Chisq_Timepoint Chisq_EULAR_6m Chisq_Timepoint.EULAR_6m P_Timepoint P_EULAR_6m P_Timepoint.EULAR_6m
IGHV7-4-1 3.2851952 1027.3249 -505.66246 3.2663050 -0.1681991 0.5045693 -1.5043291 -0.1483570 -0.2476215 3.020997e+03 2.715858e+04 8.865098e+03 0.0000000 0.0000000 0.0000000
IGHV3-49 2.3843452 1213.2192 -598.60958 3.8564428 -0.1310463 0.0870902 0.3356671 -0.2535793 -0.5576211 4.351648e+02 2.880120e+03 9.031936e+03 0.0000000 0.0000000 0.0000000
IGHJ6 4.4994184 1215.0594 -599.52969 3.1816009 -0.1150711 0.5914506 0.6796337 -0.1734701 -0.4120140 3.549184e+03 1.820147e+05 5.047856e+04 0.0000000 0.0000000 0.0000000
IGHV3OR16-12 1.4074644 419.5304 -201.76519 -0.0819321 -0.0327344 -0.0809701 -0.4644101 -0.2202112 -0.1084345 1.054081e+03 2.183599e+05 5.920973e+04 0.0000000 0.0000000 0.0000000
IGHV1-69 2.2817464 1146.2559 -565.12795 2.8746624 -0.0636055 0.9844362 0.1969424 -0.2488542 -0.5112757 1.974164e+02 4.550179e+04 1.483963e+04 0.0000000 0.0000000 0.0000000
IGHV5-10-1 3.0087428 1098.2063 -541.10315 2.9128345 -0.0895266 0.2780136 -0.5314339 -0.1938497 -0.4042939 1.701713e+03 7.684466e+04 4.225208e+04 0.0000000 0.0000000 0.0000000
IGHV3OR16-8 1.9863992 727.9288 -355.96441 1.3635681 -0.0996723 0.2473497 -0.1721759 -0.1603650 -0.2265692 6.525887e+03 5.962962e+04 5.087477e+04 0.0000000 0.0000000 0.0000000
IGHJ3P 2.3867019 1098.9404 -541.47021 3.0329762 -0.0905663 0.6639360 -1.1837380 -0.2503756 -0.2047263 2.442162e+02 5.266089e+04 2.999977e+03 0.0000000 0.0000000 0.0000000
IGHV7-27 1.9313574 327.0458 -155.52289 -0.1858293 -0.3179804 0.0917328 -0.3285713 0.0480037 0.2211036 1.312878e+04 1.374890e+04 6.739602e+03 0.0000000 0.0000000 0.0000000
CXCL9 2.2637545 1022.3385 -503.16927 3.4984207 -0.2086260 0.1620894 -1.4831678 -0.0631385 0.1264203 1.630773e+03 7.878927e+04 7.107016e+02 0.0000000 0.0000000 0.0000000
IGLV4-3 1.9003322 439.7279 -211.86395 0.5249324 -0.0841200 -0.2654372 -0.8750102 -0.1689169 -0.0451954 1.163595e+01 1.987154e+03 8.156976e+01 0.0006469 0.0000000 0.0000000
IGHV7-40 2.0988009 341.0842 -162.54208 0.0161516 -0.2289910 -0.4289696 -1.0192924 0.0224127 0.0337400 1.976403e+03 4.776836e+04 6.226033e+01 0.0000000 0.0000000 0.0000000
CXCL13 2.3526243 946.8630 -465.43149 2.7085891 -0.1004039 0.8843805 -1.3714223 -0.3199418 0.0436665 1.572207e+01 6.276206e+00 9.471960e+00 0.0000734 0.0433650 0.0087738
MS4A1 0.7951463 856.8937 -420.44683 1.8105814 0.0321673 0.6605145 -1.0688655 -0.2257699 0.0089804 1.888090e+00 6.875578e+00 8.710286e+00 0.1694175 0.0321357 0.0128406
IGHV3-23 1.9764054 1604.2499 -794.12493 5.1154581 -0.0825775 0.5813393 -0.5692484 -0.2591130 -0.4769502 1.734050e+01 6.170885e+00 8.613146e+00 0.0000312 0.0457098 0.0134797
ADAM12 2.1345465 633.6054 -308.80269 1.6569079 -0.1244486 -0.2683597 -0.5047247 -0.0327808 0.2707610 2.852569e+00 2.512715e+00 5.765926e+00 0.0912281 0.2846891 0.0559687
ADAM12.1 2.1345465 633.6054 -308.80269 1.6569079 -0.1244486 -0.2683597 -0.5047247 -0.0327808 0.2707610 2.852569e+00 2.512715e+00 5.765926e+00 0.0912281 0.2846891 0.0559687
BLK 0.4615021 524.2257 -254.11284 0.7634585 0.0148447 0.3267961 -0.5725879 -0.1495720 0.0204334 1.156196e+00 2.708080e+00 5.656018e+00 0.2822558 0.2581950 0.0591305
IL24 0.3842330 378.3689 -181.18447 0.3610892 -0.0377177 0.2145140 -0.7680385 -0.1072013 0.1365525 2.158975e+00 1.015672e+00 5.520695e+00 0.1417392 0.6017966 0.0632698
IL2RG 0.7828466 1083.3953 -533.69765 3.4904554 -0.0593842 0.2554795 -0.6133692 -0.0926156 0.0781855 7.529213e+00 3.210710e+00 4.901794e+00 0.0060706 0.2008182 0.0862162
EMILIN3 3.5375635 542.7327 -263.36633 1.0154229 0.1272457 -2.0741731 -0.8821688 0.2842837 0.0876341 1.380526e+01 8.478584e+00 3.985180e+00 0.0002028 0.0144178 0.1363418
IL36RN 4.6331400 289.3523 -136.67615 -0.3329543 -0.2082603 -1.0691994 -0.1998857 0.3315983 -0.0919322 1.144986e+00 2.820906e-01 3.932056e+00 0.2846012 0.8684500 0.1400119
PADI4 1.7980343 470.8495 -227.42475 0.2320863 0.1339818 0.3673813 0.1562910 -0.2058310 -0.1080089 8.650596e-01 2.468016e-01 3.781890e+00 0.3523267 0.8839093 0.1509291
IGHV2-70D 2.2946862 832.2096 -408.10478 1.6581318 -0.1446614 0.7038978 -0.1503531 -0.1048764 -0.3344248 1.613128e+01 2.780130e+00 3.684883e+00 0.0000591 0.2490592 0.1584302
IL36G 4.6410978 271.2917 -127.64583 -0.4583626 -0.2177700 -1.1515527 -0.7112490 0.3382520 0.0129267 8.782735e-01 8.418474e-01 3.155987e+00 0.3486750 0.6564402 0.2063888
S100B 2.2078701 1009.2901 -496.64504 2.6610443 -0.0071322 0.3670854 -0.0726711 -0.0887257 0.1557117 1.322500e-03 8.405754e-01 3.053453e+00 0.9709907 0.6568578 0.2172457
IGHV5-78 1.0112983 511.2719 -247.63593 0.2825675 -0.0327322 0.2959688 0.1957975 -0.1678010 -0.1325774 5.522931e+00 3.734270e-02 2.668774e+00 0.0187688 0.9815019 0.2633196
SAA1 4.5130513 715.7066 -349.85330 2.0021705 -0.1337644 -0.5896751 0.4246927 -0.1498307 0.1721035 8.805368e+00 7.533313e+00 2.427085e+00 0.0030035 0.0231293 0.2971428
IGHV1OR21-1 1.6964380 278.7881 -131.39406 -1.2931431 0.0931411 1.1234556 -0.6743562 -0.1972479 -0.1504572 1.850600e-02 5.809509e+00 2.126558e+00 0.8917922 0.0547622 0.3453216
HHIP 2.6145990 577.8645 -280.93223 0.7998605 -0.0255812 -0.1533630 0.5663556 0.1638981 0.0188249 4.928654e-01 1.998117e+00 2.086007e+00 0.4826519 0.3682259 0.3523946
FGF14 1.7659535 443.4933 -213.74666 -0.0846561 0.1142763 -0.0436742 0.9233238 -0.0060532 -0.1725240 2.525964e+00 1.835437e+00 2.063881e+00 0.1119864 0.3994293 0.3563149
PHOSPHO1 1.8094035 714.8166 -349.40828 1.4569268 0.0922208 0.0654148 0.0650941 -0.1375254 -0.0836019 4.454308e-01 9.091811e-01 2.060105e+00 0.5045128 0.6347078 0.3569881
IGHV1OR15-4 1.5843911 735.1463 -359.57317 1.0647778 -0.0751346 0.8342548 -0.6541352 -0.1184912 -0.2068972 6.429820e+00 5.085462e+00 1.704192e+00 0.0112220 0.0786513 0.4265200
LILRA5 1.3588557 830.4412 -407.22061 2.1304753 -0.0351868 0.3077720 -0.0538474 -0.0580985 0.0755745 1.235490e+00 5.906689e-01 1.647021e+00 0.2663421 0.7442826 0.4388881
CXCL11 1.8277524 519.0822 -251.54109 1.1846624 -0.1234842 0.0213862 -1.2730697 -0.0798814 0.1018420 6.902011e+00 4.759180e+00 1.549515e+00 0.0086099 0.0925885 0.4608155
IL2RA 0.9253515 544.2460 -264.12302 1.2868210 -0.1306894 0.0470235 -0.7623192 -0.0040900 0.1089655 9.494256e+00 2.685496e+00 1.331378e+00 0.0020612 0.2611272 0.5139193
IL16 1.2379626 1091.4116 -537.70582 3.2827704 -0.0125302 -0.1353727 -0.3379663 0.0365392 0.0946826 2.874926e-01 9.216870e-02 1.035468e+00 0.5918319 0.9549614 0.5958693
EAF2 0.6141013 730.3064 -357.15320 2.1048581 -0.0800130 0.2331011 -0.3584361 -0.0420156 0.0294982 1.111203e+01 3.195855e+00 9.641550e-01 0.0008577 0.2023154 0.6174992
FGFR2 2.1014485 429.3980 -206.69898 0.7437064 -0.0639926 -0.8606098 -1.0970974 0.1089732 0.1013839 2.215830e-02 3.846131e+00 9.559309e-01 0.8816667 0.1461583 0.6200436
VIL1 0.6183578 194.2405 -89.12023 -1.1619565 0.0525457 -0.0309993 0.3525958 -0.0219418 -0.1490192 1.105745e-01 7.325100e-02 9.109431e-01 0.7394911 0.9640371 0.6341489
IL26 0.6181075 239.0144 -111.50720 -0.7621345 -0.0074300 0.4808342 0.0683890 -0.0942476 -0.0926694 1.292617e+00 1.245948e+00 8.341110e-01 0.2555663 0.5363469 0.6589844
LILRB3 1.1933120 1004.9665 -494.48327 2.8953372 -0.0256229 0.2114848 -0.1435780 -0.0177522 0.0661034 3.230496e-01 5.747663e-01 7.845889e-01 0.5697807 0.7502242 0.6755052
SAA2 3.5694731 460.5391 -222.26957 0.5710606 -0.1215510 -0.5742976 -0.1356134 -0.0911361 0.0728808 3.267983e+00 2.481261e+00 6.864212e-01 0.0706444 0.2892019 0.7094888
CILP 3.5502875 1187.3836 -585.69181 3.6294582 0.1329390 -0.8723346 0.1801614 0.0731402 -0.0512005 5.994802e+00 2.935150e+00 6.005211e-01 0.0143481 0.2304838 0.7406252
FGF12 0.5461526 206.5879 -95.29396 -0.9119156 -0.0044117 -0.0118226 -0.1815855 -0.0142521 0.0823507 1.859240e-02 8.812260e-02 4.512377e-01 0.8915414 0.9568953 0.7980222
PPIL4 1.4755536 893.5181 -438.75906 2.3066754 -0.0007983 -0.0610189 -0.0166817 0.0393739 0.0425411 3.123808e-01 9.683130e-02 2.811167e-01 0.5762229 0.9527377 0.8688730
EMILIN2 1.2635469 1154.9780 -569.48900 3.4282169 -0.0457832 0.1907900 0.2159081 0.0340239 0.0350085 6.098077e-01 1.944400e+00 2.448315e-01 0.4348602 0.3782501 0.8847804
IL12A 0.4505427 234.8968 -109.44840 -0.5284814 0.0037258 -0.1602920 -0.2927997 -0.0219157 0.0208392 4.580000e-05 6.255961e-01 1.010101e-01 0.9946016 0.7313976 0.9507491
FGFRL1 1.7213855 788.1793 -386.08964 1.8388329 0.0397231 -0.2536428 -0.0581760 0.0118867 0.0063352 1.175615e+00 6.910527e-01 1.634450e-02 0.2782504 0.7078477 0.9918610
FGF2 1.8685417 580.1500 -282.07498 0.6276046 0.0867638 0.1511294 0.0584818 0.0023200 0.0033362 3.747805e+00 2.694922e-01 9.124000e-04 0.0528769 0.8739378 0.9995439

And the final model output 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.responder y_6_Good.responder y_0_Moderate.responder y_6_Moderate.responder y_0_Non.responder y_6_Non.responder LCI_0_Good.responder LCI_6_Good.responder LCI_0_Moderate.responder LCI_6_Moderate.responder LCI_0_Non.responder LCI_6_Non.responder UCI_0_Good.responder UCI_6_Good.responder UCI_0_Moderate.responder UCI_6_Moderate.responder UCI_0_Non.responder UCI_6_Non.responder
MS4A1 6.1140011 7.4155950 11.8354104 3.7042493 2.0995350 2.6874744 3.5888084 4.0080183 6.6660704 1.8429689 0.8519067 1.0337382 10.4159949 13.7202590 21.0134203 7.4453035 5.1743311 6.9867967
ADAM12 5.2430737 2.4848599 4.0090258 1.5607585 3.1650953 7.6145225 3.1403622 1.2959724 2.2145480 0.7229322 1.3053524 3.1283084 8.7537106 4.7643985 7.2575928 3.3695652 7.6744245 18.5342831
IGHV7-4-1 26.2142986 9.5554392 43.4180072 6.4982815 5.8239338 0.4804906 26.0574027 9.2136118 43.0508098 6.1717413 1.5777513 0.1304357 26.3721391 9.9099486 43.7883366 6.8420987 21.4978147 1.7699999
IGHV3-49 47.2968088 21.5454555 51.6005880 5.1334562 66.1622386 1.0619486 46.7008212 19.9640521 50.6848109 4.6090754 64.9879378 0.9530841 47.9004022 23.2521260 52.5329115 5.7174965 67.3577584 1.1832480
IGHV3-23 166.5770642 101.4933031 297.9117769 38.3461579 94.2743997 3.2839385 60.0895171 34.2590265 100.0447863 11.6131390 20.3600073 0.6647136 461.7763569 300.6766865 887.1169614 126.6176033 436.5255042 16.2239071
ADAM12.1 5.2430737 2.4848599 4.0090258 1.5607585 3.1650953 7.6145225 3.1403622 1.2959724 2.2145480 0.7229322 1.3053524 3.1283084 8.7537106 4.7643985 7.2575928 3.3695652 7.6744245 18.5342831
FGFRL1 6.2891940 7.9818619 4.8802191 6.6515452 5.9337534 7.8225232 3.9665551 4.5804448 2.8624987 3.5353573 2.7259706 3.5021629 9.9718674 13.9091555 8.3201918 12.5144506 12.9162905 17.4725938
IL36G 0.6323181 0.1711898 0.1999045 0.4118800 0.3104875 0.0908386 0.2589587 0.0441482 0.0571591 0.1134455 0.0590332 0.0086142 1.5439769 0.6638082 0.6991333 1.4953889 1.6330219 0.9579121
BLK 2.1456843 2.3455666 2.9750314 1.3256347 1.2103029 1.4956184 1.4464980 1.4808256 1.9625497 0.7555843 0.5990408 0.7312496 3.1828327 3.7152806 4.5098537 2.3257597 2.4452976 3.0589754
SAA1 7.4051112 3.3187365 4.1061890 0.7489560 11.3233071 14.2520087 6.5601870 1.7916903 1.4735282 0.2245448 2.9513654 2.8846055 8.3588580 6.1472745 11.4424601 2.4980988 43.4433782 70.4150894
CILP 37.6923885 83.6877254 15.7544605 54.2497822 45.1332644 73.7035225 17.1382414 35.7638670 5.8701696 20.2757303 14.4006440 22.3616830 82.8974293 195.8299247 42.2820879 145.1508196 141.4528101 242.9248834
EMILIN3 2.7605305 5.9233185 0.3468891 4.0977018 1.1425403 4.1476108 1.4213706 2.6968678 0.1309781 1.6687103 0.3480850 1.3149168 5.3613946 13.0097970 0.9187187 10.0623575 3.7502284 13.0827104
EMILIN2 30.8216346 23.4183023 37.3005036 34.7594362 38.2492888 35.8547740 21.0216378 14.6790747 24.1037093 20.5585291 20.1171615 18.3182698 45.1902545 37.3604532 57.7225502 58.7696911 72.7243797 70.1793800
IGHJ6 24.0852804 12.0754443 43.5126419 7.7044941 47.5239891 2.0111845 23.9940982 11.8003235 43.2798211 7.4573113 47.2655071 1.9469590 24.1768090 12.3569795 43.7467152 7.9598700 47.7838847 2.0775287
CXCL9 33.0631943 9.4561420 38.8811713 7.6134981 7.5026243 4.5814636 32.7185462 8.8920603 38.3109322 6.9723993 7.3925766 4.1958253 33.4114728 10.0560070 39.4598981 8.3135447 7.6143103 5.0025459
IGHV3OR16-12 0.9213345 0.7570405 0.8496743 0.1862667 0.5790641 0.2482409 0.9195141 0.7479924 0.8473013 0.1831259 0.5774469 0.2440552 0.9231585 0.7661980 0.8520539 0.1894614 0.5806858 0.2524984
IGHV1OR21-1 0.2744069 0.4798428 0.8439285 0.4518846 0.1398060 0.0991228 0.1216579 0.2010454 0.4078275 0.1688627 0.0248114 0.0127150 0.6189420 1.1452589 1.7463640 1.2092646 0.7877723 0.7727335
IGHV1-69 17.7194409 12.0978689 47.4225837 7.2742455 21.5765002 0.6854482 17.5585504 11.4485334 46.8142453 6.7197148 21.2956850 0.6339637 17.8818058 12.7840333 48.0388272 7.8745376 21.8610185 0.7411137
IL2RA 3.6212562 1.6531530 3.7956076 1.6907425 1.6896169 1.4831353 2.3730403 0.9347946 2.4166851 0.9251125 0.8220273 0.6775672 5.5260320 2.9235457 5.9613215 3.0900137 3.4728837 3.2464535
IL12A 0.5894995 0.6028262 0.5021917 0.4502679 0.4398678 0.5097215 0.3602347 0.3329725 0.2754279 0.2116138 0.1728890 0.2035655 0.9646757 1.0913796 0.9156534 0.9580714 1.1191206 1.2763261
FGF12 0.4017539 0.3912589 0.3970321 0.3549704 0.3350414 0.5347974 0.2236816 0.1904838 0.2025693 0.1521403 0.1157088 0.2134330 0.7215889 0.8036566 0.7781755 0.8282090 0.9701313 1.3400378
IGHV2-70D 5.2494944 2.2037559 10.6124683 2.3745373 4.5166864 0.2549370 2.3212359 0.8154658 4.2255012 0.7826105 1.1510792 0.0435740 11.8717752 5.9555406 26.6535205 7.2046404 17.7228956 1.4915512
IGHV5-10-1 18.4089045 10.7582697 24.3090353 4.4397139 10.8200471 0.5590466 18.3311529 10.4844850 24.1599741 4.2793595 10.7536925 0.5390150 18.4869860 11.0392039 24.4590162 4.6060770 10.8868110 0.5798227
S100B 14.3112268 13.7117200 20.6585593 11.6229532 13.3081050 32.4548503 8.5434619 7.3137170 9.2631919 5.7090449 5.5905333 12.8277145 23.9728595 25.7066640 46.0722478 23.6629843 31.6795647 82.1126249
IL24 1.4348914 1.1442883 1.7782027 0.7453430 0.6656779 1.2044938 0.9874282 0.7021362 1.1995367 0.4015591 0.3031314 0.6041081 2.0851271 1.8648742 2.6360220 1.3834479 1.4618317 2.4015657
IGHV3OR16-8 3.9101202 2.1501435 5.0074050 1.0520013 3.2916608 0.4648427 2.3751636 1.3050429 3.0420720 0.6389758 1.9992303 0.2818413 6.4370470 3.5425020 8.2424430 1.7320008 5.4196012 0.7666681
FGFR2 2.1037184 1.4329725 0.8896712 1.1653023 0.7023026 0.8789361 1.0017626 0.6883085 0.4377430 0.5050062 0.2392186 0.2972537 4.4178442 2.9832702 1.8081726 2.6889362 2.0618336 2.5988864
VIL1 0.3128734 0.4288354 0.3033234 0.3644623 0.4451425 0.2495236 0.1622754 0.2125331 0.1417720 0.1562315 0.1708399 0.0700006 0.6032326 0.8652761 0.6489651 0.8502308 1.1598690 0.8894493
FGF14 0.9188282 1.8239454 0.8795628 1.6837250 2.3132830 1.6309814 0.5207170 0.9794962 0.4569215 0.8280562 0.9997816 0.6555149 1.6213128 3.3964163 1.6931371 3.4235959 5.3524470 4.0580317
IGHJ3P 20.7589233 12.0561939 40.3226018 5.2135491 6.3549767 1.0805637 20.5194615 11.2400926 39.6667041 4.7214996 6.2515991 0.9784302 21.0011796 12.9315493 40.9893448 5.7568775 6.4600638 1.1933583
FGF2 1.8731183 3.1524727 2.1787122 3.7181892 1.9859282 3.4099108 1.1124911 1.7228338 1.2120379 1.8923511 0.8308215 1.4301651 3.1537981 5.7684521 3.9163682 7.3056903 4.7470013 8.1301745
IL16 26.6494998 24.7194233 23.2754166 26.8818705 19.0069379 31.1158859 18.2314242 15.5723576 15.0647283 15.9550285 9.9944788 15.9859716 38.9544903 39.2393950 35.9611542 45.2919881 36.1463258 60.5654999
IL26 0.4666692 0.4463221 0.7548016 0.4100952 0.4997009 0.2740781 0.2669581 0.2236428 0.4417114 0.1828588 0.1999962 0.0808602 0.8157843 0.8907210 1.2898136 0.9197155 1.2485284 0.9289961
LILRA5 8.4188674 6.8165553 11.4529494 6.5439138 7.9775228 10.1650340 5.5949975 4.1283963 7.2137931 3.7081591 4.0041069 4.9806908 12.6679821 11.2550789 18.1832287 11.5482660 15.8938988 20.7456998
CXCL11 3.2695829 1.5585520 3.3402601 0.9859544 0.9153880 0.8039158 1.7441950 0.6790263 1.7792998 0.4281204 0.3321032 0.2719739 6.1290009 3.5773052 6.2706339 2.2706374 2.5231165 2.3762600
IGHV7-40 1.0162827 0.2572275 0.6617827 0.1916110 0.3667258 0.1136483 1.0061626 0.2419271 0.6525041 0.1755609 0.3615002 0.1042599 1.0265046 0.2734956 0.6711933 0.2091284 0.3720269 0.1238822
PHOSPHO1 4.2927466 7.4651960 4.5829437 3.4921380 4.5814744 4.8246311 2.6552766 4.2229008 2.6511577 1.7891179 2.0469215 2.0844625 6.9400203 13.1968886 7.9223401 6.8162234 10.2543784 11.1669386
EAF2 8.2059389 5.0773027 10.3600725 4.9817820 5.7340498 4.2347873 6.1502435 3.5138806 7.4926749 3.2830419 3.4694770 2.4557095 10.9487426 7.3363342 14.3248044 7.5594989 9.4767386 7.3027461
IGHV5-78 1.3265313 1.0899962 1.7834260 0.5354419 1.6134343 0.5984019 0.6998779 0.4987590 0.9082586 0.2133724 0.5502143 0.1672740 2.5142746 2.3820958 3.5018754 1.3436508 4.7311933 2.1407080
CXCL13 15.0080855 8.2166772 36.3418371 2.9179926 3.8082386 2.7094417 7.8026953 3.7003900 17.2759652 1.1614251 1.2326767 0.8193941 28.8672852 18.2450459 76.4489339 7.3312353 11.7651944 8.9591494
IGHV7-27 0.8304153 0.1232289 0.9101949 0.1801515 0.5978588 0.3343183 0.8259059 0.1193206 0.9034816 0.1721466 0.5933415 0.3196921 0.8349493 0.1272651 0.9169581 0.1885286 0.6024105 0.3496137
PPIL4 10.0409869 9.9930089 9.4466156 11.9068047 9.8748753 12.6853835 6.5850767 5.9827750 5.8256622 6.6871443 4.8525500 6.0718815 15.3105912 16.6912890 15.3181809 21.2006788 20.0952409 26.5023209
IGLV4-3 1.6903446 1.0204181 1.2962755 0.2840156 0.7046333 0.3243377 1.6244732 0.7433349 1.2256389 0.1927811 0.6661908 0.2195193 1.7588870 1.4007862 1.3709831 0.4184273 0.7452940 0.4792059
IL36RN 0.7168029 0.2054575 0.2460664 0.5157539 0.5869357 0.0969078 0.2965164 0.0559107 0.0742251 0.1475546 0.1293239 0.0097421 1.7328093 0.7550042 0.8157445 1.8027365 2.6638038 0.9639695
IL2RG 32.8008833 22.9691044 42.3485840 17.0122898 17.7624423 19.8835661 23.4876538 15.2951383 29.1643407 10.6983969 10.3168550 10.8995308 45.8069570 34.4932977 61.4929918 27.0524646 30.5814475 36.2727726
IGHV1OR15-4 2.9001946 1.8477530 6.6794297 2.0902384 1.5077865 0.2776073 1.2378190 0.6861022 2.5161355 0.6909029 0.3695392 0.0480506 6.7951198 4.9762135 17.7314699 6.3237488 6.1520399 1.6038466
HHIP 2.2252305 1.9086059 1.9088433 4.3771542 3.9204878 3.7647380 1.2355029 0.9252825 0.9661138 2.0086106 1.5034156 1.3809111 4.0078019 3.9369344 3.7714844 9.5386725 10.2235365 10.2636963
PADI4 1.2612286 2.8178576 1.8211490 1.1833769 1.4745861 1.7232541 0.7342018 1.5469981 1.0106272 0.5583956 0.6053390 0.6927001 2.1665671 5.1327288 3.2817081 2.5078650 3.5920434 4.2869991
SAA2 1.7701436 0.8536404 0.9967683 0.2782154 1.5456542 1.1542213 0.6216668 0.3060970 0.3732246 0.0779862 0.4107317 0.2693852 5.0403337 2.3806243 2.6620618 0.9925330 5.8165636 4.9454337
LILRB3 18.0895994 15.5117842 22.3499024 17.2286133 15.6701735 19.9782088 12.4253100 9.8085359 14.5722850 10.2757458 8.3092618 10.3365043 26.3360515 24.5312300 34.2786417 28.8859923 29.5518835 38.6135209

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

Individual Genes

Similarly you can run the script for an individual gene:

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

or to view the lmer fit alone using glmmGene:

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

This has the advantage of increased model flexibility.

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