CRUK CI Summer School 2020

Pseudotime Analysis with HCA Data

In many situations, one is studying a process where cells change continuously. This includes, for example, many differentiation processes taking place during development: following a stimulus, cells will change from one cell-type to another. Ideally, we would like to monitor the expression levels of an individual cell over time. Unfortunately, such monitoring is not possible with scRNA-seq since the cell is lysed (destroyed) when the RNA is extracted.

Instead, we must sample at multiple time-points and obtain snapshots of the gene expression profiles. Since some of the cells will proceed faster along the differentiation than others, each snapshot may contain cells at varying points along the developmental progression. We use statistical methods to order the cells along one or more trajectories which represent the underlying developmental trajectories, this ordering is referred to as “pseudotime”.

A recent benchmarking paper by Saelens et al provides a detailed summary of the various computational methods for trajectory inference from single-cell transcriptomics. They discuss 45 tools and evaluate them across various aspects including accuracy, scalability, and usability. They provide dynverse, an open set of packages to benchmark, construct and interpret single-cell trajectories (currently they have a uniform interface for 60 methods).

1 Setting up the data

library(SingleCellExperiment)
library(scran)
library(scater)
library(batchelor)
library(cowplot)
library(pheatmap)
library(tidyverse)
library(SingleR)
library(destiny)
library(gam)
library(viridis)
library(msigdbr)
library(clusterProfiler)
library(cellAlign)

We load the SCE object we have generated previously. This object contains only the T-cells from 8 healthy donors. We will first prepare the data by identifying variable genes, integrating the data across donors and calculating principal components.

sce.tcell<-readRDS(file="~/Course_Materials/scRNAseq/pseudotime/sce.tcell.RDS")

sce.tcell

class: SingleCellExperiment 
dim: 20352 21060 
metadata(0):
assays(2): counts logcounts
rownames(20352): AL627309.1 AL669831.5 ... AC233755.1 AC240274.1
rowData names(12): ensembl_gene_id external_gene_name ... gene_sparsity Chr
colnames: NULL
colData names(24): Sample Barcode ... PC1 PC2
reducedDimNames(4): MNN UMAP TSNE PCA
altExpNames(0):

dec.tcell <- modelGeneVar(sce.tcell, block=sce.tcell$Sample.Name)
top.tcell <- getTopHVGs(dec.tcell, n=5000)

set.seed(1010001)
merged.tcell <- fastMNN(sce.tcell, batch = sce.tcell$Sample.Name, subset.row = top.tcell)
reducedDim(sce.tcell, 'MNN') <- reducedDim(merged.tcell, 'corrected')

sce.tcell <- runPCA(sce.tcell, dimred="MNN")

plotPCA(sce.tcell, colour_by="Sample.Name")

1.1 Trajectory inference with destiny

Diffusion maps were introduced by Ronald Coifman and Stephane Lafon, and the underlying idea is to assume that the data are samples from a diffusion process. The method infers the low-dimensional manifold by estimating the eigenvalues and eigenvectors for the diffusion operator related to the data. Angerer et al have applied the diffusion maps concept to the analysis of single-cell RNA-seq data to create an R package called destiny.

For ease of computation, we will perform pseudotime analysis only on one sample, and we will downsample the object to 1000 cells. We will select the sample named MantonBM1.

# pull the barcodes for MantonBM1 sample & and downsample the set to 1000 genes 
vec.bc <- colData(sce.tcell) %>%
    data.frame() %>%
    filter(Sample.Name == "MantonBM1") %>%
    group_by(Sample.Name) %>%
    sample_n(1000) %>%
    pull(Barcode)

Number of cells in the sample:

table(colData(sce.tcell)$Barcode %in% vec.bc)

FALSE  TRUE 
20060  1000 

Subset cells from the main SCE object:

tmpInd <- which(colData(sce.tcell)$Barcode %in% vec.bc)

sce.tcell.BM1 <- sce.tcell[,tmpInd]

sce.tcell.BM1

class: SingleCellExperiment 
dim: 20352 1000 
metadata(0):
assays(2): counts logcounts
rownames(20352): AL627309.1 AL669831.5 ... AC233755.1 AC240274.1
rowData names(12): ensembl_gene_id external_gene_name ... gene_sparsity Chr
colnames: NULL
colData names(24): Sample Barcode ... PC1 PC2
reducedDimNames(4): MNN UMAP TSNE PCA
altExpNames(0):

Identify top 500 highly variable genes

dec.tcell.BM1 <- modelGeneVar(sce.tcell.BM1)
top.tcell.BM1 <- getTopHVGs(dec.tcell.BM1, n=500)

We will extract normalized counts for HVG to use in pseudotime alignment

tcell_BM1_counts<-logcounts(sce.tcell.BM1)
tcell_BM1_counts<-t(as.matrix(tcell_BM1_counts[top.tcell.BM1,]))
cellLabels <- sce.tcell.BM1$Barcode
rownames(tcell_BM1_counts)<-cellLabels

tcell_BM1_counts[1:4,1:4]

                        CCL5   S100A4     NKG7 KLRB1
AAACCTGGTACAGTTC-13 3.346141 3.639412 1.719101     0
AAACCTGGTGACTCAT-13 0.000000 0.000000 0.000000     0
AAACCTGTCCCAACGG-13 0.000000 0.000000 0.000000     0
AAACCTGTCCCTAACC-13 4.385532 3.608167 1.165970     0

And finally, we can run pseudotime alignment with destiny

dm_tcell_BM1 <- DiffusionMap(tcell_BM1_counts,n_pcs = 50)

Plot diffusion component 1 vs diffusion component 2 (DC1 vs DC2).

tmp <- data.frame(DC1 = eigenvectors(dm_tcell_BM1)[, 1],
                  DC2 = eigenvectors(dm_tcell_BM1)[, 2])

ggplot(tmp, aes(x = DC1, y = DC2)) +
    geom_point() + 
    xlab("Diffusion component 1") + 
    ylab("Diffusion component 2") +
    theme_classic()

Stash diffusion components to SCE object

sce.tcell.BM1$pseudotime_destiny_1<-eigenvectors(dm_tcell_BM1)[, 1]
sce.tcell.BM1$pseudotime_destiny_2<-eigenvectors(dm_tcell_BM1)[, 2]

1.2 Find temporally expressed genes

After running destiny, an interesting next step may be to find genes that change their expression over the course of time We demonstrate one possible method for this type of analysis on the 500 most variable genes. We will regress each gene on the pseudotime variable we have generated, using a general additive model (GAM). This allows us to detect non-linear patterns in gene expression. We are going to use HVG we identified in the previous step, but this analysis can also be done using the whole transcriptome.

# Only look at the 500 most variable genes when identifying temporally expressesd genes.
# Identify the variable genes by ranking all genes by their variance.
# We will use the first diffusion components as a measure of pseudotime 
Y<-log2(counts(sce.tcell.BM1)+1)
colnames(Y)<-cellLabels
Y<-Y[top.tcell.BM1,]
# Fit GAM for each gene using pseudotime as independent variable.
t <- eigenvectors(dm_tcell_BM1)[, 1]
gam.pval <- apply(Y, 1, function(z){
  d <- data.frame(z=z, t=t)
  tmp <- gam(z ~ lo(t), data=d)
  p <- summary(tmp)[4][[1]][1,5]
  p
})

Select top 30 genes for visualization

# Identify genes with the most significant time-dependent model fit.
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:30]  

Visualize these genes in a heatmap

heatmapdata <- Y[topgenes,]
heatmapdata <- heatmapdata[,order(t, na.last = NA)]
t_ann<-as.data.frame(t)
colnames(t_ann)<-"pseudotime"
pheatmap(heatmapdata, cluster_rows = T, cluster_cols = F, color = plasma(200), show_colnames = F, annotation_col = t_ann)

Visualize how some of the temporally expressed genes change in time

Following individual genes is very helpful for identifying genes that play an important role in the differentiation process. We illustrate the procedure using the GZMA gene. We have added the pseudotime values computed with destiny to the colData slot of the SCE object. Having done that, the full plotting capabilities of the scater package can be used to investigate relationships between gene expression, cell populations and pseudotime.

plotExpression(sce.tcell.BM1, "GZMA", x = "pseudotime_destiny_1", 
               show_violin = TRUE,
               show_smooth = TRUE)

1.2.1 Pseudotime analysis for another HCA sample

# pull barcodes for MantonBM2 
vec.bc <- colData(sce.tcell) %>%
    data.frame() %>%
    filter(Sample.Name == "MantonBM2") %>%
    group_by(Sample.Name) %>%
    sample_n(1000) %>%
    pull(Barcode)

# create another object for MantonBM2
tmpInd <- which(colData(sce.tcell)$Barcode %in% vec.bc)
sce.tcell.BM2 <- sce.tcell[,tmpInd]

# Identift HVG
dec.tcell.BM2 <- modelGeneVar(sce.tcell.BM2)
top.tcell.BM2 <- getTopHVGs(dec.tcell.BM2, n=500)

# extract normalized count data for HVG 
tcell_BM2_counts<-logcounts(sce.tcell.BM2)
tcell_counts_BM2<-t(as.matrix(tcell_BM2_counts[top.tcell.BM2,]))
cellLabels <- sce.tcell.BM2$Barcode
rownames(tcell_counts_BM2)<-cellLabels

dm_tcell_BM2 <- DiffusionMap(tcell_counts_BM2,n_pcs = 50)

tmp <- data.frame(DC1 = eigenvectors(dm_tcell_BM2)[, 1],
                  DC2 = eigenvectors(dm_tcell_BM2)[, 2])

ggplot(tmp, aes(x = DC1, y = DC2)) +
    geom_point() + 
    xlab("Diffusion component 1") + 
    ylab("Diffusion component 2") +
    theme_classic()

2 Challenge 1

Obtain pseudotime for one of the Caron samples.

sce_PRET1<-readRDS("~/Course_Materials/scRNAseq/pseudotime/sce_caron_PRET1.RDS")

You need to perform: * variance remodelling * HVG identification * extract normalized counts * run destiny * visualize diffusion components

sce_PRET1

class: SingleCellExperiment 
dim: 18372 1000 
metadata(0):
assays(2): counts logcounts
rownames(18372): TSPAN6 DPM1 ... AC003043.2 AL356417.3
rowData names(12): ensembl_gene_id external_gene_name ... gene_sparsity Chr
colnames: NULL
colData names(23): Sample Barcode ... label batch
reducedDimNames(1): PCA
altExpNames(0):

# Identift HVG
dec.caron.PRET1 <- modelGeneVar(sce_PRET1)
top.caron.PRET1 <- getTopHVGs(dec.caron.PRET1, n=500)

# extract normalized count data for HVG 
caron.PRET1_counts<-logcounts(sce_PRET1)
caron.PRET1_counts<-t(as.matrix(caron.PRET1_counts[top.caron.PRET1,]))
cellLabels <- sce_PRET1$Barcode
rownames(caron.PRET1_counts)<-cellLabels

dm_caron.PRET1 <- DiffusionMap(caron.PRET1_counts,n_pcs = 50)

tmp <- data.frame(DC1 = eigenvectors(dm_caron.PRET1)[, 1],
                  DC2 = eigenvectors(dm_caron.PRET1)[, 2])

ggplot(tmp, aes(x = DC1, y = DC2)) +
    geom_point() + 
    xlab("Diffusion component 1") + 
    ylab("Diffusion component 2") +
    theme_classic()

Y<-log2(counts(sce_PRET1)+1)
colnames(Y)<-cellLabels
Y<-Y[top.caron.PRET1,]
# Fit GAM for each gene using pseudotime as independent variable.
t <- eigenvectors(dm_caron.PRET1)[, 1]
gam.pval <- apply(Y, 1, function(z){
  d <- data.frame(z=z, t=t)
  tmp <- gam(z ~ lo(t), data=d)
  p <- summary(tmp)[4][[1]][1,5]
  p
})

Select top 30 genes for visualization

# Identify genes with the most significant time-dependent model fit.
topgenes <- names(sort(gam.pval, decreasing = FALSE))[1:30]  

Visualize these genes in a heatmap

heatmapdata <- Y[topgenes,]
heatmapdata <- heatmapdata[,order(t, na.last = NA)]
t_ann<-as.data.frame(t)
colnames(t_ann)<-"pseudotime"
pheatmap(heatmapdata, cluster_rows = T, cluster_cols = F, color = plasma(200), show_colnames = F, annotation_col = t_ann)

What kind of a dynamic process might be taking place in this cancer cell?

We can quickly check in which pathways these top genes are enriched using MSigDB.

Molecular Signatures Database contains 8 major collections:

• H: hallmark gene sets
• C1: positional gene sets
• C2: curated gene sets
• C3: motif gene sets
• C4: computational gene sets
• C5: GO gene sets
• C6: oncogenic signatures
• C7: immunologic signatures

We are going to use hallmark gene sets (H) and perform a hypergeometric test with our top 30 genes for all HALLMARK sets.

msigdb_hallmark<-msigdbr(species = "Homo sapiens", category = "H")  %>% select(gs_name, gene_symbol)
em<-enricher(topgenes, TERM2GENE=msigdb_hallmark)
head(em)[,"qvalue",drop=F]

3 Ackowledgements

This notebook uses material from SVI course, OSCA Book, Broad Institute Workshop and Hemberg Group Course.

sessionInfo()

R version 4.0.2 (2020-06-22)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 18.04.4 LTS

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.7.1
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.7.1

locale:
 [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C         LC_TIME=C            LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C       
 [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C         LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 

attached base packages:
 [1] splines   parallel  stats4    stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] scRNAseq_2.2.0              knitr_1.29                  cellAlign_0.1.0             clusterProfiler_3.16.0     
 [5] msigdbr_7.1.1               viridis_0.5.1               viridisLite_0.3.0           gam_1.20                   
 [9] foreach_1.5.0               destiny_3.2.0               SingleR_1.2.4               forcats_0.5.0              
[13] stringr_1.4.0               dplyr_1.0.0                 purrr_0.3.4                 readr_1.3.1                
[17] tidyr_1.1.0                 tibble_3.0.3                tidyverse_1.3.0             pheatmap_1.0.12            
[21] cowplot_1.0.0               batchelor_1.4.0             scater_1.16.2               ggplot2_3.3.2              
[25] scran_1.16.0                SingleCellExperiment_1.10.1 SummarizedExperiment_1.18.2 DelayedArray_0.14.1        
[29] matrixStats_0.56.0          Biobase_2.48.0              GenomicRanges_1.40.0        GenomeInfoDb_1.24.2        
[33] IRanges_2.22.2              S4Vectors_0.26.1            BiocGenerics_0.34.0        

loaded via a namespace (and not attached):
  [1] rappdirs_0.3.1                ggthemes_4.2.0                bit64_0.9-7.1                 irlba_2.3.3                  
  [5] data.table_1.12.8             RCurl_1.98-1.2                generics_0.0.2                RSQLite_2.2.0                
  [9] europepmc_0.4                 proxy_0.4-24                  bit_1.1-15.2                  enrichplot_1.8.1             
 [13] xml2_1.3.2                    lubridate_1.7.9               httpuv_1.5.4                  assertthat_0.2.1             
 [17] xfun_0.15                     hms_0.5.3                     evaluate_0.14                 promises_1.1.1               
 [21] DEoptimR_1.0-8                fansi_0.4.1                   progress_1.2.2                dbplyr_1.4.4                 
 [25] readxl_1.3.1                  igraph_1.2.5                  DBI_1.1.0                     ellipsis_0.3.1               
 [29] RSpectra_0.16-0               backports_1.1.8               bookdown_0.20                 vctrs_0.3.2                  
 [33] TTR_0.23-6                    abind_1.4-5                   RcppEigen_0.3.3.7.0           withr_2.2.0                  
 [37] ggforce_0.3.2                 packrat_0.5.0                 triebeard_0.3.0               robustbase_0.93-6            
 [41] vcd_1.4-7                     xts_0.12-0                    prettyunits_1.1.1             DOSE_3.14.0                  
 [45] ExperimentHub_1.14.0          laeken_0.5.1                  crayon_1.3.4                  labeling_0.3                 
 [49] edgeR_3.30.3                  pkgconfig_2.0.3               tweenr_1.0.1                  nlme_3.1-148                 
 [53] vipor_0.4.5                   nnet_7.3-14                   rlang_0.4.7                   lifecycle_0.2.0              
 [57] downloader_0.4                BiocFileCache_1.12.0          modelr_0.1.8                  rsvd_1.0.3                   
 [61] AnnotationHub_2.20.0          cellranger_1.1.0              polyclip_1.10-0               RcppHNSW_0.2.0               
 [65] lmtest_0.9-37                 Matrix_1.2-18                 urltools_1.7.3                carData_3.0-4                
 [69] boot_1.3-25                   zoo_1.8-8                     reprex_0.3.0                  base64enc_0.1-3              
 [73] beeswarm_0.2.3                ggridges_0.5.2                knn.covertree_1.0             bitops_1.0-6                 
 [77] blob_1.2.1                    DelayedMatrixStats_1.10.1     qvalue_2.20.0                 gridGraphics_0.5-0           
 [81] scales_1.1.1                  memoise_1.1.0                 magrittr_1.5                  plyr_1.8.6                   
 [85] hexbin_1.28.1                 zlibbioc_1.34.0               compiler_4.0.2                scatterpie_0.1.4             
 [89] dqrng_0.2.1                   RColorBrewer_1.1-2            pcaMethods_1.80.0             cli_2.0.2                    
 [93] dtw_1.21-3                    XVector_0.28.0                mgcv_1.8-31                   ggplot.multistats_1.0.0      
 [97] MASS_7.3-51.6                 tidyselect_1.1.0              stringi_1.4.6                 yaml_2.2.1                   
[101] GOSemSim_2.14.0               BiocSingular_1.4.0            locfit_1.5-9.4                ggrepel_0.8.2                
[105] grid_4.0.2                    fastmatch_1.1-0               tools_4.0.2                   rio_0.5.16                   
[109] rstudioapi_0.11               foreign_0.8-80                gridExtra_2.3                 smoother_1.1                 
[113] scatterplot3d_0.3-41          farver_2.0.3                  ggraph_2.0.3                  digest_0.6.25                
[117] rvcheck_0.1.8                 BiocManager_1.30.10           shiny_1.5.0                   Rcpp_1.0.5                   
[121] car_3.0-8                     broom_0.7.0                   BiocVersion_3.11.1            later_1.1.0.1                
[125] httr_1.4.1                    AnnotationDbi_1.50.1          colorspace_1.4-1              rvest_0.3.5                  
[129] fs_1.4.2                      ranger_0.12.1                 statmod_1.4.34                graphlayouts_0.7.0           
[133] sp_1.4-2                      ggplotify_0.0.5               xtable_1.8-4                  jsonlite_1.7.0               
[137] tidygraph_1.2.0               R6_2.4.1                      pillar_1.4.6                  htmltools_0.5.0              
[141] mime_0.9                      glue_1.4.1                    fastmap_1.0.1                 VIM_6.0.0                    
[145] BiocParallel_1.22.0           BiocNeighbors_1.6.0           class_7.3-17                  interactiveDisplayBase_1.26.3
[149] codetools_0.2-16              fgsea_1.14.0                  lattice_0.20-41               curl_4.3                     
[153] ggbeeswarm_0.6.0              zip_2.0.4                     GO.db_3.11.4                  openxlsx_4.1.5               
[157] limma_3.44.3                  rmarkdown_2.3                 munsell_0.5.0                 e1071_1.7-3                  
[161] DO.db_2.9                     GenomeInfoDbData_1.2.3        iterators_1.0.12              haven_2.3.1                  
[165] reshape2_1.4.4                gtable_0.3.0