---
title: "Time Aware Trajectory Analysis (TATA) with scProkaR"
author:
- name: Ning Liu
affiliation: University of Adelaide
- name: Jinxin Zhao
affiliation: Monash University
date: "`r Sys.Date()`"
package: scProkaR
output:
BiocStyle::html_document:
toc: true
toc_float: true
toc_depth: 2
abstract: |
Time Aware Trajectory Analysis (TATA) combines expression neighbourhoods
and sampling times to construct a graph of cell states. This vignette
demonstrates the individual functions and the run_tata() wrapper using a
simulated branching response to three drugs. It compares time-informed
and topology-only runs, then shows how to inspect pseudotime, terminal-state
probabilities, branch assignments and gene-expression trends.
Comparisons with the supplied time labels assess consistency with those
labels.
vignette: >
%\VignetteIndexEntry{2. Time Aware Trajectory Analysis (TATA)}
%\VignetteEngine{knitr::rmarkdown}
%\VignetteEncoding{UTF-8}
---
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>",
crop = NULL,
fig.width = 6,
fig.height = 4.5,
dpi = 96,
fig.retina = 1
)
```
# 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 `r BiocStyle::Biocpkg("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 `r BiocStyle::Biocpkg("scater")`,
`r BiocStyle::Biocpkg("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
```{r install, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("scProkaR")
```
# Setting up
```{r load-packages, message = FALSE}
library(scProkaR)
library(SingleCellExperiment)
library(SummarizedExperiment)
```
TATA uses `r BiocStyle::CRANpkg("FNN")` for nearest-neighbour search when
available, followed by `r BiocStyle::Biocpkg("BiocNeighbors")` as a fallback.
If neither is installed, it uses `stats::dist()`. The distance-matrix fallback
requires substantially more memory for large datasets.
```{r knn-backends}
c(
FNN = requireNamespace("FNN", quietly = TRUE),
BiocNeighbors = requireNamespace("BiocNeighbors", quietly = 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.
```{r simulate-multidrug}
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
```
All three conditions are represented at each sampled time, with varying
cell numbers:
```{r simulation-design}
with(
as.data.frame(colData(sce)),
table(condition, timepoint)
)
```
The stored UMAP views show sampling time, treatment and simulated state.
Compare the locations of the earliest cells with the three response branches.
```{r plot-simulation-time}
PlotReduction(sce, reduction = "UMAP", colour_by = "timepoint_factor")
```
```{r plot-simulation-condition}
PlotReduction(sce, reduction = "UMAP", colour_by = "condition")
```
```{r plot-simulation-branch}
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.
```{r build-knn}
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
)
```
```{r knn-engine}
knn$knn_engine
```
## 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.
```{r cluster-states}
sce$tata_cluster_step <- cluster_graph_states(
cell_graph = knn$graph,
method = "louvain"
)
table(sce$tata_cluster_step)
```
```{r plot-cluster-states}
PlotReduction(sce, reduction = "UMAP", colour_by = "tata_cluster_step")
```
When `r BiocStyle::CRANpkg("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.
```{r cluster-silhouette}
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.")
}
```
## 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`.
```{r topology-weights}
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
)
```
`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.
```{r time-weights}
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
)
```
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.
```{r plot-edge-score-scatter}
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.
```{r build-graph-manual}
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
)
```
The embedding layout places each cluster node at the centre of its cells,
allowing graph connections to be compared with the cell distribution.
```{r plot-manual-graph-embedded}
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.
```{r plot-manual-graph-abstract}
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.
```{r manual-pseudotime}
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
summary(pt_manual$cell_pseudotime_scaled)
```
Inspect the median sampling time of the selected root cluster:
```{r manual-root-check}
tapply(
sce$timepoint,
sce$tata_cluster_step,
median
)[pt_manual$root_cluster]
```
# 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.
```{r run-tata}
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)
```
The updated object's annotations and embeddings can be inspected as follows:
```{r inspect-coldata}
sce_tata <- tata$sce
sce_tata$tata_cluster <- factor(sce_tata$tata_cluster)
grep("^tata_", colnames(colData(sce_tata)), value = TRUE)
reducedDimNames(sce_tata)
```
```{r plot-tata-clusters}
PlotReduction(sce_tata, reduction = "UMAP", colour_by = "tata_cluster")
```
```{r plot-tata-pseudotime}
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.
```{r plot-tata-graph-embedded}
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"
)
```
```{r plot-tata-graph-abstract}
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.
```{r plot-trajectory-all-paths}
plot_tata_trajectory_embedding(
tata,
dimred = "UMAP",
colour_by = "timepoint",
smooth_paths = TRUE
)
```
```{r plot-trajectory-top-paths}
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.
```{r run-tata-topology-only}
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.
```{r orient-pseudotime}
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.
```{r benchmark-local-order}
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)
```
**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.
```{r benchmark-adjacent-auc}
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)
```
Inspect the per-pair table to identify intervals with weaker temporal
separation:
```{r benchmark-auc-pairs}
head(auc_tata$pair_table)
```
**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.
```{r benchmark-reachability}
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
head(reach_tata$pair_table)
```
**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.
```{r benchmark-anticausal}
anticausal_tata <- anticausal_edge_mass(
cluster_graph = tata$cluster_graph,
clusters = tata$sce$tata_cluster,
timepoint = tata$sce$timepoint,
tolerance = 24
)
anticausal_tata$summary
```
```{r benchmark-anticausal-edges}
head(
anticausal_tata$edge_table[
order(anticausal_tata$edge_table$time_delta),
c("from", "to", "weight", "time_delta", "temporal_class")
]
)
```
**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.
```{r benchmark-branch-onset}
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)
```
```{r benchmark-onset-table}
onset_tata$onset_table
```
## 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.
```{r plot-calibration-tata}
plot_tata_time_calibration(
tata,
branch_col = "condition",
summary = "median_iqr"
)
```
```{r plot-calibration-topology}
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.
```{r plot-edge-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`.
```{r branch-probabilities}
branch_result <- compute_tata_branch_probabilities(
tata_result = tata,
expected_branches = 3
)
branch_result$terminal_clusters
```
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.
```{r cluster-branch-table}
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")
]
)
```
Compare entropy across the shared early population and the simulated
branches to locate cells with diffuse terminal-state predictions:
```{r plot-branch-entropy}
PlotReduction(
sce_tata,
reduction = "UMAP",
colour_by = "tata_branch_entropy"
)
```
```{r plot-branch-assignment}
PlotReduction(
sce_tata,
reduction = "UMAP",
colour_by = "tata_branch_assignment"
)
```
Separate panels show the probability assigned to each terminal branch:
```{r plot-branch-probabilities}
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.
```{r select-branch-cells}
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)
)
```
```{r plot-branch-selection}
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.
```{r plot-terminal-probabilities}
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.
```{r plot-branch-confusion}
plot_tata_branch_confusion(
tata,
truth_col = "simulated_branch",
predicted_col = "tata_branch_assignment",
normalize = "truth"
)
```
# Branch-specific gene trends
`compute_tata_gene_trends()` estimates expression trends along pseudotime
for each gene and branch, weighting cells by their branch probabilities.
A cell may therefore contribute to more than one branch.
We pick genes from the simulation's branch-specific and late feature modules
recorded in `rowData()`.
```{r select-trend-genes}
feature_modules <- rowData(sce_tata)$feature_module
trend_genes <- unique(c(
rownames(sce_tata)[feature_modules == "branch_A"][1:2],
rownames(sce_tata)[feature_modules == "branch_B"][1:2],
rownames(sce_tata)[feature_modules == "branch_C"][1:2],
rownames(sce_tata)[feature_modules == "late"][1:2]
))
trend_genes
```
```{r compute-gene-trends}
tata_trends <- compute_tata_gene_trends(
tata,
genes = trend_genes,
grid_length = 60
)
tata_trends$branches
head(tata_trends$trend_table)
```
Compare the fitted curves with the simulated feature modules to assess
whether branch-specific and shared late responses are represented.
```{r plot-gene-trends}
plot_tata_gene_trends(
tata_trends,
ncol = 3
)
```
`cluster_tata_gene_trends()` groups gene-branch curves by shape after
standardisation. It uses `stats::kmeans()`, so set a seed when clustering
for reproducible results.
```{r cluster-gene-trends}
set.seed(101)
trend_clusters <- cluster_tata_gene_trends(
tata_trends,
n_clusters = 3
)
trend_clusters
```
```{r plot-gene-trend-heatmap}
set.seed(101)
plot_tata_gene_trend_heatmap(
tata_trends,
scale_rows = TRUE,
n_clusters = 3
)
```
# Compact result overview
`plot_tata_results()` returns a named list of `ggplot2` summary plots.
Inspect or customise individual plots using their list names.
```{r tata-overview}
tata_overview <- plot_tata_results(
tata,
dimred = "UMAP",
n_probability_cells = 4
)
names(tata_overview)
```
```{r plot-overview-pseudotime}
tata_overview$pseudotime
```
```{r plot-overview-entropy}
tata_overview$entropy
```
```{r plot-overview-calibration}
tata_overview$time_calibration
```
```{r plot-overview-branch-probabilities}
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.
* 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
`r BiocStyle::Biocpkg("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.
. Available in R via
`r BiocStyle::CRANpkg("uwot")`.
* Amezquita RA *et al.* (2020). Orchestrating single-cell analysis with
Bioconductor. *Nature Methods* 17, 137--145.
# Session information
```{r session-info}
sessionInfo()
```