Cajal-Retzius cells Trajectory

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 the full dataset

WT <- readRDS("../QC.filtered.clustered.cells.RDS")
KO <- readRDS("./GmncKO.cells.RDS") %>% subset(idents = c(6:9), invert = T)

p1 <- DimPlot(object = WT,
              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"
              ),
              pt.size = 0.5)  & NoAxes()


p2 <- DimPlot(object = KO,
              group.by = "Cell.ident",
              reduction = "spring",
              cols = c( "#4cabdc", # Chp_prog
                        "#68b041", #"Dorso-Medial_pallium" 
                        "#e46b6b", #"Hem" 
                        "#e3c148", #"Medial_pallium"
                        "#a9961b",
                        "#ebcb2e",
                        "#046c9a", # Pallial
                        "#4990c9"#"Thalamic_eminence"
              ),
              pt.size = 0.5)  & NoAxes()

p1 + p2

Expression of MCC and stress response in WT and KO

MCC.genes <- list(c("Trp73", "Gmnc", "Foxj1", "Myb", "Ccno", "Ccdc67", "Deup1","Mcidas",
                    "E2f4", "E2f5", "Ahr", "Trrap", "Cdc20b", "Ccdc78", "Rfx2",
                    "Rfx3", "Foxn4", "Fank1", "Jazf1", "Ccna1", "Nek10", "Plk4",
                    "Cep63", "Cep152", "Sass6", "Pcnt", "Pcm1", "Cetn2", "Tfdp1"))

KO <- AddModuleScore(KO,
                     features = MCC.genes,
                     name = "MCC_score")

WT <- AddModuleScore(WT,
                     features = MCC.genes,
                     name = "MCC_score")

WT.CR.goterm <- read.table("../CajalRetzius_trajectory/CR_GO_res-by-waves.csv", sep = ";", header = T)

DNA_damage_GOterm <- WT.CR.goterm %>% filter(term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571",
                                                            "GO:0006974", "GO:0006977","GO:0033554",
                                                            "GO:0044773", "GO:0042771", "GO:0042770",
                                                            "GO:2001021", "GO:1902229")
                                             )

DNA_damage_genes <- DNA_damage_GOterm %>%
                    filter(query %in% c("Clust.2", "Clust.3", "Clust.4")) %>%
                    filter(term_id == "GO:0033554") %>%
                    pull(intersection) %>% strsplit("\\,") %>% unlist() %>% unique()

KO <- AddModuleScore(KO,
                     features = list(DNA_damage_genes),
                     name = "cellular_response_to_stress_score")

WT <- AddModuleScore(WT,
                     features = list(DNA_damage_genes),
                     name = "cellular_response_to_stress_score")

gradient <- rev(brewer.pal(8,"RdYlBu"))
lim <-  c(-0.5,0.8)

p1 <- ggplot(KO@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=MCC_score1), size=1, shape=16) + 
  scale_color_gradientn(colours=gradient,
                        limits = lim,
                        name='Multiciliation score')


p2 <- ggplot(WT@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=MCC_score1), size=1, shape=16) + 
  scale_color_gradientn(colours=gradient,
                        limits = lim,
                        name='Multiciliation score')

p3 <- ggplot(KO@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=cellular_response_to_stress_score1), size=1, shape=16) + 
  scale_color_gradientn(colours=gradient,
                        limits = lim,
                        name='Stress response score')

p4 <- ggplot(WT@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=cellular_response_to_stress_score1), size=1, shape=16) + 
  scale_color_gradientn(colours= gradient,
                        limits = lim,
                        name='Stress response score')


MCC.scores.plot <- p2 + p1 
Stress.score.plot <- p4 + p3

for (i in 1:2){
  MCC.scores.plot[[i]]$data <- MCC.scores.plot[[i]]$data[order(MCC.scores.plot[[i]]$data$MCC_score1),]
}

for (i in 1:2){
  Stress.score.plot[[i]]$data <- Stress.score.plot[[i]]$data[order(Stress.score.plot[[i]]$data$cellular_response_to_stress_score1),]
}

MCC.scores.plot / Stress.score.plot

# Compute differentiation states scores

AP

APgenes <- c("Rgcc", "Sparc", "Hes5","Hes1", "Slc1a3",
             "Ddah1", "Ldha", "Hmga2","Sfrp1", "Id4",
             "Creb5", "Ptn", "Lpar1", "Rcn1","Zfp36l1",
             "Sox9", "Sox2", "Nr2e1", "Ttyh1", "Trip6")

KO <- AddModuleScore(KO,
                     features = list(APgenes),
                     name = "AP_signature")

BP

BPgenes <- c("Eomes", "Igsf8", "Insm1", "Elavl2", "Elavl4",
             "Hes6","Gadd45g", "Neurog2", "Btg2", "Neurog1")

KO <- AddModuleScore(KO,
                     features = list(BPgenes),
                     name = "BP_signature")

EN

ENgenes <- c("Mfap4", "Nhlh2", "Nhlh1", "Ppp1r14a", "Nav1",
             "Neurod1", "Sorl1", "Svip", "Cxcl12", "Tenm4",
             "Dll3", "Rgmb", "Cntn2", "Vat1")

KO <- AddModuleScore(KO,
                     features = list(ENgenes),
                     name = "EN_signature")

LN

LNgenes <- c("Snhg11", "Pcsk1n", "Mapt", "Ina", "Stmn4",
             "Gap43", "Tubb2a", "Ly6h","Ptprd", "Mef2c")

KO <- AddModuleScore(KO,
                     features = list(LNgenes),
                     name = "LN_signature")

FeaturePlot(object = KO,
            features = c("AP_signature1", "BP_signature1",
                              "EN_signature1", "LN_signature1"),
            pt.size = 0.75,
            cols = rev(brewer.pal(10,"Spectral")),
            reduction = "spring",
            order = T) & NoAxes() & NoLegend()

Fit pseudotime on CR and CP differentiating neurons

Group cells in Pallial or CR lineage

KO$Lineage <- sapply(KO$Cell.ident,
                              FUN = function(x) {
                                if (x %in% c("Neuron_prob.2", "Hem")) {
                                  x = "Cajal-Retzius_neurons"
                                } else if (x %in% c("Neuron_prob.3", "Medial_pallium")) {
                                  x = "Pallial_neurons"
                                } else {
                                  x = "other"
                                  }
                              })

DimPlot(KO,
        reduction = "spring",
        group.by = "Lineage",
        pt.size = 0.5,
        cols =  c("#cc391b","#969696","#026c9a")
        ) + NoAxes()

Fit principale curve on the two lineages

Neurons.data <-  subset(KO,  subset = Lineage %in% c("Cajal-Retzius_neurons", "Pallial_neurons") & Cell.ident %in% c("Neuron_prob.2", "Neuron_prob.3"))

DimPlot(Neurons.data ,
        reduction = "spring",
        group.by = "Lineage",
        pt.size = 1,
        cols =  c("#cc391b","#026c9a")
        ) + NoAxes()

fit <- principal_curve(as.matrix(Neurons.data@meta.data[,c("Spring_1", "Spring_2")]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 1,
                       stretch=0)

## Starting curve---distance^2: 76242448164
## Iteration 1---distance^2: 48308567
## Iteration 2---distance^2: 49237931
## Iteration 3---distance^2: 50119026
## Iteration 4---distance^2: 51105768
## Iteration 5---distance^2: 51819999
## Iteration 6---distance^2: 52369574
## Iteration 7---distance^2: 52732675
## Iteration 8---distance^2: 52936995
## Iteration 9---distance^2: 53060989
## Iteration 10---distance^2: 53117718

#Pseudotime score
PseudotimeScore <- fit$lambda/max(fit$lambda)

if (cor(PseudotimeScore, Neurons.data@assays$SCT@data['Hmga2', ]) > 0) {
  Neurons.data$PseudotimeScore <- -(PseudotimeScore - max(PseudotimeScore))
}

cols <- brewer.pal(n =11, name = "Spectral")

ggplot(Neurons.data@meta.data, aes(Spring_1, Spring_2)) +
  geom_point(aes(color=PseudotimeScore), size=2, shape=16) + 
  scale_color_gradientn(colours=rev(cols), name='Pseudotime score')

Plot pan-neuronal genes along this axis

Neurons.data <- NormalizeData(Neurons.data, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

Trajectories.neurons <- Neurons.data@meta.data %>% select(Barcodes, Spring_1, Spring_2,
                                                          AP_signature1, BP_signature1, EN_signature1, LN_signature1,
                                                          Lineage, PseudotimeScore)

# Neurog2
p1 <- FeaturePlot(object = Neurons.data,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons$Neurog2 <- Neurons.data@assays$RNA@data["Neurog2", Trajectories.neurons$Barcodes]

p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()
Trajectories.neurons$Tbr1 <- Neurons.data@assays$RNA@data["Tbr1", Trajectories.neurons$Barcodes]

p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

Trajectories.neurons$Mapt <- Neurons.data@assays$RNA@data["Mapt", Trajectories.neurons$Barcodes]

p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

## Shift pseudotime score

score <- sapply(Trajectories.neurons$PseudotimeScore,
                FUN = function(x) if (x <= 0.15) {x= 0.15} else { x=x })

Trajectories.neurons$PseudotimeScore.shifted <- (score - min(score)) / (max(score) - min(score))

# Neurog2
p1 <- FeaturePlot(object = Neurons.data ,
            features = c("Neurog2"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p2 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Neurog2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Tbr1 
p3 <- FeaturePlot(object = Neurons.data ,
            features = c("Tbr1"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p4 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Tbr1)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

# Mapt 
p5 <- FeaturePlot(object = Neurons.data ,
            features = c("Mapt"),
            pt.size = 0.5,
            cols = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction = "spring",
            order = T) & NoAxes()

p6 <- ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= Mapt)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

p1 + p2 + p3 + p4 + p5 + p6 + patchwork::plot_layout(ncol = 2)

Trajectories.neurons$nUMI <- Neurons.data$nCount_RNA

ggplot(Trajectories.neurons, aes(x= PseudotimeScore.shifted, y= nUMI/10000)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(method="loess", n= 50, aes(color= Lineage)) +
        ylim(0,NA)

Fit cell cycle trajectory on progenitors

To balance the number of progenitors in both domain we will only work with Hem and Medial_pallium annotated cells. Since we are using pallial cell to contrast CR specific trajectory we think this approximation will not significantly affect our analysis.

Progenitors.data <-  subset(KO,  subset = Lineage %in% c("Cajal-Retzius_neurons", "Pallial_neurons") & Cell.ident %in% c("Hem", "Medial_pallium"))

DimPlot(Progenitors.data,
        reduction = "spring",
        group.by = "Cell.ident",
        pt.size = 1,
        cols =  c("#cc391b","#026c9a")
        ) + NoAxes()

table(Progenitors.data$Cell.ident)

## 
##            Hem Medial_pallium 
##           1114           2983

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

##             used  (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells   3468720 185.3    6151937  328.6    6151937  328.6
## Vcells 102142773 779.3  821756700 6269.6 1026242019 7829.7

Prepare data for Revelio

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, ]

Run Revelio

CCgenes <- read.table("../ChoroidPlexus_trajectory/CCgenes.csv", sep = ";", header = T)

We can now follow the tutorial form the package github page

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

## 2022-05-09 16:04:23: reading data: 5.99secs

rm("rawCounts")
gc()

##             used   (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells   3474932  185.6    6151937  328.6    6151937  328.6
## Vcells 165303030 1261.2  657405360 5015.7 1026242019 7829.7

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)

## 2022-05-09 16:04:29: assigning cell cycle phases: 34.65secs

myData <- getPCAData(dataList = myData)

## 2022-05-09 16:05:04: calculating PCA: 27.97secs

myData <- getOptimalRotation(dataList = myData)

## 2022-05-09 16:05:32: calculating optimal rotation: 17.79secs

Graphical assesment of cell cycle fitting

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,
                       Domain= Progenitors.data$Cell.ident)

p1 <- ggplot(CellCycledata, aes(DC1, DC2)) +
        geom_point(aes(color = Revelio.phase), size= 0.5) +
        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 = Domain), size = 0.5) +
        scale_color_manual(values= c("#cc391b","#026c9a"))

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

p4 <- 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)

(p1 + p2) /(p3 + p4)

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

Transfert learn cell cycle axis

Progenitors <- Progenitors.data@meta.data %>% select(Barcodes, Spring_1, Spring_2,
                                                     AP_signature1, BP_signature1, EN_signature1, LN_signature1,
                                                     Lineage)

Progenitors$PseudotimeScore <- CellCycledata$Revelio.cc
Progenitors$nUMI <- Progenitors.data$nCount_RNA

Combine Progenitors and differentiating neurons data

# Start with neurons data
Trajectories.all <- Trajectories.neurons %>% select(Barcodes, nUMI, Spring_1, Spring_2, AP_signature1, BP_signature1, EN_signature1, LN_signature1, Lineage)

Trajectories.all$Pseudotime <- Trajectories.neurons$PseudotimeScore.shifted + 0.5
Trajectories.all$Phase <- NA

# Add progenitors data
Trajectories.progenitors <- Progenitors %>%
                              select(Barcodes, nUMI, Spring_1, Spring_2, AP_signature1, BP_signature1, EN_signature1, LN_signature1, Lineage) %>% 
                              mutate(Pseudotime= Progenitors.data$Revelio.cc/2,
                                     Phase = Progenitors.data$Revelio.phase)

Trajectories.all <- rbind(Trajectories.all, Trajectories.progenitors)

Trajectories.all$Phase <- factor(Trajectories.all$Phase, levels = c("G1.S", "S", "G2", "G2.M", "M.G1"))

p1 <- ggplot(Trajectories.all, aes(Spring_1, Spring_2)) +
        geom_point(aes(color=Pseudotime), size=0.5) + 
        scale_color_gradientn(colours=rev(brewer.pal(n =11, name = "Spectral")), name='Pseudotime score')

p2 <- ggplot(Trajectories.all, aes(Spring_1, Spring_2)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a"))

p1 + p2

p1 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= nUMI/10000)) +
        geom_point(aes(color= 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)

p2 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= nUMI/10000)) +
        geom_point(aes(color= Lineage), size=0.5) +
        scale_color_manual(values= c("#cc391b", "#026c9a")) +
        geom_smooth(aes(color= Lineage), method="loess", n= 50, fill="grey") +
        ylim(0,NA)

p1 / p2

p1 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= AP_signature1)) +
  geom_point(aes(color= Lineage), size=0.5) +
  scale_color_manual(values= c("#cc391b", "#026c9a")) +
  geom_smooth(aes(color= Lineage), method="loess", n= 50, fill="grey")


p2 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= BP_signature1)) +
  geom_point(aes(color= Lineage), size=0.5) +
  scale_color_manual(values= c("#cc391b", "#026c9a")) +
  geom_smooth(aes(color= Lineage), method="loess", n= 50, fill="grey")

p3 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= EN_signature1)) +
  geom_point(aes(color= Lineage), size=0.5) +
  scale_color_manual(values= c("#cc391b", "#026c9a")) +
  geom_smooth(aes(color= Lineage), method="loess", n= 50, fill="grey")

p4 <- ggplot(Trajectories.all, aes(x= Pseudotime, y= LN_signature1)) +
  geom_point(aes(color= Lineage), size=0.5) +
  scale_color_manual(values= c("#cc391b", "#026c9a")) +
  geom_smooth(aes(color= Lineage), method="loess", n= 50, fill="grey")


p1 / p2 / p3 / p4

Subset the full Seurat object

KO <- readRDS("./GmncKO.cells.RDS") %>% subset(idents = c(6:9), invert = T)

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

spring <- as.matrix(Neuro.trajectories@meta.data %>% select("Spring_1", "Spring_2"))
  
Neuro.trajectories[["spring"]] <- CreateDimReducObject(embeddings = spring, key = "Spring_", assay = DefaultAssay(Neuro.trajectories))

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

p2 <- DimPlot(object = Neuro.trajectories,
        group.by = "Lineage",
        pt.size = 0.5,
        reduction = "spring",
        cols =  c("#cc391b", "#026c9a")) & NoAxes()


p3 <- DimPlot(object = Neuro.trajectories,
        group.by = "Phase",
        pt.size = 0.5,
        reduction = "spring",
        cols =  c(wes_palette("FantasticFox1")[1:3],"grey40",wes_palette("FantasticFox1")[5])) & NoAxes()

p1 + p2 + p3

rm(list = ls()[!ls() %in% c("Neuro.trajectories")])
gc()

##            used  (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells  3517520 187.9    6151937  328.6    6151937  328.6
## Vcells 35253774 269.0  757471775 5779.1 1026242019 7829.7

Normalization

Neuro.trajectories<- NormalizeData(Neuro.trajectories, normalization.method = "LogNormalize", scale.factor = 10000, assay = "RNA")

Neuro.trajectories <- FindVariableFeatures(Neuro.trajectories, selection.method = "disp", nfeatures = 3000, assay = "RNA")

Plot some genes along pseudotime

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

Plot.Genes.trend(Seurat.data= Neuro.trajectories,
                 group.by = "Lineage",
                 genes= c("Gas1","Sox2",
                          "Neurog2", "Btg2",
                          "Tbr1", "Mapt",
                          "Trp73", "Foxg1"))

Plot.Genes.trend(Seurat.data= Neuro.trajectories,
                 group.by = "Lineage",
                 genes= c("Gmnc", "Mcidas",
                          "Foxj1", "Trp73",
                          "Lhx1", "Cdkn1a"))

Plot.Genes.trend(Seurat.data= Neuro.trajectories,
                 group.by = "Lineage",
                 genes= c("Mki67", "Top2a",
                          "H2afx", "Cdkn1c"))

Use monocle2 to model gene expression along cycling axis

Initialize a monocle object

# Transfer metadata
meta.data <- data.frame(Barcode= Neuro.trajectories$Barcodes,
                        Lineage= Neuro.trajectories$Lineage,
                        Pseudotime= Neuro.trajectories$Pseudotime,
                        Phase= Neuro.trajectories$Phase)

Annot.data  <- new('AnnotatedDataFrame', data = meta.data)

# Transfer counts data
var.genes <- Neuro.trajectories[["RNA"]]@var.features
count.data = data.frame(gene_short_name = rownames(Neuro.trajectories[["RNA"]]@data[var.genes,]),
                        row.names = rownames(Neuro.trajectories[["RNA"]]@data[var.genes,]))

feature.data <- new('AnnotatedDataFrame', data = count.data)

# Create the CellDataSet object including variable genes only
gbm_cds <- newCellDataSet(Neuro.trajectories[["RNA"]]@counts[var.genes,],
                          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("Neuro.trajectories", "gbm_cds", "Gene.Trend", "Plot.Genes.trend")])
gc()

##            used  (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells  3587157 191.6    6151937  328.6    6151937  328.6
## Vcells 71157383 542.9  605977420 4623.3 1026242019 7829.7

Find Pan-neuronal genes

# Split pallial and subpallial cells for gene expression fitting
#Pallial cells
Pallialcells <- Neuro.trajectories@meta.data %>%
                filter(Lineage == "Pallial_neurons") %>%
                pull(Barcodes)

# Cajal-Retzius cells
CRcells <- Neuro.trajectories@meta.data %>%
                   filter(Lineage == "Cajal-Retzius_neurons") %>%
                   pull(Barcodes)

# We filter-out genes detected in less than 200 or 200 cells along Pallial or CR lineages
num.cells <- Matrix::rowSums(Neuro.trajectories@assays$RNA@counts[,Pallialcells] > 0)
Pallial.expressed <- names(x = num.cells[which(x = num.cells >= 200)])

num.cells <- Matrix::rowSums(Neuro.trajectories@assays$RNA@counts[,CRcells] > 0)
CR.expressed <- names(x = num.cells[which(x = num.cells >= 200)])

Input.genes <- rownames(gbm_cds)[rownames(gbm_cds) %in% intersect(Pallial.expressed, CR.expressed)]

Pallial.genes <- differentialGeneTest(gbm_cds[Input.genes, Pallialcells], 
                                                 fullModelFormulaStr = "~sm.ns(Pseudotime, df = 3)", 
                                                 reducedModelFormulaStr = "~1", 
                                                 cores = parallel::detectCores() - 2)

#Filter based on FDR
Pallial.genes.filtered <- Pallial.genes  %>% filter(qval < 1e-3)

CRcells.genes <- differentialGeneTest(gbm_cds[Input.genes, CRcells], 
                                                 fullModelFormulaStr = "~sm.ns(Pseudotime, df = 3)", 
                                                 reducedModelFormulaStr = "~1", 
                                                 cores = parallel::detectCores() - 2)

#Filter based on FDR
CRcells.genes.filtered <- CRcells.genes  %>% filter(qval < 1e-3)

Common.genes <- intersect(Pallial.genes.filtered$gene_short_name, CRcells.genes.filtered$gene_short_name)

# Smooth genes expression along the two trajectories
nPoints <- 300

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

new_data = do.call(rbind, new_data)

# Smooth gene expression
curve_matrix <- genSmoothCurves(gbm_cds[as.character(Common.genes),],
                                trend_formula = "~sm.ns(Pseudotime, df = 3)*Lineage",
                                relative_expr = TRUE,
                                new_data = new_data,
                                cores= parallel::detectCores() - 2)

# Extract genes with person's cor > 0.6 between the 2 trajectories

Pallial.smoothed <- scale(t(curve_matrix[,c(1:300)]))  #Pallial points
CR.smoothed <- scale(t(curve_matrix[,c(301:600)])) #CR points

mat <- cor(Pallial.smoothed, CR.smoothed, method = "pearson")

Gene.Cor <- diag(mat)
hist(Gene.Cor, breaks = 100)
abline(v=0.8,col=c("blue"))

PanNeuro.genes <- names(Gene.Cor[Gene.Cor > 0.8])

# Order rows using seriation
dst <- as.dist((1-cor(scale(t(curve_matrix[PanNeuro.genes,c(600:301)])), method = "pearson")))
row.ser <- seriate(dst, method ="MDS_angle") #MDS_angle
gene.order <- PanNeuro.genes[get_order(row.ser)]

anno.colors <- list(lineage = c(Pallial="#026c9a",CR="#cc391b"))


pheatmap::pheatmap(curve_matrix[rev(gene.order),
                                c(1:300, 301:600)], #CR
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_col = data.frame(lineage = rep(c("Pallial","CR"), each=300)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = T,
                   fontsize_row = 2,
                   color =  viridis::viridis(10),
                   breaks = seq(-2.5,2.5, length.out = 10),
                   main = "")

rm(list = ls()[!ls() %in% c("Neuro.trajectories", "gbm_cds", "Gene.Trend", "Plot.Genes.trend")])
gc()

##            used  (Mb) gc trigger   (Mb)   max used   (Mb)
## Ncells  3621897 193.5    6151937  328.6    6151937  328.6
## Vcells 71246487 543.6  484781936 3698.6 1026242019 7829.7

Test each gene trend over pseudotime score

Find genes DE along pseudomaturation axis

pseudo.maturation.diff <- differentialGeneTest(gbm_cds[fData(gbm_cds)$num_cells_expressed > 80,], 
                                                 fullModelFormulaStr = "~sm.ns(Pseudotime, df = 3)*Lineage", 
                                                 reducedModelFormulaStr = "~sm.ns(Pseudotime, df = 3)", 
                                                 cores = parallel::detectCores() - 2)

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

Direction of the DEG by calculating the area between curves (ABC)

Smooth commun genes along the two trajectories

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

new_data = list()
for (Lineage in unique(pData(gbm_cds)$Lineage)){
  new_data[[length(new_data) + 1]] = data.frame(Pseudotime = seq(min(pData(gbm_cds)$Pseudotime), max(pData(gbm_cds)$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),],
                                      trend_formula = "~sm.ns(Pseudotime, df = 3)*Lineage",
                                      relative_expr = TRUE,
                                      new_data = new_data,
                                      cores= parallel::detectCores() - 2)

Compute the ABC for each gene

# Extract matrix containing smoothed curves for each lineages
Pal_curve_matrix <- Diff.curve_matrix[, 1:nPoints] #Pallial points
CR_curve_matrix <- Diff.curve_matrix[, (nPoints + 1):(2 * nPoints)] #CR points

# Direction of the comparison : postive ABCs <=> Upregulated in CR lineage
ABCs_res <- CR_curve_matrix - Pal_curve_matrix

# Average logFC between the 2 curves
ILR_res <- log2(CR_curve_matrix/ (Pal_curve_matrix + 0.1)) 
  
ABCs_res <- apply(ABCs_res, 1, function(x, nPoints) {
                  avg_delta_x <- (x[1:(nPoints - 1)] + x[2:(nPoints)])/2
                  step <- (100/(nPoints - 1))
                  res <- round(sum(avg_delta_x * step), 3)
                  return(res)},
                  nPoints = nPoints) # Compute the area below the curve
  
ABCs_res <- cbind(ABCs_res, ILR_res[,ncol(ILR_res)])
colnames(ABCs_res)<- c("ABCs", "Endpoint_ILR")

# Import ABC values into the DE test results table
pseudo.maturation.diff.filtered <- cbind(pseudo.maturation.diff.filtered[,1:4],
                                         ABCs_res,
                                         pseudo.maturation.diff.filtered[,5:6])

Cajal-Retzius cells specific trajectory analysis

# Extract Cajal-Retzius expressed genes
CR.res <- as.data.frame(pseudo.maturation.diff.filtered[pseudo.maturation.diff.filtered$ABCs > 0,])
CR.genes <- row.names(CR.res)

CR_curve_matrix <- CR_curve_matrix[CR.genes, ]

Gene expression profiles along the trajectory

## Cluster gene by expression profiles
Pseudotime.genes.clusters <- cluster::pam(as.dist((1 - cor(Matrix::t(CR_curve_matrix),method = "pearson"))), k= 5)

CR.Gene.dynamique <- data.frame(Gene= names(Pseudotime.genes.clusters$clustering),
                                 Waves= Pseudotime.genes.clusters$clustering,
                                 Gene.Clusters = Pseudotime.genes.clusters$clustering,
                                 q.val = CR.res$qval,
                                 ABCs= CR.res$ABCs
                                 ) %>% arrange(Gene.Clusters)

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

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

# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(lineage = c(Pallial_neurons="#026c9a", Cajal_Retzius="#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(300:1,#Pal 
                                  301:600)], #CR
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=300)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

We manually correct the reordering so genes are aligned from top left to bottom rigth

pheatmap::pheatmap(Diff.curve_matrix[gene.order,
                                c(300:1,#Pal 
                                  301:600)], #CR
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=300)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

anno.colors <- list(Cell.state = c(Cycling_RG="#046c9a", Differentiating_cells="#ebcb2e"),
                    Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))

col.anno <- data.frame(Cell.state = rep(c("Cycling_RG","Differentiating_cells"), c(100,200)))
rownames(col.anno) <- 301:600

pheatmap::pheatmap(CR_curve_matrix[gene.order,],
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = col.anno,
                   annotation_colors = anno.colors,
                   gaps_col = 100,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

diff.state <- Neuro.trajectories@meta.data %>%
              filter(Lineage ==  "Cajal-Retzius_neurons") %>%
              select("AP_signature1", "BP_signature1", "EN_signature1", "LN_signature1", "Pseudotime")

AP.loess <- loess(AP_signature1 ~ Pseudotime, diff.state)
AP.smooth <- predict(AP.loess,
                     seq(0.01,1.5, length.out= 300))

BP.loess <- loess(BP_signature1 ~ Pseudotime, diff.state)
BP.smooth <- predict(BP.loess,
                     seq(0.01,1.5, length.out= 300))

EN.loess <- loess(EN_signature1 ~ Pseudotime, diff.state)
EN.smooth <- predict(EN.loess,
                     seq(0.01,1.5, length.out= 300))

LN.loess <- loess(LN_signature1 ~ Pseudotime, diff.state)
LN.smooth <- predict(LN.loess,
                     seq(0.01,1.5, length.out= 300))

Smoothed.states <- cbind(AP.smooth, BP.smooth, EN.smooth, LN.smooth)

heatmap.states <- pheatmap::pheatmap(as.data.frame(t(Smoothed.states)),
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   gaps_col = 100,
                   gaps_row = c(1,2,3),
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  rev(colorRampPalette(brewer.pal(n= 8, name = "RdBu"))(100)),
                   breaks = seq(-1,1, length.out = 100),
                   main = "")

heatmap.gene <- pheatmap::pheatmap(CR_curve_matrix[gene.order,],
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   gaps_col = 100,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

cowplot::plot_grid(heatmap.states$gtable, heatmap.gene$gtable,
                   ncol = 1,
                   align = "v",
                   rel_heights = c(1,3),
                   greedy = T)

Gene cluster trend

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

Plot.clust.trends(Neuro.trajectories,
                   Lineage = "Cajal-Retzius_neurons",
                   Which.cluster = 1:5,
                   clust.list = Pseudotime.genes.clusters$clustering,
                   Smooth.method = "gam")

GO term enrichment in gene clusters using gprofiler2

CR.gostres <- gost(query = list("Clust.1" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.1") %>% pull(Gene) %>% as.character(),
                             "Clust.2" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.2") %>% pull(Gene) %>% as.character(),
                             "Clust.3" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.3") %>% pull(Gene) %>% as.character(),
                             "Clust.4" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.4") %>% pull(Gene) %>% as.character(),
                             "Clust.5" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.5") %>% pull(Gene) %>% as.character()),
                organism = "mmusculus", ordered_query = F, 
                multi_query = F, significant = T, exclude_iea = T, 
                measure_underrepresentation = F, evcodes = T, 
                user_threshold = 0.05, correction_method = "fdr", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)

write.table(apply(CR.gostres$result,2,as.character),
            "KO_CR_GO_res-by-waves.csv", sep = ";", quote = F, row.names = F)

DNA_damage_GOterm <- CR.gostres$result[CR.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977","GO:0033554",
                                                                                 "GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]

DNA_damage_GOterm[,c(9,1,2,3,5,6,7,11)]

## # A tibble: 0 × 8
## # … with 8 variables: term_id <chr>, query <chr>, significant <lgl>,
## #   p_value <dbl>, query_size <int>, intersection_size <int>, precision <dbl>,
## #   term_name <chr>

Go term on all CR genes

CR.gostres <- gost(query = as.character(CR.Gene.dynamique$Gene),
                organism = "mmusculus", ordered_query = F, 
                multi_query = F, significant = T, exclude_iea = T, 
                measure_underrepresentation = F, evcodes = T, 
                user_threshold = 0.05, correction_method = "fdr", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)

write.table(apply(CR.gostres$result,2,as.character),
            "KOCR_GO_res.csv", sep = ";", quote = F, row.names = F)

DNA_damage_GOterm <- CR.gostres$result[CR.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977",
                                                                                 "GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]

DNA_damage_GOterm[,c(1,2,3,5,6,7,11)]

## # A tibble: 0 × 7
## # … with 7 variables: query <chr>, significant <lgl>, p_value <dbl>,
## #   query_size <int>, intersection_size <int>, precision <dbl>, term_name <chr>

Intersection with ChP dynamicaly expressed genes

ChP_dynamic_genes <- read.table("../ChoroidPlexus_trajectory/ChP.Gene.dynamique.csv", sep = ";", header = T, row.names = 1)

CR_ChP_common_genes <- CR.Gene.dynamique %>% filter(Gene %in% ChP_dynamic_genes$Gene)

gene.order2 <- gene.order[gene.order %in% CR_ChP_common_genes$Gene]

pheatmap::pheatmap(CR_curve_matrix[gene.order2,],
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   #annotation_row = CR.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   #annotation_col = col.anno,
                   #annotation_colors = anno.colors,
                   gaps_col = 100,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

CR_ChP_common.gostres <- gost(query = list("Clust.1" = CR_ChP_common_genes %>% filter(Gene.Clusters == "Clust.1") %>% pull(Gene) %>% as.character(),
                             "Clust.2" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.2") %>% pull(Gene) %>% as.character(),
                             "Clust.3" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.3") %>% pull(Gene) %>% as.character(),
                             "Clust.4" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.4") %>% pull(Gene) %>% as.character(),
                             "Clust.5" = CR.Gene.dynamique %>% filter(Gene.Clusters == "Clust.5") %>% pull(Gene) %>% as.character()),
                organism = "mmusculus", ordered_query = F, 
                multi_query = F, significant = T, exclude_iea = T, 
                measure_underrepresentation = F, evcodes = T, 
                user_threshold = 0.05, correction_method = "fdr", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)

write.table(apply(CR_ChP_common.gostres$result,2,as.character),
            "KO_CR_ChP_common_GO_res-by-waves.csv", sep = ";", quote = F, row.names = F)

DNA_damage_GOterm <- CR_ChP_common.gostres$result[CR_ChP_common.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977",
                                                                                       "GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]

DNA_damage_GOterm[,c(1,2,3,5,6,7,11)]

## # A tibble: 0 × 7
## # … with 7 variables: query <chr>, significant <lgl>, p_value <dbl>,
## #   query_size <int>, intersection_size <int>, precision <dbl>, term_name <chr>

CR_ChP_common.gostres <- gost(query = as.character(CR_ChP_common_genes$Gene),
                organism = "mmusculus", ordered_query = F, 
                multi_query = F, significant = T, exclude_iea = T, 
                measure_underrepresentation = F, evcodes = T, 
                user_threshold = 0.05, correction_method = "fdr", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)

write.table(apply(CR_ChP_common.gostres$result,2,as.character),
            "KO_CR_ChP_common_GO_res_all.csv", sep = ";", quote = F, row.names = F)

Pallial neurons trajectory analysis

# Extract Pallial neurons trajectory genes
Pal.res <- as.data.frame(pseudo.maturation.diff.filtered[pseudo.maturation.diff.filtered$ABCs < 0,])
Pal.genes <- row.names(Pal.res)

Pal_curve_matrix <- Pal_curve_matrix[Pal.genes, ]

Gene expression profiles along the trajectory

## Cluster gene by expression profiles
Pseudotime.genes.clusters <- cluster::pam(as.dist((1 - cor(Matrix::t(Pal_curve_matrix),method = "pearson"))), k= 5)

Pal.Gene.dynamique <- data.frame(Gene= names(Pseudotime.genes.clusters$clustering),
                             Waves= Pseudotime.genes.clusters$clustering,
                             Gene.Clusters = Pseudotime.genes.clusters$clustering,
                             q.val = Pal.res$pval,
                             ABCs= Pal.res$ABCs
                             ) %>% arrange(Gene.Clusters)

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

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

# Set annotation colors
pal <- wes_palette("Darjeeling1")
anno.colors <- list(lineage = c(Pallial_neurons="#026c9a", Cajal_Retzius="#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(300:1,#Pal
                                  301:600)], #CR
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = Pal.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=300)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

We manually correct the reordering so genes are aligned from top right to bottom left

pheatmap::pheatmap(Diff.curve_matrix[gene.order,
                                c(300:1,#Pal
                                  301:600)], #CR
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = Pal.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = data.frame(lineage = rep(c("Pallial_neurons","Cajal_Retzius"), each=300)),
                   annotation_colors = anno.colors,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

anno.colors <- list(Cell.state = c(Cycling_RG="#046c9a", Differentiating_cells="#ebcb2e"),
                    Gene.Clusters = c(Clust.1 =pal[1] , Clust.2=pal[2], Clust.3=pal[3], Clust.4=pal[4], Clust.5=pal[5]))

col.anno <- data.frame(Cell.state = rep(c("Differentiating_cells","Cycling_RG"), c(200,100)))
rownames(col.anno) <- 300:1

pheatmap::pheatmap(Pal_curve_matrix[gene.order,300:1],
                   scale = "row",
                   cluster_rows = F,
                   cluster_cols = F,
                   annotation_row = Pal.Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = col.anno,
                   annotation_colors = anno.colors,
                   gaps_col = 200,
                   show_colnames = F,
                   show_rownames = F,
                   fontsize_row = 8,
                   color =  viridis::viridis(9),
                   breaks = seq(-2.5,2.5, length.out = 9),
                   main = "")

Gene cluster trend

Plot.clust.trends(Neuro.trajectories,
                   Lineage = "Pallial_neurons",
                   Which.cluster = 1:5,
                   clust.list = Pseudotime.genes.clusters$clustering,
                   Smooth.method = "gam")

Pal.gostres <- gost(query = as.character(Pal.Gene.dynamique$Gene),
                organism = "mmusculus", ordered_query = F, 
                multi_query = F, significant = T, exclude_iea = T, 
                measure_underrepresentation = F, evcodes = T, 
                user_threshold = 0.05, correction_method = "fdr", 
                domain_scope = "annotated", custom_bg = NULL, 
                numeric_ns = "", sources = c("GO:MF", "GO:BP"), as_short_link = F)

write.table(apply(Pal.gostres$result, 2, as.character),
            "KO_Pal.gostres.csv", sep = ";", quote = F, row.names = F)

DNA_damage_GOterm <- Pal.gostres$result[Pal.gostres$result$term_id %in% c("GO:0008630", "GO:0030330", "GO:0031571", "GO:0006974", "GO:0006977",
                                                                                 "GO:0044773", "GO:0042771", "GO:0042770", "GO:2001021", "GO:1902229"),]

DNA_damage_GOterm[,c(1,2,3,5,6,7,11)]

## # A tibble: 7 × 7
##   query   significant p_value query_size intersection_size precision term_name  
##   <chr>   <lgl>         <dbl>      <int>             <int>     <dbl> <chr>      
## 1 query_1 TRUE        0.00461        160                 3    0.0188 mitotic G1…
## 2 query_1 TRUE        0.00649        160                 4    0.025  DNA damage…
## 3 query_1 TRUE        0.0132         160                 3    0.0188 regulation…
## 4 query_1 TRUE        0.0207         160                 3    0.0188 intrinsic …
## 5 query_1 TRUE        0.0211         160                 2    0.0125 DNA damage…
## 6 query_1 TRUE        0.0287         160                 5    0.0312 signal tra…
## 7 query_1 TRUE        0.0358         160                 4    0.025  intrinsic …

Session Info

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

## [1] "09 mai, 2022, 16,07"

#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] RCurl_1.98-1.5        sparsesvd_0.2         crayon_1.4.2         
##  [28] jsonlite_1.7.2        spatstat.data_2.1-0   survival_3.2-13      
##  [31] zoo_1.8-9             iterators_1.0.13      glue_1.5.1           
##  [34] polyclip_1.10-0       registry_0.5-1        gtable_0.3.0         
##  [37] leiden_0.3.9          future.apply_1.8.1    abind_1.4-5          
##  [40] scales_1.1.1          pheatmap_1.0.12       DBI_1.1.1            
##  [43] miniUI_0.1.1.1        Rcpp_1.0.8            viridisLite_0.4.0    
##  [46] xtable_1.8-4          reticulate_1.22       spatstat.core_2.3-1  
##  [49] htmlwidgets_1.5.4     httr_1.4.2            FNN_1.1.3            
##  [52] ellipsis_0.3.2        ica_1.0-2             farver_2.1.0         
##  [55] pkgconfig_2.0.3       sass_0.4.0            uwot_0.1.10          
##  [58] deldir_1.0-6          utf8_1.2.2            labeling_0.4.2       
##  [61] tidyselect_1.1.1      rlang_0.4.12          reshape2_1.4.4       
##  [64] later_1.3.0           munsell_0.5.0         tools_4.2.0          
##  [67] cli_3.1.0             generics_0.1.1        ggridges_0.5.3       
##  [70] evaluate_0.14         stringr_1.4.0         fastmap_1.1.0        
##  [73] yaml_2.2.1            goftest_1.2-3         knitr_1.36           
##  [76] fitdistrplus_1.1-6    purrr_0.3.4           RANN_2.6.1           
##  [79] pbapply_1.5-0         future_1.23.0         nlme_3.1-153         
##  [82] mime_0.12             slam_0.1-49           rstudioapi_0.13      
##  [85] compiler_4.2.0        plotly_4.10.0         png_0.1-7            
##  [88] spatstat.utils_2.2-0  tibble_3.1.6          bslib_0.3.1          
##  [91] stringi_1.7.6         highr_0.9             lattice_0.20-45      
##  [94] HSMMSingleCell_1.14.0 vctrs_0.3.8           pillar_1.6.4         
##  [97] lifecycle_1.0.1       spatstat.geom_2.3-0   combinat_0.0-8       
## [100] lmtest_0.9-39         jquerylib_0.1.4       RcppAnnoy_0.0.19     
## [103] bitops_1.0-7          data.table_1.14.2     httpuv_1.6.3         
## [106] patchwork_1.1.1       R6_2.5.1              promises_1.2.0.1     
## [109] TSP_1.1-11            KernSmooth_2.23-20    gridExtra_2.3        
## [112] parallelly_1.29.0     codetools_0.2-18      MASS_7.3-56          
## [115] assertthat_0.2.1      withr_2.4.3           qlcMatrix_0.9.7      
## [118] sctransform_0.3.2     mgcv_1.8-40           parallel_4.2.0       
## [121] grid_4.2.0            rpart_4.1.16          tidyr_1.1.4          
## [124] rmarkdown_2.11        Rtsne_0.15            shiny_1.7.1

---

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