15 Lab 8: Pseudotime analyses

Notes

The estimated time for this lab is around 1h30.

Aims

Understand the requirements for RNA velocity computation
Perform pseudotime and lineage analysis
Model and visualize gene expression along pseudotime
Compute RNA velocity and use it to orientate lineages

15.1 Process testis data in R

15.1.1 Import testis data from GSE112013 and pre-process it

We will work on a testis dataset, from the “The adult human testis transcriptional cell atlas” study (Guo et al., 2018). This dataset contains single-cell RNAseq data from human testis. It is available on the GEO database under the accession number GSE112013.

The pre-processing of the testis dataset has already been done for you. We have applied the different analysis steps we have learned over the past three days to the count matrix GSE112013_Combined_UMI_table.txt.gz. See the code below for the details of pre-processing steps.

👀 Data pre-processing 👇

Check the GEO page for supplementary processed data, and download the GSE112013_Combined_UMI_table.txt.gz table.

With this data, we will:

Import the dataset in R as a SingleCellExperiment object,
Remove doublets,
Filter genes,
Normalize counts,
Correct for batch effect,
Cluster cells.

library(SingleCellExperiment)
library(tidyverse)

# Download raw file
dir.create('Guo_testis')
download.file(
    'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE112nnn/GSE112013/suppl/GSE112013_Combined_UMI_table.txt.gz', 
    'Guo_testis/GSE112013_Combined_UMI_table.txt.gz'
)
system('gunzip Guo_testis/GSE112013_Combined_UMI_table.txt.gz')

# Create SingleCellExperiment object
x <- vroom::vroom('Guo_testis/GSE112013_Combined_UMI_table.txt')
counts <- as.matrix(x[, -1])
counts <- as(counts, 'dgCMatrix')
genes <- as.data.frame(x[, 1])
cells <- data.frame(cellid = colnames(x[, -1]))
testis <- SingleCellExperiment(
    assays = list(counts = counts), 
    colData = cells, 
    rowData = genes
)
testis$Barcode <- str_replace(testis$cellid, 'Donor.-', '') |> str_replace('-.', '')
testis <- testis[, !duplicated(testis$Barcode)]
testis$donor <- str_replace(testis$cellid, '-.*', '')
testis$replicate <- str_replace(testis$cellid, '.*-', '')
rownames(testis) <- rowData(testis)$Gene
set.seed(1000)

# Remove doublets
testis <- scDblFinder::scDblFinder(testis)
testis <- testis[, testis$scDblFinder.class == 'singlet']

# Recover human genomic, protein-coding gene informations
library(plyranges)
ah <- AnnotationHub::AnnotationHub()
AnnotationHub::query(ah, c('gene annotation', 'ensembl', '102', 'homo_sapiens', 'GRCh38'))
gtf <- AnnotationHub::query(ah, c('Homo_sapiens.GRCh38.102.chr.gtf'))[[1]]
genes <- gtf |> 
    filter(type == 'gene') |> 
    filter(gene_biotype == 'protein_coding') |> 
    filter(gene_source == 'ensembl_havana') 
genes <- genes[!duplicated(genes$gene_name)]

# Annotate genes in testis dataset and filter out non-relevant genes
testis <- testis[genes$gene_name[genes$gene_name %in% rownames(testis)], ]
rowRanges(testis) <- genes[match(rownames(testis), genes$gene_name)]
rowData(testis) <- rowData(testis)[, c('gene_name', 'gene_id')]
rownames(testis) <- scuttle::uniquifyFeatureNames(rowData(testis)$gene_id, rowData(testis)$gene_name)

# Get preliminary QCs per cell and per gene
testis <- scuttle::addPerCellQCMetrics(testis)
testis <- scuttle::addPerFeatureQCMetrics(testis)

# Filter out genes not expressed in at least 10 cells
testis <- testis[rowSums(assay(testis, 'counts') > 0) >= 10, ]

# Normalize counts using VST
cnts <- as(assay(testis, 'counts'), 'dgCMatrix')
colnames(cnts) <- testis$cellid
rownames(cnts) <- rownames(testis)
testis_vst <- sctransform::vst(cnts, return_cell_attr = TRUE)
corrected_cnts <- sctransform::correct(testis_vst)
assay(testis, 'corrected_counts', withDimnames = FALSE) <- corrected_cnts
assay(testis, 'logcounts', withDimnames = FALSE) <- log1p(corrected_cnts)

# Computing biological variance of each gene
testis_variance <- scran::modelGeneVar(testis)
HVGs <- scran::getTopHVGs(testis_variance, prop = 0.1)
rowData(testis)$isHVG <- rownames(testis) %in% HVGs

# Embedding dataset in PCA space, correcting for batch effect
mergedBatches <- batchelor::fastMNN(
    testis, 
    batch = testis$donor, 
    subset.row = HVGs, 
    BPPARAM = BiocParallel::MulticoreParam(workers = 12)
)
reducedDim(testis, 'corrected') <- reducedDim(mergedBatches, 'corrected')

# Embedding dataset in shared k-nearest neighbors graph for clustering 
graph <- scran::buildSNNGraph(testis, use.dimred = 'corrected')

# Cluster cells using Louvain community finding algorithm
testis_clust <- igraph::cluster_louvain(graph)$membership
table(testis_clust)
testis$cluster <- factor(testis_clust)

# Embedding dataset in TSNE space for visualization
set.seed(10)
testis <- scater::runTSNE(testis, dimred = 'corrected')

# Visualize embeddings
p <- cowplot::plot_grid(
    scater::plotReducedDim(testis, 'corrected', colour_by = 'donor'),
    scater::plotReducedDim(testis, 'corrected', colour_by = 'cluster'),
    scater::plotReducedDim(testis, 'TSNE', colour_by = 'donor'),
    scater::plotReducedDim(testis, 'TSNE', colour_by = 'cluster')
)
p
saveRDS(testis, '~/Share/Guo_testis.rds')

15.1.2 Annotate cells using HPA resources

The Human Protein Atlas (HPA) is a valuable resource for cell type annotations. It provides single-cell RNAseq data for a large number of cell types.

Question

Download HPA scRNAseq atlas (combined per cell type) from here, and import it in R as a SummarizedExperiment object.

Answer

library(SingleCellExperiment)

Loading required package: SummarizedExperiment

Loading required package: MatrixGenerics

Loading required package: matrixStats


Attaching package: 'MatrixGenerics'

The following objects are masked from 'package:matrixStats':

    colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse, colCounts, colCummaxs, colCummins, colCumprods, colCumsums, colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs, colMads, colMaxs, colMeans2,
    colMedians, colMins, colOrderStats, colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds, colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads, colWeightedMeans, colWeightedMedians,
    colWeightedSds, colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet, rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods, rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins, rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks, rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars, rowWeightedMads,
    rowWeightedMeans, rowWeightedMedians, rowWeightedSds, rowWeightedVars

Loading required package: GenomicRanges

Loading required package: stats4

Loading required package: BiocGenerics

Loading required package: generics


Attaching package: 'generics'

The following objects are masked from 'package:base':

    as.difftime, as.factor, as.ordered, intersect, is.element, setdiff, setequal, union


Attaching package: 'BiocGenerics'

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind, colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep, grepl, is.unsorted, lapply, Map, mapply, match, mget, order, paste,
    pmax, pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce, rownames, sapply, saveRDS, table, tapply, unique, unsplit, which.max, which.min

Loading required package: S4Vectors


Attaching package: 'S4Vectors'

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    expand.grid, I, unname

Loading required package: IRanges

Loading required package: Seqinfo

Loading required package: Biobase

Welcome to Bioconductor

    Vignettes contain introductory material; view with 'browseVignettes()'. To cite Bioconductor, see 'citation("Biobase")', and for packages 'citation("pkgname")'.


Attaching package: 'Biobase'

The following object is masked from 'package:MatrixGenerics':

    rowMedians

The following objects are masked from 'package:matrixStats':

    anyMissing, rowMedians

library(tidyverse)

── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.4     ✔ readr     2.1.5
✔ forcats   1.0.1     ✔ stringr   1.6.0
✔ ggplot2   4.0.1     ✔ tibble    3.3.0
✔ lubridate 1.9.4     ✔ tidyr     1.3.1
✔ purrr     1.2.0

── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ lubridate::%within%() masks IRanges::%within%()
✖ dplyr::collapse()     masks IRanges::collapse()
✖ dplyr::combine()      masks Biobase::combine(), BiocGenerics::combine()
✖ dplyr::count()        masks matrixStats::count()
✖ dplyr::desc()         masks IRanges::desc()
✖ tidyr::expand()       masks S4Vectors::expand()
✖ dplyr::filter()       masks stats::filter()
✖ dplyr::first()        masks S4Vectors::first()
✖ dplyr::lag()          masks stats::lag()
✖ ggplot2::Position()   masks BiocGenerics::Position(), base::Position()
✖ purrr::reduce()       masks GenomicRanges::reduce(), IRanges::reduce()
✖ dplyr::rename()       masks S4Vectors::rename()
✖ lubridate::second()   masks S4Vectors::second()
✖ lubridate::second<-() masks S4Vectors::second<-()
✖ dplyr::slice()        masks IRanges::slice()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

HPA <- read_tsv('~/Share/rna_single_cell_type.tsv.zip') |> 
    pivot_wider(names_from = `Cell type`, values_from = 'nTPM') |> 
    dplyr::select(-Gene) |> 
    distinct(`Gene name`, .keep_all = TRUE) |> 
    column_to_rownames('Gene name')

Rows: 1626642 Columns: 4
── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (3): Gene, Gene name, Cell type
dbl (1): nTPM

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

HPA_se <- SummarizedExperiment::SummarizedExperiment(HPA, colData = tibble(cell_type = colnames(HPA)), assays = list('logcounts' = HPA))

How many and what cell types are profiled?

Answer

length(HPA_se$cell_type)

[1] 81

HPA_se$cell_type

 [1] "Adipocytes"                      "Alveolar cells type 1"           "Alveolar cells type 2"           "Astrocytes"                      "B-cells"                         "Basal keratinocytes"            
 [7] "Basal prostatic cells"           "Basal respiratory cells"         "Basal squamous epithelial cells" "Bipolar cells"                   "Breast glandular cells"          "Breast myoepithelial cells"     
[13] "Cardiomyocytes"                  "Cholangiocytes"                  "Ciliated cells"                  "Club cells"                      "Collecting duct cells"           "Cone photoreceptor cells"       
[19] "Cytotrophoblasts"                "dendritic cells"                 "Distal enterocytes"              "Distal tubular cells"            "Ductal cells"                    "Early spermatids"               
[25] "Endometrial stromal cells"       "Endothelial cells"               "Enteroendocrine cells"           "Erythroid cells"                 "Excitatory neurons"              "Exocrine glandular cells"       
[31] "Extravillous trophoblasts"       "Fibroblasts"                     "Gastric mucus-secreting cells"   "Glandular and luminal cells"     "granulocytes"                    "Granulosa cells"                
[37] "Hepatocytes"                     "Hofbauer cells"                  "Horizontal cells"                "Inhibitory neurons"              "Intestinal goblet cells"         "Ionocytes"                      
[43] "Kupffer cells"                   "Langerhans cells"                "Late spermatids"                 "Leydig cells"                    "Lymphatic endothelial cells"     "Macrophages"                    
[49] "Melanocytes"                     "Mesothelial cells"               "Microglial cells"                "monocytes"                       "Mucus glandular cells"           "Muller glia cells"              
[55] "NK-cells"                        "Oligodendrocyte precursor cells" "Oligodendrocytes"                "Oocytes"                         "Ovarian stromal cells"           "Pancreatic endocrine cells"     
[61] "Paneth cells"                    "Peritubular cells"               "Plasma cells"                    "Prostatic glandular cells"       "Proximal enterocytes"            "Proximal tubular cells"         
[67] "Rod photoreceptor cells"         "Salivary duct cells"             "Schwann cells"                   "Secretory cells"                 "Serous glandular cells"          "Sertoli cells"                  
[73] "Skeletal myocytes"               "Smooth muscle cells"             "Spermatocytes"                   "Spermatogonia"                   "Squamous epithelial cells"       "Suprabasal keratinocytes"       
[79] "Syncytiotrophoblasts"            "T-cells"                         "Undifferentiated cells"

We can use the HPA reference to annotate cell types in the testis dataset.

Question

Use SingleR to annotate cell types in the testis dataset using HPA reference.

Answer

testis <- readRDS('~/Share/Guo_testis.rds')
predictions <- SingleR::SingleR(
    test = testis, 
    ref = HPA_se, 
    labels = colData(HPA_se)$cell_type, 
    clusters = testis$cluster
)
testis$annotation <- predictions$labels[testis$cluster]
cowplot::plot_grid(
    scater::plotReducedDim(testis, dimred = 'corrected', colour_by = 'cluster', text_by = 'cluster') + ggtitle('Testis data, PCA, graph-based clusters'), 
    scater::plotReducedDim(testis, dimred = 'corrected', colour_by = 'annotation', text_by = 'annotation') + ggtitle('PCA, Annotations transferred from HPA'),
    scater::plotReducedDim(testis, dimred = 'TSNE', colour_by = 'cluster', text_by = 'cluster') + ggtitle('Testis data, tSNE, graph-based clusters'), 
    scater::plotReducedDim(testis, dimred = 'TSNE', colour_by = 'annotation', text_by = 'annotation') + ggtitle('tSNE, Annotations transferred from HPA')
)

15.1.3 Filter the testis dataset to only germinal cells.

In this notebook, we will focus on germinal cells only.

Question

Filter the testis dataset to only retain germinal cells.

Answer

germcells <- testis[
    , 
    testis$annotation %in% c("Spermatogonia", "Spermatocytes", "Early spermatids", "Late spermatids")
]

15.2 Trajectory inference (TI) in scRNAseq

An important question in scRNAseq field of research is: how to identify a cell trajectory from high-dimensional expression data and map individual cells onto it? A large number of methods have currently emerged, each one with their own specificities, assumptions, and strengths. A nice breakdown (from 2019, so already very outdated!) is available from Saelens et al., Nat. Biotech. 2018 (doi: 10.1038/s41587-019-0071-9):

15.2.1 Trajectory

Slingshot is perhaps one of the most widely used algorithms for users who want to focus on R-based approaches.

Question

Read Slingshot documentation to understand how to identify lineages in a scRNAseq dataset in R.

Why is it recommended to infer lineages from PCA space rather than t-SNE or UMAP space, even though these spaces do “reveal” an obvious trajectory in 2D?

Infer lineages, using cluster annotations as information to build the MST. Note that you will first need to remove the 50th PCA dimension for slingshot to work (bug reported).

Answer

reducedDim(germcells, 'corrected_2') <- reducedDim(germcells, 'corrected')[, 1:49]
germcells_slingshot <- slingshot::slingshot(
    germcells, 
    reducedDim = 'corrected_2', 
    clusterLabels = germcells$annotation
)
slingshot::slingLineages(germcells_slingshot)

$Lineage1
[1] "Late spermatids"  "Early spermatids" "Spermatocytes"    "Spermatogonia"

15.2.2 Pseudotime

slingshot automatically deduces a pseudotime for each cell, based on the inferred lineage. It can be useful to visualize the trajectory of cells in a 2D space, and to order cells along the trajectory.

Question

Recover pseudotime values and principal curves from slingshot output.

Answer

# Fetching pseudotime
germcells$pseudotime <- slingshot::slingPseudotime(germcells_slingshot)[, 'Lineage1']

# Fetching principal curve in PCA space
pca_curve <- slingshot::slingCurves(germcells_slingshot, as.df = TRUE)
colnames(pca_curve) <- paste0('PC', 1:ncol(pca_curve))
head(pca_curve)

        PC1       PC2        PC3        PC4        PC5         PC6       PC7          PC8        PC9        PC10       PC11          PC12        PC13        PC14        PC15       PC16         PC17        PC18         PC19        PC20
1 0.5325495 0.1276968 -0.3338294 0.03607803 -0.2867371 -0.04204578 0.1373407 -0.009995032 0.02694292 -0.05042603 0.02989087  0.0052385866 -0.01891094 -0.02557450 -0.01074822 0.02006208 -0.004049899 0.005945587 -0.008946025 -0.01650377
2 0.5331398 0.1257388 -0.3247271 0.03485780 -0.2736647 -0.04013422 0.1309875 -0.009137282 0.02570672 -0.04767513 0.02789986  0.0036906281 -0.01921118 -0.02469666 -0.01185168 0.02123153 -0.005693288 0.008123469 -0.008691642 -0.01882372
3 0.5337301 0.1237808 -0.3156248 0.03363757 -0.2605924 -0.03822266 0.1246343 -0.008279533 0.02447052 -0.04492423 0.02590885  0.0021426696 -0.01951141 -0.02381882 -0.01295515 0.02240099 -0.007336677 0.010301352 -0.008437259 -0.02114368
4 0.5343203 0.1218225 -0.3065221 0.03241733 -0.2475201 -0.03631096 0.1182814 -0.007421829 0.02323440 -0.04217374 0.02391820  0.0005949673 -0.01981168 -0.02294110 -0.01405882 0.02357044 -0.008979922 0.012479056 -0.008182887 -0.02346336
5 0.5349071 0.1198600 -0.2974091 0.03119693 -0.2344467 -0.03439599 0.1119345 -0.006565157 0.02200032 -0.03943262 0.02193567 -0.0009469702 -0.02011281 -0.02206623 -0.01516704 0.02473980 -0.010619931 0.014652768 -0.007928743 -0.02577685
6 0.5354859 0.1178876 -0.2882723 0.02997607 -0.2213710 -0.03247363 0.1056018 -0.005710976 0.02077099 -0.03671313 0.01997189 -0.0024758689 -0.02041580 -0.02119776 -0.01628533 0.02590878 -0.012252541 0.016817534 -0.007674945 -0.02807629
           PC21        PC22         PC23        PC24          PC25        PC26        PC27         PC28         PC29         PC30        PC31         PC32         PC33        PC34          PC35        PC36         PC37         PC38
1 -0.0029292732 -0.01740418 -0.008778491 0.007464003  5.469541e-05 0.010087137 -0.01025984 -0.003594421 -0.006246169 -0.006892746 0.007830963 -0.001501482 -0.006010969 0.008765395 -5.878827e-04 0.008315953 -0.004101814 -0.002588275
2 -0.0023672892 -0.01709908 -0.009320922 0.007174715 -6.963008e-04 0.010016335 -0.01051595 -0.003669928 -0.006987168 -0.006512014 0.008072177 -0.001601186 -0.005422918 0.008652475 -2.692971e-04 0.008265890 -0.004417917 -0.002858201
3 -0.0018053053 -0.01679399 -0.009863353 0.006885427 -1.447297e-03 0.009945533 -0.01077205 -0.003745435 -0.007728167 -0.006131283 0.008313391 -0.001700890 -0.004834868 0.008539555  4.928844e-05 0.008215826 -0.004734021 -0.003128126
4 -0.0012433614 -0.01648882 -0.010405773 0.006596116 -2.198223e-03 0.009874717 -0.01102819 -0.003820893 -0.008469058 -0.005750444 0.008554561 -0.001800608 -0.004246773 0.008426616  3.678258e-04 0.008165781 -0.005050146 -0.003398065
5 -0.0006823155 -0.01618220 -0.010947931 0.006306299 -2.947566e-03 0.009803577 -0.01128523 -0.003895246 -0.009207520 -0.005367197 0.008794751 -0.001900646 -0.003657683 0.008313243  6.852795e-04 0.008116147 -0.005366749 -0.003668313
6 -0.0001230492 -0.01587233 -0.011489337 0.006015552 -3.693320e-03 0.009731523 -0.01154415 -0.003967304 -0.009940636 -0.004978854 0.009032770 -0.002001198 -0.003066598 0.008199056  1.000334e-03 0.008067239 -0.005684088 -0.003939072
         PC39        PC40          PC41        PC42         PC43         PC44         PC45        PC46         PC47          PC48          PC49 PC50 PC51
1 0.005797593 0.001019230 -1.190253e-03 0.002487841 -0.002784071 0.0006400310 0.0034053070 0.005315365 -0.002891229 -0.0001925498 -0.0022104639    1    1
2 0.006058814 0.001529594 -9.689303e-04 0.002640566 -0.002655263 0.0007108583 0.0028290951 0.004851459 -0.002537740 -0.0002736222 -0.0017384198    2    1
3 0.006320035 0.002039959 -7.476072e-04 0.002793292 -0.002526455 0.0007816857 0.0022528832 0.004387553 -0.002184250 -0.0003546946 -0.0012663757    3    1
4 0.006581281 0.002550333 -5.262806e-04 0.002946005 -0.002397629 0.0008525023 0.0016766696 0.003923651 -0.001830763 -0.0004357665 -0.0007943262    4    1
5 0.006843073 0.003060905 -3.048751e-04 0.003098456 -0.002268415 0.0009230780 0.0011004229 0.003459854 -0.001477308 -0.0005168298 -0.0003221607    5    1
6 0.007105849 0.003571608 -8.337284e-05 0.003250342 -0.002138492 0.0009931212 0.0005243531 0.002996226 -0.001123931 -0.0005980446  0.0001500939    6    1

We would like to visualize the trajectory of cells in our TSNE space, but the principal curve is only available in PCA space. We can use slingshot to embed it in the TSNE space!

Answer

tsne_curve <- slingshot::embedCurves(germcells_slingshot, 'TSNE', smoother = 'loess', span = 0.1) |> 
    slingshot::slingCurves(as.df = TRUE)
tsne_curve <- tsne_curve[order(tsne_curve$Order), ]
head(tsne_curve)

                              TSNE1    TSNE2 Order Lineage
Donor1-TTTATGCCATGACATC-1 -25.23033 11.72853     1       1
Donor1-GAACATCGTTTCCACC-1 -25.23033 11.72853     2       1
Donor1-GACCTGGCAGGGTATG-1 -25.23033 11.72853     3       1
Donor1-GCACTCTTCCCAAGTA-1 -25.23033 11.72853     4       1
Donor1-GCCAAATAGCGATAGC-1 -25.23033 11.72853     5       1
Donor1-GGCTCGATCCTAGAAC-1 -25.23033 11.72853     6       1

Plot PCA and tSNE embeddings, coloring cells by either their annotation or their pseudotime value. Overlay the principal curves to each embedding.

Answer

df <- tibble(
    PC1 = reducedDim(germcells, 'corrected')[,1], 
    PC2 = reducedDim(germcells, 'corrected')[,2], 
    TSNE1 = reducedDim(germcells, 'TSNE')[,1], 
    TSNE2 = reducedDim(germcells, 'TSNE')[,2], 
    annotation = germcells$annotation, 
    pseudotime = germcells$pseudotime
)
cowplot::plot_grid(
    df |> 
        ggplot() + 
        geom_point(aes(PC1, PC2, col = annotation)) + 
        geom_path(data = pca_curve, aes(x = PC1, y = PC2)) + 
        theme_bw() + 
        coord_fixed(),
    df |> 
        ggplot() + 
        geom_point(aes(TSNE1, TSNE2, col = annotation)) + 
        geom_path(data = tsne_curve, aes(x = TSNE1, y = TSNE2)) + 
        theme_bw() + 
        coord_fixed(),
    df |> 
        ggplot() + 
        geom_point(aes(PC1, PC2, col = pseudotime)) + 
        geom_path(data = pca_curve, aes(x = PC1, y = PC2)) + 
        theme_bw() + 
        coord_fixed(),
    df |> 
        ggplot() + 
        geom_point(aes(TSNE1, TSNE2, col = pseudotime)) + 
        geom_path(data = tsne_curve, aes(x = TSNE1, y = TSNE2)) + 
        theme_bw() + 
        coord_fixed()
)

Question

Check the distribution of pseudotime values across the different cell clusters.

What do you observe? Where you expecting this?

Answer

tibble(
    annotation = factor(germcells$annotation, c("Spermatogonia", "Spermatocytes", "Early spermatids", "Late spermatids")), 
    pseudotime = germcells$pseudotime
) |> 
    ggplot(aes(x = annotation, y = pseudotime, fill = annotation)) + 
    geom_violin(scale = 'width') + 
    geom_boxplot(outlier.shape = NULL, width = 0.1, fill = 'white') + 
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1))

Correct pseudotime values as you would expect it to be

germcells$pseudotime <- scales::rescale((-1 * slingshot::slingPseudotime(germcells_slingshot)[, 'Lineage1']), c(0, 1))

15.2.3 BONUS: Daunting snippet but that makes a cool figure for a paper: modeling pseudotime-dependent gene expression

Using pseudotime / cell, one can model gene expression along the differentiation process. This alleviates the need to study gene expression per cell, and allows one to focus on process-related effects (e.g. gene expression during a developmental trajectory).

Question

Try to do so for few markers of spermatogonia, spermatocytes and spermatids.

Answer

genes <- c('ID4', 'SYCP3', 'DMC1', 'ACR', 'PRM1', 'PGK2')
fitExprs <- logcounts(germcells[genes, ]) |> # ----------------------------------- Get norm. counts for genes of interest
    as.matrix() |> 
    t() |> 
    as_tibble() |> 
    mutate(  # ----------------------------------------------------------------- Add information for each cell
        cellID = colnames(germcells), 
        annotation = factor(germcells$annotation, c("Spermatogonia", "Spermatocytes", "Early spermatids", "Late spermatids")), 
        pseudotime = germcells$pseudotime
    ) |> 
    pivot_longer(contains(genes), names_to = 'gene', values_to = 'obs_expr') |> # - Pivot in "long" tidy format 
    mutate(gene = factor(gene, genes)) |> 
    group_by(gene) |> # ------------------------------------------------------- Group rows by genes
    nest(.key = 'data') |> # -------------------------------------------------- For each gene, extract the subtable into a column named data
    mutate(
        gamModel = map(data, ~mgcv::gam(obs_expr ~ s(pseudotime, bs = "cs"), data = .)), 
        gamFitted_expr = map(gamModel, predict) # ------------------------------ For each gene, fit the expression values ~ pseudotime with a GAM
    ) |> 
    dplyr::select(-ends_with('Model')) |> 
    unnest(c(data, ends_with('_expr'))) # -------------------------------------- Unnest all the modelled expressions
ggplot(fitExprs) + 
    geom_point(aes(x = pseudotime, y = obs_expr, col = annotation), alpha = 0.5) + 
    geom_line(aes(x = pseudotime, y = gamFitted_expr), col = 'white', size = 2, alpha = 0.5) + 
    geom_line(aes(x = pseudotime, y = gamFitted_expr), col = '#af2d0c', size = 1) +
    theme_bw() + 
    facet_grid(gene~., scales = 'free') + 
    labs(y = 'logcounts') + 
    ggtitle('Fitted models of pseudotime-dependent gene expression')

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.

15.3 Ordering trajectory with RNA velocity

As we saw earlier, TI does not necessarily know which direction is right for a given trajectory. This orientation can be sometimes refined using RNA velocity. Let’s see whether RNA velocity helps orientating our spermatocyte differentiation lineage trajectory here.

Question

Read velociraptor documentation. What do you need to compute RNA velocity scores in R?
Import spliced and unspliced counts computed with velocyto in R.
Try and compute RNA velocity.

Answer

## - Import entire GSE112013 dataset with spliced/unspliced counts 
full_GSE112013_counts <- LoomExperiment::import('~/Share/Guo_testis/Guo_testis_full-counts.loom')
full_GSE112013_counts

## - Filter `germcells` genes and cells to only retain those present in `full_GSE112013_counts`
germcells <- germcells[
    rowData(germcells)$gene_id %in% rowData(full_GSE112013_counts)$Accession, 
    germcells$Barcode %in% full_GSE112013_counts$Barcode
]

## - Reorder rows of `full_GSE112013_counts_germcells` to match those of `germcells`
full_GSE112013_counts_germcells <- full_GSE112013_counts[match(rowData(germcells)$gene_id, rowData(full_GSE112013_counts)$Accession), match(germcells$Barcode, full_GSE112013_counts$Barcode)]
dim(germcells)
dim(full_GSE112013_counts_germcells)

## - Add spliced/unspliced counts to germcells
assay(germcells, 'spliced', withDimnames=FALSE) <- as(assay(full_GSE112013_counts_germcells, 'spliced'), 'dgCMatrix')
assay(germcells, 'unspliced', withDimnames=FALSE) <- as(assay(full_GSE112013_counts_germcells, 'unspliced'), 'dgCMatrix')
germcells

## - Run velociraptor
velo_out <- velociraptor::scvelo(
    germcells, 
    assay.X = "counts", 
    subset.row = rowData(germcells)$isHVG, 
    use.dimred = "corrected"
)
velo_out

Embed the velocity field in tSNE scRNAseq embedding and plot the RNA velocity field on top of tSNE projection. Conclude.

embedded_velo <- velociraptor::embedVelocity(reducedDim(germcells, "TSNE"), velo_out)
head(embedded_velo)
grid_df <- velociraptor::gridVectors(reducedDim(germcells, "TSNE"), embedded_velo, resolution = 30, scale = TRUE)
head(grid_df)
scater::plotReducedDim(germcells, 'TSNE', colour_by = "annotation", point_alpha = 0.5) +
    geom_segment(
        data = grid_df, 
        mapping = aes(x = start.1, y = start.2, xend = end.1, yend = end.2), 
        arrow = arrow(length = unit(0.05, "inches"), type = "closed")
    )

15.4 Session info

─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.5.1 (2025-06-13)
 os       Ubuntu 24.04.2 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  en_IN.utf8
 ctype    en_IN.utf8
 tz       Europe/Paris
 date     2025-11-15
 pandoc   3.1.3 @ /usr/bin/ (via rmarkdown)
 quarto   1.5.55 @ /opt/quarto/bin/quarto

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package              * version date (UTC) lib source
 abind                  1.4-8   2024-09-12 [1] CRAN (R 4.5.1)
 beachmat               2.26.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 beeswarm               0.4.0   2021-06-01 [1] CRAN (R 4.5.1)
 Biobase              * 2.70.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 BiocGenerics         * 0.56.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 BiocNeighbors          2.4.0   2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 BiocParallel           1.44.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 BiocSingular           1.26.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 bit                    4.6.0   2025-03-06 [1] CRAN (R 4.5.0)
 bit64                  4.6.0-1 2025-01-16 [1] CRAN (R 4.5.0)
 cachem                 1.1.0   2024-05-16 [1] CRAN (R 4.5.1)
 cli                    3.6.5   2025-04-23 [1] CRAN (R 4.5.1)
 codetools              0.2-20  2024-03-31 [2] CRAN (R 4.5.1)
 cowplot                1.2.0   2025-07-07 [1] CRAN (R 4.5.1)
 crayon                 1.5.3   2024-06-20 [1] CRAN (R 4.5.1)
 DelayedArray           0.36.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 DelayedMatrixStats     1.32.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 devtools               2.4.6   2025-10-03 [1] CRAN (R 4.5.1)
 dichromat              2.0-0.1 2022-05-02 [1] CRAN (R 4.5.1)
 digest                 0.6.38  2025-11-12 [1] CRAN (R 4.5.1)
 dplyr                * 1.1.4   2023-11-17 [1] CRAN (R 4.5.1)
 ellipsis               0.3.2   2021-04-29 [1] CRAN (R 4.5.1)
 evaluate               1.0.5   2025-08-27 [1] CRAN (R 4.5.1)
 farver                 2.1.2   2024-05-13 [1] CRAN (R 4.5.1)
 fastmap                1.2.0   2024-05-15 [1] CRAN (R 4.5.1)
 forcats              * 1.0.1   2025-09-25 [1] CRAN (R 4.5.1)
 fs                     1.6.6   2025-04-12 [1] CRAN (R 4.5.1)
 generics             * 0.1.4   2025-05-09 [1] CRAN (R 4.5.1)
 GenomicRanges        * 1.62.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 ggbeeswarm             0.7.2   2023-04-29 [1] CRAN (R 4.5.1)
 ggplot2              * 4.0.1   2025-11-14 [1] CRAN (R 4.5.1)
 ggrepel                0.9.6   2024-09-07 [1] CRAN (R 4.5.1)
 glue                   1.8.0   2024-09-30 [1] CRAN (R 4.5.1)
 gridExtra              2.3     2017-09-09 [1] CRAN (R 4.5.1)
 gtable                 0.3.6   2024-10-25 [1] CRAN (R 4.5.1)
 hms                    1.1.4   2025-10-17 [1] CRAN (R 4.5.1)
 htmltools              0.5.8.1 2024-04-04 [1] CRAN (R 4.5.1)
 htmlwidgets            1.6.4   2023-12-06 [1] CRAN (R 4.5.1)
 igraph                 2.2.1   2025-10-27 [1] CRAN (R 4.5.1)
 IRanges              * 2.44.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 irlba                  2.3.5.1 2022-10-03 [1] CRAN (R 4.5.1)
 jsonlite               2.0.0   2025-03-27 [1] CRAN (R 4.5.1)
 knitr                  1.50    2025-03-16 [1] CRAN (R 4.5.1)
 labeling               0.4.3   2023-08-29 [1] CRAN (R 4.5.1)
 lattice                0.22-7  2025-04-02 [2] CRAN (R 4.5.1)
 lifecycle              1.0.4   2023-11-07 [1] CRAN (R 4.5.1)
 lubridate            * 1.9.4   2024-12-08 [1] CRAN (R 4.5.0)
 magrittr               2.0.4   2025-09-12 [1] CRAN (R 4.5.1)
 Matrix                 1.7-3   2025-03-11 [1] CRAN (R 4.5.1)
 MatrixGenerics       * 1.22.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 matrixStats          * 1.5.0   2025-01-07 [1] CRAN (R 4.5.0)
 memoise                2.0.1   2021-11-26 [1] CRAN (R 4.5.1)
 mgcv                   1.9-3   2025-04-04 [2] CRAN (R 4.5.1)
 nlme                   3.1-168 2025-03-31 [2] CRAN (R 4.5.1)
 pillar                 1.11.1  2025-09-17 [1] CRAN (R 4.5.1)
 pkgbuild               1.4.8   2025-05-26 [1] CRAN (R 4.5.1)
 pkgconfig              2.0.3   2019-09-22 [1] CRAN (R 4.5.1)
 pkgload                1.4.1   2025-09-23 [1] CRAN (R 4.5.1)
 princurve              2.1.6   2021-01-18 [1] CRAN (R 4.5.1)
 purrr                * 1.2.0   2025-11-04 [1] CRAN (R 4.5.1)
 R6                     2.6.1   2025-02-15 [1] CRAN (R 4.5.0)
 RColorBrewer           1.1-3   2022-04-03 [1] CRAN (R 4.5.1)
 Rcpp                   1.1.0   2025-07-02 [1] CRAN (R 4.5.1)
 readr                * 2.1.5   2024-01-10 [1] CRAN (R 4.5.1)
 remotes                2.5.0   2024-03-17 [1] CRAN (R 4.5.1)
 rlang                  1.1.6   2025-04-11 [1] CRAN (R 4.5.1)
 rmarkdown              2.30    2025-09-28 [1] CRAN (R 4.5.1)
 rsvd                   1.0.5   2021-04-16 [1] CRAN (R 4.5.1)
 S4Arrays               1.10.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 S4Vectors            * 0.48.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 S7                     0.2.0   2024-11-07 [1] CRAN (R 4.5.1)
 ScaledMatrix           1.18.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 scales                 1.4.0   2025-04-24 [1] CRAN (R 4.5.1)
 scater                 1.38.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 scrapper               1.4.0   2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 scuttle                1.20.0  2025-10-30 [1] Bioconductor 3.22 (R 4.5.1)
 Seqinfo              * 1.0.0   2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 sessioninfo            1.2.3   2025-02-05 [1] CRAN (R 4.5.0)
 SingleCellExperiment * 1.32.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 SingleR                2.12.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 slingshot              2.18.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 SparseArray            1.10.1  2025-10-31 [1] Bioconductor 3.22 (R 4.5.1)
 sparseMatrixStats      1.22.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 stringi                1.8.7   2025-03-27 [1] CRAN (R 4.5.1)
 stringr              * 1.6.0   2025-11-04 [1] CRAN (R 4.5.1)
 SummarizedExperiment * 1.40.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 tibble               * 3.3.0   2025-06-08 [1] CRAN (R 4.5.1)
 tidyr                * 1.3.1   2024-01-24 [1] CRAN (R 4.5.1)
 tidyselect             1.2.1   2024-03-11 [1] CRAN (R 4.5.1)
 tidyverse            * 2.0.0   2023-02-22 [1] CRAN (R 4.5.1)
 timechange             0.3.0   2024-01-18 [1] CRAN (R 4.5.1)
 TrajectoryUtils        1.18.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)
 tzdb                   0.5.0   2025-03-15 [1] CRAN (R 4.5.1)
 usethis                3.2.1   2025-09-06 [1] CRAN (R 4.5.1)
 vctrs                  0.6.5   2023-12-01 [1] CRAN (R 4.5.1)
 vipor                  0.4.7   2023-12-18 [1] CRAN (R 4.5.1)
 viridis                0.6.5   2024-01-29 [1] CRAN (R 4.5.1)
 viridisLite            0.4.2   2023-05-02 [1] CRAN (R 4.5.1)
 vroom                  1.6.6   2025-09-19 [1] CRAN (R 4.5.1)
 withr                  3.0.2   2024-10-28 [1] CRAN (R 4.5.1)
 xfun                   0.54    2025-10-30 [1] CRAN (R 4.5.1)
 XVector                0.50.0  2025-10-29 [1] Bioconductor 3.22 (R 4.5.1)

 [1] /home/rsg/R/x86_64-pc-linux-gnu-library/4.5
 [2] /opt/R/4.5.1/lib/R/library
 * ── Packages attached to the search path.

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────