# 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 (si). 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 `ggplot`2::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.