Temporal and spatial aspects of progenitors identity

Load libraries

# Load library
library(Seurat)
library(princurve)
library(caret)
library(monocle)
library(dplyr)
library(ggplot2)
library(pheatmap)
library(gridExtra)
library(ggExtra)
library(cowplot)
library(patchwork)
library(RColorBrewer)

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

Load Telley et al 2019 data and order cells along pseudo-maturation axis

To define temporally versus spatially regulated genes, we used dorso-pallial apical progenitors scRNAseq data obtained by facs sorting 1H after flashtag injection from :

Telley L, Agirman G et al. Temporal patterning of apical progenitors and their daughter neurons in the developing neocortex. Science 2019 May 10;364(6440).

The R script used for the preprocessing of the data can be find here

# Load the full annotated dataset
Telley.data <- readRDS("./Telley2019data/Telley2019.RDS")
Telley.data <- SetAllIdent(Telley.data, id = "orig.ident")

table(Telley.data@ident)

## 
## E12.1H E13.1H E14.1H E15.1H 
##    178    202    124    221

Run PCA

# Select all DEG between sampling days as input to PCA
Allmarkers <- FindAllMarkers(object = Telley.data,
                             min.pct = 0.3,
                             logfc.threshold = 0.6,
                             print.bar = F,
                             only.pos = T)

DE.genes <- unique(Allmarkers$gene)

# PCA with Seurat
Telley.data <- RunPCA(object = Telley.data,
                  pcs.compute = 10,
                  pc.genes = DE.genes,
                  do.print = F,
                  pcs.print = 1,
                  genes.print = 1,
                  weight.by.var=T)

PrintPCAParams(Telley.data)

## Parameters used in latest PCA calculation run on: 2020-11-30 11:03:06
## =============================================================================
## PCs computed    Genes used in calculation    PCs Scaled by Variance Explained
##     10                  815                                TRUE
## -----------------------------------------------------------------------------
## rev.pca 
##  FALSE
## -----------------------------------------------------------------------------
## Full gene list can be accessed using 
##  GetCalcParam(object = object, calculation = "RunPCA", parameter = "pc.genes")

# Select pcs by % variance explained
eigs <- (Telley.data@dr$pca@sdev)^2
PercentVAr <- data.frame(PercentVAr = eigs*100 / sum(eigs),
                         PCs = 1:length(eigs))

ggplot(PercentVAr,aes(x= as.factor(PCs),y=PercentVAr)) + geom_point() + geom_hline(yintercept = 5,linetype = 2, colour="red")

# Keep only PCs explaining > 5% of the variance over the 10 first
PCs <- PercentVAr %>% filter(PercentVAr > 5) %>% pull(PCs)

Fit the principal curve over the first 6 PCs

data <- as.data.frame(Telley.data@dr$pca@cell.embeddings[,PCs])

fit <- principal_curve(as.matrix(data[,PCs]),
                       smoother='lowess',
                       trace=TRUE,
                       f = .5,
                       stretch=0)

## Starting curve---distance^2: 20282698
## Iteration 1---distance^2: 26243.03
## Iteration 2---distance^2: 25451.13
## Iteration 3---distance^2: 24542.64
## Iteration 4---distance^2: 23847.34
## Iteration 5---distance^2: 23528.5
## Iteration 6---distance^2: 23464.76
## Iteration 7---distance^2: 23435.29
## Iteration 8---distance^2: 23404.87
## Iteration 9---distance^2: 23372.66
## Iteration 10---distance^2: 23359.4

PseudoMaturation.score <- fit$lambda/max(fit$lambda)
pc.line <- as.data.frame(fit$s[order(fit$lambda),])

data$Cluster <- Telley.data@meta.data$orig.ident
data$PseudoMaturation.score <- PseudoMaturation.score

Telley.data@meta.data$PseudoMaturation.score <- PseudoMaturation.score

# Plot the principal curve
p1 <- ggplot(data, aes(PC1, PC2)) + 
      geom_point(aes(color=Cluster), size=3, shape=16) +
      geom_line(data=pc.line, color='red', size=0.77) +
      scale_color_manual(values=c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896"))

# Plot pseudo-maturation score
cols <- colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)

p2 <- ggplot(data, aes(PC1, PC2)) +
     geom_point(aes(color=PseudoMaturation.score), size=3, shape=16) + 
     scale_color_gradientn(colours=rev(cols), name='Speudotime score') +
     geom_line(data=pc.line, color='red', size=0.77)

p1 + p2

# Plot density on speudotime
ggplot(data=data, aes(x=PseudoMaturation.score,fill=Cluster)) +
   scale_fill_manual(values=c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896")) +
   geom_density(alpha=.4) 

rm(list=ls()[!ls() %in% "Telley.data"])

Load E12 dataset

E12.AP.data <- readRDS("./AP.data.RDS")

Extract pallial progenitors

# Filter subpallial progenitors
AP.ident <- unique(E12.AP.data@ident)
E12.AP.data <-  SubsetData(E12.AP.data,
                          ident.use = grep("Sub", AP.ident, value = T, invert = T),
                          subset.raw = T,
                          do.clean = F)

# Reset the pseudo-DV score to [0,1]
score <- E12.AP.data@meta.data$DorsoVentral.Score
E12.AP.data@meta.data$DorsoVentral.Score <- (score - min(score)) / (max(score) - min(score))

p1 <- DimPlot(E12.AP.data,
               group.by = "Domaine",
               reduction.use = "spring",
               dim.1 = 1,
               dim.2 = 2,
               do.label=F,
               label.size = 4,
               no.legend = F,
               do.return = T,
               cols.use = tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B")))


data <- data.frame(spring1 = E12.AP.data@dr$spring@cell.embeddings[,1],
                   spring2 = E12.AP.data@dr$spring@cell.embeddings[,2],
                   DorsoVentral.Score = E12.AP.data@meta.data$DorsoVentral.Score)

cols <- colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)
p2 <- ggplot(data, aes(spring1, spring2)) + 
      geom_point(aes(color=DorsoVentral.Score), size=2, shape=16)  + 
      scale_color_gradientn(colours=rev(cols), name='Speudotime score')

p1 + p2

Filter, normalize and scale the expression matrix

# Remove non epressed genes
num.cells <- Matrix::rowSums(E12.AP.data@data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 20)])
E12.AP.data@raw.data <- E12.AP.data@raw.data[genes.use, ]
E12.AP.data@data <- E12.AP.data@data[genes.use, ]

# Normalize the counts matrix
E12.AP.data <- NormalizeData(object = E12.AP.data,
                         normalization.method = "LogNormalize", 
                         scale.factor = round(median(E12.AP.data@meta.data$nUMI)), display.progress = F)

# Scale expression matrix
E12.AP.data <- ScaleData(object = E12.AP.data,
                         vars.to.regress = c("CC.Difference", "percent.ribo" ,"nUMI"),
                         display.progress = F)


# Find most variable genes
E12.AP.data <- FindVariableGenes(object = E12.AP.data,
                                 mean.function = ExpMean,
                                 dispersion.function = LogVMR,
                                 x.low.cutoff = 0.0125,
                                 x.high.cutoff = 3,
                                 y.cutoff = 1,
                                 do.plot = F,
                                 display.progress = F)
length(E12.AP.data@var.genes)

## [1] 1129

rm(list=ls()[!ls() %in% c("E12.AP.data","Telley.data")])

Predict E12 AP pseudo-maturation score with a regression model trainned on Telley 2019 dataset

Train a random forest regression model on Telley 2019 data

# We trained the model using only genes found to be highly variable in our E12 dataset
model.data <- data.frame(PseudoMaturation.score = Telley.data@meta.data$PseudoMaturation.score)

genes <- t(Telley.data@scale.data[rownames(Telley.data@scale.data) %in% E12.AP.data@var.genes,])
genes <- as.data.frame(as.matrix(genes))
model.data <- cbind(model.data,genes)

# We did not perform cross-validatin but relied on the Out Of Bag error
trctrl <- trainControl(method = "none")

set.seed(1234)

# Train the random forest model
Telley.model <- train(PseudoMaturation.score~., 
                      data = model.data, 
                      method = "ranger",
                      trControl=trctrl,
                      importance="permutation", 
                      tuneGrid = data.frame(mtry=50, min.node.size = 5, splitrule="variance"))

We reoptimize the regression model by selecting the 100 most importante features of the first model

# Select the top 100 importante features
FeatureImportance <- sort(Telley.model$finalModel$variable.importance, decreasing = T)
SelectedFeatures <- names(FeatureImportance[1:100])

# Re-train the model with the restricted features
Telley.model <- train(PseudoMaturation.score~., 
                      data = model.data[,colnames(model.data) %in% c(SelectedFeatures,"PseudoMaturation.score")], 
                      method = "ranger",
                      trControl=trctrl,
                      importance="permutation", 
                      tuneGrid = data.frame(mtry=50, min.node.size = 5, splitrule="variance"))

# Plot Observed vs OOB predicted values from the model
plot(model.data$PseudoMaturation.score, Telley.model$finalModel$predictions,
     pch=19, xlab="observed Pseudo-maturation",
     ylab="OOB predicted Pseudo-maturation")
mtext(paste("R-squared", format(Telley.model$finalModel$r.squared,digits=2)))

Predict maturation score of E12 pallial apical progenitors

# Predict maturation score of E12 APs
E12.df <- as.data.frame(t(E12.AP.data@scale.data[rownames(E12.AP.data@scale.data) %in% SelectedFeatures,]))

E12.pred.maturation <- predict(Telley.model, E12.df)

E12.AP.data@meta.data$PseudoMaturation.score <- E12.pred.maturation

E12.plot.df <- E12.AP.data@meta.data %>% dplyr::select(Barcodes, Domaine, DorsoVentral.Score,PseudoMaturation.score)

# Plot E12 predicted maturation score by cluster
domaine <- factor(E12.plot.df$Domaine, level = c("Ventral.Pallium","lateral.Pallium.1","lateral.Pallium.2","Dorsal.Pallium"))

ggplot(E12.plot.df, aes(x= domaine, y= PseudoMaturation.score, fill=Domaine)) +
  geom_boxplot(notch=TRUE) +
  geom_jitter(color="black", size=0.8, alpha=0.9) +
  scale_fill_manual(values = tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B"))) +
  xlab("")

Manuscript Fig. 6B

rm(list = ls()[!ls() %in% c("Telley.data", "E12.AP.data")])

Decomposition of spatial vs. temoral aspect of progenitor cells identity

Smooth expression over pseudo-DV and pseudo maturation

Find all variable genes along the pseudo-DV axis in E12 AP

# Transfert metadata 
meta.data <- data.frame(barcode = rownames(E12.AP.data@meta.data),
                        Cluster = E12.AP.data@meta.data$old.ident,
                        DorsoVentral.Score =  E12.AP.data@meta.data$DorsoVentral.Score,
                        CellcyclePhase = E12.AP.data@meta.data$Phase,
                        row.names = rownames(E12.AP.data@meta.data))
                   
Annot.data  <- new('AnnotatedDataFrame', data = meta.data)

# Transfert count data
count.data = data.frame(gene_short_name = rownames(E12.AP.data@raw.data),
                  row.names = rownames(E12.AP.data@raw.data))

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

# Create the CellDataSet object
E12.gbm_cds <- newCellDataSet(as.matrix(E12.AP.data@raw.data),
                          phenoData = Annot.data,
                          featureData = feature.data,
                          lowerDetectionLimit = 1,
                          expressionFamily = negbinomial())

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

# Exclude cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
Input.genes <- E12.AP.data@var.genes[!E12.AP.data@var.genes %in% CCgenes]

# Perform the test for differential expression as a function of pseudo-DV score
DV.Axis.genes <- differentialGeneTest(E12.gbm_cds[Input.genes,], 
                                      fullModelFormulaStr = "~sm.ns(DorsoVentral.Score, df = 3)*CellcyclePhase", 
                                      reducedModelFormulaStr = "~CellcyclePhase", 
                                      cores = detectCores() -2)

# Filter genes with a FDR < 0.001
DV.Axis.genes.FDR.filtered <- DV.Axis.genes %>% filter(qval < 1e-3)

# Create a new pseudo-DV vector of 500 points
nPoints <- 200
E12.new_data <- data.frame(DorsoVentral.Score = seq(min(pData(E12.gbm_cds)$DorsoVentral.Score), max(pData(E12.gbm_cds)$DorsoVentral.Score), length.out = nPoints))

# Smooth gene expression
E12.smooth.curve.matrix <- genSmoothCurves(E12.gbm_cds[as.character(DV.Axis.genes.FDR.filtered$gene_short_name),],
                                           trend_formula = "~sm.ns(DorsoVentral.Score, df = 3)",
                                           relative_expr = TRUE,
                                           new_data = E12.new_data,
                                           cores= detectCores() - 2)

rm(list=ls()[!ls() %in% c("Telley.data", "E12.AP.data", "E12.smooth.curve.matrix", "E12.new_data", "DV.Axis.genes.FDR.filtered")])

Find all variable genes along the pseudo-maturation axis in Telley 2019 APs

# Transfert metadata 
meta.data <- data.frame(barcode = rownames(Telley.data@meta.data),
                        Cluster = Telley.data@meta.data$orig.ident,
                        PseudoMaturation.score =  Telley.data@meta.data$PseudoMaturation.score,
                        CellcyclePhase = Telley.data@meta.data$Phase,
                        row.names = rownames(Telley.data@meta.data))
                   
Annot.data  <- new('AnnotatedDataFrame', data = meta.data)

# Transfert count data
count.data = data.frame(gene_short_name = rownames(Telley.data@raw.data),
                  row.names = rownames(Telley.data@raw.data))

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

# Create the CellDataSet object
Telley.gbm_cds <- newCellDataSet(as.matrix(Telley.data@raw.data),
                                 phenoData = Annot.data,
                                 featureData = feature.data,
                                 lowerDetectionLimit = 1,
                                 expressionFamily = negbinomial())

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

# Exclude cell cycle associated, mt and ribo ribo
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
GenesToRemove <- c(grep(pattern = "(^Rpl|^Rps|^Mrp)", x = row.names(Telley.data@data), value = TRUE),
                   grep(pattern = "^mt-", x = row.names(Telley.data@data), value = TRUE),
                   "Xist", CCgenes)

Input.genes <- row.names(Telley.data@data)[!row.names(Telley.data@data) %in% GenesToRemove]

# We select all genes which are differentialy expressed between E12 and E15 AP 
E12vsE15.DEgenes <- FindMarkers(object = Telley.data,
                                ident.1 = "E12.1H",
                                ident.2 = "E15.1H",
                                min.pct = 0.25,
                                logfc.threshold = 0.4,
                                print.bar = F,
                                only.pos = F,
                                genes.use = Input.genes)

E12vsE15.DEgenes$Gene <- rownames(E12vsE15.DEgenes)

E12vsE15.DEgenes.filtered <- E12vsE15.DEgenes %>% filter(p_val_adj < 1e-3) %>% pull(Gene)

# Create a new pseudo-DV vector of 500 points
nPoints <- 200
Telley.new_data <- data.frame(PseudoMaturation.score = seq(min(pData(Telley.gbm_cds)$PseudoMaturation.score), max(pData(Telley.gbm_cds)$PseudoMaturation.score), length.out = nPoints))

# Smooth gene expression along pseudo-maturation axis
Telley.smooth.curve.matrix <- genSmoothCurves(Telley.gbm_cds[E12vsE15.DEgenes.filtered,],
                                              trend_formula = "~sm.ns(PseudoMaturation.score, df = 3)",
                                              relative_expr = TRUE,
                                              new_data = Telley.new_data,
                                              cores= detectCores() - 2)

Identify temporal genes as highly anti-correlating between the two datasets

Intersect the two smoothed expression matrices

CommonGenes <- intersect(rownames(E12.smooth.curve.matrix), rownames(Telley.smooth.curve.matrix))

E12.Smooth.Common <- t(scale(t(E12.smooth.curve.matrix[CommonGenes,])))
Telley.Smooth.Common <- t(scale(t(Telley.smooth.curve.matrix[CommonGenes,])))

Cluster temporal genes according to their smoothed expression inverse correlation between the two datasets

# Build temporal genes dendrogram
cor.mat <- cor(x= t(E12.Smooth.Common[,200:1]), y= t(Telley.Smooth.Common), method = "pearson")
dst <- dist(cor.mat, method = "euclidean")
hc <-  hclust(dst, method = "complete")

# Find 3 cluster based on hierarchical clustering
clusters <- cutree(hc, k=3)

Temporal.gene.clusters <- data.frame(Gene = names(clusters),
                             Gene.Clusters = as.numeric(clusters)) %>% arrange(Gene.Clusters)

row.names(Temporal.gene.clusters) <- Temporal.gene.clusters$Gene
Temporal.gene.clusters$Gene.Clusters <- paste0("Clust.",Temporal.gene.clusters$Gene.Clusters)

write.table(Temporal.gene.clusters, "./Progenitors/Temporal.gene.clusters.csv", sep = ";", quote = F)

pheatmap(cor.mat,
         cluster_rows= hc,
         scale = "none",
         annotation_row = Temporal.gene.clusters %>% dplyr::select(Gene.Clusters),
         annotation_colors = list(Gene.Clusters = c(Clust.1 ="#4784a2" , Clust.2="#ebcb2e",Clust.3="#cc3a1b")),
         show_colnames = T,
         show_rownames = T,
         border_color = NA,
         color = rev(brewer.pal(8,"RdBu")),
         main = "Genes inverse pearson correlation matrix")

Manuscript Fig. 6E

Plot correlation between smoothed expression profiles

Smoothed.point.cor <- cor(x= E12.Smooth.Common, y= Telley.Smooth.Common[,200:1], method = "pearson")

Domaine.E12.AP <- E12.AP.data@meta.data %>% select(DorsoVentral.Score, Domaine)
Domaine.boundaries <- aggregate(DorsoVentral.Score ~ Domaine, Domaine.E12.AP, max) %>% arrange(DorsoVentral.Score) %>% pull(DorsoVentral.Score)

Domaine.Ident <- sapply(E12.new_data$DorsoVentral.Score,
                        function(x){ if(x<Domaine.boundaries[1]){ x = "Ventral.Pallium"
                        } else if(x> Domaine.boundaries[1] & x< Domaine.boundaries[2]){ x ="lateral.Pallium.1"
                        } else if(x> Domaine.boundaries[3]){ x = "Dorsal.Pallium"
                        } else x="lateral.Pallium.2"}) 

Domaines.Clust <- data.frame(Domaine.Ident)

pheatmap(t(Smoothed.point.cor),
         cluster_rows = F,
         cluster_cols = F,
         annotation_col = Domaines.Clust,
         gaps_col = cumsum(as.numeric(table(Domaines.Clust$Domaine.Ident))[c(4,2,3,1)]),
         annotation_colors = list(Domaine.Ident = c(Ventral.Pallium = "#e46b6b",lateral.Pallium.1 = "#e3c148", lateral.Pallium.2 = "#b7d174",Dorsal.Pallium = "#68b041")),
         show_colnames = F,
         show_rownames = F,
         color = colorRampPalette(rev(brewer.pal(n = 10, name = "RdBu")))(100),
         main = "Pearson's Correlation matrix")

Manuscript Fig. 6D

Extract E12 spatial domain genes which do not show E12/E15 difference in Telley data

# Select genes variable in E12 dataset which do not show significative difference between E12 and E15 AP in Telley et al dataset
E12.Domain.specific.genes <- rownames(E12.smooth.curve.matrix)[!rownames(E12.smooth.curve.matrix) %in% CommonGenes]

E12.smooth.specific <- E12.smooth.curve.matrix[E12.Domain.specific.genes,]

# Build temporal genes dendrogram
cor.mat <- cor(t(E12.smooth.specific), method = "spearman")
dst <- as.dist(1 - cor.mat)
hc <-  hclust(dst, method = "ward.D") #"ward.D"

# Seriate the dendrogram
dend <- as.dendrogram(hc)
dend <- dendextend::seriate_dendrogram(dend, dst, method="GW")

# Find 2 cluster based on hierarchical clustering
clusters <- cutree(hc, k=4)

rownames(DV.Axis.genes.FDR.filtered)  <- DV.Axis.genes.FDR.filtered$gene_short_name
DV.Axis.genes.FDR.filtered <- DV.Axis.genes.FDR.filtered[names(clusters),]

Spatial.gene.clusters <- data.frame(Gene = names(clusters),
                                    pval= DV.Axis.genes.FDR.filtered$pval,
                                    qval= DV.Axis.genes.FDR.filtered$qval,
                                    num_cells_expressed= DV.Axis.genes.FDR.filtered$num_cells_expressed,
                                    Gene.Clusters = as.numeric(clusters)) %>% arrange(Gene.Clusters)

row.names(Spatial.gene.clusters) <- Spatial.gene.clusters$Gene
Spatial.gene.clusters$Gene.Clusters <- paste0("Clust.",Spatial.gene.clusters$Gene.Clusters)

write.table(Spatial.gene.clusters, "./Progenitors/Spatial.gene.clusters.csv", sep = ";", quote = F)

anno.col <- list(Domaine.Ident = c(Ventral.Pallium = "#e46b6b",lateral.Pallium.1 = "#e3c148", lateral.Pallium.2 = "#b7d174",Dorsal.Pallium = "#68b041"),
                 Gene.Clusters = c(Clust.1 ="#4784a2" , Clust.2="#ebcb2e", Clust.3="#cc3a1b",Clust.4="#4cabdc",Clust.5="#3ca73f"))

pheatmap(E12.smooth.specific,
         cluster_rows= as.hclust(dend),
         scale = "row",
         cluster_cols = F,
         annotation_col = Domaines.Clust,
         gaps_col = cumsum(as.numeric(table(Domaines.Clust$Domaine.Ident))[c(4,2,3,1)]),
         annotation_row = Spatial.gene.clusters %>% dplyr::select(Gene.Clusters),
         annotation_colors = anno.col,
         show_colnames = F,
         show_rownames = F,
         color = rev(brewer.pal(11,"RdBu")),
         breaks = seq(-2.5,2.5, length.out = 11),
         main = "Spatial smoothed expression along pseudoDV axis")

Manuscript Fig. 6G

Plot the representative trends

dataE12 <- data.frame(Cluster = paste0("Cluster",as.character(E12.AP.data@ident)),
                   DorsoVentral.Score =  E12.AP.data@meta.data$DorsoVentral.Score,
                   Zbtb20 = E12.AP.data@data["Zbtb20",],
                   Lrrn1 = E12.AP.data@data["Lrrn1",],
                   Lmo4 = E12.AP.data@data["Lmo4",],
                   Mpped2 =E12.AP.data@data["Mpped2",])

dataTelley <- data.frame(Cluster = paste0("Cluster",as.character(Telley.data@ident)),
                         PseudoMaturation.score =  Telley.data@meta.data$PseudoMaturation.score,
                         Zbtb20 = Telley.data@data["Zbtb20",],
                         Lrrn1 = Telley.data@data["Lrrn1",],
                         Lmo4 = Telley.data@data["Lmo4",],
                         Mpped2 = Telley.data@data["Mpped2",])

ggplot(data=dataE12, aes(x=DorsoVentral.Score, y=Zbtb20, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B"))) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Zbtb20") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataE12, aes(x=DorsoVentral.Score, y=Lrrn1, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B"))) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Lrrn1") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataE12, aes(x=DorsoVentral.Score, y=Lmo4, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B"))) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Lmo4") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataE12, aes(x=DorsoVentral.Score, y=Mpped2, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B"))) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Mpped2") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataTelley, aes(x=PseudoMaturation.score, y=Zbtb20, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896")) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey",span = 0.5) + 
    ggtitle("Zbtb20") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataTelley, aes(x=PseudoMaturation.score, y=Lrrn1, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896")) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Lrrn1") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataTelley, aes(x=PseudoMaturation.score, y=Lmo4, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896")) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Lmo4") +
    ylim(0,NA) +
    theme(legend.position="none")

ggplot(data=dataTelley, aes(x=PseudoMaturation.score, y=Mpped2, color=Cluster)) +
    geom_point() +
    scale_color_manual(values= c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896")) + 
    geom_smooth(method="loess", n= 30, color="red", fill="grey") + 
    ggtitle("Mpped2") +
    ylim(0,NA) +
    theme(legend.position="none")

Session Info

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

## [1] "30 novembre, 2020, 11,09"

#Packages used
sessionInfo()

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.5 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/atlas/libblas.so.3.10.3
## LAPACK: /usr/lib/x86_64-linux-gnu/atlas/liblapack.so.3.10.3
## 
## 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    parallel  stats     graphics  grDevices utils    
##  [8] datasets  methods   base     
## 
## other attached packages:
##  [1] RColorBrewer_1.1-2  patchwork_0.0.1     ggExtra_0.9        
##  [4] gridExtra_2.3       pheatmap_1.0.12     dplyr_0.8.3        
##  [7] monocle_2.14.0      DDRTree_0.1.5       irlba_2.3.3        
## [10] VGAM_1.1-2          Biobase_2.46.0      BiocGenerics_0.32.0
## [13] caret_6.0-84        lattice_0.20-41     princurve_2.1.4    
## [16] Seurat_2.3.4        Matrix_1.2-17       cowplot_1.0.0      
## [19] ggplot2_3.2.1      
## 
## loaded via a namespace (and not attached):
##   [1] snow_0.4-3           backports_1.1.5      Hmisc_4.3-0         
##   [4] plyr_1.8.4           igraph_1.2.5         lazyeval_0.2.2      
##   [7] densityClust_0.3     fastICA_1.2-2        digest_0.6.25       
##  [10] foreach_1.4.7        htmltools_0.5.0      viridis_0.5.1       
##  [13] lars_1.2             gdata_2.18.0         magrittr_1.5        
##  [16] checkmate_1.9.4      cluster_2.1.0        gclus_1.3.2         
##  [19] mixtools_1.1.0       ROCR_1.0-7           limma_3.42.0        
##  [22] recipes_0.1.7        gower_0.2.1          matrixStats_0.55.0  
##  [25] docopt_0.6.1         R.utils_2.9.0        colorspace_1.4-1    
##  [28] ggrepel_0.8.1        xfun_0.18            sparsesvd_0.2       
##  [31] crayon_1.3.4         jsonlite_1.7.0       zeallot_0.1.0       
##  [34] survival_2.44-1.1    zoo_1.8-6            iterators_1.0.12    
##  [37] ape_5.3              glue_1.4.1           registry_0.5-1      
##  [40] gtable_0.3.0         ipred_0.9-9          kernlab_0.9-29      
##  [43] prabclus_2.3-1       DEoptimR_1.0-8       scales_1.1.0        
##  [46] bibtex_0.4.2         miniUI_0.1.1.1       Rcpp_1.0.5          
##  [49] metap_1.1            dtw_1.21-3           xtable_1.8-4        
##  [52] viridisLite_0.3.0    htmlTable_1.13.2     reticulate_1.13     
##  [55] foreign_0.8-72       bit_4.0.4            proxy_0.4-23        
##  [58] mclust_5.4.5         SDMTools_1.1-221.1   Formula_1.2-3       
##  [61] tsne_0.1-3           lava_1.6.6           prodlim_2019.11.13  
##  [64] htmlwidgets_1.5.1    httr_1.4.1           FNN_1.1.3           
##  [67] gplots_3.0.1.1       fpc_2.2-3            acepack_1.4.1       
##  [70] modeltools_0.2-22    ica_1.0-2            farver_2.0.1        
##  [73] pkgconfig_2.0.3      R.methodsS3_1.7.1    flexmix_2.3-15      
##  [76] nnet_7.3-14          labeling_0.3         later_1.0.0         
##  [79] tidyselect_0.2.5     rlang_0.4.7          reshape2_1.4.3      
##  [82] munsell_0.5.0        tools_3.6.3          generics_0.0.2      
##  [85] ranger_0.11.2        ggridges_0.5.1       fastmap_1.0.1       
##  [88] evaluate_0.14        stringr_1.4.0        yaml_2.2.1          
##  [91] npsurv_0.4-0         ModelMetrics_1.2.2   knitr_1.26          
##  [94] bit64_4.0.2          fitdistrplus_1.0-14  robustbase_0.93-5   
##  [97] caTools_1.17.1.2     purrr_0.3.3          RANN_2.6.1          
## [100] dendextend_1.12.0    pbapply_1.4-2        nlme_3.1-141        
## [103] mime_0.7             slam_0.1-46          R.oo_1.23.0         
## [106] hdf5r_1.3.2.9000     compiler_3.6.3       rstudioapi_0.11     
## [109] png_0.1-7            e1071_1.7-2          lsei_1.2-0          
## [112] tibble_2.1.3         stringi_1.4.6        highr_0.8           
## [115] HSMMSingleCell_1.6.0 vctrs_0.2.0          pillar_1.4.2        
## [118] lifecycle_0.1.0      combinat_0.0-8       Rdpack_0.11-0       
## [121] lmtest_0.9-37        data.table_1.12.6    bitops_1.0-6        
## [124] seriation_1.2-9      gbRd_0.4-11          httpuv_1.5.2        
## [127] R6_2.4.1             latticeExtra_0.6-28  TSP_1.1-10          
## [130] promises_1.1.0       KernSmooth_2.23-15   codetools_0.2-16    
## [133] MASS_7.3-53          gtools_3.8.1         assertthat_0.2.1    
## [136] withr_2.1.2          qlcMatrix_0.9.7      diptest_0.75-7      
## [139] doSNOW_1.0.18        grid_3.6.3           rpart_4.1-15        
## [142] timeDate_3043.102    tidyr_1.0.0          class_7.3-17        
## [145] rmarkdown_2.5        segmented_1.0-0      Rtsne_0.15          
## [148] shiny_1.4.0          lubridate_1.7.4      base64enc_0.1-3

---

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