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).
Here is the code from the main vignette:
library(dreamlet)
library(muscat)
library(ExperimentHub)
library(scater)
# Download data, specifying EH2259 for the Kang, et al study
eh <- ExperimentHub()
sce <- eh[["EH2259"]]
# only keep singlet cells with sufficient reads
sce <- sce[rowSums2(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
)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 CD4 T cells THEMIS 1.16752 6.92840
## 2 CD4 T cells PTPLB 0.68874 7.12379
## 3 CD14+ Monocytes SLC36A4 1.31412 5.39654
## splines..ns.Age..3.3 AveExpr F P.Value adj.P.Val
## <numeric> <numeric> <numeric> <numeric> <numeric>
## 1 1.111439 3.91103 18.3347 1.47128e-05 0.213217
## 2 -0.149684 3.77629 11.9246 1.94126e-04 0.999985
## 3 0.282197 4.21266 12.0068 2.86562e-04 0.999985We 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)
# quadratic (i.e. Q)
# 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
## AveExpr F P.Value adj.P.Val
## <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.0853298Due 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 FALSEWhle 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.
## R version 4.3.2 Patched (2023-11-13 r85521)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.18-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
## [9] LC_ADDRESS=C LC_TELEPHONE=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.64.1 zenith_1.4.2
## [5] muscData_1.16.0 scater_1.30.1
## [7] scuttle_1.12.0 ExperimentHub_2.10.0
## [9] AnnotationHub_3.10.0 BiocFileCache_2.10.1
## [11] dbplyr_2.4.0 muscat_1.16.0
## [13] dreamlet_1.0.3 SingleCellExperiment_1.24.0
## [15] SummarizedExperiment_1.32.0 Biobase_2.62.0
## [17] GenomicRanges_1.54.1 GenomeInfoDb_1.38.6
## [19] IRanges_2.36.0 S4Vectors_0.40.2
## [21] BiocGenerics_0.48.1 MatrixGenerics_1.14.0
## [23] matrixStats_1.2.0 variancePartition_1.32.5
## [25] BiocParallel_1.36.0 limma_3.58.1
## [27] ggplot2_3.5.0 BiocStyle_2.30.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.46.0 numDeriv_2016.8-1.1
## [7] tools_4.3.2 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.62.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.1 gtools_3.9.5
## [33] dplyr_1.1.4 Matrix_1.6-5
## [35] metadat_1.2-0 ggbeeswarm_0.7.2
## [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.0.16
## [43] mathjaxr_1.6-0 gplots_3.1.3.1
## [45] SparseArray_1.2.4 grid_4.3.2
## [47] blob_1.2.4 promises_1.2.1
## [49] crayon_1.5.2 lattice_0.22-5
## [51] beachmat_2.18.1 msigdbr_7.5.1
## [53] annotate_1.80.0 KEGGREST_1.42.0
## [55] magick_2.8.3 pillar_1.9.0
## [57] knitr_1.45 ComplexHeatmap_2.18.0
## [59] rjson_0.2.21 boot_1.3-30
## [61] estimability_1.5 corpcor_1.6.10
## [63] future.apply_1.11.1 codetools_0.2-19
## [65] glue_1.7.0 data.table_1.15.0
## [67] vctrs_0.6.5 png_0.1-8
## [69] Rdpack_2.6 gtable_0.3.4
## [71] assertthat_0.2.1 cachem_1.0.8
## [73] xfun_0.42 mime_0.12
## [75] rbibutils_2.2.16 S4Arrays_1.2.0
## [77] Rfast_2.1.0 coda_0.19-4.1
## [79] survival_3.5-8 iterators_1.0.14
## [81] statmod_1.5.0 ellipsis_0.3.2
## [83] interactiveDisplayBase_1.40.0 TH.data_1.1-2
## [85] nlme_3.1-164 pbkrtest_0.5.2
## [87] bit64_4.0.5 filelock_1.0.3
## [89] progress_1.2.3 EnvStats_2.8.1
## [91] bslib_0.6.1 TMB_1.9.10
## [93] irlba_2.3.5.1 vipor_0.4.7
## [95] KernSmooth_2.23-22 colorspace_2.1-0
## [97] rmeta_3.0 DBI_1.2.2
## [99] DESeq2_1.42.0 tidyselect_1.2.0
## [101] emmeans_1.10.0 curl_5.2.0
## [103] bit_4.0.5 compiler_4.3.2
## [105] graph_1.80.0 BiocNeighbors_1.20.2
## [107] DelayedArray_0.28.0 bookdown_0.37
## [109] scales_1.3.0 caTools_1.18.2
## [111] remaCor_0.0.18 rappdirs_0.3.3
## [113] stringr_1.5.1 digest_0.6.34
## [115] minqa_1.2.6 rmarkdown_2.25
## [117] aod_1.3.3 XVector_0.42.0
## [119] RhpcBLASctl_0.23-42 htmltools_0.5.7
## [121] pkgconfig_2.0.3 lme4_1.1-35.1
## [123] sparseMatrixStats_1.14.0 highr_0.10
## [125] mashr_0.2.79 fastmap_1.1.1
## [127] rlang_1.1.3 GlobalOptions_0.1.2
## [129] shiny_1.8.0 DelayedMatrixStats_1.24.0
## [131] farver_2.1.1 jquerylib_0.1.4
## [133] zoo_1.8-12 jsonlite_1.8.8
## [135] BiocSingular_1.18.0 RCurl_1.98-1.14
## [137] magrittr_2.0.3 GenomeInfoDbData_1.2.11
## [139] munsell_0.5.0 Rcpp_1.0.12
## [141] babelgene_22.9 viridis_0.6.5
## [143] EnrichmentBrowser_2.32.0 RcppZiggurat_0.1.6
## [145] stringi_1.8.3 zlibbioc_1.48.0
## [147] MASS_7.3-60.0.1 plyr_1.8.9
## [149] listenv_0.9.1 parallel_4.3.2
## [151] ggrepel_0.9.5 Biostrings_2.70.2
## [153] splines_4.3.2 hms_1.1.3
## [155] circlize_0.4.16 locfit_1.5-9.8
## [157] reshape2_1.4.4 ScaledMatrix_1.10.0
## [159] BiocVersion_3.18.1 XML_3.99-0.16.1
## [161] evaluate_0.23 RcppParallel_5.1.7
## [163] BiocManager_1.30.22 httpuv_1.6.14
## [165] nloptr_2.0.3 foreach_1.5.2
## [167] tidyr_1.3.1 purrr_1.0.2
## [169] future_1.33.1 clue_0.3-65
## [171] scattermore_1.2 ashr_2.2-63
## [173] rsvd_1.0.5 broom_1.0.5
## [175] xtable_1.8-4 fANCOVA_0.6-1
## [177] later_1.3.2 viridisLite_0.4.2
## [179] truncnorm_1.0-9 tibble_3.2.1
## [181] lmerTest_3.1-3 glmmTMB_1.1.8
## [183] memoise_2.0.1 beeswarm_0.4.0
## [185] cluster_2.1.6 globals_0.16.2
## [187] GSEABase_1.64.0