Skip to contents

Calculates dynamic features for lineages

Source: R/RunDynamicFeatures.R
RunDynamicFeatures.Rd

Calculates dynamic features for lineages

Usage

RunDynamicFeatures(
  srt,
  lineages,
  features = NULL,
  suffix = lineages,
  n_candidates = 1000,
  minfreq = 5,
  family = NULL,
  layer = "counts",
  assay = NULL,
  libsize = NULL,
  fit_method = c("gam", "pretsa"),
  knot = 0,
  max_knot_allowed = 10,
  pretsa_backend = c("cpp", "r"),
  padjust_method = "fdr",
  cores = 1,
  verbose = TRUE,
  seed = 11
)

Arguments

srt

A Seurat object.

lineages

Lineage names for which dynamic features should be calculated.

features

Features to use. If NULL, n_candidates must be provided.

suffix

Suffix to append to the output layer names for each lineage. Default is the lineage names.

n_candidates

A number of candidate features to select when features is NULL.

minfreq

Minimum frequency threshold for candidate features. Features with a frequency less than minfreq will be excluded.

family

A character or character vector specifying the family of distributions to use for the GAM. If family is set to NULL, the appropriate family will be automatically determined based on the data. If length(family) is 1, the same family will be used for all features. Otherwise, family must have the same length as features.

layer

Assay layer to use.

assay

Assay to use. NULL uses the default assay.

libsize

A numeric or numeric vector specifying the library size correction factors for each cell. If NULL, the library size correction factors will be calculated based on the expression matrix. If length(libsize) is 1, the same value will be used for all cells. Otherwise, libsize must have the same length as the number of cells in srt.

fit_method

The method used for fitting features. Either "gam" (generalized additive models) or "pretsa" (Pattern recognition in Temporal and Spatial Analyses).

knot

For fit_method = "pretsa": B-spline knots. 0 or "auto".

max_knot_allowed

For fit_method = "pretsa" when knot = "auto": max knots.

pretsa_backend

PreTSA fitting backend. "cpp" batches model selection and fitting; "r" retains the reference implementation.

padjust_method

The method used for p-value adjustment.

cores

The number of worker processes to use for parallelization. Default is 1.

verbose

Whether to print the message. Default is TRUE.

seed

Optional integer seed. When supplied, every input receives a deterministic independent L'Ecuyer-CMRG random-number stream, making results reproducible across worker counts and scheduling order. The caller's random number state is restored when the call finishes.

Value

Returns the modified Seurat object with the calculated dynamic features stored in the tools slot.

References

Zhuang, H., Ji, Z. PreTSA: computationally efficient modeling of temporal and spatial gene expression patterns. Genome Biol (2026). https://doi.org/10.1186/s13059-026-03994-3

See also

DynamicHeatmap, DynamicPlot, RunDynamicEnrichment

Examples

data(pancreas_sub)
pancreas_sub <- RunStandardWorkflow(pancreas_sub)
#> ℹ [2026-08-30 04:59:41] Start standard processing workflow...
#> ℹ [2026-08-30 04:59:41] Checking a list of <Seurat>...
#> ! [2026-08-30 04:59:42] Data 1/1 of the `srt_list` is "unknown"
#> Warning: Data 1/1 of the `srt_list` is "unknown"
#> ℹ [2026-08-30 04:59:42] Perform `NormalizeData()` with `normalization.method = 'LogNormalize'` on 1/1 of `srt_list`...
#> ℹ [2026-08-30 04:59:42] Perform `FindVariableFeatures()` on 1/1 of `srt_list`...
#> ℹ [2026-08-30 04:59:42] Use the separate HVF from `srt_list`
#> ℹ [2026-08-30 04:59:42] Number of available HVF: 2000
#> ℹ [2026-08-30 04:59:42] Finished check
#> ℹ [2026-08-30 04:59:42] Perform `ScaleData()`
#> ℹ [2026-08-30 04:59:42] Perform pca linear dimension reduction
#> ℹ [2026-08-30 04:59:43] Use stored estimated dimensions 1:23 for Standardpca
#> ℹ [2026-08-30 04:59:43] Perform `Seurat::FindClusters()` with `cluster_algorithm = 'louvain'` and `cluster_resolution = 0.6`
#> ℹ [2026-08-30 04:59:43] Reorder clusters...
#> ℹ [2026-08-30 04:59:43] Skip `log1p()` because `layer = data` is not "counts"
#> ℹ [2026-08-30 04:59:43] Perform umap nonlinear dimension reduction
#> ✔ [2026-08-30 04:59:51] Standard processing workflow completed
pancreas_sub <- RunSlingshot(
  pancreas_sub,
  group.by = "SubCellType",
  reduction = "UMAP"
)
#> Warning: Removed 7 rows containing missing values or values outside the scale range
#> (`geom_path()`).
#> Warning: Removed 7 rows containing missing values or values outside the scale range
#> (`geom_path()`).


pancreas_sub <- RunDynamicFeatures(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  n_candidates = 200,
  fit_method = "gam"
)
#> ℹ [2026-08-30 04:59:53] Start find dynamic features
#> ℹ [2026-08-30 04:59:53] Data type is raw counts
#> ℹ [2026-08-30 04:59:54] Number of candidate features (union): 225
#> ℹ [2026-08-30 04:59:54] Data type is raw counts
#> ℹ [2026-08-30 04:59:54] Calculating dynamic features for "Lineage1"...
#> ℹ [2026-08-30 04:59:54] Using 1 core
#> ⠙ [2026-08-30 04:59:54] Running for Ghrl [1/225]              0% | ETA:  9s
#> ⠹ [2026-08-30 04:59:54] Running for Cck [9/225]              4% | ETA:  6s
#> ⠸ [2026-08-30 04:59:54] Running for Smc4 [115/225] ■■■■■       51% | ETA:  3s
#> ⠼ [2026-08-30 04:59:54] Running for Pam [222/225] ■■■■■■■■■   99% | ETA:  0s
#> ✔ [2026-08-30 04:59:54] Completed 225 tasks in 6.3s
#> 
#> ℹ [2026-08-30 04:59:54] Building results
#> ℹ [2026-08-30 05:00:01] Calculating dynamic features for "Lineage2"...
#> ℹ [2026-08-30 05:00:01] Using 1 core
#> ⠙ [2026-08-30 05:00:01] Running for Tmem171 [99/225] ■■■■        44% | ETA:  4s
#> ⠹ [2026-08-30 05:00:01] Running for H2-Q2 [198/225] ■■■■■■■■    88% | ETA:  1s
#> ✔ [2026-08-30 05:00:01] Completed 225 tasks in 6.6s
#> 
#> ℹ [2026-08-30 05:00:01] Building results
#> ✔ [2026-08-30 05:00:07] Find dynamic features done

names(
  pancreas_sub@tools$DynamicFeatures_Lineage1
)
#> [1] "DynamicFeatures" "raw_matrix"      "fitted_matrix"   "upr_matrix"     
#> [5] "lwr_matrix"      "libsize"         "lineages"        "family"         
head(
  pancreas_sub@tools$DynamicFeatures_Lineage1$DynamicFeatures
)
#>      features exp_ncells      r.sq  dev.expl peaktime valleytime pvalue padjust
#> Ghrl     Ghrl        152 0.3596130 0.7120520 16.14970 8.92943499      0       0
#> Ins1     Ins1        194 0.6676277 0.8370525 22.37234 2.94027875      0       0
#> Ins2     Ins2        128 0.7983418 0.9320678 22.37234 0.03525422      0       0
#> Nnat     Nnat        238 0.7666325 0.8470533 22.37234 5.58327966      0       0
#> Iapp     Iapp        310 0.6968821 0.8604451 22.37234 0.03525422      0       0
#> Pyy       Pyy        421 0.4148165 0.7831469 19.11470 9.09236858      0       0
ht <- DynamicHeatmap(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  cell_annotation = "SubCellType",
  n_split = 3,
  reverse_ht = "Lineage1"
)
#> ℹ [2026-08-30 05:00:07] [1] 176 features from Lineage1,Lineage2 passed the threshold (exp_ncells>[1] 20 & r.sq>[1] 0.2 & dev.expl>[1] 0.2 & padjust<[1] 0.05): 
#> ℹ                       Ghrl,Ins1,Ins2,Nnat,Iapp,Pyy,Lrpprc,Chgb,Cck,Slc38a5...
#> ℹ [2026-08-30 05:00:08] 
#> ℹ                       The size of the heatmap is fixed because certain elements are not scalable.
#> ℹ                       The width and height of the heatmap are determined by the size of the current viewport.
#> ℹ                       If you want to have more control over the size, you can manually set the parameters 'width' and 'height'.
ht$plot


DynamicPlot(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  features = c("Arxes1", "Ncoa2"),
  group.by = "SubCellType",
  compare_lineages = TRUE,
  compare_features = FALSE
)
#> ℹ [2026-08-30 05:00:10] Start find dynamic features
#> ℹ [2026-08-30 05:00:10] Data type is raw counts
#> ℹ [2026-08-30 05:00:11] Number of candidate features (union): 2
#> ℹ [2026-08-30 05:00:11] Data type is raw counts
#> ℹ [2026-08-30 05:00:11] Calculating dynamic features for "Lineage1"...
#> ℹ [2026-08-30 05:00:11] Using 1 core
#> ⠙ [2026-08-30 05:00:11] Running for Arxes1 [1/2] ■■■■■       50% | ETA:  0s
#> ✔ [2026-08-30 05:00:11] Completed 2 tasks in 123ms
#> 
#> ℹ [2026-08-30 05:00:11] Building results
#> ✔ [2026-08-30 05:00:11] Find dynamic features done
#> ℹ [2026-08-30 05:00:11] Start find dynamic features
#> ℹ [2026-08-30 05:00:11] Data type is raw counts
#> ℹ [2026-08-30 05:00:12] Number of candidate features (union): 2
#> ℹ [2026-08-30 05:00:12] Data type is raw counts
#> ℹ [2026-08-30 05:00:12] Calculating dynamic features for "Lineage2"...
#> ℹ [2026-08-30 05:00:12] Using 1 core
#> ℹ [2026-08-30 05:00:12] Building results
#> ✔ [2026-08-30 05:00:12] Find dynamic features done


pancreas_sub <- RunDynamicFeatures(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  n_candidates = 200,
  fit_method = "pretsa"
)
#> ℹ [2026-08-30 05:00:13] Start find dynamic features
#> ℹ [2026-08-30 05:00:13] Data type is raw counts
#> ℹ [2026-08-30 05:00:14] Number of candidate features (union): 225
#> ℹ [2026-08-30 05:00:15] Data type is raw counts
#> ℹ [2026-08-30 05:00:15] Calculating dynamic features for "Lineage1"...
#> ℹ [2026-08-30 05:00:15] Calculating dynamic features for "Lineage2"...
#> ✔ [2026-08-30 05:00:15] Find dynamic features done
head(
  pancreas_sub@tools$DynamicFeatures_Lineage1$DynamicFeatures
)
#>      features exp_ncells      r.sq  dev.expl peaktime  valleytime        pvalue
#> Ghrl     Ghrl        152 0.1768403 0.1768403 16.10426 22.37234415  5.167474e-28
#> Ins1     Ins1        194 0.7009764 0.7009764 22.37234 12.66106258 1.329080e-174
#> Ins2     Ins2        128 0.8568251 0.8568251 22.37234 12.57868745 2.192379e-281
#> Nnat     Nnat        238 0.8100130 0.8100130 22.37234 11.55400320 2.306570e-240
#> Iapp     Iapp        310 0.8724484 0.8724484 22.37234  0.03525422 3.839975e-298
#> Pyy       Pyy        421 0.6850789 0.6850789 22.37234  5.07400185 4.288961e-167
#>            padjust
#> Ghrl  7.549881e-28
#> Ins1 1.246012e-173
#> Ins2 1.644284e-279
#> Nnat 7.413974e-239
#> Iapp 4.319972e-296
#> Pyy  3.216721e-166
ht <- DynamicHeatmap(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  cell_annotation = "SubCellType",
  n_split = 3,
  reverse_ht = "Lineage1"
)
#> ℹ [2026-08-30 05:00:15] [1] 164 features from Lineage1,Lineage2 passed the threshold (exp_ncells>[1] 20 & r.sq>[1] 0.2 & dev.expl>[1] 0.2 & padjust<[1] 0.05): 
#> ℹ                       Ins1,Ins2,Nnat,Iapp,Pyy,Lrpprc,Chgb,Cck,Slc38a5,Npy...
#> ℹ [2026-08-30 05:00:15] 
#> ℹ                       The size of the heatmap is fixed because certain elements are not scalable.
#> ℹ                       The width and height of the heatmap are determined by the size of the current viewport.
#> ℹ                       If you want to have more control over the size, you can manually set the parameters 'width' and 'height'.
ht$plot


DynamicPlot(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  features = c("Arxes1", "Ncoa2"),
  group.by = "SubCellType",
  compare_lineages = TRUE,
  compare_features = FALSE
)
#> ℹ [2026-08-30 05:00:17] Start find dynamic features
#> ℹ [2026-08-30 05:00:17] Data type is raw counts
#> ℹ [2026-08-30 05:00:18] Number of candidate features (union): 2
#> ℹ [2026-08-30 05:00:18] Data type is raw counts
#> ℹ [2026-08-30 05:00:18] Calculating dynamic features for "Lineage1"...
#> ℹ [2026-08-30 05:00:18] Using 1 core
#> ⠙ [2026-08-30 05:00:18] Running for Arxes1 [1/2] ■■■■■       50% | ETA:  0s
#> ✔ [2026-08-30 05:00:18] Completed 2 tasks in 121ms
#> 
#> ℹ [2026-08-30 05:00:18] Building results
#> ✔ [2026-08-30 05:00:18] Find dynamic features done
#> ℹ [2026-08-30 05:00:18] Start find dynamic features
#> ℹ [2026-08-30 05:00:19] Data type is raw counts
#> ℹ [2026-08-30 05:00:19] Number of candidate features (union): 2
#> ℹ [2026-08-30 05:00:19] Data type is raw counts
#> ℹ [2026-08-30 05:00:19] Calculating dynamic features for "Lineage2"...
#> ℹ [2026-08-30 05:00:19] Using 1 core
#> ⠙ [2026-08-30 05:00:19] Running for Arxes1 [1/2] ■■■■■       50% | ETA:  0s
#> ✔ [2026-08-30 05:00:19] Completed 2 tasks in 123ms
#> 
#> ℹ [2026-08-30 05:00:19] Building results
#> ✔ [2026-08-30 05:00:19] Find dynamic features done