Time Aware Trajectory Analysis (TATA) with scProkaR

Introduction

scProkaR provides tools for microbial single-cell RNA sequencing (scRNA-seq) analysis using SingleCellExperiment objects. This vignette describes Time Aware Trajectory Analysis (TATA), the package’s trajectory component. Data import, quality control, integration and pseudobulk differential expression are covered in the core workflow vignette.

The problem TATA addresses

Perturbation time courses measure bacterial populations at defined intervals after exposure to a drug or stress. Trajectory analysis can describe changes in cell states and gene expression across these observations.

Expression similarity alone does not determine the direction of a transition. Sparse measurements can also make neighbourhoods uncertain. In a branching response, the earliest population may occupy an interior region of the embedding, so position alone is insufficient for choosing a root.

TATA uses observed sampling times together with expression neighbourhoods to score candidate connections between clusters. The inferred graph describes relationships consistent with these observations; temporal order alone does not establish a lineage relationship.

For each candidate connection, TATA computes a topology weight from cell-to-cell connectivity and a temporal score from the sampling times. These scores determine edge weights, pruning and direction. Connections with insufficient directional evidence can remain bidirectional if they pass the pruning threshold.

Relationship to Bioconductor single-cell infrastructure

TATA accepts a SingleCellExperiment with an embedding stored in reducedDims(). run_tata() returns a list containing the updated object, graphs, score tables and other results.

  • The input embedding may come from scater, scran, IntegrateBacData() or another analysis.
  • The updated object’s colData() contains clusters, pseudotime, terminal probabilities, branch entropy and branch assignments. The "TATA" embedding is stored in reducedDims(), with graph results in metadata().
  • Plotting tools can use these annotations directly. Downstream statistical analyses must still account for independent samples and the uncertainty in inferred branches.

TATA draws on PAGA’s graph abstraction and Palantir’s terminal-state probability model. Its implementation combines these ideas with observed sampling times, as described below.

Installation

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")
BiocManager::install("scProkaR")

Setting up

library(scProkaR)
library(SingleCellExperiment)
library(SummarizedExperiment)

TATA uses FNN for nearest-neighbour search when available, followed by BiocNeighbors as a fallback. If neither is installed, it uses stats::dist(). The distance-matrix fallback requires substantially more memory for large datasets.

c(
    FNN = requireNamespace("FNN", quietly = TRUE),
    BiocNeighbors = requireNamespace("BiocNeighbors", quietly = TRUE)
)
#>           FNN BiocNeighbors 
#>          TRUE          TRUE

A simulated multi-drug branching time course

simulate_tata_multidrug_sce() generates a SingleCellExperiment with three simulated drug responses sampled at six times between 0 and 120 hours. The branch activation onsets are 24 hours for drug A, 48 hours for drug C, and 96 hours for drug B.

The earliest cells occupy a shared region near the centre of the stored embedding. This design provides a case in which sampling time can help identify the root independently of the embedding’s endpoints.

set.seed(101)

sce <- simulate_tata_multidrug_sce(
    n_cells = 3000,
    n_features = 80,
    n_pcs = 15
)

sce$timepoint_factor <- factor(
    sce$timepoint,
    levels = sort(unique(sce$timepoint))
)

sce
#> class: SingleCellExperiment 
#> dim: 80 3000 
#> metadata(1): simulation
#> assays(2): counts logcounts
#> rownames(80): gene_1 gene_2 ... gene_79 gene_80
#> rowData names(1): feature_module
#> colnames(3000): bacterium_1 bacterium_2 ... bacterium_2999
#>   bacterium_3000
#> colData names(8): cell_id timepoint ... branch_activation
#>   timepoint_factor
#> reducedDimNames(3): PCA UMAP truth
#> mainExpName: NULL
#> altExpNames(0):

All three conditions are represented at each sampled time, with varying cell numbers:

with(
    as.data.frame(colData(sce)),
    table(condition, timepoint)
)
#>          timepoint
#> condition   0  24  48  72  96 120
#>    drug_A 163 165 162 181 162 165
#>    drug_B 200 146 200 165 200 138
#>    drug_C 185 123 192 138 192 123

The stored UMAP views show sampling time, treatment and simulated state. Compare the locations of the earliest cells with the three response branches.

PlotReduction(sce, reduction = "UMAP", colour_by = "timepoint_factor")

PlotReduction(sce, reduction = "UMAP", colour_by = "condition")

PlotReduction(sce, reduction = "UMAP", colour_by = "simulated_branch")

The TATA pipeline, step by step

Use run_tata() for the complete workflow. The individual functions below allow each intermediate graph and score table to be inspected, or clustering and edge scoring to be supplied separately.

Build the cell–cell graph

TATA operates on a low-dimensional representation. Here we use the stored PCA and connect each cell to its k nearest neighbours.

knn <- build_knn_graph(
    x = sce,
    dimred = "PCA",
    k = 30
)

c(
    n_cells = igraph::vcount(knn$graph),
    n_edges = igraph::ecount(knn$graph),
    k = knn$k
)
#> n_cells n_edges       k 
#>    3000   54721      30
knn$knn_engine
#> [1] "FNN"

Cluster cells into coarse states

Clustering partitions the cell graph into states that become nodes of the abstract graph. The partition affects the inferred transitions, so inspect cluster sizes and time distributions when choosing analysis settings.

sce$tata_cluster_step <- cluster_graph_states(
    cell_graph = knn$graph,
    method = "louvain"
)

table(sce$tata_cluster_step)
#> 
#>  C1  C2  C3  C4  C5  C6  C7  C8  C9 C10 C11 C12 C13 C14 C15 C16 C17 
#> 133 363 171 126 445 265 147 161 147 200  81  75 142 120  95 148 181
PlotReduction(sce, reduction = "UMAP", colour_by = "tata_cluster_step")

When cluster is installed, silhouette widths can be calculated on a subsample to inspect cluster separation. This diagnostic is optional and does not determine temporal validity.

if (requireNamespace("cluster", quietly = TRUE)) {
    set.seed(101)
    idx <- sample(seq_len(ncol(sce)), 600)
    pcs <- reducedDim(sce, "PCA")[idx, seq_len(5)]
    sil <- cluster::silhouette(
        as.integer(sce$tata_cluster_step[idx]),
        stats::dist(pcs)
    )
    round(mean(sil[, "sil_width"]), 3)
} else {
    message("Install 'cluster' for the optional silhouette diagnostic.")
}
#> [1] 0.42

Score cluster pairs on topology and on time

Each candidate edge receives a topology weight and a temporal score.

compute_topology_weights() compares observed connections between each cluster pair with the number expected from their sizes. Their ratio is reported as topology_weight.

topology_table <- compute_topology_weights(
    adjacency = knn$adjacency,
    clusters = sce$tata_cluster_step
)

head(
    topology_table[
        order(-topology_table$topology_weight),
        c(
            "cluster_a", "cluster_b", "observed_edges",
            "expected_edges", "topology_weight"
        )
    ],
    6
)
#>     cluster_a cluster_b observed_edges expected_edges topology_weight
#> 47         C4        C6            420       439.0706       0.9565660
#> 27         C2       C13            514       640.2096       0.8028620
#> 106        C9       C15            112       144.6020       0.7745396
#> 66         C5       C13            615       802.9965       0.7658813
#> 134       C15       C16            105       146.8176       0.7151732
#> 72         C6        C8            365       548.7031       0.6652049

compute_time_weights() compares sampling times across edges joining two clusters. flow_score ranges from -1 to 1 and summarises their net temporal direction. time_consistency combines directional agreement, differences in cluster median times, temporal overlap and a penalty for large time gaps. It is a score, not a statistical confidence level.

time_table <- compute_time_weights(
    adjacency = knn$adjacency,
    clusters = sce$tata_cluster_step,
    timepoint = sce$timepoint
)

head(
    time_table[
        order(-time_table$time_consistency),
        c(
            "cluster_a", "cluster_b", "n_time_edges",
            "flow_score", "time_weight", "time_consistency"
        )
    ],
    6
)
#>     cluster_a cluster_b n_time_edges flow_score time_weight time_consistency
#> 72         C6        C8          365  0.3260274   0.9760209        0.9520417
#> 122       C12       C13            3 -0.6666667   0.1519838        0.6960323
#> 19         C2        C5           12  0.0000000   0.1867981        0.6264037
#> 45         C3       C17          185  0.1621622   0.8093046        0.6186093
#> 38         C3       C10           99 -0.2222222   0.1923928        0.6152143
#> 125       C12       C16           10  0.2000000   0.8071993        0.6143985

The scatter plot compares topology weights with flow_score, colouring points by time_consistency. Values near zero indicate little net direction across the connecting cell pairs. A near-zero flow can reflect equal times or opposing directions; inspect the time distributions to distinguish them.

edge_metric_table <- merge(
    topology_table,
    time_table,
    by = c("cluster_a", "cluster_b"),
    all.x = TRUE,
    sort = FALSE
)

ggplot2::ggplot(
    edge_metric_table,
    ggplot2::aes(
        x = topology_weight,
        y = flow_score,
        size = observed_edges,
        colour = time_consistency
    )
) +
    ggplot2::geom_hline(
        yintercept = c(-0.10, 0.10),
        linetype = 2,
        colour = "grey60"
    ) +
    ggplot2::geom_point(alpha = 0.8) +
    ggplot2::scale_colour_viridis_c(option = "plasma") +
    ggplot2::theme_classic() +
    ggplot2::labs(
        title = "Topology and time consistency for cluster pairs",
        x = "topology weight",
        y = "temporal flow score",
        colour = "time consistency",
        size = "observed edges"
    )

Build the abstract graph

build_tata_graph() combines the weights as topology_weight^alpha * time_factor^beta. With the default time_weight_mode = "directional_confidence", time_factor is time_consistency; the "ordered" mode uses time_weight. The function prunes edges below prune_threshold and assigns direction using signed_time_score and direction_threshold. Retained edges within the directional threshold are represented by connections in both directions.

tata_manual <- build_tata_graph(
    topology_table = topology_table,
    time_table = time_table,
    clusters = sce$tata_cluster_step,
    embedding = reducedDim(sce, "PCA"),
    timepoint = sce$timepoint,
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    time_weight_mode = "directional_confidence"
)

head(
    tata_manual$edge_table[
        order(-tata_manual$edge_table$tata_weight),
        c(
            "cluster_a", "cluster_b", "topology_weight", "flow_score",
            "time_consistency", "tata_weight", "direction_label", "kept"
        )
    ],
    8
)
#>     cluster_a cluster_b topology_weight  flow_score time_consistency
#> 72         C6        C8       0.6652049  0.32602740        0.9520417
#> 27         C2       C13       0.8028620 -0.28210117        0.5877679
#> 106        C9       C15       0.7745396  0.22321429        0.5906241
#> 134       C15       C16       0.7151732 -0.09523810        0.5765472
#> 10         C1       C11       0.5957840 -0.03076923        0.6043222
#> 18         C2        C4       0.5895465  0.21408451        0.6101644
#> 123       C12       C14       0.5846232 -0.09433962        0.5956485
#> 45         C3       C17       0.4667605  0.16216216        0.6186093
#>     tata_weight direction_label kept
#> 72    0.6333028        C6 -> C8 TRUE
#> 27    0.4718965       C13 -> C2 TRUE
#> 106   0.4574618       C9 -> C15 TRUE
#> 134   0.4123311      C16 -> C15 TRUE
#> 10    0.3600455       C11 -> C1 TRUE
#> 18    0.3597203        C2 -> C4 TRUE
#> 123   0.3482300      C14 -> C12 TRUE
#> 45    0.2887424       C3 -> C17 TRUE

The embedding layout places each cluster node at the centre of its cells, allowing graph connections to be compared with the cell distribution.

plot_tata_cluster_graph(
    list(
        cluster_graph = tata_manual$graph,
        sce = sce,
        parameters = list(cluster_col = "tata_cluster_step")
    ),
    dimred = "UMAP",
    layout_mode = "embedding",
    show_cells = TRUE,
    colour_by = "tata_cluster_step",
    node_colour_by = "cluster",
    main = "Manual TATA graph aligned to UMAP"
)

The same graph in a PAGA-style abstract layout, with nodes coloured by cluster median time, shows the temporal gradient directly.

plot_tata_cluster_graph(
    list(
        cluster_graph = tata_manual$graph,
        sce = sce,
        parameters = list(cluster_col = "tata_cluster_step")
    ),
    layout_mode = "graph",
    node_colour_by = "median_time",
    main = "Manual TATA abstract graph"
)

Compute pseudotime

compute_tata_pseudotime() assigns pseudotime to clusters and then cells. With sampling times supplied, root selection prioritises clusters containing the largest fraction of earliest-time cells, using connectivity and size to resolve ties. Sampling times also inform cluster pseudotime and within-cluster cell ordering. Without times, cluster distances use inverse edge weights.

sce$tata_cluster_step <- factor(sce$tata_cluster_step)

pt_manual <- compute_tata_pseudotime(
    cluster_graph = tata_manual$graph,
    sce = sce,
    cluster_col = "tata_cluster_step",
    time_col = "timepoint"
)

pt_manual$root_cluster
#> [1] "C5"
summary(pt_manual$cell_pseudotime_scaled)
#>    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
#>  0.0000  0.2196  0.3790  0.4551  0.7846  1.0000

Inspect the median sampling time of the selected root cluster:

tapply(
    sce$timepoint,
    sce$tata_cluster_step,
    median
)[pt_manual$root_cluster]
#> C5 
#>  0

The run_tata() wrapper

run_tata() combines graph construction, clustering, optional temporal refinement, edge scoring, graph abstraction, pseudotime and branch probabilities. Here, expected_branches = 3 supplies the known number of simulated terminal states; branch count is therefore an input to this example.

set.seed(101)

tata <- run_tata(
    sce = sce,
    dimred = "PCA",
    time_col = "timepoint",
    cluster_col = "tata_cluster",
    k = 30,
    cluster_method = "louvain",
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    time_weight_mode = "directional_confidence",
    do_pseudotime = TRUE,
    expected_branches = 3
)

names(tata)
#>  [1] "sce"                          "cell_graph"                  
#>  [3] "adjacency_matrix"             "cluster_graph"               
#>  [5] "cluster_edge_table"           "cluster_vertex_table"        
#>  [7] "cluster_pseudotime"           "cell_pseudotime"             
#>  [9] "cell_pseudotime_scaled"       "cell_space"                  
#> [11] "terminal_clusters"            "cluster_branch_probabilities"
#> [13] "branch_probabilities"         "transition_matrix"           
#> [15] "parameters"                   "topology_table"              
#> [17] "time_table"

The updated object’s annotations and embeddings can be inspected as follows:

sce_tata <- tata$sce
sce_tata$tata_cluster <- factor(sce_tata$tata_cluster)

grep("^tata_", colnames(colData(sce_tata)), value = TRUE)
#>  [1] "tata_cluster_step"           "tata_cluster"               
#>  [3] "tata_pseudotime"             "tata_pseudotime_scaled"     
#>  [5] "tata_prob_C16"               "tata_prob_C2"               
#>  [7] "tata_prob_C15"               "tata_branch_entropy"        
#>  [9] "tata_branch_plasticity"      "tata_branch_commitment"     
#> [11] "tata_max_branch_probability" "tata_branch_assignment"
reducedDimNames(sce_tata)
#> [1] "PCA"   "UMAP"  "truth" "TATA"
PlotReduction(sce_tata, reduction = "UMAP", colour_by = "tata_cluster")

PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_pseudotime_scaled"
)

Compare the inferred ordering with the known progression from the shared early population along each simulated branch.

plot_tata_cluster_graph(
    tata,
    dimred = "UMAP",
    layout_mode = "embedding",
    show_cells = TRUE,
    colour_by = "timepoint",
    node_colour_by = "cluster",
    main = "TATA abstract graph aligned to UMAP"
)

plot_tata_cluster_graph(
    tata,
    layout_mode = "graph",
    node_colour_by = "median_time",
    main = "TATA abstract graph"
)

plot_tata_trajectory_embedding() traces root-to-terminal paths through cluster centroids. Use max_paths to limit the number of paths displayed when the graph has several terminal states.

plot_tata_trajectory_embedding(
    tata,
    dimred = "UMAP",
    colour_by = "timepoint",
    smooth_paths = TRUE
)

plot_tata_trajectory_embedding(
    tata,
    dimred = "UMAP",
    colour_by = "timepoint",
    max_paths = 3,
    smooth_paths = TRUE
)

Comparing time-informed and topology-only runs

The following run sets use_time = FALSE and keeps the same embedding, random seed and neighbourhood settings. It disables temporal edge scoring and temporal cluster refinement, and changes how roots and pseudotime are computed. This compares the two workflow settings rather than isolating the effect of edge direction alone.

set.seed(101)

tata_topology <- run_tata(
    sce = sce,
    dimred = "PCA",
    use_time = FALSE,
    time_col = NULL,
    cluster_col = "tata_cluster_topology",
    k = 30,
    cluster_method = "louvain",
    alpha = 1,
    beta = 1,
    direction_threshold = 0.10,
    prune_threshold = 0.001,
    do_pseudotime = TRUE,
    expected_branches = 3
)

Benchmarking helpers

The following helpers assess temporal order, graph direction and branch onset. orient_pseudotime_to_time() first aligns the sign of each pseudotime vector with observed sampling time. This affects pseudotime-based metrics; it does not alter graph edges.

TATA already uses these sampling times during inference. Agreement with them is therefore a consistency check, not independent validation of the trajectory. The simulation also supplies branch labels and activation onsets for comparison.

pt_tata <- orient_pseudotime_to_time(
    pseudotime = tata$cell_pseudotime,
    timepoint = tata$sce$timepoint
)

pt_topology <- orient_pseudotime_to_time(
    pseudotime = tata_topology$cell_pseudotime,
    timepoint = tata_topology$sce$timepoint
)

Local temporal order. local_temporal_order_metrics() compares the direction of pseudotime and sampling-time differences across neighbouring cells. It measures agreement within the cell graph.

local_tata <- local_temporal_order_metrics(
    pseudotime = pt_tata,
    adjacency = tata$adjacency_matrix,
    timepoint = tata$sce$timepoint
)

local_topology <- local_temporal_order_metrics(
    pseudotime = pt_topology,
    adjacency = tata_topology$adjacency_matrix,
    timepoint = tata_topology$sce$timepoint
)

rbind(TATA = local_tata, topology_only = local_topology)
#>               n_compared_edges local_time_concordance
#> TATA                     10656             0.96846847
#> topology_only            10656             0.04588964
#>               local_temporal_inversion_rate local_temporal_direction_score
#> TATA                             0.03153153                      0.9684685
#> topology_only                    0.04635886                      0.9536411
#>               pseudotime_tie_rate
#> TATA                    0.0000000
#> topology_only           0.9077515

Adjacent-timepoint AUC. adjacent_timepoint_auc() measures how often pseudotime orders a cell from the next sampled time after one from the preceding time. It compares adjacent sampled times rather than all pairs. Passing group computes the comparisons within each treatment group.

auc_tata <- adjacent_timepoint_auc(
    pseudotime = pt_tata,
    timepoint = tata$sce$timepoint,
    group = tata$sce$condition
)

auc_topology <- adjacent_timepoint_auc(
    pseudotime = pt_topology,
    timepoint = tata_topology$sce$timepoint,
    group = tata_topology$sce$condition
)

rbind(TATA = auc_tata$summary, topology_only = auc_topology$summary)
#>               adjacent_timepoint_auc n_valid_comparisons averaging
#> TATA                       0.9937871                  15  weighted
#> topology_only              0.7767164                  15  weighted

Inspect the per-pair table to identify intervals with weaker temporal separation:

head(auc_tata$pair_table)
#>    group time_earlier time_later n_earlier n_later       auc weight
#> 1 drug_B            0         24       200     146 0.9855822  29200
#> 2 drug_B           24         48       146     200 0.9933219  29200
#> 3 drug_B           48         72       200     165 0.9961212  33000
#> 4 drug_B           72         96       165     200 1.0000000  33000
#> 5 drug_B           96        120       200     138 0.9954348  27600
#> 6 drug_C            0         24       185     123 0.9915623  22755

Directed reachability. directed_reachability_concordance() checks whether an earlier cluster can reach a later cluster through directed edges. Here, pairs must differ by at least 24 hours in median sampling time. The pair table distinguishes forward, backward, bidirectional and absent connections; inspect it alongside the overall score.

reach_tata <- directed_reachability_concordance(
    cluster_graph = tata$cluster_graph,
    clusters = tata$sce$tata_cluster,
    timepoint = tata$sce$timepoint,
    min_time_gap = 24
)

reach_tata$summary
#>   n_ordered_cluster_pairs directed_reachability_concordance weighting
#> 1                      76                         0.4744113     cells
head(reach_tata$pair_table)
#>   cluster_earlier cluster_later median_time_earlier median_time_later time_gap
#> 1              C1            C5                   0                48       48
#> 2              C1            C6                   0                48       48
#> 3              C1            C7                   0                48       48
#> 4              C1            C8                   0                72       72
#> 5              C1            C9                   0                72       72
#> 6              C1           C10                   0                72       72
#>   size_earlier size_later forward_reachable backward_reachable
#> 1          452        276              TRUE              FALSE
#> 2          452        217             FALSE              FALSE
#> 3          452        148              TRUE              FALSE
#> 4          452        193              TRUE              FALSE
#> 5          452        146             FALSE              FALSE
#> 6          452         95              TRUE              FALSE
#>   reachability_score weight
#> 1                  1 124752
#> 2                  0  98084
#> 3                  1  66896
#> 4                  1  87236
#> 5                  0  65992
#> 6                  1  42940

Anti-causal edge mass. anticausal_edge_mass() reports the weight of edges directed backwards in median sampling time beyond a specified tolerance. It requires cluster times and a graph, but no pseudotime. A high value indicates disagreement with the sampled time order and warrants inspection.

anticausal_tata <- anticausal_edge_mass(
    cluster_graph = tata$cluster_graph,
    clusters = tata$sce$tata_cluster,
    timepoint = tata$sce$timepoint,
    tolerance = 24
)

anticausal_tata$summary
#>   total_edge_weight anticausal_edge_mass anticausal_edge_mass_fraction
#> 1          7.067314                    0                             0
#>   causal_edge_score
#> 1                 1
head(
    anticausal_tata$edge_table[
        order(anticausal_tata$edge_table$time_delta),
        c("from", "to", "weight", "time_delta", "temporal_class")
    ]
)
#>    from  to      weight time_delta temporal_class
#> 3    C3  C1 0.007905095        -24        neutral
#> 7    C6  C4 0.200614029        -24        neutral
#> 17  C12 C10 0.457461773        -24        neutral
#> 1    C1  C2 0.055786518          0        neutral
#> 2    C2  C1 0.055786518          0        neutral
#> 5    C3  C4 0.043642082          0        neutral

Branch onset ordering. branch_onset_metrics() compares inferred pseudotime at branch activation with the known simulated onsets: 24 hours for A, 48 hours for C, and 96 hours for B. The comparison assesses whether the onset order is retained in each workflow.

true_onset <- c(drug_A = 24, drug_B = 96, drug_C = 48)

onset_tata <- branch_onset_metrics(
    pseudotime = pt_tata,
    branch = tata$sce$simulated_branch,
    branch_activation = tata$sce$branch_activation,
    true_onset = true_onset
)

onset_topology <- branch_onset_metrics(
    pseudotime = pt_topology,
    branch = tata_topology$sce$simulated_branch,
    branch_activation = tata_topology$sce$branch_activation,
    true_onset = true_onset
)

rbind(TATA = onset_tata$summary, topology_only = onset_topology$summary)
#>               n_branches_used kendall_tau spearman_rho normalized_onset_mae
#> TATA                        3   1.0000000          1.0           0.05853579
#> topology_only               3   0.3333333          0.5           0.15880215
#>               branch_onset_order_score branch_onset_accuracy
#> TATA                         1.0000000             0.9414642
#> topology_only                0.6666667             0.8411978
#>               activation_threshold
#> TATA                          0.25
#> topology_only                 0.25
onset_tata$onset_table
#>   branch true_onset inferred_onset n_committed_cells
#> 1 drug_A         24      0.7616570               658
#> 2 drug_B         96      0.9892122               126
#> 3 drug_C         48      0.8774691               443

Temporal consistency plots

plot_tata_time_calibration() displays pseudotime against sampling time within each treatment group. Inspect their ordering and spread; equal slopes across treatments are not required. Differences may reflect response rates, branch structure or the inferred ordering.

plot_tata_time_calibration(
    tata,
    branch_col = "condition",
    summary = "median_iqr"
)

plot_tata_time_calibration(
    tata_topology,
    time_col = "timepoint",
    branch_col = "condition",
    summary = "median_iqr"
)

plot_tata_edge_causality() displays topology weight against signed temporal evidence and marks retained edges. The sign refers to the ordering of each cluster pair; either sign can support a directed connection. The plot describes temporal evidence rather than establishing biological causality.

plot_tata_edge_causality(tata)

Branch probabilities and plasticity

TATA estimates terminal-state probabilities using an absorbing Markov model on the cluster graph, following the approach used in Palantir. compute_tata_branch_probabilities() returns these estimates directly; run_tata() includes them when do_branch_probs = TRUE.

branch_result <- compute_tata_branch_probabilities(
    tata_result = tata,
    expected_branches = 3
)

branch_result$terminal_clusters
#> [1] "C16" "C2"  "C15"

The terminal-state probabilities are summarised by the following quantities:

  • Entropy is normalised Shannon entropy over terminal-state probabilities. Higher values indicate more diffuse model predictions; entropy is used here as a model-based measure of plasticity.
  • Commitment is calculated as one minus normalised entropy. The maximum terminal probability is reported separately as tata_max_branch_probability.
  • Assignment summarises the probabilities with a branch label for subsetting. It is an inferred annotation, not observed cell fate.
prob_cols <- grep(
    "^tata_prob_",
    colnames(tata$cluster_branch_probabilities),
    value = TRUE
)

head(
    tata$cluster_branch_probabilities[
        ,
        c("cluster", prob_cols, "tata_branch_entropy",
            "tata_branch_assignment")
    ]
)
#>    cluster tata_prob_C16 tata_prob_C2 tata_prob_C15 tata_branch_entropy
#> C1      C1             0  1.000000000     0.0000000           0.0000000
#> C2      C2             0  1.000000000     0.0000000           0.0000000
#> C3      C3             0  0.005414414     0.9945856           0.0306349
#> C4      C4             0  0.000000000     1.0000000           0.0000000
#> C5      C5             0  0.000000000     1.0000000           0.0000000
#> C6      C6             0  0.000000000     1.0000000           0.0000000
#>    tata_branch_assignment
#> C1                     C2
#> C2                     C2
#> C3                    C15
#> C4                    C15
#> C5                    C15
#> C6                    C15

Compare entropy across the shared early population and the simulated branches to locate cells with diffuse terminal-state predictions:

PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_branch_entropy"
)

PlotReduction(
    sce_tata,
    reduction = "UMAP",
    colour_by = "tata_branch_assignment"
)

Separate panels show the probability assigned to each terminal branch:

plot_tata_branch_probabilities(
    tata,
    dimred = "UMAP",
    ncol = 2
)

Selecting cells for a branch

select_tata_branch_cells() selects cells for a terminal branch. The "soft" mode applies a probability threshold; the "hard" mode selects cells with the corresponding branch assignment. Their relative sizes depend on the threshold and assignments. At a threshold of 0.60, cells with diffuse probabilities may be excluded even if they precede the branch.

branch_one <- tata$terminal_clusters[1]

soft_selection <- select_tata_branch_cells(
    tata,
    branch = branch_one,
    selection_mode = "soft",
    probability_threshold = 0.60
)

hard_selection <- select_tata_branch_cells(
    tata,
    branch = branch_one,
    selection_mode = "hard"
)

c(
    branch = branch_one,
    soft = sum(soft_selection$selected),
    hard = sum(hard_selection$selected)
)
#> branch   soft   hard 
#>  "C16"  "294"  "294"
plot_tata_branch_selection(
    tata,
    branch = branch_one,
    dimred = "UMAP",
    selection_mode = "soft",
    probability_threshold = 0.60
)

plot_tata_terminal_probabilities() displays terminal-state probabilities for individual cells, selecting ambiguous cells by default. Similar bar heights indicate diffuse predictions; one dominant bar indicates a concentrated prediction.

plot_tata_terminal_probabilities(
    tata,
    n_cells = 4
)

The simulation provides labels for comparing inferred branch assignments. Inspect agreement separately for the shared founder population and the branch-specific states, whose labels have different interpretations.

plot_tata_branch_confusion(
    tata,
    truth_col = "simulated_branch",
    predicted_col = "tata_branch_assignment",
    normalize = "truth"
)

Compact result overview

plot_tata_results() returns a named list of ggplot2 summary plots. Inspect or customise individual plots using their list names.

tata_overview <- plot_tata_results(
    tata,
    dimred = "UMAP",
    n_probability_cells = 4
)

names(tata_overview)
#> [1] "pseudotime"           "entropy"              "time_calibration"    
#> [4] "edge_causality"       "branch_probabilities" "representative_cells"
tata_overview$pseudotime

tata_overview$entropy

tata_overview$time_calibration

tata_overview$branch_probabilities

Practical guidance

  • Inspect neighbourhood size. Check graph connectivity, cluster sizes and time distributions across plausible values of k.
  • Inspect directional evidence. Increasing direction_threshold leaves more retained connections bidirectional. The threshold is a score cutoff, not a calibrated confidence level.
  • Check temporal disagreement. Use anticausal_edge_mass() with a tolerance suited to the sampling schedule, then inspect the flagged edges.
  • Retain probability information. Branch assignments simplify uncertain predictions. Downstream inference still requires an appropriate model and independent biological replicates.
  • Use comparable time labels. Confirm that units and the origin of time are consistent across conditions before combining their observations.

Further reading

  • Robinson MD, McCarthy DJ, Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. https://bioconductor.org/packages/edgeR/
  • Haghverdi L, Lun ATL, Morgan MD, Marioni JC (2018). Batch effects in single-cell RNA-sequencing data are corrected by matching mutual nearest neighbors. Nature Biotechnology 36, 421–427. Implemented in batchelor.
  • Korsunsky I et al. (2019). Fast, sensitive and accurate integration of single-cell data with Harmony. Nature Methods 16, 1289–1296.
  • Luecken MD et al. (2022). Benchmarking atlas-level data integration in single-cell genomics (scIB). Nature Methods 19, 41–50. The integration metrics in BenchmarkIntegration() are inspired by this work and use independent implementations.
  • Wolf FA et al. (2019). PAGA: graph abstraction reconciles clustering with trajectory inference. Genome Biology 20, 59. TATA’s abstract cluster graph follows this design.
  • Setty M et al. (2019). Characterization of cell fate probabilities in single-cell data with Palantir. Nature Biotechnology 37, 451–460. The source of the absorbing-Markov branch-probability layer.
  • McInnes L, Healy J, Melville J (2018). UMAP: Uniform Manifold Approximation and Projection for dimension reduction. https://arxiv.org/abs/1802.03426. Available in R via uwot.
  • Amezquita RA et al. (2020). Orchestrating single-cell analysis with Bioconductor. Nature Methods 17, 137–145. https://bioconductor.org/books/release/OSCA/

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] scProkaR_0.99.2             BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] rlang_1.3.0               magrittr_2.0.5           
#>  [3] RcppAnnoy_0.0.23          otel_0.2.0               
#>  [5] compiler_4.6.1            DelayedMatrixStats_1.35.0
#>  [7] vctrs_0.7.3               pkgconfig_2.0.3          
#>  [9] fastmap_1.2.0             XVector_0.53.0           
#> [11] scuttle_1.23.2            labeling_0.4.3           
#> [13] rmarkdown_2.32            xfun_0.61                
#> [15] bluster_1.23.1            cachem_1.1.0             
#> [17] beachmat_2.29.3           jsonlite_2.0.0           
#> [19] DelayedArray_0.39.8       BiocParallel_1.47.0      
#> [21] irlba_2.4.1               parallel_4.6.1           
#> [23] cluster_2.1.8.3           R6_2.6.1                 
#> [25] bslib_0.12.0              RColorBrewer_1.1-3       
#> [27] limma_3.99.1              parallelly_1.48.0        
#> [29] jquerylib_0.1.4           Rcpp_1.1.2               
#> [31] knitr_1.52                future.apply_1.20.2      
#> [33] FNN_1.1.4.1               splines_4.6.1            
#> [35] batchelor_1.29.0          Matrix_1.7-6             
#> [37] igraph_2.3.4              tidyselect_1.2.1         
#> [39] abind_1.4-8               yaml_2.3.12              
#> [41] codetools_0.2-20          listenv_1.1.0            
#> [43] lattice_0.23-1            tibble_3.3.1             
#> [45] withr_3.0.3               S7_0.2.2                 
#> [47] evaluate_1.0.5            future_1.76.0            
#> [49] pillar_1.11.1             BiocManager_1.30.27      
#> [51] sp_2.2-3                  sparseMatrixStats_1.25.0 
#> [53] scales_1.4.0              scrapper_1.7.9           
#> [55] globals_0.19.1            RhpcBLASctl_0.23-42      
#> [57] glue_1.8.1                maketools_1.3.2          
#> [59] tools_4.6.1               BiocNeighbors_2.7.3      
#> [61] sys_3.4.3                 ScaledMatrix_1.21.0      
#> [63] RSpectra_0.16-2           locfit_1.5-9.12          
#> [65] RANN_2.6.3                buildtools_1.0.0         
#> [67] dotCall64_1.2             cowplot_1.2.0            
#> [69] grid_4.6.1                edgeR_4.99.6             
#> [71] BiocSingular_1.29.2       cli_3.6.6                
#> [73] rsvd_1.0.5                spam_2.11-4              
#> [75] viridisLite_0.4.3         S4Arrays_1.13.2          
#> [77] dplyr_1.2.1               uwot_0.2.5               
#> [79] ResidualMatrix_1.23.0     gtable_0.3.6             
#> [81] sass_0.4.10               digest_0.6.39            
#> [83] progressr_1.0.0           SparseArray_1.13.4       
#> [85] SeuratObject_5.4.0        farver_2.1.2             
#> [87] htmltools_0.5.9           lifecycle_1.0.5          
#> [89] statmod_1.5.2             harmony_2.0.5