# Load library
library(Seurat)
library(princurve)
library(monocle)
library(FateID)
library(velocyto.R)
library(parallel)
library(seriation)
library(dplyr)
library(reshape2)
library(ggplot2)
library(cowplot)
library(ggExtra)
library(patchwork)
library(RColorBrewer)
library(wesanderson)
#Set ggplot theme as classic
theme_set(theme_classic())
colors <- c("#969696",
tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B")),
"#ec756d", "#c773a7", "#7293c8", "#b79f0b", "#3ca73f","#31b6bd",
"#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")
DimPlot(Allcells.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 2,
no.legend = F,
cols.use = colors)
We perform K-means clustering on the 4 cell state scores :
SP
Pal
SP.BP
Pal.BP
# Calculate Pallial and Sub-pallial BP scores
Pal.BP.genes <- c("Eomes", "Neurog2", "Neurog1", "Prmt8", "Nrp1")
genes.list <- list(Pal.BP.genes)
enrich.name <- "Pal.BP_signature"
Allcells.data <- AddModuleScore(Allcells.data, genes.list = genes.list, genes.pool = NULL, n.bin = 5,
seed.use = 1, ctrl.size = length(genes.list), use.k = FALSE, enrich.name = enrich.name ,
random.seed = 1)
SP.BP.genes <- c("Dlx1", "Dlx2", "Dlx5","Ascl1", "Gsx2")
genes.list <- list(SP.BP.genes)
enrich.name <- "SP.BP_signature"
Allcells.data <- AddModuleScore(Allcells.data, genes.list = genes.list, genes.pool = NULL, n.bin = 5,
seed.use = 1, ctrl.size = length(genes.list), use.k = FALSE, enrich.name = enrich.name ,
random.seed = 1)
set.seed(100)
# Run K-means clustering
cl <- kmeans(cbind(Allcells.data@meta.data$SP_signature1,
Allcells.data@meta.data$Pal_signature1,
Allcells.data@meta.data$SP.BP_signature1,
Allcells.data@meta.data$Pal.BP_signature1), 4)
Allcells.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
col.pal <- wes_palette("GrandBudapest1", 4, type = "discrete")
p1 <- ggplot(Allcells.data@meta.data, aes(x=SP_signature1, y=Pal_signature1, colour = kmeanClust)) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey")
DimPlot(Allcells.data,
group.by = "kmeanClust",
reduction.use = "spring",
cols.use = col.pal,
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
no.legend = F)
We then extract the Pallial cells branch. We also excludes CR cells cluster form the trajectory inference.
# Remove the sub-pallial cells branch
MeanKclust.SPscore <- aggregate(SP_signature1 ~ kmeanClust, Allcells.data@meta.data, mean)
SPclust <- MeanKclust.SPscore %>% filter(SP_signature1 == max(SP_signature1)) %>% pull(kmeanClust)
SP.cells <- Allcells.data@meta.data %>% filter(kmeanClust == SPclust) %>% pull(Barcodes)
# Remove cells not use for trajectory inference
Excluded.clusters <- Allcells.data@meta.data %>%
filter(Cluster.ident %in% grep("*Sub|GABA|LN.Glut.13|LN.Glut.14|LN.Glut.1$", unique(as.character(Allcells.data@ident)), value = T)) %>%
pull(Barcodes)
# We further keep only pallial apical progenitor clusters
MeanKclust.APscore <- aggregate(AP_signature1 ~ kmeanClust, Allcells.data@meta.data, mean)
APclust <- MeanKclust.APscore %>% filter(AP_signature1 == max(AP_signature1)) %>% pull(kmeanClust)
All.AP <- Allcells.data@meta.data %>% filter(kmeanClust == APclust) %>% pull(Barcodes)
Valide.AP <- Allcells.data@meta.data %>% filter(Cluster.ident %in% grep("Dorsal.Pallium|lateral.Pallium.1|lateral.Pallium.2|Ventral.Pallium",
unique(as.character(Allcells.data@ident)), value = T)) %>% pull(Barcodes)
filtered.AP <- All.AP[!All.AP %in% Valide.AP]
# Remove all invalide cells + 3 pallial outlier cells
Cells.to.use <- rownames(Allcells.data@meta.data)[!rownames(Allcells.data@meta.data) %in% unique(c(SP.cells, Excluded.clusters, filtered.AP, c("ATTTCTGCACGGCCAT" , "CAAGTTGCAAGCCTAT", "CTACATTGTAGCTGCC")))]
Allcells.data <- SubsetData(Allcells.data, cells.use = Cells.to.use, subset.raw = T, do.clean = F)
colors <- c("#969696",
tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B")),
"#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")
DimPlot(Allcells.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 2,
no.legend = T,
cols.use = colors)
# Filter genes
num.cells <- Matrix::rowSums(Allcells.data@raw.data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 20)])
Allcells.data@raw.data <- Allcells.data@raw.data[genes.use, ]
# Normalization and variable genes selection
Allcells.data <- NormalizeData(object = Allcells.data,
normalization.method = "LogNormalize",
scale.factor = round(median(Allcells.data@meta.data$nUMI)),
display.progress = F)
Allcells.data <- FindVariableGenes(object = Allcells.data,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.03,
x.high.cutoff = 4,
y.cutoff = 2, do.plot = F,
display.progress = F)
# Take the cluster id from Seurat analysis
cluster.label <- Allcells.data@ident
Cluster.ident <- as.character(Allcells.data@ident)
# Set color Palette for layout
palette <- c(scales::hue_pal()(length(unique(Cluster.ident))))
colorsident <- cbind(ident = sort(unique(Cluster.ident)),
colors = c("#969696",
tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B")),
"#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9") )
# Create annotation data.frame
Cells.Color.df <- data.frame(sample_name = row.names(Allcells.data@meta.data),
primary_type_label = as.character(Cluster.ident),
primary_type_color = as.character(colorsident[match(Cluster.ident, colorsident[,1]),2]))
cell.colors <- Cells.Color.df$primary_type_color
names(cell.colors) <- Cells.Color.df$sample_name
rm(Cells.Color.df,colorsident,palette)
# Load Loom file containing spliced VS unspliced transcripts count matrices
LoomPath <- "./velocyto/E12-WT.loom"
ldat <- read.loom.matrices(LoomPath)
## reading loom file via hdf5r...
BarcodesVelocity <- stringi::stri_sub(ldat$spliced@Dimnames[[2]],40,55)
ldat$spliced@Dimnames[[2]] <- BarcodesVelocity
ldat$unspliced@Dimnames[[2]] <- BarcodesVelocity
ldat$ambiguous@Dimnames[[2]] <- BarcodesVelocity
# Filter matrix
ldat$spliced <- ldat$spliced[,rownames(Allcells.data@meta.data)]
ldat$unspliced <- ldat$unspliced[,rownames(Allcells.data@meta.data)]
ldat$spliced <- ldat$spliced[!duplicated(ldat[["unspliced"]]@Dimnames[[1]]),]
ldat$unspliced <- ldat$unspliced[!duplicated(ldat[["unspliced"]]@Dimnames[[1]]),]
# Clean workspace
rm(list=ls()[!ls() %in% c("Allcells.data", "ldat", "cell.colors", "cluster.label")])
spliced <- ldat$spliced
unspliced <- ldat$unspliced
# Filter genes by cluster expression
Filtered.spliced <- filter.genes.by.cluster.expression(spliced,cluster.label,min.max.cluster.average = 0.3)
Filtered.unspliced <- filter.genes.by.cluster.expression(unspliced,cluster.label,min.max.cluster.average = 0.3)
# Filter Cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
Filtered.spliced <- Filtered.spliced[!rownames(Filtered.spliced) %in% CCgenes,]
Filtered.unspliced <- Filtered.unspliced[!rownames(Filtered.unspliced) %in% CCgenes,]
# Velocity estimation
fit.quantile <- 0.04
rvel.cd <- gene.relative.velocity.estimates(Filtered.spliced,
Filtered.unspliced,
deltaT=1,
kCells = 8,
kGenes = 10,
cell.dist=cell.dist,
fit.quantile=fit.quantile,
n.cores = detectCores() -2)
## calculating cell knn ... done
## calculating convolved matrices ... done
## gene kNN ... scaling gene weights ... convolving matrices ... done
## fitting gamma coefficients ... done. succesfful fit for 2423 genes
## filtered out 3 out of 2423 genes due to low nmat-emat correlation
## filtered out 10 out of 2420 genes due to low nmat-emat slope
## calculating RNA velocity shift ... re-estimating gamma of individual genes ... done
## calculating RNA velocity shift ... done
## calculating extrapolated cell state ... done
# Velocity on embedding
emb <- Allcells.data@dr$spring@cell.embeddings/50
Spring.velo <- show.velocity.on.embedding.cor(emb,
rvel.cd,
n=50,
scale='sqrt',
cex=0.8,
arrow.scale=0.7,
show.grid.flow=T,
min.grid.cell.mass=0.5,
grid.n=40,
arrow.lwd=1,
do.par=F,
cell.border.alpha = 0,
cell.colors=cell.colors,
expression.scaling = T,
return.details= T,
n.core= detectCores() -2)
## delta projections ... sqrt knn ... transition probs ... done
## calculating arrows ... done
## grid estimates ... grid.sd= 0.288576 min.arrow.size= 0.005771519 max.grid.arrow.length= 0.04106169 done
## expression shifts .... done
Spring.plot <- as.data.frame(Allcells.data@dr$spring@cell.embeddings/50) %>%
mutate(x0 = Spring.velo$arrows[, "x0"],
x1 = Spring.velo$arrows[, "x1"],
y0 = Spring.velo$arrows[, "y0"],
y1 = Spring.velo$arrows[, "y1"]) %>%
mutate(x2 = x0 + (x1 - x0),
y2 = y0 + (y1 - y0)) %>%
mutate(Cluster = Allcells.data@meta.data$Cluster.ident)
colors <- c("#969696",
tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B")),
"#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")
ggplot(Spring.plot) +
geom_point(aes(x = spring1, y = spring2, colour = Cluster)) +
scale_color_manual(values = colors) +
geom_segment(aes(x = x0, xend = x2, y = y0, yend = y2),
arrow = arrow(length = unit(3, "points"), type = "closed"),
colour = "grey20", alpha = 0.8) +
theme(legend.position="none")
global.vel.arrow <- Spring.velo$garrows %>%
as.data.frame() %>%
mutate(x2 = x0 + (x1 - x0),
y2 = y0 + (y1 - y0))
ggplot(Spring.plot) +
geom_point(aes(x = spring1, y = spring2, colour = Cluster)) +
scale_color_manual(values = colors) +
geom_segment(data = global.vel.arrow,
aes(x = x0, xend = x2, y = y0, yend = y2),
size = 1,
arrow = arrow(length = unit(3, "points"), type = "open"),
colour = "grey20", alpha = 0.8) +
theme(legend.position="none")
# FateID requieres the terminal clusters ID to be set as integers
TerminalFates <- grep("16|24|22|19|26|20|21", unique(as.character(Allcells.data@ident)), value = T)
Allcells.data@meta.data$NewClusterID <- sapply(as.character(Allcells.data@ident),
function(x) if(x == TerminalFates[1]){x=1}
else if(x== TerminalFates[2]){x=2}
else if(x== TerminalFates[3]){x=3}
else if(x== TerminalFates[4]){x=4}
else if(x== TerminalFates[5]){x=5}
else if(x== TerminalFates[6]){x=6}
else if(x== TerminalFates[7]){x=7}
else{x=x})
Allcells.data <- SetAllIdent(Allcells.data, id = "NewClusterID")
colors <- c("#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b","#cc3a1b",
"#969696",
tolower(c("#68B041", "#E3C148", "#B7D174", "#E46B6B")),
"#046c9a", "#4784a2" , "#4990c9")
DimPlot(Allcells.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
no.legend = F,
cols.use = colors)
We restricted the analysis to the most variable genes as dertermined by the Seurat function “FindVariableGenes” excluding cell cycle phase genes
# Exclude cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
Input.genes <- Allcells.data@var.genes[!Allcells.data@var.genes %in% CCgenes]
Norm.Mat <- as.data.frame(as.matrix(Allcells.data@data[Input.genes,]))
# Set a cluster assignment factor for each cells
ClusterIdent <- as.character(Allcells.data@ident)
names(ClusterIdent) <- names(Allcells.data@ident)
Attractors <- 1:7
Infered.Fate.bias <- fateBias(Norm.Mat, ClusterIdent, Attractors,
z = 1-cor(Norm.Mat, method = "spearman"),
minnr=10,
minnrh=30,
adapt=TRUE,
confidence=0.95,
nbfactor=5,
use.dist=FALSE,
seed=1234,
nbtree=NULL)
## minnr: 10
## minnrh: 30
## test set size iteration 1 : 10 10 10 10 10 10 10
## randomforest iteration 1 of 50 cells
## test set size iteration 2 : 10 10 10 10 10 10 10
## randomforest iteration 2 of 62 cells
## test set size iteration 3 : 10 10 10 10 10 10 10
## randomforest iteration 3 of 53 cells
## test set size iteration 4 : 10 10 10 10 10 10 10
## randomforest iteration 4 of 55 cells
## test set size iteration 5 : 10 10 10 10 10 10 10
## randomforest iteration 5 of 49 cells
## test set size iteration 6 : 10 10 10 10 10 10 10
## randomforest iteration 6 of 48 cells
## test set size iteration 7 : 10 10 10 10 10 10 10
## randomforest iteration 7 of 54 cells
## test set size iteration 8 : 10 10 10 10 10 10 10
## randomforest iteration 8 of 49 cells
## test set size iteration 9 : 10 10 10 10 10 10 10
## randomforest iteration 9 of 50 cells
## test set size iteration 10 : 10 10 10 10 10 10 10
## randomforest iteration 10 of 51 cells
## test set size iteration 11 : 10 10 10 10 10 10 10
## randomforest iteration 11 of 52 cells
## test set size iteration 12 : 10 10 10 10 10 10 10
## randomforest iteration 12 of 44 cells
## test set size iteration 13 : 10 10 10 10 10 10 10
## randomforest iteration 13 of 45 cells
## test set size iteration 14 : 10 10 10 10 10 10 10
## randomforest iteration 14 of 42 cells
## test set size iteration 15 : 10 10 10 10 10 10 10
## randomforest iteration 15 of 43 cells
## test set size iteration 16 : 10 10 10 10 10 10 10
## randomforest iteration 16 of 40 cells
## test set size iteration 17 : 10 10 10 10 10 10 10
## randomforest iteration 17 of 40 cells
## test set size iteration 18 : 10 10 10 10 10 10 10
## randomforest iteration 18 of 37 cells
## test set size iteration 19 : 10 10 10 10 10 10 10
## randomforest iteration 19 of 42 cells
## test set size iteration 20 : 10 10 10 10 10 10 10
## randomforest iteration 20 of 43 cells
## test set size iteration 21 : 10 10 10 10 10 10 10
## randomforest iteration 21 of 48 cells
## test set size iteration 22 : 10 10 10 10 10 10 10
## randomforest iteration 22 of 47 cells
## test set size iteration 23 : 10 10 10 10 10 10 10
## randomforest iteration 23 of 50 cells
## test set size iteration 24 : 10 10 10 10 10 10 10
## randomforest iteration 24 of 52 cells
## test set size iteration 25 : 10 10 10 10 10 10 10
## randomforest iteration 25 of 53 cells
## test set size iteration 26 : 10 10 10 10 10 10 10
## randomforest iteration 26 of 57 cells
## test set size iteration 27 : 10 10 10 10 10 10 10
## randomforest iteration 27 of 54 cells
## test set size iteration 28 : 10 10 10 10 10 10 10
## randomforest iteration 28 of 49 cells
## test set size iteration 29 : 10 10 10 10 10 10 10
## randomforest iteration 29 of 54 cells
## test set size iteration 30 : 10 10 10 10 10 10 10
## randomforest iteration 30 of 51 cells
## test set size iteration 31 : 10 10 10 10 10 10 10
## randomforest iteration 31 of 50 cells
## test set size iteration 32 : 10 10 10 10 10 10 10
## randomforest iteration 32 of 55 cells
## test set size iteration 33 : 5 5 5 5 5 5 10
## randomforest iteration 33 of 36 cells
## test set size iteration 34 : 10 10 10 10 10 10 10
## randomforest iteration 34 of 56 cells
## test set size iteration 35 : 10 10 10 10 10 10 10
## randomforest iteration 35 of 55 cells
## test set size iteration 36 : 10 10 10 10 10 10 10
## randomforest iteration 36 of 55 cells
## test set size iteration 37 : 10 10 10 10 10 10 10
## randomforest iteration 37 of 54 cells
## test set size iteration 38 : 10 10 10 10 10 10 10
## randomforest iteration 38 of 56 cells
## test set size iteration 39 : 10 10 10 10 10 10 10
## randomforest iteration 39 of 58 cells
## test set size iteration 40 : 10 10 10 10 10 10 10
## randomforest iteration 40 of 47 cells
## test set size iteration 41 : 10 10 10 10 10 10 10
## randomforest iteration 41 of 44 cells
## test set size iteration 42 : 10 10 10 10 10 10 10
## randomforest iteration 42 of 30 cells
## test set size iteration 43 : 10 10 10 10 10 10 10
## randomforest iteration 43 of 18 cells
## test set size iteration 44 : 10 10 10 10 10 10 10
## randomforest iteration 44 of 2 cells
Allcells.data@meta.data$FateID.iteration <- "0"
Allcells.data <- SetAllIdent(Allcells.data, id = "FateID.iteration")
for (i in seq(0, length(Infered.Fate.bias$rfl), by = 5)[-1]) {
iter <- seq(i-4,i)
Barcodes <- c()
for (j in iter) {
Barcodes <- c(Barcodes, names(Infered.Fate.bias$rfl[[j]]$test$predicted))
}
Allcells.data <- SetIdent(Allcells.data, cells.use = Barcodes, ident.use = paste0("iter ",iter[1]," to ", iter[4]))
}
DimPlot(Allcells.data,
reduction.use = "spring",
dim.1 = 1,dim.2 = 2,
do.label=T,
label.size = 3,
no.legend = F,
cols.use = c("#dfdfdf", colors))
probs <- Infered.Fate.bias$probs[,seq(length(Attractors))]
Allcells.data@meta.data$prob.Nr4a2 <- probs$t1
Allcells.data@meta.data$prob.Foxp2c <- probs$t2
Allcells.data@meta.data$prob.Ppp1r14c <- probs$t3
Allcells.data@meta.data$prob.Fezf1 <- probs$t4
Allcells.data@meta.data$prob.Foxp2a <- probs$t5
Allcells.data@meta.data$prob.Foxp2b <- probs$t6
Allcells.data@meta.data$prob.Nfix <- probs$t7
FeaturePlot(object = Allcells.data,
features.plot = c("prob.Nr4a2", "prob.Nfix",
"prob.Ppp1r14c","prob.Fezf1",
"prob.Foxp2a", "prob.Foxp2b", "prob.Foxp2c"),
cols.use = rev(RColorBrewer::brewer.pal(n = 11, name = "Spectral")),
reduction.use = "spring",
no.legend = T)
Manuscript Fig. 7A and S7A
New.data <- data.frame(barcode=Allcells.data@meta.data$Barcodes,
cluster= Allcells.data@meta.data$Cluster.ident,
spring1= Allcells.data@dr$spring@cell.embeddings[,1],
spring2= Allcells.data@dr$spring@cell.embeddings[,2],
Nr4a2.biased= Allcells.data@meta.data$prob.Nr4a2,
Nfix.biased= Allcells.data@meta.data$prob.Nfix,
Ppp1r14c.biased = Allcells.data@meta.data$prob.Ppp1r14c,
Fezf1.biased = Allcells.data@meta.data$prob.Fezf1)
New.data <- New.data %>% filter(!cluster %in% grep("26|20|21", unique(as.character(Allcells.data@ident)), value = T))
New.data$lineage.bias <- colnames(New.data[,5:8])[apply(New.data[,5:8],1,which.max)]
ggplot(New.data, aes(spring1, spring2, colour = lineage.bias)) +
scale_color_manual(values=c("#e7823a","#cc391b","#026c9a","#d14c8d")) +
geom_point()
Manuscript Fig. 7B
Allcells.data@meta.data$Nr4a2.biase <- ifelse(Allcells.data@meta.data$prob.Nr4a2 > 0.5 & abs(Allcells.data@meta.data$prob.Nfix- Allcells.data@meta.data$prob.Nr4a2) > 0.25, "Nr4a2.lineage", "Other.lineages" )
table(Allcells.data@meta.data$Nr4a2.biase)
##
## Nr4a2.lineage Other.lineages
## 617 1682
Allcells.data@meta.data$Nfix.biase <- ifelse(Allcells.data@meta.data$prob.Nfix > 0.5 & abs(Allcells.data@meta.data$prob.Nfix- Allcells.data@meta.data$prob.Nr4a2) > 0.25, "Nfix.lineage", "Other.lineages" )
table(Allcells.data@meta.data$Nfix.biase)
##
## Nfix.lineage Other.lineages
## 520 1779
Nr4a2 <- rownames(subset(Allcells.data@meta.data, Allcells.data@meta.data$Nr4a2.biase == "Nr4a2.lineage"))
Nfix <- rownames(subset(Allcells.data@meta.data, Allcells.data@meta.data$Nfix.biase == "Nfix.lineage"))
Allcells.data@meta.data$lineage.bias <- sapply(Allcells.data@meta.data$Barcodes,
function(x) if(x %in% Nr4a2){x="Nr4a2.lineage"}
else if(x %in% Nfix ){x= "Nfix.lineage"}
else{x= "Not.assigned"})
DimPlot(Allcells.data,
group.by = "lineage.bias",
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
cols.use = c("#cc391b","#969696","#026c9a"),
no.legend = F)
# Fate biase in progenitor domaines
data <- melt(data.frame(Clusters = Allcells.data@meta.data$Cluster.ident,
Nfix.score = Allcells.data@meta.data$prob.Nfix,
Nr4a2.score = Allcells.data@meta.data$prob.Nr4a2))
colnames(data) <- c("Clusters", "lineage", "Score")
data %>%
filter(Clusters %in% c("AP.Ventral.Pallium", "AP.lateral.Pallium.1", "AP.lateral.Pallium.2", "AP.Dorsal.Pallium")) %>%
mutate(Clusters = factor(Clusters, levels =c("AP.Ventral.Pallium", "AP.lateral.Pallium.1", "AP.lateral.Pallium.2", "AP.Dorsal.Pallium"))) %>%
ggplot(aes(x=factor(lineage, levels = c("Nr4a2.score", "Nfix.score") ), y=Score, fill= Clusters)) +
geom_boxplot(notch=TRUE) +
geom_point(position=position_jitterdodge(jitter.width = 0.2), size = .5, aes(group=Clusters)) +
scale_fill_manual(values= c("#68b041", "#e3c148", "#b7d174", "#e46b6b")) +
xlab("")
# Subset dataset to keep only cells whith confident differences in fate probabilities
VP.DP.lineage.cells <- c(rownames(subset(Allcells.data@meta.data, Allcells.data@meta.data$Nr4a2.biase == "Nr4a2.lineage")),
rownames(subset(Allcells.data@meta.data, Allcells.data@meta.data$Nfix.biase == "Nfix.lineage")))
# We further only retain VP and DP progenitor clusters
Allcells.data <- SubsetData(Allcells.data,
cells.use = VP.DP.lineage.cells,
ident.remove = c("AP.lateral.Pallium.1", "AP.lateral.Pallium.2") ,
subset.raw = T,
do.clean = F)
DimPlot(Allcells.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
cols.use = c( "#dfdfdf","#68b041", "#e46b6b","#cc3a1b","#cc8778","#e6bb9b" ,"#046c9a","#4784a2"),
no.legend = T)
Allcells.data@meta.data$Lineage <- ifelse(Allcells.data@meta.data$Nr4a2.biase == "Nr4a2.lineage", "Nr4a2.lineage", "Nfix.lineage")
DimPlot(Allcells.data,
group.by = "Lineage",
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
cols.use = c("#cc391b","#026c9a"),
label.size = 4,
no.legend = F)
Manuscript Fig. 7C
We decided to use spring space dimensionality reduction to fit the principale curve because it has been calculated on all cells together. Thus reflecting pan neuronal differentiation axis of variation.
Nr4a2.lineage.cells <- Allcells.data@meta.data %>%
filter(Lineage == "Nr4a2.lineage") %>%
pull(Barcodes)
## Fit a principale curve on the data in the spring space
data.Nr4a2 <- as.data.frame(Allcells.data@dr$spring@cell.embeddings[Nr4a2.lineage.cells,1:2])
fit <- principal_curve(as.matrix(data.Nr4a2[,1:2]),
smoother='smooth_spline',
trace=TRUE,
stretch=2)
## Starting curve---distance^2: 4868119869
## Iteration 1---distance^2: 734019.9
## Iteration 2---distance^2: 770558.7
## Iteration 3---distance^2: 787358.8
## Iteration 4---distance^2: 794176.1
## Iteration 5---distance^2: 797375.2
## Iteration 6---distance^2: 799004.2
## Iteration 7---distance^2: 800140.8
## Iteration 8---distance^2: 800778.4
Nr4a2.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) #The principal curve smoothed
data.Nr4a2$PseudotimeScore <- fit$lambda/max(fit$lambda)
data.Nr4a2$Cluster <- Allcells.data@meta.data %>%
filter(Lineage == "Nr4a2.lineage") %>%
pull(NewClusterID)
# Direction of the maturation score using Nes expression (reverte if positive correlation)
if (cor(data.Nr4a2$PseudotimeScore, Allcells.data@data['Hmga2', Nr4a2.lineage.cells]) > 0) {
data.Nr4a2$PseudotimeScore <- -(data.Nr4a2$PseudotimeScore - max(data.Nr4a2$PseudotimeScore))
}
Nfix.lineage.cells <- Allcells.data@meta.data %>%
filter(Lineage == "Nfix.lineage") %>%
pull(Barcodes)
## Fit a principale curve on the data in the spring space
data.Nfix <- as.data.frame(Allcells.data@dr$spring@cell.embeddings[Nfix.lineage.cells,1:2])
fit <- principal_curve(as.matrix(data.Nfix[,1:2]),
smoother='smooth_spline',
trace=TRUE,
stretch=2)
## Starting curve---distance^2: 741666427
## Iteration 1---distance^2: 263692.7
## Iteration 2---distance^2: 257894.2
## Iteration 3---distance^2: 257508.9
## Iteration 4---distance^2: 257452.1
Nfix.pc.line <- as.data.frame(fit$s[order(fit$lambda),]) #The principal curve smoothed
data.Nfix$PseudotimeScore <- fit$lambda/max(fit$lambda)
data.Nfix$Cluster <- Allcells.data@meta.data %>%
filter(Lineage == "Nfix.lineage") %>%
pull(NewClusterID)
# Direction of the maturation score using Nes expression (reverte if positive correlation)
if (cor(data.Nfix$PseudotimeScore, Allcells.data@data['Hmga2', Nfix.lineage.cells]) > 0) {
data.Nfix$PseudotimeScore <- -(data.Nfix$PseudotimeScore - max(data.Nfix$PseudotimeScore))
}
Pseudotime <- rbind(data.Nr4a2 %>% dplyr::select(PseudotimeScore),
data.Nfix %>% dplyr::select(PseudotimeScore))
Allcells.data@meta.data$Pseudotime <- Pseudotime[rownames(Allcells.data@meta.data),]
#Plot pseudotime score along the two principal curves
data <- as.data.frame(Allcells.data@dr$spring@cell.embeddings[,1:2])
data$PseudotimeScore <- Allcells.data@meta.data$Pseudotime
cols <- colorRampPalette(brewer.pal(n =11, name = "Spectral"))(100)
ggplot(data, aes(spring1, spring2)) +
geom_point(aes(color=PseudotimeScore), size=2, shape=16) +
scale_color_gradientn(colours=rev(cols), name='Speudotime score') +
geom_line(data=Nfix.pc.line, color="#cc391b", size=0.77) +
geom_line(data=Nr4a2.pc.line, color="#026c9a", size=0.77)
# Filter genes
num.cells <- Matrix::rowSums(Allcells.data@raw.data > 0)
genes.use <- names(x = num.cells[num.cells >= 10])
Allcells.data@raw.data <- Allcells.data@raw.data[genes.use, ]
# Normalization and variable gene selection
Allcells.data <- NormalizeData(object = Allcells.data,
normalization.method = "LogNormalize",
scale.factor = round(median(Allcells.data@meta.data$nUMI)),
display.progress = F)
Allcells.data <- FindVariableGenes(object = Allcells.data,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.03,
x.high.cutoff = 4,
y.cutoff = 1.3, do.plot = F,
display.progress = F)
length(Allcells.data@var.genes)
## [1] 926
# Transfert metadata
meta.data <- data.frame(barcode = rownames(Allcells.data@meta.data),
Cluster = Allcells.data@meta.data$Cluster.ident,
Lineage = Allcells.data@meta.data$Lineage,
Pseudotime.Score = Allcells.data@meta.data$Pseudotime,
row.names = rownames(Allcells.data@meta.data))
Annot.data <- new('AnnotatedDataFrame', data = meta.data)
# Transfert count data
count.data = data.frame(gene_short_name = rownames(Allcells.data@raw.data),
row.names = rownames(Allcells.data@raw.data))
feature.data <- new('AnnotatedDataFrame', data = count.data)
# Create the CellDataSet object
gbm_cds <- newCellDataSet(as.matrix(Allcells.data@raw.data),
phenoData = Annot.data,
featureData = feature.data,
lowerDetectionLimit = 1,
expressionFamily = negbinomial())
# Exclude cell cycle associated genes
CCgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
Input.genes <- Allcells.data@var.genes[!Allcells.data@var.genes %in% CCgenes]
Speudotime.lineages.diff <- differentialGeneTest(gbm_cds[Input.genes,],
fullModelFormulaStr = "~sm.ns(Pseudotime.Score, df = 3)*Lineage",
reducedModelFormulaStr = "~sm.ns(Pseudotime.Score, df = 3)",
cores = detectCores() - 2)
# Filter genes based on FDR
Speudotime.lineages.diff.filtered <- Speudotime.lineages.diff %>% filter(qval < 5e-2)
We find direction of the DEG by calculating the area between curves (ABC) for branch-dependent genes by adapting the monocle package function calABCs
. Genes specific ABC is computed on smoothed expression value over 100 points along the pseudotime
# Create a new pseudo-DV vector of 500 points
nPoints <- 100
new_data = list()
for (Lineage in unique(pData(gbm_cds)$Lineage)){
new_data[[length(new_data) + 1]] = data.frame(Pseudotime.Score = seq(min(pData(gbm_cds)$Pseudotime.Score), max(pData(gbm_cds)$Pseudotime.Score), length.out = nPoints), Lineage=Lineage)
}
new_data = do.call(rbind, new_data)
# Smooth gene expression
curve_matrix <- genSmoothCurves(gbm_cds[as.character(Speudotime.lineages.diff.filtered$gene_short_name),],
trend_formula = "~sm.ns(Pseudotime.Score, df = 3)*Lineage",
relative_expr = TRUE,
new_data = new_data,
cores= detectCores() - 2)
# Extract matrix containing smoothed curves for each lineages
Nr4a2_curve_matrix <- curve_matrix[, 1:nPoints] # the first 100 points correspond to Nr4a2 cells
Nfix_curve_matrix <- curve_matrix[, (nPoints + 1):(2 * nPoints)]
ABCs_res <- Nr4a2_curve_matrix - Nfix_curve_matrix # Direction of the comparison : postive ABCs <=> Upregulated in Nr4a2 lineage
ILR_res <- log2(Nr4a2_curve_matrix/ (Nfix_curve_matrix + 0.1)) # Average logFC between the 2 curves
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
Speudotime.lineages.diff.filtered <- cbind(Speudotime.lineages.diff.filtered[,1:4],
ABCs_res,
Speudotime.lineages.diff.filtered[,5:6])
# Extract Nr4a2 neurons trajectory genes
Nr4a2.res <- as.data.frame(Speudotime.lineages.diff.filtered[Speudotime.lineages.diff.filtered$ABCs > 0,])
Nr4a2.genes <- row.names(Nr4a2.res)
Nr4a2_curve_matrix <- Nr4a2_curve_matrix[rownames(Nr4a2_curve_matrix) %in% Nr4a2.genes, ]
# Groupe genes in 6 clusters using partitioning round medoids
Nr4a2.genes.clusters <- cluster::pam(as.dist((1 - cor(t(Nr4a2_curve_matrix), method = "pearson"))), k=6)
# Create a dataframe storing DEG test and clustering results
Nr4a2.Gene.dynamique <- data.frame(Gene= names(Nr4a2.genes.clusters$clustering),
pval=Nr4a2.res$pval,
qval=Nr4a2.res$qval,
ABCs=Nr4a2.res$ABCs,
Gene.Clusters= paste0("Clust.", Nr4a2.genes.clusters$clustering)) %>% arrange(Gene.Clusters)
row.names(Nr4a2.Gene.dynamique) <- Nr4a2.Gene.dynamique$Gene
write.table(Nr4a2.Gene.dynamique, "./Nr4a2.Gene.dynamique.csv", sep=";")
# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(Nr4a2_curve_matrix)), method = "pearson")))
row.ser <- seriate(dst, method ="R2E") #"R2E"
gene.order <- rownames(Nr4a2_curve_matrix[get_order(row.ser),])
anno.colors <- list(lineage = c(Nfix="#cc391b", Nr4a2="#026c9a"),
Gene.Clusters = c(Clust.1 ="#b79f0b" , Clust.2="#e46B6b", Clust.3="#e7823a", Clust.4="#cc8778", Clust.5="#68b041", Clust.6="#5ab793"))
pheatmap::pheatmap(curve_matrix[gene.order,
c(200:101, #Nfix points
1:100)], #Nr4a2 points
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = Nr4a2.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = data.frame(lineage = rep(c("Nr4a2", "Nfix"), each=100)),
annotation_colors = anno.colors,
show_colnames = F,
show_rownames = T,
fontsize_row = 2,
color = rev(brewer.pal(11,"RdBu")),
breaks = seq(-2.5,2.5, length.out = 11),
main = "Nr4a2 cells enriched genes expression along pseudotime")
Manuscript Fig. 7D
# Extract Nfix neurons trajectory genes
Nfix.res <- as.data.frame(Speudotime.lineages.diff.filtered[Speudotime.lineages.diff.filtered$ABCs < 0,])
Nfix.genes <- row.names(Nfix.res)
Nfix_curve_matrix <- Nfix_curve_matrix[rownames(Nfix_curve_matrix) %in% Nfix.genes, ]
# Groupe genes in 6 clusters by partitioning round medoids
Nfix.genes.clusters <- cluster::pam(as.dist((1 - cor(t(Nfix_curve_matrix), method = "pearson"))), k=6)
# Create a dataframe storing DEG test and clustering results
Nfix.Gene.dynamique <- data.frame(Gene= names(Nfix.genes.clusters$clustering),
pval=Nfix.res$pval,
qval=Nfix.res$qval,
ABCs=Nfix.res$ABCs,
Gene.Clusters= paste0("Clust.", Nfix.genes.clusters$clustering)) %>% arrange(Gene.Clusters)
row.names(Nfix.Gene.dynamique) <- Nfix.Gene.dynamique$Gene
write.table(Nfix.Gene.dynamique, "./Nfix.Gene.dynamique.csv", sep=";")
# Order the rows using seriation
dst <- as.dist((1-cor(scale(t(Nfix_curve_matrix)), method = "pearson")))
row.ser <- seriate(dst, method ="R2E")
gene.order <- rownames(Nfix_curve_matrix[get_order(row.ser),])
pheatmap::pheatmap(curve_matrix[gene.order,
c(200:101, #Nfix points
1:100)], #Nr4a2 points
scale = "row",
cluster_rows = F,
cluster_cols = F,
annotation_row = Nfix.Gene.dynamique %>% dplyr::select(Gene.Clusters),
annotation_col = data.frame(lineage = rep(c("Nr4a2", "Nfix"), each=100)),
annotation_colors = anno.colors,
show_colnames = F,
show_rownames = T,
fontsize_row = 2,
color = rev(brewer.pal(11,"RdBu")),
breaks = seq(-2.5,2.5, length.out = 11),
main = "Nfix cells enriched genes expression along pseudotime")
Manuscript Fig. 7C
Plot.Genes.trend(Allcells.data,
genes = c("Nfib", "Dbx1",
"Neurod6", "Pbx3",
"Pcp4", "Zfhx3"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. 7E
Plot.Genes.trend(Allcells.data,
genes = c("Neurod4", "Neurog2",
"Zfhx3", "Tubb3",
"Hey1", "Mapt"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. 7E
Plot.Genes.trend(Allcells.data,
genes = c("Tbr1", "Emx1",
"Neurod1", "Hey1",
"Neurod2", "Flrt3"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. S7B
Plot.Genes.trend(Allcells.data,
genes = c("Nfia", "Dmrt3",
"Cnr1", "Slc30a10"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. S7B
Plot.Genes.trend(Allcells.data,
genes = c("Meis2", "Gm29260",
"Nr2f2", "Dleu7",
"Zbtb20", "Svil"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. S7B
Plot.Genes.trend(Allcells.data,
genes = c("Epha3", "Pdlim4",
"Tox", "Tfap2c"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. S7B
Plot.Genes.trend(Allcells.data,
genes = c("Neurog1", "Eomes",
"Btg2", "Foxg1",
"Emx2"),
color.by = "lineage",
Use.scale.data = F)
Manuscript Fig. S7B
set.seed(100)
# Run Kmeans clustering
cl <- kmeans(cbind(Allcells.data@meta.data$AP_signature1,
Allcells.data@meta.data$BP_signature1,
Allcells.data@meta.data$EN_signature1,
Allcells.data@meta.data$LN_signature1), 4)
Allcells.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
col.pal <- wes_palette("GrandBudapest1", 4, type = "discrete")
p1 <- ggplot(Allcells.data@meta.data, aes(x=AP_signature1, y=BP_signature1, colour = kmeanClust)) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey")
DimPlot(Allcells.data,
group.by = "kmeanClust",
reduction.use = "spring",
cols.use = col.pal,
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
no.legend = F)
# Find the kmeans cluster with the highest mean BP_signature1 score
MeanKclust.BPscore <- aggregate(BP_signature1 ~ kmeanClust, Allcells.data@meta.data, mean)
BPclust <- MeanKclust.BPscore %>% filter(BP_signature1 == max(BP_signature1)) %>% pull(kmeanClust)
# Extract barcodes and filter the seurat object
Glut.cells <- Allcells.data@meta.data %>% filter(kmeanClust == BPclust) %>% pull(Barcodes)
BP.data <- SubsetData(Allcells.data, cells.use = Glut.cells , subset.raw = T, do.clean = F)
DimPlot(BP.data,
group.by = "Lineage",
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
cols.use = c("#cc391b","#026c9a"),
label.size = 4,
no.legend = F)
# Filter genes
num.cells <- Matrix::rowSums(BP.data@raw.data > 0)
genes.use <- names(x = num.cells[num.cells >= 10])
BP.data@raw.data <- BP.data@raw.data[genes.use, ]
# Normalization and variable gene selection
BP.data <- NormalizeData(object = BP.data,
normalization.method = "LogNormalize",
scale.factor = round(median(BP.data@meta.data$nUMI)),
display.progress = F)
BP.data <- FindVariableGenes(object = BP.data,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.03,
x.high.cutoff = 4,
y.cutoff = 1.3, do.plot = F,
display.progress = F)
# Exclude cell cycle associated, mt and ribo ribo
CCgenes <- as.character(read.table("/home/matthieu/Bureau/R/FiguresCodes/Progenitors/CellCycleGenes.csv", sep = "\t", header = F)[,1])
GenesToRemove <- c(grep(pattern = "(^Rpl|^Rps|^Mrp)", x = row.names(BP.data@data), value = TRUE),
grep(pattern = "^mt-", x = row.names(BP.data@data), value = TRUE),
"Xist", CCgenes)
Input.genes <- BP.data@var.genes[!BP.data@var.genes %in% GenesToRemove]
# Set ident
BP.data <- SetAllIdent(BP.data, id = "Lineage")
NfixvsNr4a2.DEgenes <- FindMarkers(object = BP.data,
test.use = "roc",
ident.1 = "Nfix.lineage",
ident.2 = "Nr4a2.lineage",
min.pct = 0,
logfc.threshold = 0,
print.bar = F,
only.pos = F,
genes.use = Input.genes)
NfixvsNr4a2.DEgenes$Gene <- rownames(NfixvsNr4a2.DEgenes)
write.table(NfixvsNr4a2.DEgenes, "./BP_ROC_res.csv", sep=";")
# Set a new data frame which will store all variable needed to generate the plot
Vplot.data <- data.frame(gene = NfixvsNr4a2.DEgenes$Gene,
AUC = NfixvsNr4a2.DEgenes$myAUC,
LogFC = NfixvsNr4a2.DEgenes$avg_logFC)
temporal.module <- row.names(read.table("./Progenitors/Temporal.gene.clusters.csv", sep = ";", header = T, row.names = 1))
spatial.module <- row.names(read.table("./Progenitors/Spatial.gene.clusters.csv", sep = ";", header = T, row.names = 1))
Vplot.data <- dplyr::mutate(Vplot.data, color = dplyr::case_when(Vplot.data$gene %in% temporal.module ~ "Temporal",
Vplot.data$gene %in% spatial.module ~ "Spatial",
TRUE ~ "NA"))
ggplot(Vplot.data, aes(x = LogFC, y = AUC, color = color)) +
geom_point(size = 2, alpha = 0.8, na.rm = T) +
scale_color_manual(name = "AP module",
values = c("Temporal" = "#cc391b",
"Spatial" = "#026c9a",
"NA" = "#969696")) +
xlab(expression(log("Fold change"))) +
ylab("AUC") +
geom_vline(xintercept = c(-0.4, 0.4), colour = "#969696") +
geom_hline(yintercept = c(0.65, 0.35), colour = "#969696") +
ggrepel::geom_text_repel(aes(label=ifelse(gene %in% c(temporal.module,spatial.module) & abs(Vplot.data$LogFC) > 0.4 & (Vplot.data$AUC > 0.65 | Vplot.data$AUC < 0.35), as.character(gene),'')),
box.padding = 0.35,
point.padding = 0.5,
segment.color = 'grey50',
colour = "black")
## [1] "30 novembre, 2020, 11,27"
## 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] wesanderson_0.3.6 RColorBrewer_1.1-2 patchwork_0.0.1
## [4] ggExtra_0.9 reshape2_1.4.3 dplyr_0.8.3
## [7] seriation_1.2-9 velocyto.R_0.6 FateID_0.1.9
## [10] monocle_2.14.0 DDRTree_0.1.5 irlba_2.3.3
## [13] VGAM_1.1-2 Biobase_2.46.0 BiocGenerics_0.32.0
## [16] princurve_2.1.4 Seurat_2.3.4 Matrix_1.2-17
## [19] cowplot_1.0.0 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 lle_1.1 fastICA_1.2-2
## [10] digest_0.6.25 foreach_1.4.7 htmltools_0.5.0
## [13] viridis_0.5.1 lars_1.2 gdata_2.18.0
## [16] magrittr_1.5 checkmate_1.9.4 cluster_2.1.0
## [19] mixtools_1.1.0 ROCR_1.0-7 limma_3.42.0
## [22] matrixStats_0.55.0 R.utils_2.9.0 docopt_0.6.1
## [25] askpass_1.1 colorspace_1.4-1 ggrepel_0.8.1
## [28] xfun_0.18 sparsesvd_0.2 crayon_1.3.4
## [31] jsonlite_1.7.0 zeallot_0.1.0 survival_2.44-1.1
## [34] zoo_1.8-6 iterators_1.0.12 ape_5.3
## [37] glue_1.4.1 registry_0.5-1 gtable_0.3.0
## [40] kernlab_0.9-29 prabclus_2.3-1 DEoptimR_1.0-8
## [43] scales_1.1.0 pheatmap_1.0.12 som_0.3-5.1
## [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 umap_0.2.3.1 htmlwidgets_1.5.1
## [64] httr_1.4.1 FNN_1.1.3 gplots_3.0.1.1
## [67] fpc_2.2-3 acepack_1.4.1 modeltools_0.2-22
## [70] ica_1.0-2 farver_2.0.1 pkgconfig_2.0.3
## [73] R.methodsS3_1.7.1 flexmix_2.3-15 nnet_7.3-14
## [76] locfit_1.5-9.1 labeling_0.3 later_1.0.0
## [79] tidyselect_0.2.5 rlang_0.4.7 munsell_0.5.0
## [82] tools_3.6.3 ggridges_0.5.1 fastmap_1.0.1
## [85] evaluate_0.14 stringr_1.4.0 yaml_2.2.1
## [88] npsurv_0.4-0 knitr_1.26 bit64_4.0.2
## [91] fitdistrplus_1.0-14 robustbase_0.93-5 caTools_1.17.1.2
## [94] randomForest_4.6-14 purrr_0.3.3 RANN_2.6.1
## [97] pbapply_1.4-2 nlme_3.1-141 mime_0.7
## [100] slam_0.1-46 R.oo_1.23.0 hdf5r_1.3.2.9000
## [103] compiler_3.6.3 rstudioapi_0.11 png_0.1-7
## [106] lsei_1.2-0 tibble_2.1.3 stringi_1.4.6
## [109] highr_0.8 lattice_0.20-41 HSMMSingleCell_1.6.0
## [112] vctrs_0.2.0 pillar_1.4.2 lifecycle_0.1.0
## [115] combinat_0.0-8 Rdpack_0.11-0 lmtest_0.9-37
## [118] data.table_1.12.6 bitops_1.0-6 gbRd_0.4-11
## [121] httpuv_1.5.2 pcaMethods_1.78.0 R6_2.4.1
## [124] latticeExtra_0.6-28 promises_1.1.0 TSP_1.1-10
## [127] KernSmooth_2.23-15 gridExtra_2.3 codetools_0.2-16
## [130] MASS_7.3-53 gtools_3.8.1 assertthat_0.2.1
## [133] openssl_1.4.1 withr_2.1.2 qlcMatrix_0.9.7
## [136] mgcv_1.8-33 diptest_0.75-7 doSNOW_1.0.18
## [139] grid_3.6.3 rpart_4.1-15 tidyr_1.0.0
## [142] class_7.3-17 rmarkdown_2.5 segmented_1.0-0
## [145] Rtsne_0.15 shiny_1.4.0 snowfall_1.84-6.1
## [148] scatterplot3d_0.3-41 base64enc_0.1-3
Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩