# 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]))
}
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)
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()
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)
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")
# 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")
Plot.Genes.trend(Allcells.data,
genes = c("Nfib", "Dbx1",
"Neurod6", "Pbx3",
"Pcp4", "Zfhx3"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Neurod4", "Neurog2",
"Zfhx3", "Tubb3",
"Hey1", "Mapt"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Tbr1", "Emx1",
"Neurod1", "Hey1",
"Neurod2", "Flrt3"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Nfia", "Dmrt3",
"Cnr1", "Slc30a10"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Meis2", "Gm29260",
"Nr2f2", "Dleu7",
"Zbtb20", "Svil"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Epha3", "Pdlim4",
"Tox", "Tfap2c"),
color.by = "lineage",
Use.scale.data = F)
Plot.Genes.trend(Allcells.data,
genes = c("Neurog1", "Eomes",
"Btg2", "Foxg1",
"Emx2"),
color.by = "lineage",
Use.scale.data = F)
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↩