Fitting models with seqwrap

Authors

Affiliations

Chidimma Echebiri

University of Inland Norway

Daniel Hammarström

University of Inland Norway

The seqwrap package is designed to iterate over multiple targets in a high-dimensional data set to fit a user-defined regression model using a user-selected model fitting function. Each target may represent, e.g., a gene, an exon, a protein, a methylation site, or any other outcome from a high-throughput omics platform. Each target is measured across a set of samples, with a metadata set used to describe study covariates. In the basic case, data and metadata are combined in each iteration to fit the selected model.

Dungan et al. (Dungan et al. 2022) evaluated a senolytic treatment (dasatinib and quercetin) in mouse skeletal muscle subjected to mechanical overload through synergist ablation. The study included sham surgery and a vehicle control, allowing analysis of the effect of senolytic treatment and the interaction between mechanical overloading and senolytic treatment. RNA-seq was performed on 20 animals (5 per group). The data are available in the seqwrap package under dungan_counts. We will load the data and filter low-expression genes.

# Load packages
library(dplyr)
library(ggplot2)
library(seqwrap)
library(edgeR)
library(gt)

# This vignette is pre-rendered by the package maintainer (not built on CRAN),
# so n_cores can be tuned to the rendering machine.

n_cores <- parallel::detectCores()

# Load data 
metadata <- dungan_counts$metadata
counts <- dungan_counts$countdata

The metadata contains the sample ID corresponding to columns in the counts data set. We are using seq_sample_id as the column name here. The first rows of the prepared metadata are displayed in Table 1.

Table 1: Selected rows from the Dungan et al. metadata data frame.

seq_sample_id	treatment	surgery
OS14d5	Senolytic	Overload
OS14d6	Senolytic	Overload
OSS4	Senolytic	Sham
OSS5	Senolytic	Sham
OV14D5	Vehicle	Overload
OV14D6	Vehicle	Overload
OVS4	Vehicle	Sham
OVS5	Vehicle	Sham

The count data matrix can be manipulated as a data frame, with the first column indicating the targets (in this case, gene symbols). Importantly, additional column names have names corresponding to the seq_sample_id variable in the metadata. As an example, the first rows of the count data frame are shown for the same samples as in the metadata data frame (see Table 2). seqwrap expects that targets are unique so before going further we should validate this in the data. When doing so we find that five targets have duplicate rows in the data set, for simplicity, we will remove them.

Table 2: Selected rows from the Dungan et al. counts data frame.

gene_name	OS14d5	OS14d6	OSS4	OSS5	OV14D5	OV14D6	OVS4	OVS5
Gnai3	858	773	280	295	700	357	380	236
Cdc45	63	71	45	43	61	39	53	54
H19	63437	69156	4539	4694	75553	63084	3565	4981
Narf	1935	2033	2573	2415	2275	1807	2999	2543
Cav2	730	701	491	400	584	441	548	320
Klf6	1300	1232	356	327	820	458	530	518
Scmh1	672	635	738	547	599	464	826	637
Cox5a	5384	6501	7794	6355	7189	5354	8651	8478
Tbx2	91	103	59	55	131	107	59	35
Tbx4	43	77	25	33	65	45	37	27

# Find targets with more than a unique row in the data
rm_targets <- counts |> 
  summarise(.by = gene_name, 
            n = n()) |> 
  filter(n > 1) |> 
  pull(gene_name)

# Filter out non-uniqe rows
counts <- counts |> 
  filter(!(gene_name %in% rm_targets))

The data suggests a negative binomial model were the linear predictor can be used to assess the difference between senolytic treatment and control (treat; \(\beta_1g\)), mechanical overloading compared to sham surgery (mov; \(\beta_2g\)), and their interaction (\(\beta_3g\)) for gene g in 1 to 12273 and row i in 1 to 20 of the meta data (Equation 1). We want to control for differences in sequencing depth by including an offset (s_i). We will use the TMM method (Robinson and Oshlack 2010) to calculate the effective library size.

\[ \begin{align} y_{g[i]} &\sim \operatorname{NB2}(\mu_{g[i]}, \theta_g) \\ \operatorname{log}(\mu_{g[i]}) &= \beta_0 + \beta_{1g}\text{treat}_i + \beta_{2g}\text{mov}_i + \beta_{3g}(\text{treat}_i \times \text{mov}_i) + \operatorname{log}(s_i) \end{align} \tag{1}\]

Below, we are using the edgeR package to calculate the library size normalization factor. Effective library sizes are then expressed as relative to the median effective library and log transformed.

# Use EdgeR to calculate the TMM
y <- edgeR::DGEList(counts[,-1])
y <- edgeR::calcNormFactors(y)

calcNormFactors has been renamed to normLibSizes

# Combine the data into metadata
metadata <- metadata |> 
  inner_join(
    y$samples |> 
  tibble::rownames_to_column("seq_sample_id")
  ) |> 
   mutate(
     efflibsize = (lib.size * norm.factors) / median(lib.size * norm.factors),
     efflibsize = log(efflibsize)
   ) 

Joining with `by = join_by(seq_sample_id)`

seqwrap_compose is used to combine the data, metadata and the selected modelling function with its arguments into a swobject. A negative binomial model can be fitted using e.g., glmmTMB from the glmmTMB package (glmmTMB::glmmTMB) or MASS::glm.nb, here we will use glmmTMB as it can be extended in subsequent modelling steps.

The specified glmmTMB function (modelfun slot) needs a formula and a family function, these are supplied in the arguments slot. The count data are supplied in the data slot and the metadata in the metadata slot. Notice that we use the namespace operator (::) to specify selected functions to avoid difficulties in parallel processing. The samplename indicates the column name in the metadata that gives sample names.

m1.1 <- seqwrap_compose(
  modelfun = glmmTMB::glmmTMB,
  data = counts, 
  metadata = metadata, 
  arguments = list(formula = y ~ treatment * surgery + offset(efflibsize), 
                family = glmmTMB::nbinom2),
  samplename = "seq_sample_id"
  )

Next, we will run the model using seqwrap. We recommend using subset to specify a smaller number of targets when prototyping the model. Here, we will assess whether the model output matches expectations. By specifying return_models = TRUE, all model objects will be returned in the results object, which can be inspected if needed. This is not the default, because retaining every fitted model is the largest memory cost a run can incur; seqwrap warns when more than ten models would be kept, as it does here for our 24 prototyping targets. Here we are using a single core for fitting the models, specified in the cores argument.

m1.1_temp <- seqwrap(
  m1.1, 
  cores = 1, 
  subset = 1:10, 
  return_models = TRUE, 
  verbose = FALSE
)

Following iterative fitting using seqwrap, we can summarize the models using seqwrap_summarise, which produces a list of two summaries. Below, we are saving the summaries in an object for later access. Running the seqwrap_summarise function also gives us an overview of the model and successfully summarized targets.

In the summaries slot, model parameter values are returned (see Table 3). These are produced using broom.mixed::tidy which is a general function to retrieve model parameters from many regression model packages in R.

In the evaluations slot, the convergence diagnostics reported by the fitting algorithm are returned (see Table 4). These are read directly from the fitted model object, so the columns are the same whichever engine was used: converged and code record what the optimiser reported, singular flags a degenerate fit, and message carries any convergence message. The two flags are separate on purpose: glmmTMB regularly returns convergence code 0 together with a non-positive-definite Hessian, so screening a large set of fits usually means requiring converged & !singular rather than converged alone.

Residual-based diagnostics are available through residual_diagnostics(), which simulates scaled residuals using the DHARMa package and can be passed as eval_fun. It is roughly 400 times more expensive per target than reading convergence codes, so it is best applied to a subset of targets.

m1.1_temp_sum <- seqwrap_summarise(m1.1_temp, verbose = FALSE)

Table 3: The first ten rows from summarised preliminary models (model summaries).

target	effect	component	term	estimate	std.error	statistic	p.value
Cav2	fixed	cond	(Intercept)	6.072	0.047	128.043	0
Cav2	fixed	cond	treatmentSenolytic	0.072	0.067	1.079	0.281
Cav2	fixed	cond	surgeryOverload	0.193	0.066	2.939	0.003
Cav2	fixed	cond	treatmentSenolytic:surgeryOverload	0.047	0.092	0.509	0.611
Cdc45	fixed	cond	(Intercept)	3.953	0.068	58.036	0
Cdc45	fixed	cond	treatmentSenolytic	−0.052	0.096	−0.547	0.584
Cdc45	fixed	cond	surgeryOverload	0.112	0.09	1.241	0.215
Cdc45	fixed	cond	treatmentSenolytic:surgeryOverload	0.033	0.127	0.264	0.792
Cox5a	fixed	cond	(Intercept)	9.066	0.034	265.863	0
Cox5a	fixed	cond	treatmentSenolytic	−0.129	0.048	−2.667	0.008

Table 4: The first ten rows from summarised preliminary models (model evaluations). Values are convergence diagnostics read from the fitted model objects.

target	engine	converged	code	message	singular	iterations
Cav2	glmmTMB	TRUE	0	relative convergence (4)	FALSE	31
Cdc45	glmmTMB	TRUE	0	relative convergence (4)	FALSE	30
Cox5a	glmmTMB	TRUE	0	relative convergence (4)	FALSE	43
Gnai3	glmmTMB	TRUE	0	relative convergence (4)	FALSE	27
H19	glmmTMB	TRUE	0	relative convergence (4)	FALSE	31
Klf6	glmmTMB	TRUE	0	relative convergence (4)	FALSE	27
Narf	glmmTMB	TRUE	0	relative convergence (4)	FALSE	40
Scmh1	glmmTMB	TRUE	0	relative convergence (4)	FALSE	34
Tbx2	glmmTMB	TRUE	0	relative convergence (4)	FALSE	29
Tbx4	glmmTMB	TRUE	0	relative convergence (4)	FALSE	21

Specifying summary and evaluation function

We can specify custom functions for evaluating the model object. seqwrap allows for two slots (summaries and evaluations). Each function should take the model object as an argument, and for compatibility with seqwrap_summarise, return a data frame. However, for our next step, we would like to retrieve the dispersion parameter from each model. seqwrap provides dispersion_evaluation() for this purpose: it reads the dispersion parameter on the log scale from a glmmTMB model object, together with its standard error and the log of the mean observed counts. We will place this function in the evaluation slot. Using our preliminary models, we can validate the function for one model.

dispersion_evaluation(m1.1_temp@models[[1]])

# A tibble: 1 × 3
  dispersion dispersion.se log_mu
       <dbl>         <dbl>  <dbl>
1       4.76         0.394   6.22

The output is a data frame with three columns. The dispersion parameter is expected to be on the log scale. To validate this we can directly access the estimated dispersion parameter from the model object using sigma(), which return the dispersion parameter on the natural scale (see ?sigma.glmmTMB for details).

log(sigma(m1.1_temp@models[[1]]))

[1] 4.760462

Running the full model

When we are happy with preliminary performance on a subset of targets, we can initialize seqwrap to iterate over all targets. We are updating the previous model with the new evaluation function directly in seqwrap. Notice that we also remove the subset and set return_models to FALSE. Saving all model objects will be too memory-consuming.

We should expect some errors from modelling many thousands of targets. As seen in the note from seqwrap, some models are associated with warnings from the modelling algorithm and/or errors from the evaluation function.

m1.2_results <- seqwrap(
  m1.1,
  eval_fun = dispersion_evaluation,
  cores = n_cores,
  return_models = FALSE,
  verbose = FALSE
)

Errors and warnings can be inspected through the errors slot, which records one row per condition raised, holding the target, the stage it arose at (fit, summary or evaluation), whether it was an error or a warning, and the message. Targets that completed cleanly do not appear, so the table stays small even for very large runs. seqwrap_errors() filters it. From the inspection below, warnings are at least to some extent related to model convergence.

# The errors slot holds one row per condition raised, with columns
# target, stage, type, class and message. seqwrap_errors() filters it.
warnings <- seqwrap_errors(m1.2_results, stage = "fit", type = "warning")

# Which warnings dominate, and how many targets each affects
warnings |>
  count(message, sort = TRUE) |>
  slice_head(n = 3)

# A tibble: 2 × 2
  message                                                                      n
  <chr>                                                                    <int>
1 Model convergence problem; false convergence (8). See vignette('trouble…    58
2 Model convergence problem; non-positive-definite Hessian matrix. See vi…    53

Using target data for prior information in glmmTMB

A major benefit of canonical R packages for RNA-seq data modelling is the sharing of information between targets through shrinkage and empirical Bayes strategies. We suggest using priors in glmmTMB as convenient way of flexibly regularize models using an empirical Bayes strategy. In glmmTMB, the priors argument can be used to set priors on parameters. We will use the estimated parameters from the first model to set up the priors.

First, we need to extract the parameters. From the evaluation function we get the average log count and the dispersion parameter estimate. From the summaries, we will get the other parameters (summarized in Figure 1).

m1.2_sum <- seqwrap_summarise(m1.2_results, verbose = FALSE)

Figure 1: Estimates from modelling the Dungan data. The dispersion trend is visualized using the default ggplot2::geom_smooth`.

Except for the dispersion values, we can easily summarize parameter values by their mean and standard deviation. A more elaborate summary of the dispersion values could include an estimate of the mean (log) dispersion, conditional on the average log counts. As the dispersion prior will be slightly different for each target (based on the observed counts), priors have to be combined into target-specific data frames and stored in a list. seqwrap will make target-specific list elements available in each iteration, making it possible to build the prior data frame for glmmTMB.

seqwrap_priors() does this work. Given the results from the first run and the count data, it summarises the fixed-effect estimates across targets into normal priors, centered on zero by default so that estimates are shrunk towards no effect, fits a loess trend of the log dispersion against the log mean count, and gives every target a normal prior for the dispersion centered on the trend at its own mean count. Random-effect standard deviations, when the model has any, receive a gamma prior. The result is a list with one prior data frame per target, in the order of the rows in counts, which prints a summary of the priors it holds.

Priors_list <- seqwrap_priors(m1.2_results, data = counts)

# The priors built for the first target
Priors_list[[1]]

                 prior      class                               coef
1     normal(0, 0.137)      fixef                 treatmentSenolytic
2     normal(0, 0.587)      fixef                    surgeryOverload
3     normal(0, 0.164)      fixef treatmentSenolytic:surgeryOverload
4 normal(5.506, 2.296) fixef_disp                                  1

Plotting the object draws each prior over the estimates it summarises, which is a quick check that the priors are sensible. The same panels can be requested while building the priors with plot = TRUE.

plot(Priors_list)

The location and scale used for each shared prior, the per-target dispersion priors, and the fitted loess trend are kept as attributes (common, dispersion and trend) for inspection, and the per-target estimates the panels are drawn from are kept in the data attribute for displays of your own (attr(Priors_list, "data")). Priors can be centered on the mean estimate across targets instead (center = "mean"), restricted to selected terms (terms), made robust to targets whose fits went astray (robust = TRUE), or built without a dispersion trend (trend = "constant"). Importantly, we urge the user to evaluate these choices through simulations to guide informed decisions.

The list is named by target identifier, so seqwrap matches each element to its target by name. A list built by hand can be used in the same way: name its elements by target identifier, or leave it unnamed to match by position, in which case seqwrap warns that element i is taken to belong to row i of the data.

Again, we will combine data, metadata, and arguments using seqwrap_compose. We are now also adding the priors to the arguments. The variables used in the priors data frame are made available from the target-specific list referenced in targetdata = Prior_list. Using arguments = alist(... avoids evaluation of the list so that the data frame is not evaluated before being used in each iteration. Priors have been built corresponding what is expected by glmmTMB (see Priors in glmmTMB).

m1.3 <- seqwrap_compose(
  modelfun = glmmTMB::glmmTMB,
  data = counts, 
  metadata = metadata, 
  eval_fun = dispersion_evaluation, 
  arguments =  alist(
    formula = y ~ treatment * surgery + offset(efflibsize), 
    family = glmmTMB::nbinom2,
      priors = data.frame(
        prior = prior,
        class = class,
        coef = coef)), 
  targetdata = Priors_list
)


m1.3_temp <- seqwrap(
  m1.3, 
  subset = 1:10, 
  return_models = TRUE, 
  verbose = FALSE
)


summary(m1.3_temp@models[[1]])

 Family: nbinom2  ( log )
Formula:          y ~ treatment * surgery + offset(efflibsize)
Data: list(y = c(730L, 701L, 617L, 541L, 697L, 491L, 400L, 374L, 442L,  
460L, 584L, 441L, 548L, 635L, 526L, 548L, 320L, 331L, 318L, 370L 
), treatment = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,  
2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), levels = c("Vehicle",  
"Senolytic"), class = "factor"), surgery = structure(c(2L, 2L,  
2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L,  
1L, 1L), levels = c("Sham", "Overload"), class = "factor"), efflibsize = c(0.166155699457707,  
0.196470816938507, 0.0682347925124099, 0.0387654466787985, 0.0392563629230484,  
-0.0785846626255147, -0.225337792488665, -0.000892800687716758,  
0.000892004305606193, -0.0490532038569703, 0.144822834668174,  
-0.0984341203861923, -0.0169305794837407, 0.049199853171865,  
0.100216656584421, 0.0604222716589449, -0.167592705247281, -0.28188788975239,  
-0.182898895921145, -0.191409863872739))

      AIC       BIC    logLik -2*log(L)  df.resid 
    224.5     229.5    -107.3     214.5        15 


Dispersion parameter for nbinom2 family ():  119 

Conditional model:
                                   Estimate Std. Error z value Pr(>|z|)    
(Intercept)                         6.07647    0.04363  139.29  < 2e-16 ***
treatmentSenolytic                  0.06364    0.05630    1.13  0.25828    
surgeryOverload                     0.19124    0.05990    3.19  0.00141 ** 
treatmentSenolytic:surgeryOverload  0.04895    0.07667    0.64  0.52319    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Priors:
fixef(treatmentSenolytic) ~ normal(0, 0.137) 
fixef(surgeryOverload) ~ normal(0, 0.587) 
fixef(treatmentSenolytic:surgeryOverload) ~ normal(0, 0.164) 
fixef_disp(1) ~ normal(5.499, 2.296) 

Running the prototyping seems to work, we will proceed with the full data set.

m1.3_results <- seqwrap(
  m1.3,
  cores = n_cores,
  return_models = FALSE,
  verbose = FALSE
)

m1.3_sum <- seqwrap_summarise(m1.3_results, verbose = FALSE)

Regularization increases the number of models without convergence issues: 0 of 12273 targets now raise a fitting warning, compared with 111 for the unregularized model.

The summaries data frame is easily manipulated in subsequent analysis steps to, e.g., adjust p-values and create basic visualizations of the results (Figure 2).

Figure 2: A visualization of results from the regularized model of the Dungan et al. data.

The above example shows how to build models with seqwrap. Importantly, we make no claims about the validity of the final model but urge the user to assess assumptions in their data and models as part of a principled workflow.

References

Dungan, Cory M., Vandre C. Figueiredo, Yuan Wen, et al. 2022. “Senolytic Treatment Rescues Blunted Muscle Hypertrophy in Old Mice.” GeroScience 44 (4): 1925–40. https://doi.org/10.1007/s11357-022-00542-2.

Robinson, Mark D., and Alicia Oshlack. 2010. “A Scaling Normalization Method for Differential Expression Analysis of RNA-seq Data.” Genome Biology 11 (3): R25. https://doi.org/10.1186/gb-2010-11-3-r25.