Analysis of pallial and hem progenitors cycling dynamics using revelio

In this part of the analysis we apply Revelio algorithm to explore cell cycle dynamic of pallial and hem domain radial glial cells

Load libraries

library(Seurat)
library(Revelio)
library(orthologsBioMART)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)

#Set ggplot theme as classic
theme_set(theme_classic())

Load and filter progenitors data

Hem.data <- readRDS("../QC.filtered.clustered.cells.RDS")

DimPlot(object = Hem.data,
        group.by = "Cell_ident",
        reduction = "spring",
        cols = c("#ebcb2e", #"CR"
            "#e7823a", #"ChP"
            "#4cabdc", # Chp_prog
            "#68b041", #"Dorso-Medial_pallium" 
            "#e46b6b", #"Hem" 
            "#e3c148", #"Medial_pallium"
            "#046c9a", # Pallial
            "#4990c9"#"Thalamic_eminence"
                     ))

We first test the tool on the hem progenitor subset

Idents(Hem.data) <- Hem.data$Cell_ident

Progenitors.data <-  subset(Hem.data, idents = c("Dorso-Medial_pallium", "Hem","Medial_pallium"))

DimPlot(Progenitors.data,
        reduction = "spring",
        pt.size = 1,
        cols =  c("#68b041", "#e3c148", "#e46b6b")) + NoAxes()

Prepare data for revelio input

Find mous ortologues to provided human cell cycle genes

Cellcyclegenes <- revelioTestData_cyclicGenes
head(Cellcyclegenes)

## # A tibble: 6 × 5
##   G1.S    S        G2        G2.M    M.G1  
##   <fct>   <fct>    <fct>     <fct>   <fct> 
## 1 ABCA7   ABCC2    ALKBH1    ADH4    AFAP1 
## 2 ACD     ABCC5    ANLN      AHI1    AGFG1 
## 3 ACYP1   ABHD10   AP3D1     AKIRIN2 AGPAT3
## 4 ADAMTS1 ACPP     ARHGAP11B ANKRD40 AKAP13
## 5 ADCK2   ADAM22   ARHGAP19  ANLN    AMD1  
## 6 ADCY6   ANKRD18A ARL4A     ANP32B  ANP32E

We use orthologsBioMART library to map human to mouse mouse orthologs

G1.S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G1.S), 
                          to_attributes = "external_gene_name")

S <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$S), 
                          to_attributes = "external_gene_name")

G2 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G2), 
                          to_attributes = "external_gene_name")

G2.M <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$G2.M), 
                          to_attributes = "external_gene_name")

M.G1 <- findOrthologsHsMm(from_filters = "hgnc_symbol",
                          from_values = as.character(Cellcyclegenes$M.G1), 
                          to_attributes = "external_gene_name")


gene.list <- list(G1.S$external_gene_name,
                  S$external_gene_name,
                  G2$external_gene_name,
                  G2.M$external_gene_name,
                  M.G1$external_gene_name)

CCgenes <- t(plyr::ldply(gene.list, rbind))

colnames(CCgenes) <- colnames(Cellcyclegenes)

Export counts matrix

rawCounts <- as.matrix(Progenitors.data[["RNA"]]@counts)

# Filter genes expressed by less than 10 cells
num.cells <- Matrix::rowSums(rawCounts > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])
rawCounts <- rawCounts[genes.use, ]

rm(list = ls()[!ls() %in% c("rawCounts", "CCgenes", "Progenitors.data")])
gc()

##             used   (Mb) gc trigger   (Mb)  max used   (Mb)
## Ncells   5647931  301.7    8262065  441.3   7820564  417.7
## Vcells 354504046 2704.7 1023099440 7805.7 975369557 7441.5

Run Revelio

We can now follow the tutorial form the package github page

myData <- createRevelioObject(rawData = rawCounts,
                              cyclicGenes = CCgenes,
                              lowernGeneCutoff = 0,
                              uppernUMICutoff = Inf,
                              ccPhaseAssignBasedOnIndividualBatches = F)

## 2021-12-03 14:38:29: reading data: 11.8secs

rm("rawCounts")
gc()

##             used   (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells   5745695  306.9    8872374  473.9    8872374  473.9
## Vcells 354807483 2707.0 1030600065 7862.9 1030598716 7862.9

The getCellCyclePhaseAssignInformation filter “outlier” cells for cell cycle phase assignation. We modified the function to keep all cells as we observed this does not affect the global cell cycle fitting procedure

source("../Functions/functions_InitializationCCPhaseAssignFiltering.R")

myData <- getCellCyclePhaseAssign_allcells(myData)

## 2021-12-03 14:39:03: assigning cell cycle phases: 1.22mins

myData <- getPCAData(dataList = myData)

## 2021-12-03 14:40:55: calculating PCA: 56.98secs

myData <- getOptimalRotation(dataList = myData)

## 2021-12-03 14:43:25: calculating optimal rotation: 42.48secs

gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   5806337  310.1    8872374   473.9    8872374   473.9
## Vcells 619112417 4723.5 1781168111 13589.3 1781163515 13589.2

Graphical exploration of the infered cell cycle axis

CellCycledata <- cbind(as.data.frame(t(myData@transformedData$dc$data[1:2,])),
                       nUMI= myData@cellInfo$nUMI,
                       Revelio.phase = factor(myData@cellInfo$ccPhase, levels = c("G1.S", "S", "G2", "G2.M", "M.G1")),
                       Revelio.cc= myData@cellInfo$ccPercentageUniformlySpaced,
                       Seurat.cc= Progenitors.data@meta.data[myData@cellInfo$cellID,"CC.Difference"])

Cells distribution in the DC1-DC2 space

ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color = Revelio.phase)) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p1 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color = Revelio.phase)) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p2 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color=Revelio.cc), size=2, shape=16) + 
        scale_color_gradientn(colours=rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
                              name='Revelio_cc')


p1 + p2 

Cellular transcript content along the two computed axis

• Mean UMI +/- 95CI per .05 interval along the Revelio axis

meanUMI <- CellCycledata %>%
            mutate(Revelio.interval = cut(CellCycledata$Revelio.cc, seq(0, max(CellCycledata$Revelio.cc) + 0.05, 0.05), dig.lab = 2, right = FALSE)) %>%
            select(Revelio.interval, nUMI) %>%
            group_by(Revelio.interval) %>%
            summarise(n=n(),
                      meanUMI= mean(nUMI),
                      sd= sd(nUMI)) %>%
            mutate(se= sd/sqrt(n))  %>%
            mutate(ic= se * qt((1-0.05)/2 + .5, n-1))

meanUMI$max.Revelio.phase <- CellCycledata %>%
            mutate(Revelio.interval = cut(CellCycledata$Revelio.cc, seq(0, max(CellCycledata$Revelio.cc) + 0.05, 0.05), dig.lab = 2, right = FALSE)) %>%
            select(Revelio.interval, Revelio.phase) %>%
            group_by(Revelio.interval, Revelio.phase) %>%
            summarise(nb = n()) %>%
            filter(nb == max(nb)) %>%
            pull(Revelio.phase)

p1 <- ggplot(meanUMI[1:20,], aes(x=Revelio.interval, y=meanUMI/10000, fill=max.Revelio.phase)) +
        geom_bar(stat = "identity", width = 0.90) +
        geom_errorbar(aes(ymin=(meanUMI-ic)/10000, ymax=(meanUMI+ic)/10000), width=.2) +
        theme(axis.text.x = element_text(angle = 45, hjust=1))+
        scale_fill_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))


p2 <- ggplot(CellCycledata, aes(x= Revelio.cc, y= nUMI/10000)) +
        geom_point(aes(color= Revelio.phase), size=0.5) +
        scale_color_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) +
        geom_smooth(method="loess", n= 50, fill="grey") +
        ylim(0,NA)

p3 <- ggplot(CellCycledata, aes(x=Revelio.cc, y = Revelio.phase, fill= Revelio.phase)) +
        ggridges::geom_density_ridges() +
        scale_fill_manual(values= c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5]))

p1 / p2 / p3

Save cycling trajectories into Seurat metadata slot

Progenitors.data$Revelio.DC1 <- CellCycledata$DC1
Progenitors.data$Revelio.DC2 <- CellCycledata$DC2

Progenitors.data$Revelio.phase <- CellCycledata$Revelio.phase

Progenitors.data$Revelio.cc <- CellCycledata$Revelio.cc

p1 <- FeaturePlot(object = Progenitors.data,
            features = "Revelio.cc",
            pt.size = 1,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- DimPlot(object = Progenitors.data,
        group.by = "Revelio.phase",
        pt.size = 1,
        reduction = "spring",
        cols =  c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()

p1 + p2

Save Progenitors seurat obj

saveRDS(Progenitors.data,"Progenitors.RDS")

Session Info

#date
format(Sys.time(), "%d %B, %Y, %H,%M")

## [1] "24 février, 2022, 14,09"

#Packages used
sessionInfo()

## R version 4.1.2 (2021-11-01)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
## 
## locale:
##  [1] LC_CTYPE=fr_FR.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=fr_FR.UTF-8        LC_COLLATE=fr_FR.UTF-8    
##  [5] LC_MONETARY=fr_FR.UTF-8    LC_MESSAGES=fr_FR.UTF-8   
##  [7] LC_PAPER=fr_FR.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=fr_FR.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] wesanderson_0.3.6      cowplot_1.1.1          ggExtra_0.9           
##  [4] ggplot2_3.3.5          RColorBrewer_1.1-2     dplyr_1.0.7           
##  [7] Matrix_1.4-0           orthologsBioMART_0.1.0 data.table_1.14.2     
## [10] biomaRt_2.50.1         Revelio_0.1.0          SeuratObject_4.0.4    
## [13] Seurat_4.0.5          
## 
## loaded via a namespace (and not attached):
##   [1] BiocFileCache_2.2.0    plyr_1.8.6             igraph_1.2.11         
##   [4] lazyeval_0.2.2         splines_4.1.2          listenv_0.8.0         
##   [7] scattermore_0.7        GenomeInfoDb_1.30.0    digest_0.6.29         
##  [10] htmltools_0.5.2        fansi_0.5.0            magrittr_2.0.2        
##  [13] memoise_2.0.1          tensor_1.5             cluster_2.1.2         
##  [16] ROCR_1.0-11            globals_0.14.0         Biostrings_2.62.0     
##  [19] matrixStats_0.61.0     spatstat.sparse_2.0-0  prettyunits_1.1.1     
##  [22] colorspace_2.0-2       rappdirs_0.3.3         blob_1.2.2            
##  [25] ggrepel_0.9.1          xfun_0.28              crayon_1.4.2          
##  [28] RCurl_1.98-1.5         jsonlite_1.7.2         spatstat.data_2.1-0   
##  [31] survival_3.2-13        zoo_1.8-9              glue_1.5.1            
##  [34] polyclip_1.10-0        gtable_0.3.0           zlibbioc_1.40.0       
##  [37] XVector_0.34.0         leiden_0.3.9           future.apply_1.8.1    
##  [40] BiocGenerics_0.40.0    abind_1.4-5            scales_1.1.1          
##  [43] DBI_1.1.1              miniUI_0.1.1.1         Rcpp_1.0.8            
##  [46] viridisLite_0.4.0      xtable_1.8-4           progress_1.2.2        
##  [49] reticulate_1.22        spatstat.core_2.3-1    bit_4.0.4             
##  [52] stats4_4.1.2           htmlwidgets_1.5.4      httr_1.4.2            
##  [55] ellipsis_0.3.2         ica_1.0-2              farver_2.1.0          
##  [58] pkgconfig_2.0.3        XML_3.99-0.8           dbplyr_2.1.1          
##  [61] sass_0.4.0             uwot_0.1.10            deldir_1.0-6          
##  [64] utf8_1.2.2             labeling_0.4.2         tidyselect_1.1.1      
##  [67] rlang_0.4.12           reshape2_1.4.4         later_1.3.0           
##  [70] AnnotationDbi_1.56.2   munsell_0.5.0          tools_4.1.2           
##  [73] cachem_1.0.6           cli_3.1.0              generics_0.1.1        
##  [76] RSQLite_2.2.8          ggridges_0.5.3         evaluate_0.14         
##  [79] stringr_1.4.0          fastmap_1.1.0          yaml_2.2.1            
##  [82] goftest_1.2-3          knitr_1.36             bit64_4.0.5           
##  [85] fitdistrplus_1.1-6     purrr_0.3.4            RANN_2.6.1            
##  [88] KEGGREST_1.34.0        pbapply_1.5-0          future_1.23.0         
##  [91] nlme_3.1-153           mime_0.12              xml2_1.3.3            
##  [94] rstudioapi_0.13        compiler_4.1.2         filelock_1.0.2        
##  [97] curl_4.3.2             plotly_4.10.0          png_0.1-7             
## [100] spatstat.utils_2.2-0   tibble_3.1.6           bslib_0.3.1           
## [103] stringi_1.7.6          highr_0.9              lattice_0.20-45       
## [106] vctrs_0.3.8            pillar_1.6.4           lifecycle_1.0.1       
## [109] spatstat.geom_2.3-0    lmtest_0.9-39          jquerylib_0.1.4       
## [112] RcppAnnoy_0.0.19       bitops_1.0-7           irlba_2.3.3           
## [115] httpuv_1.6.3           patchwork_1.1.1        R6_2.5.1              
## [118] promises_1.2.0.1       KernSmooth_2.23-20     gridExtra_2.3         
## [121] IRanges_2.28.0         parallelly_1.29.0      codetools_0.2-18      
## [124] MASS_7.3-55            assertthat_0.2.1       withr_2.4.3           
## [127] sctransform_0.3.2      S4Vectors_0.32.3       GenomeInfoDbData_1.2.7
## [130] mgcv_1.8-38            parallel_4.1.2         hms_1.1.1             
## [133] grid_4.1.2             rpart_4.1.16           tidyr_1.1.4           
## [136] rmarkdown_2.11         Rtsne_0.15             Biobase_2.54.0        
## [139] shiny_1.7.1

---

1. Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩︎