Introduction

A trajectory inference tool orders the cells of a single-cell RNA-seq experiment along a pseudotime and splits them into branching paths, such as the paths from a progenitor to two cell types. scMaSigPro finds the genes whose expression changes along the pseudotime of a path and the genes whose expression differs between paths (Srivastava et al. 2024). It brings the method of the maSigPro package for time-course experiments (Conesa et al. 2006; Nueda et al. 2014) to single cells in two steps:

The cells of each path are pooled into bins along pseudotime, and their counts are added up into pseudo-bulk counts. The bins play the part of the time points of a time-course experiment.
A polynomial generalised linear model in the binned pseudotime, with one set of terms per path, is fitted to the pseudo-bulk counts of each gene, by default with a negative binomial distribution and an offset for the sequencing depth of each bin. The genes whose model is better than a flat line, and the terms that make it so, tell how a gene changes along each path.

This vignette runs these steps on splat.sim, a small data set bundled with the package. The other vignettes go further:

Trajectory input: pseudotime and path labels from slingshot, Monocle3 or another tool.
Shared root: trajectories whose paths share their first cells.
Model and tuning: the model, its settings and how they change the results.
maSigPro and the ScMaSigPro class: the object, its accessors and the equivalence with maSigPro.
Frequently asked questions.

Installation

scMaSigPro is installed from GitHub with BiocManager, which also installs its Bioconductor dependencies:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("BioBam/scMaSigPro")

The data

splat.sim is a SingleCellExperiment with the counts of 100 genes in 200 cells, simulated with splatter (Zappia et al. 2017) along two paths that start from the same cells. Each cell has a simulated step along its path, which serves as its pseudotime, and the label of its path.

library(scMaSigPro)
library(SingleCellExperiment)
library(ggplot2)

data("splat.sim", package = "scMaSigPro")
splat.sim
#> class: SingleCellExperiment 
#> dim: 100 200 
#> metadata(0):
#> assays(1): counts
#> rownames(100): Gene1 Gene2 ... Gene99 Gene100
#> rowData names(10): Gene BaseGeneMean ... SigmaFacPath1 SigmaFacPath2
#> colnames(200): Cell1 Cell2 ... Cell199 Cell200
#> colData names(6): Cell Batch ... Step sizeFactor
#> reducedDimNames(1): PCA
#> mainExpName: NULL
#> altExpNames(0):
head(colData(splat.sim)[, c("Step", "Group")])
#> DataFrame with 6 rows and 2 columns
#>            Step    Group
#>       <numeric> <factor>
#> Cell1        42    Path2
#> Cell2        41    Path1
#> Cell3        23    Path1
#> Cell4        42    Path1
#> Cell5        79    Path2
#> Cell6         1    Path1

The first two principal components show the two paths, Path1 and Path2, leaving the same starting cells:

pca <- as.data.frame(reducedDim(splat.sim, "PCA"))
pca$Step <- colData(splat.sim)$Step
pca$Path <- colData(splat.sim)$Group
ggplot(pca, aes(x = PC1, y = PC2, colour = Step, shape = Path)) +
    geom_point() +
    scale_colour_viridis_c() +
    theme_minimal()

Principal components of splat.sim: the cells of Path1 and Path2 leave the same starting cells, coloured from low to high step.

The workflow

Create the object

scMaSigPro keeps the data and every result in an object of class ScMaSigPro. as_scmp() creates it from a SingleCellExperiment whose cell metadata holds the pseudotime and the path of each cell; here they are the Step and Group columns. create_scmp() does the same from a count matrix and a table of cell metadata.

scmp_ob <- as_scmp(
    object = splat.sim,
    verbose = FALSE,
    additional_params = list(
        labels_exist = TRUE,
        exist_ptime_col = "Step",
        exist_path_col = "Group"
    )
)
scmp_ob
#> Class: ScMaSigProClass
#> nCells: 200
#> nFeatures: 100
#> Pseudotime Range: 1 100
#> Branching Paths: Path2, Path1

Bin the cells along pseudotime

sc.squeeze() divides the pseudotime of each path into bins of equal width and adds up the counts of the cells of each bin. The number of bins comes from a histogram rule, by default that of Sturges, log₂(n) + 1 bins for the n cells of a path, rounded up.

scmp_ob <- sc.squeeze(scmp_ob)
bins <- cDense(scmp_ob)
table(bins$Path)
#> 
#> Path1 Path2 
#>     8     8

The 16 bins hold between 7 and 20 cells. plotBinTile() shows the number of cells of each bin along the binned pseudotime of each path:

plotBinTile(scmp_ob)

Tile plot of the bins of Path1 and Path2 along binned pseudotime, coloured by the number of cells of each bin.

Set up the polynomial model

sc.set.poly() builds the design of the model. With two paths and a polynomial of degree 1, the model of the expected count of a gene in a bin is

\[ \log(\mu) = \beta_0 + \beta_1 D + \beta_2 t + \beta_3 t D + \text{offset}, \]

where \(t\) is the binned pseudotime and \(D\) is 1 for the bins of Path2 and 0 for those of Path1, the reference path. \(\beta_0\) and \(\beta_2\) describe the reference path; \(\beta_1\) and \(\beta_3\) are the differences of Path2 from it, at the start and in slope.

scmp_ob <- sc.set.poly(scmp_ob, poly_degree = 1)
showPoly(scmp_ob)
#> [1] "beta0 + beta1*Path2vsPath1 + beta2*scmp_binned_pseudotime + beta3*scmp_binned_pseudotimexPath2"

A higher degree adds the squares, cubes and so on of \(t\) and their products with \(D\), for profiles that rise and fall; the Model and tuning vignette compares degrees.

Find the genes that change

sc.p.vector() fits the full model of each gene and tests it against a model with the intercept only, a flat profile. The genes whose adjusted p-value is at most p_value have a non-flat profile.

scmp_ob <- sc.p.vector(scmp_ob, p_value = 0.05, verbose = FALSE)
res <- scmpResults(scmp_ob)
sum(res$non_flat)
#> [1] 50

50 of the 100 genes have a non-flat profile.

Select the significant terms

sc.t.fit() then removes, one at a time, the terms of the model of each of these genes that are not significant (backward stepwise selection). What is left says how the gene changes: a gene with only \(\beta_3\) besides the intercept, for example, has the same expression in both paths at the start but a different slope.

scmp_ob <- sc.t.fit(scmp_ob, verbose = FALSE)
head(showSol(scmp_ob)[, c("p-value", "R-squared")])
#>              p-value R-squared
#> Gene2   1.601390e-09 0.7260152
#> Gene3   3.985890e-20 0.8731628
#> Gene4   1.230465e-66 0.9553623
#> Gene9   7.977339e-35 0.9308984
#> Gene10 7.184117e-269 0.9894360
#> Gene11  4.309154e-03 0.5235111

List the significant genes

sc.filter() keeps the genes whose selected model explains enough of their variation, measured by its R-squared, and sorts them into lists. With vars = "groups" there is one list per path: the list of the reference path holds the genes that change along it, and the list Path2vsPath1 the genes whose profile in Path2 differs from that in Path1.

scmp_ob <- sc.filter(scmp_ob, rsq = 0.7, vars = "groups")
res <- scmpResults(scmp_ob)
table(Path1 = res$sig_Path1, Path2vsPath1 = res$sig_Path2vsPath1)
#>        Path2vsPath1
#> Path1   FALSE TRUE
#>   FALSE    65   13
#>   TRUE      1   21

The genes in both lists change along Path1 and differ in Path2; those only in Path2vsPath1 are flat in Path1 but change in Path2; those only in Path1 change along it in the same way in both paths.

plotIntersect() shows how the lists overlap; it draws an UpSet plot with the UpSetR package (Conway et al. 2017):

plotIntersect(scmp_ob)

UpSet plot of the gene lists Path1 and Path2vsPath1: the size of each list and of their intersection.

Plot gene trends

plotTrend() draws the pseudo-bulk counts of a gene, divided by the size factor of their bin, on the log scale against the binned pseudotime, with the fitted curve of each path. Here is one gene of each kind:

genes <- c(both[1], only_diff[1], only_ref[1])
plots <- lapply(genes, function(g) plotTrend(scmp_ob, g, lines = TRUE))
patchwork::wrap_plots(plots, ncol = 2)

Trends of three genes: one that changes in Path1 and differs in Path2, one flat in Path1 that changes in Path2, and one that changes in the same way in both paths.

Cluster the trends

sc.cluster.trend() groups the genes of the lists by the coefficients of their models, here with hierarchical clustering into four clusters; k-means and the model-based clustering of mclust (Scrucca et al. 2016) are the other methods (cluster_method). plotTrendCluster() draws the trends of each cluster:

scmp_ob <- sc.cluster.trend(scmp_ob, geneSet = "union", k = 4)
plotTrendCluster(scmp_ob, verbose = FALSE)

Trends of the four clusters of genes, one panel per cluster, with the median profile of each path.

Get the results

scmpResults() returns every result with one row per gene: the p-value of the full model and whether the profile is non-flat (sc.p.vector()), the p-value and R-squared of the selected model and whether the gene has influential observations (sc.t.fit()), one column per gene list (sc.filter()) and the cluster (sc.cluster.trend()).

head(scmpResults(scmp_ob))
#> DataFrame with 6 rows and 9 columns
#>           p_value adj_p_value  non_flat model_p_value r_squared influential
#>         <numeric>   <numeric> <logical>     <numeric> <numeric>   <logical>
#> Gene1 2.21646e-01 3.46322e-01     FALSE            NA        NA       FALSE
#> Gene2 1.30835e-13 4.84576e-13      TRUE   1.60139e-09  0.726015       FALSE
#> Gene3 2.51249e-24 1.14204e-23      TRUE   3.98589e-20  0.873163       FALSE
#> Gene4 1.96112e-56 1.96112e-55      TRUE   1.23046e-66  0.955362       FALSE
#> Gene5 6.16240e-02 1.04447e-01     FALSE            NA        NA       FALSE
#> Gene6 5.80848e-01 7.21780e-01     FALSE            NA        NA       FALSE
#>       sig_Path2vsPath1 sig_Path1   cluster
#>              <logical> <logical> <integer>
#> Gene1            FALSE     FALSE        NA
#> Gene2             TRUE     FALSE         1
#> Gene3             TRUE      TRUE         2
#> Gene4             TRUE     FALSE         2
#> Gene5            FALSE     FALSE        NA
#> Gene6            FALSE     FALSE        NA

scmpToRowData() writes these columns, with the prefix scmp_, into the rowData() of a SingleCellExperiment, so that they travel with the data:

splat.sim <- scmpToRowData(splat.sim, scmp_ob)
head(rowData(splat.sim)[, c("scmp_non_flat", "scmp_cluster")])
#> DataFrame with 6 rows and 2 columns
#>       scmp_non_flat scmp_cluster
#>           <logical>    <integer>
#> Gene1         FALSE           NA
#> Gene2          TRUE            1
#> Gene3          TRUE            2
#> Gene4          TRUE            2
#> Gene5         FALSE           NA
#> Gene6         FALSE           NA

Compare with the simulated truth

splatter stores in the rowData() how it simulated each gene. DEFacPath1 and DEFacPath2 are the factors by which the expression of a gene changes from the start to the end of each path: a gene whose factor is not 1 on a path changes along that path. In splat.sim every such gene has different factors on the two paths, so it also differs between them.

truth <- rowData(splat.sim)
changes <- truth$DEFacPath1 != 1 | truth$DEFacPath2 != 1
names(changes) <- rownames(truth)
found <- scmpResults(scmp_ob)$non_flat
table(non_flat = found, simulated = changes[rownames(scmpResults(scmp_ob))])
#>         simulated
#> non_flat FALSE TRUE
#>    FALSE    37   13
#>    TRUE      7   43

Of the 50 genes with a non-flat profile, 43 were simulated to change (a precision of 0.86), and they are 43 of the 56 genes simulated to change (a recall of 0.77). The Model and tuning vignette shows how the settings move these numbers.

The trends of the plotted genes follow the factors of the simulation:

as.data.frame(truth[genes, c("DEFacPath1", "DEFacPath2")])
#>        DEFacPath1 DEFacPath2
#> Gene3   0.3643974  1.0000000
#> Gene2   1.0000000  0.8224767
#> Gene92  0.4996190  0.4264288

The bundled object scmp.ob holds the result of this workflow, so the examples of the help pages and the other vignettes can start from it:

data("scmp.ob", package = "scMaSigPro")
identical(scmpResults(scmp.ob), scmpResults(scmp_ob))
#> [1] FALSE

Session information

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] ggplot2_4.0.3               SingleCellExperiment_1.35.2
#>  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [5] GenomicRanges_1.65.4        Seqinfo_1.3.2              
#>  [7] IRanges_2.47.5              S4Vectors_0.51.10          
#>  [9] BiocGenerics_0.59.12        generics_0.1.4             
#> [11] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [13] maSigPro_1.85.0             scMaSigPro_0.6.0           
#> [15] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6        xfun_0.62           bslib_0.12.0       
#>  [4] lattice_0.23-1      vctrs_0.7.3         tools_4.6.1        
#>  [7] parallel_4.6.1      Matrix_1.7-6        RColorBrewer_1.1-3 
#> [10] S7_0.2.2            lifecycle_1.0.5     compiler_4.6.1     
#> [13] farver_2.1.2        codetools_0.2-20    htmltools_0.5.9    
#> [16] sys_3.4.3           buildtools_1.0.0    sass_0.4.10        
#> [19] yaml_2.3.12         venn_1.13           jquerylib_0.1.4    
#> [22] MASS_7.3-66         admisc_0.41         BiocParallel_1.47.0
#> [25] DelayedArray_0.39.8 cachem_1.1.0        abind_1.4-8        
#> [28] mclust_6.1.3        nlme_3.1-171        digest_0.6.39      
#> [31] maketools_1.3.2     labeling_0.4.3      splines_4.6.1      
#> [34] fastmap_1.2.0       grid_4.6.1          cli_3.6.6          
#> [37] SparseArray_1.13.4  patchwork_1.3.2     S4Arrays_1.13.2    
#> [40] withr_3.0.3         scales_1.4.0        rmarkdown_2.32     
#> [43] XVector_0.53.0      otel_0.2.0          gridExtra_2.3.1    
#> [46] evaluate_1.0.5      knitr_1.52          UpSetR_1.4.1       
#> [49] viridisLite_0.4.3   mgcv_1.9-4          rlang_1.3.0        
#> [52] Rcpp_1.1.2          glue_1.8.1          BiocManager_1.30.27
#> [55] jsonlite_2.0.0      plyr_1.8.9          R6_2.6.1

References

Conesa, Ana, María José Nueda, Alberto Ferrer, and Manuel Talón. 2006. “maSigPro: A Method to Identify Significantly Differential Expression Profiles in Time-Course Microarray Experiments.” Bioinformatics 22 (9): 1096–102. https://doi.org/10.1093/bioinformatics/btl056.

Conway, Jake R., Alexander Lex, and Nils Gehlenborg. 2017. “UpSetR: An R Package for the Visualization of Intersecting Sets and Their Properties.” Bioinformatics 33 (18): 2938–40. https://doi.org/10.1093/bioinformatics/btx364.

Nueda, María José, Sonia Tarazona, and Ana Conesa. 2014. “Next maSigPro: Updating maSigPro Bioconductor Package for RNA-seq Time Series.” Bioinformatics 30 (18): 2598–602. https://doi.org/10.1093/bioinformatics/btu333.

Scrucca, Luca, Michael Fop, T. Brendan Murphy, and Adrian E. Raftery. 2016. “mclust 5: Clustering, Classification and Density Estimation Using Gaussian Finite Mixture Models.” The R Journal 8 (1): 289–317. https://doi.org/10.32614/RJ-2016-021.

Srivastava, Priyansh, Marta Benegas Coll, Stefan Götz, María José Nueda, and Ana Conesa. 2024. “scMaSigPro: Differential Expression Analysis Along Single-Cell Trajectories.” Bioinformatics 40 (7): btae443. https://doi.org/10.1093/bioinformatics/btae443.

Zappia, Luke, Belinda Phipson, and Alicia Oshlack. 2017. “Splatter: Simulation of Single-Cell RNA Sequencing Data.” Genome Biology 18: 174. https://doi.org/10.1186/s13059-017-1305-0.

Quick start