Skip to contents

Calculates dynamic features for lineages in a single-cell RNA-seq dataset.

Usage

RunDynamicFeatures(
  srt,
  lineages,
  features = NULL,
  suffix = lineages,
  n_candidates = 1000,
  minfreq = 5,
  family = NULL,
  layer = "counts",
  assay = NULL,
  libsize = NULL,
  cores = 1,
  verbose = TRUE,
  seed = 11
)

Arguments

srt

A Seurat object.

lineages

A character vector specifying the lineage names for which dynamic features should be calculated.

features

A character vector specifying the features (genes or metadata variables) for which dynamic features should be calculated. If NULL, n_candidates must be provided.

suffix

A character vector specifying the suffix to append to the output layer names for each lineage. Defaults to the lineage names.

n_candidates

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

minfreq

An integer specifying the minimum frequency threshold for candidate features. Features with a frequency less than minfreq will be excluded. Defaults to 5.

family

A character or character vector specifying the family of distributions to use for the generalized additive models (GAMs). 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

A character vector specifying the layer in the Seurat object to use. Default is "counts".

assay

A character vector specifying the assay in the Seurat object to use. Default is NULL.

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. Defaults to NULL.

cores

The number of cores to use for parallelization with foreach::foreach. Default is 1.

verbose

Whether to print progress messages. Default is TRUE.

seed

An integer specifying the seed for random number generation. Defaults to 11.

Value

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

Examples

data(pancreas_sub)
pancreas_sub <- standard_scop(pancreas_sub)
#>  [2025-09-20 13:40:38] Start standard scop workflow...
#>  [2025-09-20 13:40:39] Checking a list of <Seurat> object...
#> ! [2025-09-20 13:40:39] Data 1/1 of the `srt_list` is "unknown"
#>  [2025-09-20 13:40:39] Perform `NormalizeData()` with `normalization.method = 'LogNormalize'` on the data 1/1 of the `srt_list`...
#>  [2025-09-20 13:40:41] Perform `Seurat::FindVariableFeatures()` on the data 1/1 of the `srt_list`...
#>  [2025-09-20 13:40:42] Use the separate HVF from srt_list
#>  [2025-09-20 13:40:42] Number of available HVF: 2000
#>  [2025-09-20 13:40:42] Finished check
#>  [2025-09-20 13:40:42] Perform `Seurat::ScaleData()`
#> Warning: Different features in new layer data than already exists for scale.data
#>  [2025-09-20 13:40:42] Perform pca linear dimension reduction
#> StandardPC_ 1 
#> Positive:  Aplp1, Cpe, Gnas, Fam183b, Map1b, Hmgn3, Pcsk1n, Chga, Tuba1a, Bex2 
#> 	   Syt13, Isl1, 1700086L19Rik, Pax6, Chgb, Scgn, Rbp4, Scg3, Gch1, Camk2n1 
#> 	   Cryba2, Pcsk2, Pyy, Tspan7, Mafb, Hist3h2ba, Dbpht2, Abcc8, Rap1b, Slc38a5 
#> Negative:  Spp1, Anxa2, Sparc, Dbi, 1700011H14Rik, Wfdc2, Gsta3, Adamts1, Clu, Mgst1 
#> 	   Bicc1, Ldha, Vim, Cldn3, Cyr61, Rps2, Mt1, Ptn, Phgdh, Nudt19 
#> 	   Smtnl2, Smco4, Habp2, Mt2, Col18a1, Rpl12, Galk1, Cldn10, Acot1, Ccnd1 
#> StandardPC_ 2 
#> Positive:  Rbp4, Tagln2, Tuba1b, Fkbp2, Pyy, Pcsk2, Iapp, Tmem27, Meis2, Tubb4b 
#> 	   Pcsk1n, Dbpht2, Rap1b, Dynll1, Tubb2a, Sdf2l1, Scgn, 1700086L19Rik, Scg2, Abcc8 
#> 	   Atp1b1, Hspa5, Fam183b, Papss2, Slc38a5, Scg3, Mageh1, Tspan7, Ppp1r1a, Ociad2 
#> Negative:  Neurog3, Btbd17, Gadd45a, Ppp1r14a, Neurod2, Sox4, Smarcd2, Mdk, Pax4, Btg2 
#> 	   Sult2b1, Hes6, Grasp, Igfbpl1, Gpx2, Cbfa2t3, Foxa3, Shf, Mfng, Tmsb4x 
#> 	   Amotl2, Gdpd1, Cdc14b, Epb42, Rcor2, Cotl1, Upk3bl, Rbfox3, Cldn6, Cer1 
#> StandardPC_ 3 
#> Positive:  Nusap1, Top2a, Birc5, Aurkb, Cdca8, Pbk, Mki67, Tpx2, Plk1, Ccnb1 
#> 	   2810417H13Rik, Incenp, Cenpf, Ccna2, Prc1, Racgap1, Cdk1, Aurka, Cdca3, Hmmr 
#> 	   Spc24, Kif23, Sgol1, Cenpe, Cdc20, Hist1h1b, Cdca2, Mxd3, Kif22, Ska1 
#> Negative:  Anxa5, Pdzk1ip1, Acot1, Tpm1, Anxa2, Dcdc2a, Capg, Sparc, Ttr, Pamr1 
#> 	   Clu, Cxcl12, Ndrg2, Hnf1aos1, Gas6, Gsta3, Krt18, Ces1d, Atp1b1, Muc1 
#> 	   Hhex, Acadm, Spp1, Enpp2, Bcl2l14, Sat1, Smtnl2, 1700011H14Rik, Tgm2, Fam159a 
#> StandardPC_ 4 
#> Positive:  Glud1, Tm4sf4, Akr1c19, Cldn4, Runx1t1, Fev, Pou3f4, Gm43861, Pgrmc1, Arx 
#> 	   Cd200, Lrpprc, Hmgn3, Ppp1r14c, Pam, Etv1, Tsc22d1, Slc25a5, Akap17b, Pgf 
#> 	   Fam43a, Emb, Jun, Krt8, Dnajc12, Mid1ip1, Ids, Rgs17, Uchl1, Alcam 
#> Negative:  Ins2, Ins1, Ppp1r1a, Nnat, Calr, Sytl4, Sdf2l1, Iapp, Pdia6, Mapt 
#> 	   G6pc2, C2cd4b, Npy, Gng12, P2ry1, Ero1lb, Adra2a, Papss2, Arhgap36, Fam151a 
#> 	   Dlk1, Creld2, Gip, Tmem215, Gm27033, Cntfr, Prss53, C2cd4a, Lyve1, Ociad2 
#> StandardPC_ 5 
#> Positive:  Pdx1, Nkx6-1, Npepl1, Cldn4, Cryba2, Fev, Jun, Chgb, Gng12, Adra2a 
#> 	   Mnx1, Sytl4, Pdk3, Gm27033, Nnat, Chga, Ins2, 1110012L19Rik, Enho, Krt7 
#> 	   Mlxipl, Tmsb10, Flrt1, Pax4, Tubb3, Prrg2, Gars, Frzb, BC023829, Gm2694 
#> Negative:  Irx2, Irx1, Gcg, Ctxn2, Tmem27, Ctsz, Tmsb15l, Nap1l5, Pou6f2, Gria2 
#> 	   Ghrl, Peg10, Smarca1, Arx, Lrpap1, Rgs4, Ttr, Gast, Tmsb15b2, Serpina1b 
#> 	   Slc16a10, Wnk3, Ly6e, Auts2, Sct, Arg1, Dusp10, Sphkap, Dock11, Edn3 
#>  [2025-09-20 13:40:43] Perform `Seurat::FindClusters()` with louvain and `cluster_resolution` = 0.6
#>  [2025-09-20 13:40:43] Reorder clusters...
#> ! [2025-09-20 13:40:43] Using `Seurat::AggregateExpression()` to calculate pseudo-bulk data for <Assay5>
#>  [2025-09-20 13:40:44] Perform umap nonlinear dimension reduction
#>  [2025-09-20 13:40:44] Non-linear dimensionality reduction (umap) using (Standardpca) dims (1-50) as input
#>  [2025-09-20 13:40:44] UMAP will return its model
#>  [2025-09-20 13:40:48] Non-linear dimensionality reduction (umap) using (Standardpca) dims (1-50) as input
#>  [2025-09-20 13:40:48] UMAP will return its model
#>  [2025-09-20 13:40:52] Run scop standard workflow done
pancreas_sub <- RunSlingshot(
  pancreas_sub,
  group.by = "SubCellType",
  reduction = "UMAP"
)
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> Warning: Removed 8 rows containing missing values or values outside the scale range
#> (`geom_path()`).
#> Warning: Removed 8 rows containing missing values or values outside the scale range
#> (`geom_path()`).


pancreas_sub <- RunDynamicFeatures(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  n_candidates = 200
)
#>  [2025-09-20 13:40:55] Start find dynamic features
#>  [2025-09-20 13:40:55] Installing: mgcv...
#>  
#>   
#>  No downloads are needed
#>  1 pkg + 3 deps: kept 3 [642ms]
#> ! [2025-09-20 13:40:55] Failed to install: mgcv. Please install manually
#>  [2025-09-20 13:40:56] Data type is raw counts
#>  [2025-09-20 13:40:58] Number of candidate features (union): 225
#>  [2025-09-20 13:40:58] Data type is raw counts
#>  [2025-09-20 13:40:58] Calculating dynamic features for "Lineage1"...
#>  [2025-09-20 13:40:58] Using 1 core
#> ⠙ [2025-09-20 13:40:58] Running [1/225] ETA:  8s
#> ⠹ [2025-09-20 13:40:58] Running [60/225] ETA:  6s
#> ⠸ [2025-09-20 13:40:58] Running [153/225] ETA:  2s
#>  [2025-09-20 13:40:58] Completed 225 tasks in 7.8s
#> 
#>  [2025-09-20 13:41:06] Building results
#>  [2025-09-20 13:41:06] Calculating dynamic features for "Lineage2"...
#>  [2025-09-20 13:41:06] Using 1 core
#> ⠙ [2025-09-20 13:41:06] Running [9/225] ETA:  6s
#> ⠹ [2025-09-20 13:41:06] Running [103/225] ETA:  4s
#> ⠸ [2025-09-20 13:41:06] Running [193/225] ETA:  1s
#>  [2025-09-20 13:41:06] Completed 225 tasks in 7.6s
#> 
#>  [2025-09-20 13:41:14] Building results
#>  [2025-09-20 13:41:14] 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
#> Gcg       Gcg        178 0.07474585 0.8463493 29.93042   2.028632      0
#> Iapp     Iapp        308 0.26805682 0.7766207 29.34657   1.834261      0
#> Pyy       Pyy        475 0.49836581 0.7682591 27.10419  14.503163      0
#> Gast     Gast        103 0.04420017 0.6562081 27.70426   0.140757      0
#> Rbp4     Rbp4        437 0.39150230 0.7319113 26.71387  14.634662      0
#> Chgb     Chgb        302 0.48290058 0.8015278 23.69410   0.140757      0
#>      padjust
#> Gcg        0
#> Iapp       0
#> Pyy        0
#> Gast       0
#> Rbp4       0
#> Chgb       0
ht <- DynamicHeatmap(
  pancreas_sub,
  lineages = c("Lineage1", "Lineage2"),
  cell_annotation = "SubCellType",
  n_split = 6,
  reverse_ht = "Lineage1"
)
#>  [2025-09-20 13:41:14] 175 features from Lineage1,Lineage2 passed the threshold (exp_ncells>20 & r.sq>0.2 & dev.expl>0.2 & padjust<0.05): 
#>                        Iapp,Pyy,Rbp4,Chgb,Slc38a5,Lrpprc,Cck,Chga,2810417H13Rik,Cdc20...
#>  [2025-09-20 13:41:14] e1071 installed successfully
#> 'magick' package is suggested to install to give better rasterization.
#> 
#> Set `ht_opt$message = FALSE` to turn off this message.
#>  [2025-09-20 13:41: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
)
#>  [2025-09-20 13:41:18] MatrixGenerics installed successfully
#>  [2025-09-20 13:41:18] Start find dynamic features
#>  [2025-09-20 13:41:18] Installing: mgcv...
#>  
#>   
#>  No downloads are needed
#>  1 pkg + 3 deps: kept 3 [648ms]
#> ! [2025-09-20 13:41:19] Failed to install: mgcv. Please install manually
#>  [2025-09-20 13:41:19] Data type is raw counts
#>  [2025-09-20 13:41:20] Number of candidate features (union): 2
#>  [2025-09-20 13:41:20] Data type is raw counts
#>  [2025-09-20 13:41:20] Calculating dynamic features for "Lineage1"...
#>  [2025-09-20 13:41:20] Using 1 core
#> ⠙ [2025-09-20 13:41:20] Running [1/2] ETA:  0s
#>  [2025-09-20 13:41:20] Completed 2 tasks in 136ms
#> 
#>  [2025-09-20 13:41:20] Building results
#>  [2025-09-20 13:41:21] Find dynamic features done
#>  [2025-09-20 13:41:21] Start find dynamic features
#>  [2025-09-20 13:41:21] Installing: mgcv...
#>  
#>   
#>  No downloads are needed
#>  1 pkg + 3 deps: kept 3 [637ms]
#> ! [2025-09-20 13:41:21] Failed to install: mgcv. Please install manually
#>  [2025-09-20 13:41:22] Data type is raw counts
#>  [2025-09-20 13:41:22] Number of candidate features (union): 2
#>  [2025-09-20 13:41:23] Data type is raw counts
#>  [2025-09-20 13:41:23] Calculating dynamic features for "Lineage2"...
#>  [2025-09-20 13:41:23] Using 1 core
#> ⠙ [2025-09-20 13:41:23] Running [1/2] ETA:  0s
#>  [2025-09-20 13:41:23] Completed 2 tasks in 146ms
#> 
#>  [2025-09-20 13:41:23] Building results
#>  [2025-09-20 13:41:23] Find dynamic features done
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.
#> Warning: No shared levels found between `names(values)` of the manual scale and the
#> data's fill values.