Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Pseudotemporal ordering

🧠 Key takeaways

Trajectory inference methods require the start of the biological process to be known (approximately).

The nature of the biological process defines which TI algorithms can be used. dynguidelines helps selecting the appropriate TI method.

⚙️ Environment setup
Steps
yml

  1. Install conda:

    • Before creating the environment, ensure that conda is installed on your system.

  2. Save the yml content:

    • Copy the content from the yml tab into a file named environment.yml.

  3. Create the environment:

    • Open a terminal or command prompt.

    • Run the following command:

      conda env create -f environment.yml

  4. Activate the environment:

    • After the environment is created, activate it using:

      conda activate <environment_name>

    • Replace <environment_name> with the name specified in the environment.yml file. In the yml file it will look like this:

      name: <environment_name>

  5. Verify the installation:

    • Check that the environment was created successfully by running:

      conda env list
🗄️ Get data and notebooks

This book uses lamindb to store, share, and load datasets and notebooks using the theislab/sc-best-practices instance. We acknowledge free hosting from Lamin Labs.

  1. Install lamindb

    • Install the lamindb Python package:

    pip install lamindb

  2. Optionally create a lamin account

  3. Verify your setup

    • Run the lamin connect command:

    import lamindb as ln

ln.Artifact.connect("theislab/sc-best-practices").df()

    You should now see up to 100 of the stored datasets.

  4. Accessing datasets (Artifacts)

    • Search for the datasets on the Artifacts page

    • Load an Artifact and the corresponding object:

    import lamindb as ln
af = ln.Artifact.connect("theislab/sc-best-practices").get(key="key_of_dataset", is_latest=True)
obj = af.load()

    The object is now accessible in memory and is ready for analysis. Adapt the lamindb.Artifact.connect("theislab/sc-best-practices").get("SOMEIDXXXX") suffix to get respective versions.

  5. Accessing notebooks (Transforms)

    lamin load <notebook url>

    which will download the notebook to the current working directory. Analogously to Artifacts, you can adapt the suffix ID to get older versions.

Motivation

Single-cell sequencing assays provide high-resolution measurements of biological tissues Islam et al., 2011, Hwang et al., 2018. Consequently, such technologies can help decipher and understand cellular heterogeneity Briggs et al., 2018, Sikkema et al., 2022 and the dynamics of a biological process Jardine et al., 2021, He et al., 2022.
Corresponding studies include quantifying cellular fates and identifying genes that drive the process. However, as cells are destroyed when sequenced in classical single-cell RNA sequencing (scRNA-seq) protocols, it is impossible to track their development and, for example, gene expression profile over time. Although recent technological advances enable the sequential recording of the transcriptome Chen et al., 2022, they are experimentally challenging and currently fail to scale to larger datasets. Consequently, the underlying dynamic process needs to be estimated from the measured snapshot data instead.

Although samples are traditionally taken from a single experimental time point, a variety of cell types can be observed. This diversity stems from the asynchronous nature of biological processes. As such, a range of the developmental process can be observed. Reconstructing the developmental landscape is the goal of the field, which is referred to as trajectory inference (TI). This task is achieved by ordering the observed cellular states according to the developmental process. States are aligned along the developmental direction by mapping discrete annotations to a continuous domain - the so-called pseudotime.

Pseudotimes rank cells relative to each other according to their respective stage in the developmental process. Less mature cells are assigned small values, while mature cells are assigned large values. Studying a bone marrow sample, for example, hematopoietic stem cells are assigned a low pseudotime, and erythroid cells a high pseudotime. The assignment, in the case of single-cell RNA sequencing data, is based on the transcriptomic profile of a single cell. Additionally, the construction usually requires the specification of an initial, or equivalently, a root cell, where the overall process begins.

Pseudotime construction

Pseudotime construction generally follows a common workflow: As a first step, the ultra-high-dimensional single-cell data is projected onto a lower-dimensional representation. This procedure is justified by the observation that dynamical processes progress on a low-dimensional manifold Wagner et al., 2016. In practice, pseudotime methods may rely on principal components (for example, Palantir Setty et al., 2019) or diffusion components (for example, diffusion pseudotime (DPT) Haghverdi et al., 2016). Following, pseudotimes are constructed based on one of the following principles.

  1. Observations are first clustered and, following, connections between these clusters are identified. The clusters can be ordered, and thereby, a pseudotime can be constructed. Henceforth, we will refer to this approach as the cluster approach. Classical cluster algorithms include kk-means Lloyd, 1982, MacQueen, 1967, Leiden Traag et al., 2019, or hierarchical clustering Müllner, 2011. Clusters may be connected based on similarity, or by constructing a minimum spanning tree (MST) Pettie & Ramachandran, 2002.

  2. The graph approach first finds connections between the lower-dimensional representation of the observations. This procedure defines a graph based on which clusters are defined, and thus an ordering is established. PAGA Wolf et al., 2019, for example, partitions the graph into Leiden clusters and estimates connections between them. Intuitively, this approach preserves the global topology of the data while analyzing it at a lower resolution. Consequently, the computational efficiency is increased.

  3. Manifold-learning based approaches proceed similar to the cluster approach. However, connections between clusters are defined by using principal curves or graphs to estimate the underlying trajectories. Principal curves find a one-dimensional curve that connects cellular observations in the higher-dimensional space. A notable representation of this approach is Slingshot Street et al., 2018.

  4. Probabilistic frameworks assign transition probabilities to ordered cell-cell pairs. Each transition probability quantifies the likelihood that the reference cell is the ancestor of the other cell. These probabilities define random processes that are used to define a pseudotime. DPT (Diffusion Pseudotime), for example, is defined as the difference between consecutive states of a random walk. Contrastingly, Palantir Setty et al., 2019 models trajectories as Markov chains themselves. While both approaches rely on a probabilistic framework, they require a root cell to be specified. The pseudotime itself is computed with respect to this cell.

Trajectory inference (TI) is a well-studied field providing a rich set of methods. To apply the appropriate method for analyzing a single-cell dataset, it is essential to understand the biological process itself. This understanding especially includes the nature of the process, i.e., if it, for example, is linear, cyclic, or branching. Similarly, orthogonal processes within one dataset limits the TI methods, which are applicable. To help identify appropriate tools, dynguidelines Deconinck et al., 2021 provides an exhaustive overview of algorithms and their characteristics.

Orthogonal process in single-cell transcriptomics

Orthogonal processes describe two independent biological programs who are uncorrelated. A cell cycle and a disease would be orthogonal processes.

Down-stream tasks and outlook

Although TI and pseudotime can already provide valuable insights, they typically serve as a stepping stone for more fine-grained analysis. Identifying terminal states, for example, is a classical biological question that can be studied. Similarly, lineage bifurcation and driver genes can be identified based on TI and pseudotime. Which question can be answered, and how the answer is found, is usually method-specific. Palantir Setty et al., 2019, for example, identifies terminal states as absorbing states of its constructed Markov chain.

The success of trajectory inference is well-documented, and consequently, many methods have been proposed. However, with the advances of sequencing technologies, new sources of information have become available. ATAC-seq Buenrostro et al., 2015, CITE-seq Stoeckius et al., 2017, and DOGMA-seq Mimitou et al., 2021, for example, measure additional modalities beyond the transcriptome. Lineage tracing Weinreb et al., 2020 and metabolic labeling Erhard et al., 2019, Battich et al., 2020, Qiu et al., 2020, Erhard et al., 2022 even provide the (likely) future state of a given cell. Consequently, future TI tools will be able to include more information to estimate trajectories and pseudotime more accurately and robustly, allowing for the answer to novel questions. For example, RNA velocity Manno et al., 2018, Bergen et al., 2020, Bergen et al., 2021 is one technique that uses unspliced and spliced mRNA to infer directed, dynamic information beyond classical, static snapshot data.

Inferring pseudotime for adult human bone marrow

To demonstrate how a pseudotime can be constructed and how different pseudotimes can be compared, we analyze a dataset of adult human bone marrow Setty et al., 2019.

Environment setup

import lamindb as ln
import scanpy as sc

assert ln.setup.settings.instance.slug == "theislab/sc-best-practices"

ln.track("5tukHKyz8kFR")
→ connected lamindb: theislab/sc-best-practices
→ loaded Transform('5tukHKyz8kFR0001', key='pseudotemporal.ipynb'), re-started Run('Hfx6ZbLSPyO3msPK') at 2026-01-14 17:24:35 UTC
→ notebook imports: lamindb==2.0.1 scanpy==1.11.5

Data loading

af = ln.Artifact.connect("theislab/sc-best-practices").get(
    key="trajectory/pseudotemporal.h5ad", is_latest=True
)
adata = af.load()
adata
... synchronizing pseudotemporal.h5ad: 100.0%
AnnData object with n_obs × n_vars = 5780 × 27876 obs: 'clusters', 'palantir_pseudotime', 'palantir_diff_potential' var: 'palantir' uns: 'clusters_colors', 'palantir_branch_probs_cell_types' obsm: 'MAGIC_imputed_data', 'X_tsne', 'palantir_branch_probs' layers: 'spliced', 'unspliced'

To construct pseudotimes, the data must be preprocessed. Here, we filter out genes expressed in only a small number of cells (here, at least 20). Notably, the construction of the pseudotime later on is robust to the exact choice of the threshold. Following this initial gene filtering, the cell size is normalized, and the counts are log1p-transformed to reduce the effect of outliers. As usual, we also identify and annotate highly variable genes. Finally, a nearest neighbor graph is constructed based on which we will define the pseudotime. The number of principal components is chosen based on the explained variance.

sc.pp.filter_genes(adata, min_counts=20)
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata)

sc.tl.pca(adata)
sc.pp.neighbors(adata, n_pcs=10)
/Users/seohyon/miniconda3/envs/pseudotemporal/lib/python3.13/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

sc.pl.scatter(adata, basis="tsne", color="clusters")
<Figure size 852x480 with 1 Axes>

The two-dimensional t-SNE representation colored by cell type annotations shows that the cell types cluster together well.

Additionally, the developmental hierarchy is visible.

Pseudotime construction

To calculate diffusion pseudotime (DPT), the corresponding diffusion maps must first be calculated.

sc.tl.diffmap(adata)

The differentiation hierarchy in bone marrow is well understood. However, it is only known that the developmental process starts in the form of hematopoietic stem cells, but not which cell, exactly, in the corresponding cluster in our dataset. To identify a putative initial cell, we study the individual diffusion components. We identify the stem cell with the most extreme diffusion component in one dimension (in our case, dimension 3).

Diffusion map and diffusion components

Diffusion map shows us the non-linear distance between one state of a cell to another. It is important for this distance to be non-linear, since cells follow a continuous path of gradual changes, rather than “jumping” from one state to another. Unlike PCA, diffusion map ignores the Euclidean distance and tries to find the curved path the data follows.

Diffusion components is the eigenvector of a transition matrix that represents the “preferred directions” of movement in the diffusion map. Lower-order diffusion components capture more significant biological transitions.

# Setting root cell as described above
root_ixs = adata.obsm["X_diffmap"][:, 3].argmin()
sc.pl.scatter(
    adata,
    basis="diffmap",
    color=["clusters"],
    components=[2, 3],
)

adata.uns["iroot"] = root_ixs
<Figure size 852x480 with 1 Axes>
sc.tl.dpt(adata)

Different pseudotime methods give different results. Sometimes, one pseudotime captures the underlying developmental processes more accurately than others. Here, we compare the just-computed diffusion pseudotime (DPT) with the pre-computed Palantir pseudotime (see here for the corresponding tutorial). One option for comparing different pseudotimes is to color the low-dimensional embedding of the data (here, t-SNE). Here, DPT is extremely high in the cluster of CLPs compared to all other cell types. Contrastingly, the Palantir pseudotime increases continuously with developmental maturity.

sc.pl.scatter(
    adata,
    basis="tsne",
    color=["dpt_pseudotime", "palantir_pseudotime"],
    color_map="gnuplot2",
)
<Figure size 1495.76x480 with 4 Axes>

Instead of coloring the lower-dimensional representation of the data, we can study the distribution of pseudotime values assigned to each cell type cluster. This representation again shows that the CLP cluster forms an outlier in the case of DPT. Additionally, clusters such as HSC_1 and HSC_2 include several cells with increased pseudotime. These inflated values contrast with our prior biological knowledge that these clusters form the beginning of the developmental process.

sc.pl.violin(
    adata,
    keys=["dpt_pseudotime", "palantir_pseudotime"],
    groupby="clusters",
    rotation=45,
    order=[
        "HSC_1",
        "HSC_2",
        "Precursors",
        "Ery_1",
        "Ery_2",
        "Mono_1",
        "Mono_2",
        "CLP",
        "DCs",
        "Mega",
    ],
)
<Figure size 1609.2x480 with 2 Axes>

Considering these observations and prior knowledge about the development in bone marrow, we would conclude to continue working with the Palantir pseudotime.

Contributors

We gratefully acknowledge the contributions of:

Authors

  • Philipp Weiler

  • Seo H. Kim

Reviewers

  • Lukas Heumos

References
  1. Islam, S., Kjällquist, U., Moliner, A., Zajac, P., Fan, J.-B., Lönnerberg, P., & Linnarsson, S. (2011). Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. Genome Research, 21(7), 1160–1167. 10.1101/gr.110882.110
  2. Hwang, B., Lee, J. H., & Bang, D. (2018). Single-cell RNA sequencing technologies and bioinformatics pipelines. Experimental &\mathsemicolon Molecular Medicine, 50(8), 1–14. 10.1038/s12276-018-0071-8
  3. Briggs, J. A., Weinreb, C., Wagner, D. E., Megason, S., Peshkin, L., Kirschner, M. W., & Klein, A. M. (2018). The dynamics of gene expression in vertebrate embryogenesis at single-cell resolution. Science, 360(6392). 10.1126/science.aar5780
  4. Sikkema, L., Strobl, D. C., Zappia, L., Madissoon, E., Markov, N. S., Zaragosi, L.-E., Ansari, M., Arguel, M.-J., Apperloo, L., Becavin, C., Berg, M., Chichelnitskiy, E., Chung, M.-I., Collin, A., Gay, A. C. A., Kashani, B. H., Jain, M., Kapellos, T., Kole, T. M., … Theis, F. J. (2022). An integrated cell atlas of the human lung in health and disease. bioRxiv, 2022.03.10.483747. 10.1101/2022.03.10.483747
  5. Jardine, L., Webb, S., Goh, I., Londoño, M. Q., Reynolds, G., Mather, M., Olabi, B., Stephenson, E., Botting, R. A., Horsfall, D., Engelbert, J., Maunder, D., Mende, N., Murnane, C., Dann, E., McGrath, J., King, H., Kucinski, I., Queen, R., … Haniffa, M. (2021). Blood and immune development in human fetal bone marrow and Down syndrome. Nature, 598(7880), 327–331. 10.1038/s41586-021-03929-x
  6. He, P., Lim, K., Sun, D., Pett, J. P., Jeng, Q., Polanski, K., Dong, Z., Bolt, L., Richardson, L., Mamanova, L., Dabrowska, M., Wilbrey-Clark, A., Madissoon, E., Tuong, Z. K., Dann, E., Suo, C., Goh, I., Yoshida, M., Nikolić, M. Z., … Rawlins, E. L. (2022). A human fetal lung cell atlas uncovers proximal-distal gradients of differentiation and key regulators of epithelial fates. bioRxiv. 10.1101/2022.01.11.474933
  7. Chen, W., Guillaume-Gentil, O., Rainer, P. Y., Gäbelein, C. G., Saelens, W., Gardeux, V., Klaeger, A., Dainese, R., Zachara, M., Zambelli, T., Vorholt, J. A., & Deplancke, B. (2022). Live-seq enables temporal transcriptomic recording of single cells. Nature, 608(7924), 733–740. 10.1038/s41586-022-05046-9
  8. Wagner, A., Regev, A., & Yosef, N. (2016). Revealing the vectors of cellular identity with single-cell genomics. Nature Biotechnology, 34(11), 1145–1160. 10.1038/nbt.3711
  9. Setty, M., Kiseliovas, V., Levine, J., Gayoso, A., Mazutis, L., & Pe’er, D. (2019). Characterization of cell fate probabilities in single-cell data with Palantir. Nature Biotechnology, 37(4), 451–460. 10.1038/s41587-019-0068-4
  10. Haghverdi, L., Büttner, M., Wolf, F. A., Buettner, F., & Theis, F. J. (2016). Diffusion pseudotime robustly reconstructs lineage branching. Nature Methods, 13(10), 845–848. 10.1038/nmeth.3971
  11. Lloyd, S. (1982). Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2), 129–137. 10.1109/tit.1982.1056489
  12. MacQueen, J. B. (1967). Some Methods for Classification and Analysis of MultiVariate Observations. In L. M. L. Cam & J. Neyman (Eds.), Proc. of the fifth Berkeley Symposium on Mathematical Statistics and Probability (Vol. 1, pp. 281–297). University of California Press.
  13. Traag, V. A., Waltman, L., & van Eck, N. J. (2019). From Louvain to Leiden: guaranteeing well-connected communities. Scientific Reports, 9(1). 10.1038/s41598-019-41695-z
  14. Müllner, D. (2011). Modern hierarchical, agglomerative clustering algorithms. arXiv. 10.48550/ARXIV.1109.2378
  15. Pettie, S., & Ramachandran, V. (2002). An optimal minimum spanning tree algorithm. Journal of the ACM, 49(1), 16–34. 10.1145/505241.505243
Identifying cellular structure
Data integration
Inferring trajectories
RNA velocity