Comparison between trajectories in the Gmnc WT/KO

Load libraries

library(Seurat)
library(princurve)
library(Revelio)
library(monocle)
library(gprofiler2)
library(seriation)
library(Matrix)
library(dplyr)
library(RColorBrewer)
library(ggplot2)
library(ggExtra)
library(cowplot)
library(wesanderson)

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

Load integrated datasets

WT.KO.integrated <- readRDS("WT_KO.integrated.RDS")
DefaultAssay(WT.KO.integrated) <- "RNA"

DimPlot(WT.KO.integrated,
        reduction = "integrated_spring",
        group.by = "Lineage",
        pt.size = 1,
        cols =   c("#cc391b","#e7823a","#969696","#026c9a")
        ) + NoAxes()

CPx.data <-  subset(WT.KO.integrated,
                        subset = Lineage %in% c("Choroid_Plexus") & orig.ident %in% c("Hem1", "Gmnc_KO"))

DimPlot(CPx.data,
        group.by = "orig.ident",
        reduction = "integrated_spring",
        pt.size = 1,
        cols =  c("#cc391b","#026c9a")
        ) + NoAxes()

rm(WT.KO.integrated)
gc()

##             used   (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells   3373160  180.2    6316384  337.4    5080289  271.4
## Vcells 149511464 1140.7 1017816426 7765.4 1158667588 8840.0

Pseudotime in WT

WT.data <-  subset(CPx.data,
                   subset = orig.ident =="Hem1")

Exclude septal cells

FeaturePlot(object = WT.data,
            features = c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "integrated_spring",
            order = T) & NoAxes() & NoLegend()

WT.data <- AddModuleScore(WT.data,
                           features = list(c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1")),
                           ctrl = 10,
                           name = "Septum")

FeaturePlot(object = WT.data,
            features = c("Septum1"),
            pt.size = 1,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "integrated_spring",
            order = T) & NoAxes()

WT.data$Septal.prog <- WT.data$Septum1 > 0

p1 <- DimPlot(WT.data,
        reduction = "integrated_spring",
        group.by = "Septal.prog",
        pt.size = 1) + NoAxes()

p2 <- FeaturePlot(object = WT.data,
            features = c("Fgf17"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "integrated_spring",
            order = T) & NoAxes() & NoLegend()

p1 + p2

WT.data<- subset(WT.data,
                   subset = Septal.prog == FALSE & WT.data$Spring_1 > 1300)

Fit principal curve

Trajectories.ChP <- WT.data@meta.data %>%
                    dplyr::select("Barcodes", "Spring_1", "Spring_2")

fit <- principal_curve(as.matrix(Trajectories.ChP[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 0.8, 
                       stretch=2)

## Starting curve---distance^2: 32583625527
## Iteration 1---distance^2: 22649223
## Iteration 2---distance^2: 20882914
## Iteration 3---distance^2: 20559041
## Iteration 4---distance^2: 20476272
## Iteration 5---distance^2: 20454340
## Iteration 6---distance^2: 20457685

#The principal curve smoothed
ChP.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) 

#Pseudotime score
Trajectories.ChP$Pseudotime <- fit$lambda/max(fit$lambda)

#Inverse the score if positive correlation with progenitor marker
if (cor(Trajectories.ChP$Pseudotime, WT.data@assays$SCT@data['Hmga2', Trajectories.ChP$Barcodes]) > 0) {
  Trajectories.ChP$Pseudotime <- -(Trajectories.ChP$Pseudotime - max(Trajectories.ChP$Pseudotime))
}

WT.data$Pseudotime <- Trajectories.ChP$Pseudotime

FeaturePlot(object = WT.data,
            features = "Pseudotime",
            pt.size = 2,
            cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
            reduction = "integrated_spring",
            order = T) & NoAxes()

Pseudotime in KO

KO.data <-  subset(CPx.data,
                   subset = orig.ident =="Gmnc_KO")

Exclude septal cells

FeaturePlot(object = KO.data,
            features = c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "integrated_spring",
            order = T) & NoAxes() & NoLegend()

KO.data <- AddModuleScore(KO.data,
                          features = list(c("Fgf8", "Fgf17", "Adamts15", "Fgfbp1")),
                          ctrl = 10,
                          name = "Septum")

FeaturePlot(object = KO.data,
            features = c("Septum1"),
            pt.size = 1,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "integrated_spring",
            order = T) & NoAxes()

KO.data$Septal.prog <- KO.data$Septum1 > 0

p1 <- DimPlot(KO.data,
              reduction = "integrated_spring",
              group.by = "Septal.prog",
              pt.size = 1) + NoAxes()

p2 <- FeaturePlot(object = KO.data,
                  features = c("Fgf17"),
                  pt.size = 0.5,
                  cols = c("grey90", brewer.pal(9,"YlGnBu")),
                  reduction = "integrated_spring",
                  order = T) & NoAxes() & NoLegend()

p1 + p2

KO.data<- subset(KO.data,
                 subset = Septal.prog == FALSE & KO.data$Spring_1 > 2000)

Fit principal curve

Trajectories.ChP <- KO.data@meta.data %>%
  dplyr::select("Barcodes", "Spring_1", "Spring_2")

fit <- principal_curve(as.matrix(Trajectories.ChP[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 0.8, 
                       stretch=2)

## Starting curve---distance^2: 5080053369
## Iteration 1---distance^2: 3624561
## Iteration 2---distance^2: 3590262
## Iteration 3---distance^2: 3588718

#The principal curve smoothed
ChP.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) 

#Pseudotime score
Trajectories.ChP$Pseudotime <- fit$lambda/max(fit$lambda)

#Inverse the score if positive correlation with progenitor marker
if (cor(Trajectories.ChP$Pseudotime, KO.data@assays$SCT@data['Hmga2', Trajectories.ChP$Barcodes]) > 0) {
  Trajectories.ChP$Pseudotime <- -(Trajectories.ChP$Pseudotime - max(Trajectories.ChP$Pseudotime))
}

KO.data$Pseudotime <- Trajectories.ChP$Pseudotime

FeaturePlot(object = KO.data,
            features = "Pseudotime",
            pt.size = 2,
            cols = rev(colorRampPalette(brewer.pal(n =10, name = "Spectral"))(100)),
            reduction = "integrated_spring",
            order = T) & NoAxes()

Subset the full integrated dataset

Trajectories <- rbind(WT.data@meta.data %>% select(Barcodes, Pseudotime), KO.data@meta.data %>% select(Barcodes, Pseudotime))

WT.KO.integrated <- readRDS("WT_KO.integrated.RDS")

meta.data <- WT.KO.integrated@meta.data[Trajectories$Barcodes,]
meta.data$Pseudotime <- Trajectories$Pseudotime

WT.KO.integrated <- CreateSeuratObject(counts = WT.KO.integrated@assays$RNA@counts[, Trajectories$Barcodes],
                                       meta.data = meta.data)

spring <- as.matrix(WT.KO.integrated@meta.data %>% select("Integrated_Spring_1", "Integrated_Spring_2"))
  
WT.KO.integrated[["integrated_spring"]] <- CreateDimReducObject(embeddings = spring, key = "Spring_", assay = DefaultAssay(WT.KO.integrated))

p1 <- FeaturePlot(object = WT.KO.integrated,
            features = "Pseudotime",
            pt.size = 0.5,
            cols = rev(colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)),
            reduction = "integrated_spring",
            order = T) & NoAxes()

p2 <- DimPlot(object = WT.KO.integrated,
        group.by = "orig.ident",
        pt.size = 0.5,
        reduction = "integrated_spring",
        cols =  c("#cc391b", "#026c9a")) & NoAxes()


p1 + p2

Normalization

WT.KO.integrated <- NormalizeData(WT.KO.integrated, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

WT.KO.integrated <- FindVariableFeatures(WT.KO.integrated, selection.method = "disp", nfeatures = 6500, assay = "RNA")

Plot some genes along pseudotime

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

Plot.Genes.trend(Seurat.data= WT.KO.integrated,
                 group.by = "Genotype",
                 genes= c("Nasp","Ttr","Htr2c", "Gmnc", "Trp73", "Foxj1", "Pifo", "Ccdc67"))

Use monocle2 to model gene expression along cycling axis

rm(list = ls()[!ls() %in% c("WT.KO.integrated")])
gc()

##            used  (Mb) gc trigger   (Mb)   max used    (Mb)
## Ncells  3498548 186.9    6316384  337.4    6316384   337.4
## Vcells 28981446 221.2  944420894 7205.4 1335753510 10191.0

Initialize a monocle object

# Transfer metadata
Annot.data  <- new('AnnotatedDataFrame',
                   data = data.frame(Barcode= WT.KO.integrated$Barcodes,
                                     Lineage= WT.KO.integrated$Lineage,
                                     Pseudotime= WT.KO.integrated$Pseudotime,
                                     Genotype= WT.KO.integrated$orig.ident))

# Transfer counts data
feature.data <- new('AnnotatedDataFrame',
                    data = data.frame(gene_short_name = rownames(WT.KO.integrated[["RNA"]]@counts),
                                      row.names = rownames(WT.KO.integrated[["RNA"]]@counts)))

# Create the CellDataSet object including variable genes only
gbm_cds <- newCellDataSet(WT.KO.integrated[["RNA"]]@counts,
                              phenoData = Annot.data,
                              featureData = feature.data,
                              lowerDetectionLimit = 0,
                              expressionFamily = negbinomial())

gbm_cds <- estimateSizeFactors(gbm_cds)
gbm_cds <- estimateDispersions(gbm_cds)
gbm_cds <- detectGenes(gbm_cds, min_expr = 0.1)

rm(list = ls()[!ls() %in% c("WT.KO.integrated", "gbm_cds")])
gc()

##            used  (Mb) gc trigger   (Mb)   max used    (Mb)
## Ncells  3566447 190.5    6316384  337.4    6316384   337.4
## Vcells 29460924 224.8  755536716 5764.3 1335753510 10191.0

Plot WT CR dynamic genes along WT and KO trajectories

# Load WT CPx variable genes along pseudotime
WT.CPx.genes <- read.table("../ChoroidPlexus_trajectory/ChP.Gene.dynamique.csv", sep = ";", header = T)

# Create a new pseudotime vector of 300 points
nPoints <- 100

new_data = list()
for (Genotype in unique(WT.KO.integrated$orig.ident)){
  new_data[[length(new_data) + 1]] = data.frame(Pseudotime = seq(min(WT.KO.integrated$Pseudotime), max(WT.KO.integrated$Pseudotime), length.out = nPoints), Genotype= Genotype)
}

new_data = do.call(rbind, new_data)

# Smooth gene expression
Diff.curve_matrix <- genSmoothCurves(gbm_cds[WT.CPx.genes$Gene, ],
                                      trend_formula = "~sm.ns(Pseudotime, df = 3)*Genotype",
                                      relative_expr = TRUE,
                                      new_data = new_data,
                                      cores= parallel::detectCores() - 2)

# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(Diff.curve_matrix[,1:100])), method = "pearson")))
row.ser <- seriation::seriate(dst, method ="R2E") #"R2E" #TSP #"GW" "GW_ward"
gene.order <- rownames(Diff.curve_matrix[,1:100][get_order(row.ser),])

# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(lineage = c(Gmnc_KO="#026c9a", Gmnc_WT="#cc391b"),
                    Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))


pheatmap::pheatmap(Diff.curve_matrix[gene.order,
                                c(200:101,#KO 
                                  1:100)], #WT
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = WT.CPx.genes  %>% dplyr::select(Gene.Clusters),
                   annotation_col = data.frame(lineage = rep(c("Gmnc_KO","Gmnc_WT"), each=100)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   border_color = NA,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

TF only

TFs <- read.table("TF.csv", sep = ";")[,1]

gene.order <- gene.order[gene.order %in% TFs]

heatmap.gene <- pheatmap::pheatmap(Diff.curve_matrix[gene.order,
                                c(200:101,#KO 
                                  1:100)], #WT
                       scale = "row",
                       cluster_rows = F,
                       cluster_cols = F,
                       annotation_row = WT.CPx.genes %>% dplyr::select(Gene.Clusters),
                       annotation_col = data.frame(lineage = rep(c("Gmnc_KO","Gmnc_WT"), each=100)),
                       annotation_colors = anno.colors,
                       show_colnames = F,
                       show_rownames = F,
                       fontsize_row = 8,
                       border_color = NA,
                       color =  viridis::viridis(9),
                       breaks = seq(-2.5,2.5, length.out = 9),
                       main = "WT CPx dynamic genes along GmncWT trajectories")

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

Plot.Genes.trend(Seurat.data= WT.KO.integrated,
                 group.by = "Genotype",
                 genes= c("Gmnc", "Trp73", "E2f7", "Foxj1", "Irx5", "Carhsp1", "Foxa2", "Sox9", "Pou3f2", "Myb", 
                          "Plagl1", "Prdm16", "Aebp1"))

KO dynamic genes along pseudotime

Find DEG in the KO

KO.pData <- pData(gbm_cds) %>% subset(Genotype == "Gmnc_KO")

pseudo.maturation.diff <- differentialGeneTest(gbm_cds[WT.KO.integrated[["RNA"]]@var.features, KO.pData$Barcode], 
                                                 fullModelFormulaStr = "~sm.ns(Pseudotime, df = 3)",
                                                 cores = parallel::detectCores() - 2)

# Filter genes based on FDR
pseudo.maturation.diff.filtered <- pseudo.maturation.diff %>% filter(qval < 1e-30)

# Create a new pseudo-DV vector of 300 points
nPoints <- 100

new_data = list()
for (Lineage in unique(KO.pData$Lineage)){
  new_data[[length(new_data) + 1]] = data.frame(Pseudotime = seq(min(KO.pData$Pseudotime), max(KO.pData$Pseudotime), length.out = nPoints), Lineage=Lineage)
}

new_data = do.call(rbind, new_data)

# Smooth gene expression
Diff.curve_matrix <- genSmoothCurves(gbm_cds[as.character(pseudo.maturation.diff.filtered$gene_short_name), KO.pData$Barcode],
                                      trend_formula = "~sm.ns(Pseudotime, df = 3)",
                                      relative_expr = TRUE,
                                      new_data = new_data,
                                      cores= parallel::detectCores() - 2)

Pseudotime.genes.clusters <- cluster::pam(as.dist((1 - cor(Matrix::t(Diff.curve_matrix),method = "pearson"))), k= 5)

KO_CPx_Gene.dynamique <- data.frame(Gene= names(Pseudotime.genes.clusters$clustering),
                                     Waves= Pseudotime.genes.clusters$clustering,
                                     Gene.Clusters = Pseudotime.genes.clusters$clustering,
                                     q.val = pseudo.maturation.diff.filtered$qval
                                     ) %>% arrange(Gene.Clusters)

row.names(KO_CPx_Gene.dynamique) <- KO_CPx_Gene.dynamique$Gene
KO_CPx_Gene.dynamique$Gene.Clusters <- paste0("Clust.", KO_CPx_Gene.dynamique$Gene.Clusters)

write.table(KO_CPx_Gene.dynamique, "KO_CPx_dynamic_genes.csv", sep = ";", quote = F, row.names = F)

# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(Diff.curve_matrix)), method = "pearson")))
row.ser <- seriation::seriate(dst, method ="R2E") #"R2E" #TSP #"GW" "GW_ward"
gene.order <- rev(rownames(Diff.curve_matrix[get_order(row.ser),]))

# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))


pheatmap::pheatmap(Diff.curve_matrix[gene.order,],
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = KO_CPx_Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   border_color = NA,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

Session Info

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

## [1] "16 juin, 2022, 16,15"

#Packages used
sessionInfo()

## R version 4.2.0 (2022-04-22)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.4 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] splines   stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] wesanderson_0.3.6   cowplot_1.1.1       ggExtra_0.9        
##  [4] RColorBrewer_1.1-2  dplyr_1.0.7         seriation_1.3.1    
##  [7] gprofiler2_0.2.1    monocle_2.22.0      DDRTree_0.1.5      
## [10] irlba_2.3.3         VGAM_1.1-5          ggplot2_3.3.5      
## [13] Biobase_2.54.0      BiocGenerics_0.40.0 Matrix_1.4-1       
## [16] Revelio_0.1.0       princurve_2.1.6     SeuratObject_4.0.4 
## [19] Seurat_4.0.5       
## 
## loaded via a namespace (and not attached):
##   [1] plyr_1.8.6            igraph_1.2.11         lazyeval_0.2.2       
##   [4] densityClust_0.3      listenv_0.8.0         scattermore_0.7      
##   [7] fastICA_1.2-3         digest_0.6.29         foreach_1.5.1        
##  [10] htmltools_0.5.2       viridis_0.6.2         fansi_0.5.0          
##  [13] magrittr_2.0.2        tensor_1.5            cluster_2.1.3        
##  [16] ROCR_1.0-11           limma_3.50.0          globals_0.14.0       
##  [19] matrixStats_0.61.0    docopt_0.7.1          spatstat.sparse_2.0-0
##  [22] colorspace_2.0-2      ggrepel_0.9.1         xfun_0.28            
##  [25] sparsesvd_0.2         crayon_1.4.2          jsonlite_1.7.2       
##  [28] spatstat.data_2.1-0   survival_3.2-13       zoo_1.8-9            
##  [31] iterators_1.0.13      glue_1.5.1            polyclip_1.10-0      
##  [34] registry_0.5-1        gtable_0.3.0          leiden_0.3.9         
##  [37] future.apply_1.8.1    abind_1.4-5           scales_1.1.1         
##  [40] pheatmap_1.0.12       DBI_1.1.1             miniUI_0.1.1.1       
##  [43] Rcpp_1.0.8            viridisLite_0.4.0     xtable_1.8-4         
##  [46] reticulate_1.22       spatstat.core_2.3-1   htmlwidgets_1.5.4    
##  [49] httr_1.4.2            FNN_1.1.3             ellipsis_0.3.2       
##  [52] ica_1.0-2             farver_2.1.0          pkgconfig_2.0.3      
##  [55] sass_0.4.0            uwot_0.1.10           deldir_1.0-6         
##  [58] utf8_1.2.2            labeling_0.4.2        tidyselect_1.1.1     
##  [61] rlang_0.4.12          reshape2_1.4.4        later_1.3.0          
##  [64] munsell_0.5.0         tools_4.2.0           generics_0.1.1       
##  [67] ggridges_0.5.3        evaluate_0.14         stringr_1.4.0        
##  [70] fastmap_1.1.0         yaml_2.2.1            goftest_1.2-3        
##  [73] knitr_1.36            fitdistrplus_1.1-6    purrr_0.3.4          
##  [76] RANN_2.6.1            pbapply_1.5-0         future_1.23.0        
##  [79] nlme_3.1-153          mime_0.12             slam_0.1-49          
##  [82] compiler_4.2.0        plotly_4.10.0         png_0.1-7            
##  [85] spatstat.utils_2.2-0  tibble_3.1.6          bslib_0.3.1          
##  [88] stringi_1.7.6         highr_0.9             lattice_0.20-45      
##  [91] HSMMSingleCell_1.14.0 vctrs_0.3.8           pillar_1.6.4         
##  [94] lifecycle_1.0.1       spatstat.geom_2.3-0   combinat_0.0-8       
##  [97] lmtest_0.9-39         jquerylib_0.1.4       RcppAnnoy_0.0.19     
## [100] data.table_1.14.2     httpuv_1.6.3          patchwork_1.1.1      
## [103] R6_2.5.1              promises_1.2.0.1      TSP_1.1-11           
## [106] KernSmooth_2.23-20    gridExtra_2.3         parallelly_1.29.0    
## [109] codetools_0.2-18      MASS_7.3-57           assertthat_0.2.1     
## [112] withr_2.4.3           qlcMatrix_0.9.7       sctransform_0.3.2    
## [115] mgcv_1.8-40           parallel_4.2.0        grid_4.2.0           
## [118] rpart_4.1.16          tidyr_1.1.4           rmarkdown_2.11       
## [121] Rtsne_0.15            shiny_1.7.1

---

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