Transcriptomic age from single-cell RNA-seq (pseudobulk) with tAge

Overview

The clocks are trained on bulk RNA-seq. To apply them to single-cell data, aggregate cells into pseudobulk samples, then run the standard bulk pipeline. tAge aggregates by coverage — each pseudobulk sample pools enough cells to reach a target read depth — which is more robust than a fixed cells-per-sample cut-off.

This walkthrough uses Tabula Muris Senis (GSM4505404). Requires Seurat (and SeuratDisk or reticulate+anndata to read .h5ad).

1. Load single-cell data into a Seurat object

library(tAge)
library(Seurat)
library(reticulate)

Sys.setenv(RETICULATE_PYTHON = ".venv/bin/python")

# Read an .h5ad via anndata and build a Seurat object
ad    <- import("anndata")
adata <- ad$read_h5ad("GSM4505404_tabula-muris-senis-droplet-official-raw-obj.h5ad")

counts <- Matrix::t(adata$X)
rownames(counts) <- adata$var_names$to_list()
colnames(counts) <- adata$obs_names$to_list()

seurat_obj <- CreateSeuratObject(
  counts    = counts,
  meta.data = as.data.frame(adata$obs),
  project   = "TabulaMurisSenis"
)

tAge also provides load_h5ad_to_seurat() as a convenience wrapper.

2. (Optional) subset

selected_tissues <- c("Tongue", "Kidney", "Liver", "Heart")
seurat_subset <- subset(seurat_obj, subset = tissue %in% selected_tissues)

3. Coverage-based pseudobulk aggregation

Aggregate within each unique combination of grouping columns, pooling cells until each pseudobulk sample reaches coverage_threshold reads.

eset <- aggregate_on_obs_columns(
  seurat_obj       = seurat_subset,
  obs_column_names = c("age", "mouse.id", "tissue"),
  coverage_threshold = 1e7,
  verbose = TRUE
)

# Simpler variant: aggregate the whole object by coverage only
# eset <- aggregate_pseudobulk(seurat_subset, coverage_threshold = 1e6)

4. Outlier samples

flag_outliers() flags samples whose log-CPM profile is far from the median profile of their tissue (robust z of 1 - r above 5); look at them, then remove them explicitly.

eset <- flag_outliers(eset, split_by = "tissue")
Biobase::pData(eset)[eset$tage_outlier, c("tissue", "tage_outlier_z")]
eset_clean <- eset[, !eset$tage_outlier]
plot_eset_density(eset_clean)

5. Predict per tissue

tAge_by_group() runs preprocessing and prediction separately within each level of split_by (recommended for multi-tissue data).

# Download clocks (see the bulk vignette / ?download_clocks)
clocks <- download_clocks(
  list_clocks(type = "EN", outcome = "Chronological",
              species = "Multispecies", tissue = "Multi-Tissue"),
  dest_dir = "clocks"
)
model_paths <- list(
  scaled_diff = clocks$path[clocks$scaling == "Scaled"],
  yugene_diff = clocks$path[clocks$scaling == "YuGene"]
)

results <- tAge_by_group(
  eset        = eset_clean,
  split_by    = "tissue",
  model_paths = model_paths,
  species     = "mouse",
  mode        = "EN"
)

6. Compare age groups

The statistics are the same as for bulk data — pseudobulk samples are just samples. Here the comparison is old versus young within each tissue, so split_by fits one model per tissue.

stats <- tage_compare_groups(
  results,
  value_columns   = c("scaled_diff_EN_tAge", "yugene_diff_EN_tAge"),
  group_column    = "age",
  reference_group = "3m",           # the reference group (youngest age you kept)
  split_by        = "tissue",
  p_adjust        = "BH",
  p_adjust_scope  = "across_columns"
)

stats[, c("value_column", "split", "group2", "estimate",
          "ci_low", "ci_high", "p_adjusted", "label")]

Age is also a number here, so the trend across the whole range is a more natural summary than a set of pairwise contrasts:

results$age_numeric <- as.numeric(gsub("m", "", results$age))

tage_regress_continuous(
  results,
  value_columns = "scaled_diff_EN_tAge",
  predictor     = "age_numeric",
  split_by      = "tissue"
)

A caveat specific to pseudobulk: several metacells from the same animal are not independent. Pooling one metacell per (animal × tissue), as in step 3, keeps the rows independent. If you build several metacells per animal, add the animal as a covariate or aggregate the predictions before testing.

7. Figures

clocks_meta <- data.frame(
  filename = c("scaled_diff_EN_tAge", "yugene_diff_EN_tAge"),
  name     = c("EN Scaled", "EN YuGene"),
  outcome  = "Chronological",
  stringsAsFactors = FALSE
)

tage_clock_forest(
  results,
  clocks_meta     = clocks_meta,
  group_column    = "age",
  reference_group = "3m",
  split_by        = "tissue",
  title           = "Predicted age vs the 3-month group"
)

For one clock at a time, tage_boxplot() shows the pseudobulk samples themselves and annotates the brackets from the same engine:

tage_boxplot(
  results,
  x_var        = "age",
  y_var        = "scaled_diff_EN_tAge",
  subgroup_var = "tissue",
  theme_type   = "bw",
  title        = "Predicted tAge by age and tissue",
  xlab         = "Chronological age",
  ylab         = "Predicted tAge (months)"
)

With module clocks, tage_module_heatmap() shows which cellular subsystems carry the signal:

module_cols <- grep("^module_", names(results), value = TRUE)  # one column per module clock
tage_module_heatmap(
  results,
  module_columns  = module_cols,
  group_column    = "age",
  reference_group = "3m",
  split_by        = "tissue"
)

See the bulk vignette (vignette("tage-bulk")) for how to interpret the output and for the role of reference groups (centring).

Gladyshev Lab