logo



1 A bit of theory

1.1 Why ?

Because …

  • … cells in a defined environment (body, organoid, …) are :
    1. possibly of different kinds (cell types)
    2. able to temporaly evolve from a state to another (ex : stem to differentiated)
    3. not (all) synchronous …
  • … and a single cell experiment is an observation of numerous cells of various types in a common context :
    1. Various types : states ?
    2. Common context : lineage ?
    3. Numerous : internal variation in each state ?
logo

1.2 How ?

1.2.1 Main steps

The basic principle can be divided into 3 main steps :

  • Inferring a schematic representation of the transitional behavior of some cells : deciphering a topology from the data
    • Which data ?
      • Cell by feature matrix ? (spoiler : NO)
      • A reduced space ?
        • Which method ?
        • Number of dimensions ?
    • Which topology ?
      • Linear ?
      • Circular ?
      • Poly-cyclic (ie, nested cycles)?
      • Bifurcations ?
      • Multi-furcations ?
      • Simple / complex tree ?
      • A combination of some of all above ?
      • Single one / multiple disconnected ?
      • Oriented or not (rooting) ?
logo


  • Positioning each cell in the defined topology : defining cells pseudotimes
    • Directly / Proxy (clusters) ?
    • Direction : single / multiple roots ? Which one(s) ?

1.2.2 Inferring a topology

  • Most methods rely on identifying a minimal spanning tree (MST), with specific variations (to allow cycles, by example) and tweaks (branch pruning).

  • The key idea in MST is to define the simplest path that connects all cells (vertices) using ponderated edges (their distance in the reduced space)

logo


1.2.3 Defining cells pseudotime

Once a topology has been created, pseudotime of a cell consists in the rank of the cell projection to its closest topology backbone point

logo


1.2.4 Historically

  • First methods (circa 2014)
    • Strong priors :
      • Reduced space (PCA)
      • Clustering results (~ cell states) or inferred cell types
      • Expected topology type
        • Single one
        • Simple : linear (direct or cyclic), bifurcation tree
  • Evolution :
    • Lower priors (just space +/- clusters +/- root)
      • No mandatory input topology type
    • Better-fited types of reduced space (tSNE, uMAP, LLE, Diffusion maps, Poincarre maps, GPLVM, …)

1.2.5 Methods : Welcome to the dynverse !

Inferring any kind of topology is a complex task :

“It is perhaps surprising that of the 59 methods in existence today [at the time of the publication], almost all methods have a unique combination of [these prior + reduced space + topology inferring method] characteristics” (dynverse publication, 2018)

  • Multiple (dozens) of methods are available, either through R or python packages, or independent binaries.

  • Wouter Saelens and Robrecht Cannoodt did a tremendous work setting up the dynverse, a framework that :

    • re-implements/wraps 55 existing methods (dynmethods)
    • harmonizes them in a common way of setting inputs and parameters, getting results (dyno, dynparam)
    • is compatible with any novel method for new developers (dynwrap)
    • contains metrics to compare methods
    • allows a comparison of numerous methods (dyneval)
    • is available as reference map to help select the best method based on user context (dynguidelines)
dynguidelines::guidelines_shiny()
logo


2 Practice

2.1 Prerequisites

2.1.1 This Rmarkdown

Copy this Rmd file (and its dependences) from the common ‘Courses’ dir to your working directory :

## Working directory
user.id <- Sys.getenv('USER')
root.dir <- '/shared/projects/sincellte_2022'
script.dir <- paste(c(root.dir, 'Courses', 'Secondary_analyses',
                      'scripts'), collapse = '/')
work.dir <- paste(c(root.dir, user.id, 'Courses', 'Secondary_analysis',
  'trajectory_analysis'), collapse = '/')
dir.create(path = paste(c(work.dir, 'trajectory_rmd_resources'), collapse = '/'), 
           recursive = TRUE, showWarnings = FALSE)
## Copying the Rmd
file.copy(from = paste0(script.dir, '/Trajectory_TD_deliverable.Rmd'),
          to = work.dir, overwrite = TRUE)
## Copying the dependencies
file.copy(from = list.files(
  path = paste0(script.dir, '/trajectory_rmd_resources/'), full.names = TRUE), 
  to = paste0(work.dir, '/trajectory_rmd_resources/'), overwrite = TRUE)
## Moving to your working directory
setwd(work.dir)
## If you want to visualize your files in Rstudio, you can click the little greyish/white arrow on the right of the path written at the top of your Rstudio Console.

2.1.2 The dataset

  • This practice is inspired by this session from this Galaxy training
  • The dataset you will use :
    • Comes from [Bacon et al, Frontiers in Immunology, 2018].
    • Is available as a scanpy h5ad object (tool-specific format based upon the HDF5 standard)
    • Is hosted at Zenodo #4726927
    • Consists in a drop-seq preparation of thymic cells in differenciation from 7 neonatal mice, selected by mass flow cytometry using 22 immunologic markers corresponding to multiple T-cell types and subtypes (see embedded Excel file)
    • For this practical session, will focus on these subtypes :
      • Double-negative (DN) : Early immature cells
      • Double-positive Module 1 (DP-M1) : Intermediate state 1
      • Double-positive Module 2 (DP-M2) : Intermediate state 2
      • Double-positive Module 3 (DP-M3) : Intermediate state 3
      • Double-positive Module 4 (DP-M4) : Intermediate state 4
      • Double-positive Late (DP-L) : Immature late state
      • Mature T-cells (T-mat)

Download markers_celltypes.xlsx

2.1.3 The environment

The R packages needed for this tutorial are :

Additional packages needed to build this HTML report :

  • CRAN :
    • gifski To generate GIF animations of markers expression
    • knitr To set the Rmarkdown hook to compute CPU time of chunks
    • Hmisc To set the button-hidden TOC

2.2 Preprocessing steps

This is just a description of the preparation I performed on the data, building the the object you will work on. This is not part of the practical session.

Unfold

2.2.1 Setting our variables

## Remote input file
remote.h5ad.file <- 'https://zenodo.org/record/4726927/files/Final_cell_annotated_object.h5ad?download=1'
## Working directory
user.id <- Sys.getenv('USER')
root.dir <- '/shared/projects/sincellte_2022'
data.dir <- paste(c(root.dir, 'Courses', 'Secondary_analyses', 'input', 'trajectory_analysis'), collapse = '/')
work.dir <- paste(c(root.dir, user.id, 'Courses', 'Secondary_analysis', 'trajectory_analysis'), collapse = '/')
## Seurat assay name
assay <- 'RNA'
## Contaminating cell types
contam.types <- c('Macrophages', 'RBC')
## RNG seed (for irlba PCA)
my.seed <- 1337L
## HVGs to keep
nfeatures <- 3000
## 1st reduction method
reduction <- 'pca'
## Number of PCs to produce
pca.dims <- 100
## Number of PCs to retain for 2nd reduction
maxdim <- 30
## Use scaled data instead of unscaled
use.scaled <- TRUE
## Cell type to use as root to orient the trajectory
root.group <- 'DN'

2.2.2 Loading source data

  • First, the dataset was retrieved from Zenodo
dir.create(path = work.dir, recursive = TRUE, showWarnings = FALSE)
local.root.file <- paste0(work.dir, '/Original_scanpy_object.h5')
local.h5ad.file <- paste0(local.root.file, 'ad')
download.file(url = h5ad.file, destfile = local.h5ad.file,
              method = 'wget', quiet = TRUE)
SeuratDisk::Convert(source = local.h5ad.file, dest = 'h5seurat',
                    assay = assay, verbose = FALSE, overwrite = TRUE)
  • The converted object was loaded into R as a SeuratObject using the SeuratDisk::LoadH5Seurat() function. All optional slots (neighbors, graphs, reductions, images), won’t be loaded as we want to reprocess them.
local.h5s.file <- paste0(local.root.file, 'seurat')
sobj <- SeuratDisk::LoadH5Seurat(file = local.h5s.file, assays = assay,
                                 reductions = FALSE, graphs = FALSE,
                                 neighbors = FALSE, images = FALSE)

2.2.3 Prepping the object to our context

  • The downloaded object contains cell-type annotation from the original analysis. From the publication, we learnt that two cell types are not related to T cells, thus considered as contaminants : macrophages and red blood cells (RBC). As their presence in the analysis space can modify the shape of the cells position, thus impact trajectory analysis, we decided to remove them :
sobj <- sobj[, !(sobj$cell_type %in% contam.types)]
  • We then follow the Seurat4 canonical pipeline, applying log-normalization to the count data …
sobj <- Seurat::NormalizeData(object = sobj, assay = assay)
  • … then scaling. As this object is a merge of 7 experiments (labelled as ‘batch’ in the object metadata), we want to regress this effect.
sobj <- Seurat::ScaleData(object = sobj, assay = assay,
                          vars.to.regress = c('batch'))
  • … then identify HVGs. As we focus on subtypes of a single, main cell type (T-cells) in a differentiation context, we anticipate that we may have limited variation of gene expression (at least compared to, by example, multiple cell types of a differentiated organ). Thus, we increase the default number of features to output from 2,000 to 3,000.
sobj <- Seurat::FindVariableFeatures(object = sobj, assay = assay,
                                     verbose = FALSE, nfeatures = nfeatures)
  • … then reducing dimensions (through PCA). We request 100 PCs as output.
sobj <- Seurat::RunPCA(object = sobj, assay = assay, npcs = pca.dims,
                       reduction.name = reduction, seed.use = my.seed,
                       verbose = FALSE)
  • … then building additional reductions for exploration / visualization (t-SNE, uMAP). After manual exploration, we decided to retain 30 PCs as input.
## TSNE
reduction2t <- paste0(reduction, '.', maxdim, '_tsne')
sobj <- Seurat::RunTSNE(object = sobj, dims = 1:maxdim, assay = assay,
                        reduction = reduction, seed.use = my.seed,
                        verbose = FALSE, reduction.name = reduction2t)
## uMAP
reduction2u <- paste0(reduction, '.', maxdim, '_umap')
sobj <- Seurat::RunUMAP(object = sobj, dims = 1:maxdim, assay = assay,
                        reduction = reduction, seed.use = my.seed,
                        verbose = FALSE, reduction.name = reduction2u)
Seurat::DimPlot(object = sobj, reduction = reduction2u, group.by = 'cell_type',
                seed = my.seed)
  • Lastly, we save the processed object (as a RDS file).
saveRDS(object = sobj, file = paste0(work.dir,'/Thymus_neonat_mm_Seurat_ENSEMBL.RDS'),
        compress = 'bzip2')

2.3 Trajectory analysis using TInGa

We will use a tool named TInGa to perform this analysis, due to its convenience to our educational purpose :

  • It’s in full R
  • It’s very fast
  • Its parameters are limited (as compared to other renowned tools like Monocle3)
  • It is compatible with most topologies (see here)
  • It’s integrated in the versatile dynverse (as a plugin)

2.3.1 Setting parameters

## Working directory
user.id <- Sys.getenv('USER')
root.dir <- '/shared/projects/sincellte_2022'
data.dir <- paste(c(root.dir, 'Courses', 'Secondary_analyses', 'input',
                    'trajectory_analysis'), collapse = '/')
work.dir <- paste(c(root.dir, user.id, 'Courses', 'Secondary_analysis',
                    'trajectory_analysis'), collapse = '/')
if(!dir.exists(work.dir)) dir.create(path = work.dir, recursive = TRUE)
## Local input file
## The file created if you performed the preprocessing by yourself
# local.Seurat.file <- paste0(work.dir, '/Thymus_neonat_mm_Seurat_ENSEMBL.RDS')
## The preprocessed file I prepared for you
local.Seurat.file <- paste0(data.dir, '/Thymus_neonat_mm_Seurat_ENSEMBL.RDS')
## Seurat assay name
assay <- 'RNA'
## RNG seed (for irlba PCA)
my.seed <- 1337L
## 1st reduction method
reduction <- 'pca'
## 2nd reduction method
reduction2 <- 'pca.30_umap'
## Use scaled data instead of unscaled
use.scaled <- TRUE
## Cell type to use as root to orient the trajectory
root.group <- 'DN'
## Marker genes to plot
trans.markers <- c('Il2ra', 'Npm1', 'Nusap1', 'Ube2c',
                   'Pclaf','Cd8b1', 'Itm2a', 'Cd52', 'H2-D1')

2.3.2 Loading source data

We will simply load the pre-processed Seurat object, contained in a RDS archive, as input :

sobj <- readRDS(file = local.Seurat.file)

Time for this code chunk to run : 1.305 s

2.3.3 Performing trajectory analysis

2.3.3.1 The reduced space

First, to impregnate us with the cells topology on which TInGa will work, we can display our 2D dimension reduction (uMAP or t-SNE) :

## With batches
db <- Seurat::DimPlot(object = sobj, reduction = reduction2, group.by = 'batch',
                      pt.size = 2) + Seurat::DarkTheme()
## With cell types
dct <- Seurat::DimPlot(object = sobj, reduction = reduction2, group.by = 'cell_type',
                       pt.size = 2) + Seurat::DarkTheme()
## With sex
ds <- Seurat::DimPlot(object = sobj, reduction = reduction2, group.by = 'sex',
                      pt.size = 2) + Seurat::DarkTheme()
print(patchwork::wrap_plots(list(db, dct, ds)))

2.3.3.2 Extracting data

Using TInGa first requires to build a specific object from the count and normalized matrices, that we will extract from our Seurat object :

## Matrices are required to be formatted as [cells] x [features], thus a transposition is required
## Normalized matrix
my.exp <- t(as.matrix(if(use.scaled) sobj@assays[[assay]]@scale.data
                      else sobj@assays[[assay]]@data)) 
### Ordering to keep synch in the [use.scaled] context
my.exp <- my.exp[, order(colnames(my.exp))]
## Count matrix
my.counts <- t(as.matrix(if(use.scaled) sobj@assays[[assay]]@counts[rownames(sobj@assays[[assay]]@counts)%in% colnames(my.exp),]
                         else sobj@assays[[assay]]@counts))
### Ordering to keep synch in the [use.scaled] context
my.counts <- my.counts[, order(colnames(my.counts))]

2.3.3.3 Building a dynwrap object

We can now build our object (a dynwrap object, which actually consists in a list with mandatory entries) :

my.data <- dynwrap::wrap_expression(
  id = sobj@project.name,
  expression = my.exp,
  counts = my.counts
)
is(my.data)
## [1] "dynwrap::with_expression"

Time for this code chunk to run : 2.55 s

2.3.3.4 Import prior

We can import our reduced space (t-SNE or uMAP) as a prior, which will be used by TInGa to build the trajectory :

my.data <- dynwrap::add_prior_information(
  dataset = my.data,
  dimred = sobj@reductions[[reduction2]]@cell.embeddings
)

2.3.3.5 Running TInGa

Now, we can run TInGa !

## Performing TInGa analysis
set.seed(my.seed)
my.traj <- dynwrap::infer_trajectory(dataset = my.data, method = TInGa::gng_param2(),
                                     seed = my.seed, give_priors = 'dimred',
                                     parameters = list(max_nodes = 8, lambda = 200))
## Adding pseudotime
my.traj <- dynwrap::add_pseudotime(trajectory = my.traj)
## Warning in calculate_pseudotime(trajectory): Trajectory is not rooted. Add a
## root to the trajectory using dynwrap::add_root(). This will result in an error
## in future releases.
## root cell or milestone not provided, trying first outgoing milestone_id
## Using '1' as root
## Plotting with pseudotime
dynplot::plot_dimred(trajectory = my.traj, label_milestones = TRUE,
                     color_cells = 'pseudotime') + Seurat::DarkTheme()
## Registered S3 method overwritten by 'cli':
##   method     from         
##   print.boxx spatstat.geom