# 1 Introduction

Typical analysis using regression models assumes a linear affect of the covariate on the response. Here we consider testing non-linear effects in the case of 1) continuous and 2) ordered categorical variables.

We demonstrate this feature on a lightly modified analysis of PBMCs from 8 individuals stimulated with interferon-β (Kang, et al, 2018, Nature Biotech).

# 2 Standard processing

Here is the code from the main vignette:

library(dreamlet)
library(muscat)
library(ExperimentHub)
library(scater)

eh <- ExperimentHub()
sce <- eh[["EH2259"]]

# only keep singlet cells with sufficient reads
sce <- sce[rowSums(counts(sce) > 0) > 0, ]
sce <- sce[, colData(sce)$multiplets == "singlet"] # compute QC metrics qc <- perCellQCMetrics(sce) # remove cells with few or many detected genes ol <- isOutlier(metric = qc$detected, nmads = 2, log = TRUE)
sce <- sce[, !ol]

# set variable indicating stimulated (stim) or control (ctrl)
sce$StimStatus <- sce$stim

sce$id <- paste0(sce$StimStatus, sce$ind) # Create pseudobulk pb <- aggregateToPseudoBulk(sce, assay = "counts", cluster_id = "cell", sample_id = "id", verbose = FALSE ) # 3 Continuous variable Consider the continuous variable Age. Typical analysis only considers linear effects using a single regression coefficient, but we also want to consider the non-linear effects of age. We can peform a basis expansion using splines instead use 3 coefficients to model the age effect. # Simulate age between 18 and 65 pb$Age <- runif(ncol(pb), 18, 65)

# formula included non-linear effects of Age
# by using a natural spline of degree 3
# This corresponds to using 3 coefficients instead of 1
form <- ~ splines::ns(Age, 3)

# Normalize and apply voom/voomWithDreamWeights
res.proc <- processAssays(pb, form, min.count = 5)

# Differential expression analysis within each assay
res.dl <- dreamlet(res.proc, form)

# The spline has degree 3, so there are 3 coefficients
# estimated for Age effects
coefNames(res.dl)
## [1] "(Intercept)"          "splines::ns(Age, 3)1" "splines::ns(Age, 3)2"
## [4] "splines::ns(Age, 3)3"
# Jointly test effects of the 3 spline components
# The test of the 3 coefficients is performed with an F-statistic
topTable(res.dl, coef = coefNames(res.dl)[2:4], number = 3)
## DataFrame with 3 rows and 9 columns
##               assay          ID splines..ns.Age..3.1 splines..ns.Age..3.2
##         <character> <character>            <numeric>            <numeric>
## 1    Megakaryocytes       PSME2             -0.90314             -3.56179
## 2       CD8 T cells      NDUFS8              2.18670              5.87197
## 3 FCGR3A+ Monocytes       XRCC6              1.67359              5.10846
##   splines..ns.Age..3.3   AveExpr         F     P.Value  adj.P.Val
##              <numeric> <numeric> <numeric>   <numeric>  <numeric>
## 1             -5.22202   9.36461   22.2955 2.03009e-07 0.00294201
## 2              1.35281   6.50647   17.9279 8.42425e-06 0.06104213
## 3              1.86802   6.36049   12.6251 1.49791e-04 0.64665432

# 4 Ordered categorical

We can also test non-linear effects in the case of categorical variables with a natural ordering to the categories. Consider time course data with 4 time points. Each time point is a category and has a natural ordering from first to last.

We have multiple options to model the time course.

• Continuous: Modeling time point as a continuous variable uses a single regression coefficient to model the linear effects of the time course. This is simple, models the order of the time points, but ignores non-linear effects

Model using as.numeric(TimePoint)

• Categorical: Including time point as a typical categorical variable uses estimated the mean response value for each category. So it estimates 4 coefficients. While this can be useful for comparing two categories, it ignores the order of the time points.

Model using factor(TimePoint)

• Ordered categorical: Here, the trend across ordered time points is modled using orthogonal polynomials. The trend is decomposed into independent linear, quadratic, etc., effects that can be tested either jointly or by themselves.

Model using:

ord <- c("time_1", "time_2", "time_3", "time_4")
ordered(factor(TimePoint), ord)

Here we simulated 4 time points, and perform differential expression analysis.

# Consider data generated across 4 time points
# While there are no time points in the real data
# we can add some for demonstration purposes
pb$TimePoint <- ordered(paste0("time_", rep(1:4, 4))) # examine the ordering pb$TimePoint
##  [1] time_1 time_2 time_3 time_4 time_1 time_2 time_3 time_4 time_1 time_2
## [11] time_3 time_4 time_1 time_2 time_3 time_4
## Levels: time_1 < time_2 < time_3 < time_4
# Use formula including time point
form <- ~TimePoint

# Normalize and apply voom/voomWithDreamWeights
res.proc <- processAssays(pb, form, min.count = 5)

# Differential expression analysis within each assay
res.dl <- dreamlet(res.proc, form)

# Examine the coefficient estimated
# for TimePoint it estimates
# linear (i.e. L)
# and cubic (i.e. C) effects
coefNames(res.dl)
## [1] "(Intercept)" "TimePoint.L" "TimePoint.Q" "TimePoint.C"
# Test only linear effect
topTable(res.dl, coef = "TimePoint.L", number = 3)
## DataFrame with 3 rows and 9 columns
##         assay          ID     logFC   AveExpr         t     P.Value adj.P.Val
##   <character> <character> <numeric> <numeric> <numeric>   <numeric> <numeric>
## 1 CD4 T cells        DCXR -0.671645   6.52867  -5.42880 4.78666e-05  0.393058
## 2 CD4 T cells        GGA2 -0.890112   4.98987  -5.26372 6.69370e-05  0.393058
## 3 CD8 T cells        FTH1 -0.759875  14.55380  -4.99199 8.13672e-05  0.393058
##           B     z.std
##   <numeric> <numeric>
## 1  1.686147  -5.42880
## 2  0.936906  -5.26372
## 3  1.627753  -4.99199
# Test linear, quadratic and cubic effcts
coefs <- c("TimePoint.L", "TimePoint.Q", "TimePoint.C")
topTable(res.dl, coef = coefs, number = 3)
## DataFrame with 3 rows and 9 columns
##         assay                   ID TimePoint.L TimePoint.Q TimePoint.C
##   <character>          <character>   <numeric>   <numeric>   <numeric>
## 1 CD8 T cells                 CD52    0.903032   -1.528456   -0.984065
## 2 CD8 T cells CCL5_ENSG00000161570    0.626201   -0.986502   -0.645468
## 3 CD8 T cells                  CD2    0.887578   -1.152415   -0.195616
##   <numeric> <numeric>   <numeric> <numeric>
## 1   8.69458   16.0946 1.90868e-05 0.0853298
## 2  11.88811   16.0433 1.94982e-05 0.0853298
## 3   9.67181   15.7793 2.17777e-05 0.0853298

## 4.1 Sample filtering

Due to variation in cell and read count for each sample, processAssays() filters out some sample. This filtering is summarized here:

details(res.dl)
##               assay n_retain    formula formDropsTerms n_genes n_errors
## 1           B cells       16 ~TimePoint          FALSE    1961        0
## 2   CD14+ Monocytes       16 ~TimePoint          FALSE    3087        0
## 3       CD4 T cells       16 ~TimePoint          FALSE    5262        0
## 4       CD8 T cells       16 ~TimePoint          FALSE    1030        0
## 5   Dendritic cells       13 ~TimePoint          FALSE     164        0
## 6 FCGR3A+ Monocytes       16 ~TimePoint          FALSE    1160        0
## 7    Megakaryocytes       13 ~TimePoint          FALSE     172        0
## 8          NK cells       16 ~TimePoint          FALSE    1656        0
##   error_initial
## 1         FALSE
## 2         FALSE
## 3         FALSE
## 4         FALSE
## 5         FALSE
## 6         FALSE
## 7         FALSE
## 8         FALSE

Whle all 16 samples are detained in B cells, only 9 are retained for megakaryocytes. This can result in a time point being dropped, and so the polynomial expansion for some cell types can have a lower degree. The combined results will then have NA values for these coefficients. For example, for TIMP1 in Megakaryocytes above there is not enought data to fit the cubic term, so TimePoint.C is NA.

# 5 Session Info

## R Under development (unstable) (2024-01-16 r85808)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
##
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
##  [3] LC_TIME=en_GB              LC_COLLATE=C
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: America/New_York
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods
## [8] base
##
## other attached packages:
##  [1] GO.db_3.18.0                org.Hs.eg.db_3.18.0
##  [3] AnnotationDbi_1.65.2        zenith_1.5.2
##  [5] muscData_1.17.0             scater_1.31.2
##  [7] scuttle_1.13.0              ExperimentHub_2.11.1
##  [9] AnnotationHub_3.11.1        BiocFileCache_2.11.1
## [11] dbplyr_2.4.0                muscat_1.17.0
## [13] dreamlet_1.1.18             SingleCellExperiment_1.25.0
## [15] SummarizedExperiment_1.33.3 Biobase_2.63.0
## [17] GenomicRanges_1.55.3        GenomeInfoDb_1.39.6
## [19] IRanges_2.37.1              S4Vectors_0.41.3
## [21] BiocGenerics_0.49.1         MatrixGenerics_1.15.0
## [23] matrixStats_1.2.0           variancePartition_1.33.11
## [25] BiocParallel_1.37.0         limma_3.59.3
## [27] ggplot2_3.4.4               BiocStyle_2.31.0
##
## loaded via a namespace (and not attached):
##   [1] bitops_1.0-7              httr_1.4.7
##   [3] RColorBrewer_1.1-3        doParallel_1.0.17
##   [5] Rgraphviz_2.47.0          numDeriv_2016.8-1.1
##   [7] tools_4.4.0               sctransform_0.4.1
##   [9] backports_1.4.1           utf8_1.2.4
##  [11] R6_2.5.1                  metafor_4.4-0
##  [13] mgcv_1.9-1                GetoptLong_1.0.5
##  [15] withr_3.0.0               prettyunits_1.2.0
##  [17] gridExtra_2.3             cli_3.6.2
##  [19] sandwich_3.1-0            labeling_0.4.3
##  [21] sass_0.4.8                KEGGgraph_1.63.0
##  [23] SQUAREM_2021.1            mvtnorm_1.2-4
##  [25] blme_1.0-5                mixsqp_0.3-54
##  [27] parallelly_1.37.0         invgamma_1.1
##  [29] RSQLite_2.3.5             generics_0.1.3
##  [31] shape_1.4.6               gtools_3.9.5
##  [33] dplyr_1.1.4               Matrix_1.6-5
##  [37] fansi_1.0.6               abind_1.4-5
##  [39] lifecycle_1.0.4           multcomp_1.4-25
##  [41] yaml_2.3.8                edgeR_4.1.16
##  [43] mathjaxr_1.6-0            gplots_3.1.3.1
##  [45] SparseArray_1.3.4         grid_4.4.0
##  [47] blob_1.2.4                crayon_1.5.2
##  [49] lattice_0.22-5            beachmat_2.19.1
##  [51] msigdbr_7.5.1             annotate_1.81.2
##  [53] KEGGREST_1.43.0           magick_2.8.2
##  [55] pillar_1.9.0              knitr_1.45
##  [57] ComplexHeatmap_2.19.0     rjson_0.2.21
##  [59] boot_1.3-28.1             estimability_1.4.1
##  [61] corpcor_1.6.10            future.apply_1.11.1
##  [63] codetools_0.2-19          glue_1.7.0
##  [65] data.table_1.15.0         vctrs_0.6.5
##  [67] png_0.1-8                 Rdpack_2.6
##  [69] gtable_0.3.4              assertthat_0.2.1
##  [71] cachem_1.0.8              xfun_0.42
##  [73] mime_0.12                 rbibutils_2.2.16
##  [75] S4Arrays_1.3.3            Rfast_2.1.0
##  [77] coda_0.19-4.1             survival_3.5-8
##  [79] iterators_1.0.14          statmod_1.5.0
##  [81] TH.data_1.1-2             nlme_3.1-164
##  [83] pbkrtest_0.5.2            bit64_4.0.5
##  [85] filelock_1.0.3            progress_1.2.3
##  [87] EnvStats_2.8.1            bslib_0.6.1
##  [89] TMB_1.9.10                irlba_2.3.5.1
##  [91] vipor_0.4.7               KernSmooth_2.23-22
##  [93] colorspace_2.1-0          rmeta_3.0
##  [95] DBI_1.2.1                 DESeq2_1.43.2
##  [97] tidyselect_1.2.0          emmeans_1.10.0
##  [99] curl_5.2.0                bit_4.0.5
## [101] compiler_4.4.0            graph_1.81.0
## [103] BiocNeighbors_1.21.2      DelayedArray_0.29.2
## [105] bookdown_0.37             scales_1.3.0
## [107] caTools_1.18.2            remaCor_0.0.18
## [109] rappdirs_0.3.3            stringr_1.5.1
## [111] digest_0.6.34             minqa_1.2.6
## [113] rmarkdown_2.25            aod_1.3.3
## [115] XVector_0.43.1            RhpcBLASctl_0.23-42
## [117] htmltools_0.5.7           pkgconfig_2.0.3
## [119] lme4_1.1-35.1             sparseMatrixStats_1.15.0
## [121] highr_0.10                mashr_0.2.79
## [123] fastmap_1.1.1             rlang_1.1.3
## [125] GlobalOptions_0.1.2       DelayedMatrixStats_1.25.1
## [127] farver_2.1.1              jquerylib_0.1.4
## [129] zoo_1.8-12                jsonlite_1.8.8
## [131] BiocSingular_1.19.0       RCurl_1.98-1.14
## [133] magrittr_2.0.3            GenomeInfoDbData_1.2.11
## [135] munsell_0.5.0             Rcpp_1.0.12
## [137] babelgene_22.9            viridis_0.6.5
## [139] EnrichmentBrowser_2.33.1  RcppZiggurat_0.1.6
## [141] stringi_1.8.3             zlibbioc_1.49.0
## [143] MASS_7.3-60.2             plyr_1.8.9
## [145] listenv_0.9.1             parallel_4.4.0
## [147] ggrepel_0.9.5             Biostrings_2.71.2
## [149] splines_4.4.0             hms_1.1.3
## [151] circlize_0.4.15           locfit_1.5-9.8
## [153] reshape2_1.4.4            ScaledMatrix_1.11.0
## [155] BiocVersion_3.19.1        XML_3.99-0.16.1
## [157] evaluate_0.23             RcppParallel_5.1.7
## [159] BiocManager_1.30.22       nloptr_2.0.3
## [161] foreach_1.5.2             tidyr_1.3.1
## [163] purrr_1.0.2               future_1.33.1
## [165] clue_0.3-65               scattermore_1.2
## [167] ashr_2.2-63               rsvd_1.0.5
## [169] broom_1.0.5               xtable_1.8-4
## [171] fANCOVA_0.6-1             viridisLite_0.4.2
## [173] truncnorm_1.0-9           tibble_3.2.1
## [175] lmerTest_3.1-3            glmmTMB_1.1.8
## [177] memoise_2.0.1             beeswarm_0.4.0
## [179] cluster_2.1.6             globals_0.16.2
## [181] GSEABase_1.65.1