Spatial and maturation signatures along corticogenesis

Analysis of the pallial progenitor domains along corticogenesis from the Mouse Brain Atlas - Development released by the Linnarsson’s lab on July 2020 (BioRxiv preprint).

load libraries

library(loomR)
library(Seurat)
library(scrattch.hicat)
library(reticulate)

library(dplyr)
library(ggExtra)

library(ggplot2)
library(patchwork)
library(RColorBrewer)

#Set ggplot theme classic
theme_set(theme_classic())

#Set python env for UMAP 
use_python("/home/matthieu/.local/share/r-miniconda/envs/r-reticulate/bin/python")

Extraction of annotated radial glial cell from E11 to E15.5

Reading loom file containing cluster metadata

lfile <- connect(filename = "./DevelopingMouseBrainAtlas/dev_all.agg.loom",
                 skip.validate =T,
                 mode="r+")

We extract all cluster annotated as “radial glia”

ClustersAnno <- lfile$col.attrs$Clusters[][lfile$col.attrs$Class[] %in% c("Radial glia")]

lfile$close_all()

Reading loom file containing all metadata and count table for the 292,495 high-quality cells

lfile <- connect(filename = "./DevelopingMouseBrainAtlas/dev_all.loom",
                 skip.validate =T,
                 mode="r+")

Plot extracted cells on the global tsne

tSNE <- t(lfile$col.attrs$BTSNE[,])

tSNE.df <- lfile$get.attribute.df(MARGIN = 2, attributes = c("Clusters", "Age"), col.names = "CellID")
tSNE.df$GlialCell <- ifelse(tSNE.df$Clusters %in% ClustersAnno & tSNE.df$Age %in% c("e11.0","e12.0","e12.5","e13.0","e13.5","e14.0", "e14.5", "e15.0", "e15.5"), "Radial glia", "OtherCells")

tSNE.df$tSNE.1 <- tSNE[,1]
tSNE.df$tSNE.2 <- tSNE[,2]

plot <-ggplot(tSNE.df, aes(x=tSNE.1, y=tSNE.2, colour = GlialCell)) +
          geom_point( size=0.3, shape=16) +
          scale_color_manual(values= c("grey90","#e7823a")) +
          theme_classic()+
          theme(axis.text.x=element_blank(), axis.text.y=element_blank(),
                axis.ticks.x=element_blank(), axis.ticks.y=element_blank(),
                axis.title.x=element_blank(),axis.title.y=element_blank()) +
          ggtitle(label = paste(names(table(tSNE.df$GlialCell)[2]), table(tSNE.df$GlialCell)[2]))

plot$data <- plot$data[order(plot$data$GlialCell,decreasing = F),]

plot

rm(tSNE.df,tSNE)

Extract Radial glia cells count matrix and metadata

Retrieve the count matrix

cells <- lfile$col.attrs$Clusters[] %in% ClustersAnno & lfile$col.attrs$Age[] %in% c("e11.0","e12.0","e12.5","e13.0","e13.5","e14.0", "e14.5", "e15.0", "e15.5")
data.subset <- lfile[["matrix"]][ cells, ]
data.subset <- t(data.subset)

dim(data.subset)

## [1] 31053 22650

Retrieve gene names

gene.names <- lfile[["row_attrs/Gene"]][]
rownames(data.subset) <- gene.names

Extract meta.data information

attrs <- c("CellID","Age", "Batch", "Chemistry", "ChipID","Clusters", "Region", "Tissue", "PseudoTissue")
attr.df <- lfile$get.attribute.df(MARGIN = 2, attributes = attrs, col.names = "CellID")
attr.df <- attr.df[cells,]
rownames(attr.df) <- attr.df$CellID

colnames(data.subset) <- attr.df$CellID

Closing loom object

lfile$close_all()

Perfom analysis with seurat v2.3.4

Initiate a Seurat object

Raw.data <- CreateSeuratObject(raw.data = data.subset,
                               min.cells = 5,
                               min.genes = 0,
                               project = "RadialGlia_DevelopingMouseBrainAtlas")

rm(list=ls()[!ls() %in% c("Raw.data","attr.df")])
gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   3263133  174.3    6057934   323.6    6057934   323.6
## Vcells 879352160 6709.0 2453505816 18718.8 2991256468 22821.5

We exclude genes with duplicated row names

rownames(Raw.data@raw.data)[duplicated(rownames(Raw.data@raw.data))]

## [1] "Atp5o"   "Gm2464"  "Gcat"    "Nudt8"   "Aldoa"   "Ints5"   "Zc3h11a"
## [8] "Pick1"

Raw.data@raw.data <- Raw.data@raw.data[!duplicated(rownames(Raw.data@raw.data)),]

Low quality cells filtering and data normalization

Compute ribosomal and mitochondrial transcript fractions per cell

mito.genes <- grep(pattern = "^mt-", x = rownames(x = Raw.data@data), value = TRUE)
percent.mito <- Matrix::colSums(Raw.data@raw.data[mito.genes, ])/Matrix::colSums(Raw.data@raw.data)
Raw.data <- AddMetaData(object = Raw.data, metadata = percent.mito, col.name = "percent.mito")

ribo.genes <- grep(pattern = "(^Rpl|^Rps|^Mrp)", x = rownames(x = Raw.data@data), value = TRUE)
percent.ribo <- Matrix::colSums(Raw.data@raw.data[ribo.genes, ])/Matrix::colSums(Raw.data@raw.data)
Raw.data <- AddMetaData(object = Raw.data, metadata = percent.ribo, col.name = "percent.ribo")

Raw.data@meta.data <- cbind(Raw.data@meta.data, attr.df)
rm(mito.genes, percent.mito,ribo.genes,percent.ribo)

Cell.QC.Stat <- Raw.data@meta.data

Filtering cells based on number of genes and transcripts detected

Remove cells with to few genes detected or with to many UMI counts

# Set low and hight thresholds on the number of detected genes
min.Genes.thr <- median(log10(Cell.QC.Stat$nGene)) - 3*mad(log10(Cell.QC.Stat$nGene))
max.Genes.thr <- median(log10(Cell.QC.Stat$nGene)) + 3*mad(log10(Cell.QC.Stat$nGene))

# Set hight threshold on the number of transcripts
max.nUMI.thr <- median(log10(Cell.QC.Stat$nUMI)) + 3*mad(log10(Cell.QC.Stat$nUMI))

# Gene/UMI scatter plot before filtering
p <- ggplot(Cell.QC.Stat, aes(x=log10(nUMI), y=log10(nGene))) +
      geom_point() +
      geom_smooth(method="lm") +
      geom_hline(aes(yintercept = min.Genes.thr), colour = "green", linetype = 2) +
      geom_hline(aes(yintercept = max.Genes.thr), colour = "green", linetype = 2) +
      geom_vline(aes(xintercept = max.nUMI.thr), colour = "red", linetype = 2)

ggMarginal(p, type = "histogram", fill="lightgrey")

# Filter cells base on both metrics
Cell.QC.Stat <- Cell.QC.Stat %>% filter(log10(nGene) > min.Genes.thr) %>% filter(log10(nUMI) < max.nUMI.thr)

p <- ggplot(Cell.QC.Stat, aes(x=log10(nUMI), y=log10(nGene))) +
  geom_point() +
  geom_smooth(method="lm") +
  geom_hline(aes(yintercept = min.Genes.thr), colour = "green", linetype = 2) +
  geom_hline(aes(yintercept = max.Genes.thr), colour = "green", linetype = 2) +
  geom_vline(aes(xintercept = max.nUMI.thr), colour = "red", linetype = 2) +
  annotate(geom = "text", label = paste0(dim(Cell.QC.Stat)[1], " QC passed cells"), x = 4, y = 3.8)

ggMarginal(p, type = "histogram", fill="lightgrey")

Keep only valid cells

Radial.Glia.data <- SubsetData(Raw.data, cells.use = Cell.QC.Stat$CellID , subset.raw = T,  do.clean = F)
rm(list=ls()[!ls() %in% c("Radial.Glia.data")])

# Plot final QC metrics
VlnPlot(object = Radial.Glia.data,
        features.plot = c("nGene","nUMI", "percent.mito", "percent.ribo"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        nCol = 2, point.size.use = 0, group.by = "Age")

Normalization and scaling

Gene filtering

Radial.Glia.data@raw.data <- Radial.Glia.data@raw.data[!duplicated(rownames(Radial.Glia.data@raw.data)),]

num.cells <- Matrix::rowSums(Radial.Glia.data@raw.data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])
Radial.Glia.data@raw.data <- Radial.Glia.data@raw.data[genes.use, ]
Radial.Glia.data@data <- Radial.Glia.data@data[genes.use, ]

Normalization

Radial.Glia.data <- NormalizeData(object = Radial.Glia.data,
                                  normalization.method = "LogNormalize", 
                                  scale.factor = round(median(Radial.Glia.data@meta.data$nUMI)),
                                  display.progress = F)

Radial.Glia.data <- FindVariableGenes(object = Radial.Glia.data,
                                       mean.function = ExpMean,
                                       dispersion.function = LogVMR,
                                       x.low.cutoff = 0.0125,
                                       x.high.cutoff = 3,
                                       y.cutoff = 1,
                                       do.plot = F, display.progress = T)

Assign Cell-Cycle Score

s.genes <- c("Mcm5", "Pcna", "Tym5", "Fen1", "Mcm2", "Mcm4", "Rrm1", "Ung",
             "Gins2", "Mcm6", "Cdca7", "Dtl", "Prim1", "Uhrf1", "Mlf1ip",
             "Hells", "Rfc2", "Rap2", "Nasp", "Rad51ap1", "Gmnn", "Wdr76", "Slbp",
             "Ccne2", "Ubr7", "Pold3", "Msh2", "Atad2", "Rad51", "Rrm2", "Cdc45",
             "Cdc6", "Exo1", "Tipin", "Dscc1", "Blm", " Casp8ap2", "Usp1", "Clspn",
             "Pola1", "Chaf1b", "Brip1", "E2f8")

g2m.genes <- c("Hmgb2", "Ddk1","Nusap1", "Ube2c", "Birc5", "Tpx2", "Top2a", "Ndc80",
               "Cks2", "Nuf2", "Cks1b", "Mki67", "Tmpo", " Cenpk", "Tacc3", "Fam64a",
               "Smc4", "Ccnb2", "Ckap2l", "Ckap2", "Aurkb", "Bub1", "Kif11", "Anp32e",
               "Tubb4b", "Gtse1", "kif20b", "Hjurp", "Cdca3", "Hn1", "Cdc20", "Ttk", "Cdc25c",
               "kif2c", "Rangap1", "Ncapd2", "Dlgap5", "Cdca2", "Cdca8", "Ect2", "Kif23", "Hmmr",
               "Aurka", "Psrc1", "Anln", "Lbr", "Ckap5", "Cenpe", "Ctcf", "Nek2", "G2e3", "Gas2l3", "Cbx5", "Cenpa")

Radial.Glia.data<- CellCycleScoring(object = Radial.Glia.data,
                                    s.genes = s.genes,
                                    g2m.genes = g2m.genes,
                                    set.ident = F)

Radial.Glia.data@meta.data$CC.Difference <- Radial.Glia.data@meta.data$S.Score - Radial.Glia.data@meta.data$G2M.Score

VlnPlot(object = Radial.Glia.data,
        features.plot = c("S.Score", "G2M.Score","CC.Difference" ,"percent.ribo"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        nCol = 2, point.size.use = 0, group.by = "Age")

Data scaling and cell cycle regression

Radial.Glia.data <- ScaleData(object = Radial.Glia.data,
                     vars.to.regress = c("CC.Difference","percent.mito","nGene", "nUMI"),
                     display.progress = T)

## 
## Time Elapsed:  1.55461797714233 mins

Dimensionality reduction

PCA and UMAP

Radial.Glia.data <- RunPCA(object = Radial.Glia.data,
                          pcs.compute = 30,
                          do.print = TRUE,
                          pcs.print = F)

rm(list=ls()[!ls() %in% c("Radial.Glia.data")])
gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   3109520  166.1    6057934   323.6    6057934   323.6
## Vcells 888117436 6775.8 2683674431 20474.9 2991256468 22821.5

Radial.Glia.data <- RunUMAP(Radial.Glia.data, dims.use = 1:30, n_neighbors = 8, max.dim = 2)

Features plot

Age, Region, Tissue, Pseudo-tissue

Age

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        group.by = c("Age"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.axes = T,
        no.legend = F)

Region

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        group.by = c("Region"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.axes = T,
        no.legend = F)

Tissue

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        group.by = c("Tissue"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.axes = T,
        no.legend = F)

PseudoTissue

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        group.by = c("PseudoTissue"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.axes = T,
        no.legend = F)

Major radial glia population makers

plot <- FeaturePlot(object = Radial.Glia.data,
                    features.plot =  c("Foxg1","Gsx2","Otx2", "Olig2", "Dlx2",
                                      "Tfap2c", "Dmrta2", "Fezf2", "Pax6","Emx1",
                                      "Shh", "Rspo2","Sfrp2", "Fgf8", "Qrfpr"),
                    cols.use = c("grey90", brewer.pal(9,"YlGnBu")),
                    reduction.use = "umap",
                    no.legend = T,
                    overlay = F,
                    dark.theme = F,
                    do.return =T,
                    no.axes = T)

for (i in 1:length(plot)){
  plot[[i]]$data <- plot[[i]]$data[order(plot[[i]]$data$gene),]
}

cowplot::plot_grid(plotlist = plot[1:15], ncol = 5)

Extract Forbrain cells

Forbrain.cells <- Radial.Glia.data@meta.data %>% filter(Region == "Forebrain") %>% pull(CellID)

Radial.Glia.data <- SubsetData(Radial.Glia.data, cells.use = Forbrain.cells , subset.raw = T,  do.clean = F)
rm(list=ls()[!ls() %in% c("Radial.Glia.data")])

Normalization and scaling

Gene filtering

num.cells <- Matrix::rowSums(Radial.Glia.data@raw.data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])
Radial.Glia.data@raw.data <- Radial.Glia.data@raw.data[genes.use, ]
Radial.Glia.data@data <- Radial.Glia.data@data[genes.use, ]

Normalization

Radial.Glia.data <- NormalizeData(object = Radial.Glia.data,
                         normalization.method = "LogNormalize", 
                         scale.factor = round(median(Radial.Glia.data@meta.data$nUMI)),
                         display.progress = F)

Radial.Glia.data <- FindVariableGenes(object = Radial.Glia.data,
                                       mean.function = ExpMean,
                                       dispersion.function = LogVMR,
                                       x.low.cutoff = 0.0125,
                                       x.high.cutoff = 3,
                                       y.cutoff = 1,
                                       do.plot = F, display.progress = T)

Data scaling and cell cycle regression

Radial.Glia.data <- ScaleData(object = Radial.Glia.data,
                     vars.to.regress = c("CC.Difference","percent.mito","nGene", "nUMI"),
                     display.progress = T)

## 
## Time Elapsed:  1.63568867842356 mins

Dimensionality reduction

Radial.Glia.data <- RunPCA(object = Radial.Glia.data,
                          pcs.compute = 30,
                          do.print = TRUE,
                          pcs.print = F)

Radial.Glia.data <- RunUMAP(Radial.Glia.data, dims.use = 1:30, n_neighbors = 8, max.dim = 2)

colors <- c("#b7d174", "#83c3b8", "#3e69ac",
            "#e46b6b",  "#009fda", "#c773a7",
            "#ec756d", "#7293c8", "#b79f0b",
            "#ebcb2e", "#cc3a1b", "#cc8778")

DimPlot(Radial.Glia.data,
        group.by = "Age",
        reduction.use = "umap",
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.legend = F,
        no.axes = T,
        cols.use = rep(c("#bdd8ef", "#6db8e2", "#357ebc", "#1c4896"),each=2))

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        group.by = c("Tissue"),
        cols.use = c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d","#4cabdc", "#5ab793", "#e7823a"),
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.axes = T,
        no.legend = F)

rm(list=ls()[!ls() %in% c("Radial.Glia.data")])
gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   3159501  168.8    6057934   323.6    6057934   323.6
## Vcells 587530028 4482.5 2146939545 16379.9 2991256468 22821.5

Perform iterative clustering

Prepare the dataset for clustering with scrattch.hicat

Gene filtering

We decided to exclude cell cycle, ribosomal and mitochondrial associated genes, as well as Xist for the clustering step.

# Exclude cell cycle associated genes
ccgenes <- as.character(read.table("./Progenitors/CellCycleGenes.csv", sep = ";")[,1])

# Exclude genes detected in less than 10 cells
num.cells <- Matrix::rowSums(Radial.Glia.data@data > 0)
genes.use <- names(x = num.cells[which(x = num.cells >= 10)])

GenesToRemove <- c(grep(pattern = "(^Rpl|^Rps|^Mrp)", x = genes.use, value = TRUE),
                   grep(pattern = "^mt-", x = genes.use, value = TRUE),
                   "Xist", ccgenes)

genes.use <- genes.use[!genes.use %in% GenesToRemove]

Normalization

dgeMatrix_count <- as.matrix(Radial.Glia.data@raw.data)[rownames(Radial.Glia.data@raw.data) %in% genes.use,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1)

norm.dat <- Matrix(norm.dat, sparse = TRUE)
Data.matrix <- list(raw.dat= dgeMatrix_count, norm.dat= norm.dat)
attach(Data.matrix)

rm(list = ls()[!ls() %in% c("Radial.Glia.data", "Data.matrix")])
gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   3166012  169.1    6057934   323.6    6057934   323.6
## Vcells 899707492 6864.3 2146939545 16379.9 2991256468 22821.5

Filter unwanted sources of variation

gene.counts <- log2(colSums(as.matrix(Data.matrix$norm.dat) > 0))
nUMI <- log2(colSums(Data.matrix$raw.dat))
perctMito <- Radial.Glia.data@meta.data$percent.mito
perctRibo <- Radial.Glia.data@meta.data$percent.ribo
ccphasediff <- Radial.Glia.data@meta.data$CC.Difference

rm.eigen <- as.matrix(cbind(gene.counts,
                            nUMI,
                            perctMito,
                            perctRibo,
                            ccphasediff))

row.names(rm.eigen) <- names(gene.counts)

colnames(rm.eigen) <- c("log2nGenes",
                        "log2nUMI",
                        "perctMito",
                        "perctRibo",
                        "ccphasediff")

Set Iterative clustering parameters

# Parameters for iterative clustering
de.param <- de_param(padj.th     = 0.01, 
                     lfc.th      = 1,
                     low.th      = 1, 
                     q1.th       = 0.25, 
                     q2.th       = NULL,
                     q.diff.th   = 0.7,
                     de.score.th = 300,
                     min.cells = 50)

Run iterative clustering

iter.result <- iter_clust(norm.dat, 
                          counts = raw.dat,
                          dim.method = "pca",
                          max.dim = 30,
                          de.param = de.param,
                          type = "undirectional",
                          rm.eigen = rm.eigen,
                          k.param = 15,
                          rm.th = 0.7,
                          vg.padj.th = 0.5,
                          method = "louvain",
                          prefix = "test-iter_clust",
                          verbose = F)

## [1] "test-iter_clust"
##   Finding nearest neighbors...DONE ~ 2.976 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.774 s
##   Build undirected graph from the weighted links...DONE ~ 0.765 s
##   Run louvain clustering on the graph ...DONE ~ 0.587 s
##   Return a community class
##   -Modularity value: 0.8365338 
##   -Number of clusters: 22[1] "test-iter_clust.1"
##   Finding nearest neighbors...DONE ~ 0.006 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.022 s
##   Build undirected graph from the weighted links...DONE ~ 0.026 s
##   Run louvain clustering on the graph ...DONE ~ 0.01 s
##   Return a community class
##   -Modularity value: 0.6294141 
##   -Number of clusters: 7[1] "test-iter_clust.2"
##   Finding nearest neighbors...DONE ~ 0.009 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.035 s
##   Build undirected graph from the weighted links...DONE ~ 0.033 s
##   Run louvain clustering on the graph ...DONE ~ 0.02 s
##   Return a community class
##   -Modularity value: 0.7529057 
##   -Number of clusters: 11[1] "test-iter_clust.3"
##   Finding nearest neighbors...DONE ~ 0.017 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.042 s
##   Build undirected graph from the weighted links...DONE ~ 0.04 s
##   Run louvain clustering on the graph ...DONE ~ 0.021 s
##   Return a community class
##   -Modularity value: 0.7050087 
##   -Number of clusters: 13[1] "test-iter_clust.4"
##   Finding nearest neighbors...DONE ~ 0.384 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.232 s
##   Build undirected graph from the weighted links...DONE ~ 0.229 s
##   Run louvain clustering on the graph ...DONE ~ 0.156 s
##   Return a community class
##   -Modularity value: 0.8127554 
##   -Number of clusters: 19[1] "test-iter_clust.4.1"
##   Finding nearest neighbors...DONE ~ 0.005 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.021 s
##   Build undirected graph from the weighted links...DONE ~ 0.019 s
##   Run louvain clustering on the graph ...DONE ~ 0.01 s
##   Return a community class
##   -Modularity value: 0.671172 
##   -Number of clusters: 7[1] "test-iter_clust.4.2"
##   Finding nearest neighbors...DONE ~ 0.257 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.195 s
##   Build undirected graph from the weighted links...DONE ~ 0.18 s
##   Run louvain clustering on the graph ...DONE ~ 0.124 s
##   Return a community class
##   -Modularity value: 0.8213536 
##   -Number of clusters: 17[1] "test-iter_clust.4.2.1"
##   Finding nearest neighbors...DONE ~ 0.242 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.185 s
##   Build undirected graph from the weighted links...DONE ~ 0.173 s
##   Run louvain clustering on the graph ...DONE ~ 0.113 s
##   Return a community class
##   -Modularity value: 0.8166865 
##   -Number of clusters: 18[1] "test-iter_clust.4.2.2"
##   Finding nearest neighbors...DONE ~ 0.003 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.017 s
##   Build undirected graph from the weighted links...DONE ~ 0.028 s
##   Run louvain clustering on the graph ...DONE ~ 0.01 s
##   Return a community class
##   -Modularity value: 0.6669176 
##   -Number of clusters: 6[1] "test-iter_clust.4.3"
##   Finding nearest neighbors...DONE ~ 0.005 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.019 s
##   Build undirected graph from the weighted links...DONE ~ 0.018 s
##   Run louvain clustering on the graph ...DONE ~ 0.008 s
##   Return a community class
##   -Modularity value: 0.6211829 
##   -Number of clusters: 6[1] "test-iter_clust.5"
##   Finding nearest neighbors...DONE ~ 0.503 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.261 s
##   Build undirected graph from the weighted links...DONE ~ 0.249 s
##   Run louvain clustering on the graph ...DONE ~ 0.196 s
##   Return a community class
##   -Modularity value: 0.8039505 
##   -Number of clusters: 16[1] "test-iter_clust.5.1"
##   Finding nearest neighbors...DONE ~ 0.297 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.164 s
##   Build undirected graph from the weighted links...DONE ~ 0.163 s
##   Run louvain clustering on the graph ...DONE ~ 0.098 s
##   Return a community class
##   -Modularity value: 0.7740656 
##   -Number of clusters: 17[1] "test-iter_clust.5.1.1"
##   Finding nearest neighbors...DONE ~ 0.318 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.157 s
##   Build undirected graph from the weighted links...DONE ~ 0.151 s
##   Run louvain clustering on the graph ...DONE ~ 0.099 s
##   Return a community class
##   -Modularity value: 0.7643717 
##   -Number of clusters: 13[1] "test-iter_clust.5.1.2"
##   Finding nearest neighbors...DONE ~ 0.001 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.01 s
##   Build undirected graph from the weighted links...DONE ~ 0.01 s
##   Run louvain clustering on the graph ...DONE ~ 0.006 s
##   Return a community class
##   -Modularity value: 0.643927 
##   -Number of clusters: 6[1] "test-iter_clust.5.2"
##   Finding nearest neighbors...DONE ~ 0.13 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.076 s
##   Build undirected graph from the weighted links...DONE ~ 0.074 s
##   Run louvain clustering on the graph ...DONE ~ 0.047 s
##   Return a community class
##   -Modularity value: 0.6347574 
##   -Number of clusters: 9[1] "test-iter_clust.5.3"
##   Finding nearest neighbors...DONE ~ 0.017 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.039 s
##   Build undirected graph from the weighted links...DONE ~ 0.037 s
##   Run louvain clustering on the graph ...DONE ~ 0.017 s
##   Return a community class
##   -Modularity value: 0.7105503 
##   -Number of clusters: 12[1] "test-iter_clust.6"
##   Finding nearest neighbors...DONE ~ 0.028 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.054 s
##   Build undirected graph from the weighted links...DONE ~ 0.051 s
##   Run louvain clustering on the graph ...DONE ~ 0.027 s
##   Return a community class
##   -Modularity value: 0.7669241 
##   -Number of clusters: 12[1] "test-iter_clust.6.1"
##   Finding nearest neighbors...DONE ~ 0.004 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.025 s
##   Build undirected graph from the weighted links...DONE ~ 0.023 s
##   Run louvain clustering on the graph ...DONE ~ 0.01 s
##   Return a community class
##   -Modularity value: 0.7934622 
##   -Number of clusters: 9[1] "test-iter_clust.6.2"
##   Finding nearest neighbors...DONE ~ 0.015 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.035 s
##   Build undirected graph from the weighted links...DONE ~ 0.034 s
##   Run louvain clustering on the graph ...DONE ~ 0.015 s
##   Return a community class
##   -Modularity value: 0.7163692 
##   -Number of clusters: 11[1] "test-iter_clust.7"
##   Finding nearest neighbors...DONE ~ 0.016 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.036 s
##   Build undirected graph from the weighted links...DONE ~ 0.033 s
##   Run louvain clustering on the graph ...DONE ~ 0.017 s
##   Return a community class
##   -Modularity value: 0.696045 
##   -Number of clusters: 10[1] "test-iter_clust.8"
##   Finding nearest neighbors...DONE ~ 0.03 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.045 s
##   Build undirected graph from the weighted links...DONE ~ 0.043 s
##   Run louvain clustering on the graph ...DONE ~ 0.025 s
##   Return a community class
##   -Modularity value: 0.7293799 
##   -Number of clusters: 12[1] "test-iter_clust.9"
##   Finding nearest neighbors...DONE ~ 0.092 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.081 s
##   Build undirected graph from the weighted links...DONE ~ 0.075 s
##   Run louvain clustering on the graph ...DONE ~ 0.044 s
##   Return a community class
##   -Modularity value: 0.7855973 
##   -Number of clusters: 13[1] "test-iter_clust.9.1"
##   Finding nearest neighbors...DONE ~ 0.001 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.011 s
##   Build undirected graph from the weighted links...DONE ~ 0.011 s
##   Run louvain clustering on the graph ...DONE ~ 0.004 s
##   Return a community class
##   -Modularity value: 0.5786325 
##   -Number of clusters: 8[1] "test-iter_clust.9.2"
##   Finding nearest neighbors...DONE ~ 0.008 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.026 s
##   Build undirected graph from the weighted links...DONE ~ 0.023 s
##   Run louvain clustering on the graph ...DONE ~ 0.01 s
##   Return a community class
##   -Modularity value: 0.7372565 
##   -Number of clusters: 11[1] "test-iter_clust.9.3"
##   Finding nearest neighbors...DONE ~ 0.037 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.054 s
##   Build undirected graph from the weighted links...DONE ~ 0.052 s
##   Run louvain clustering on the graph ...DONE ~ 0.026 s
##   Return a community class
##   -Modularity value: 0.7374588 
##   -Number of clusters: 10[1] "test-iter_clust.10"
##   Finding nearest neighbors...DONE ~ 0.085 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.077 s
##   Build undirected graph from the weighted links...DONE ~ 0.072 s
##   Run louvain clustering on the graph ...DONE ~ 0.046 s
##   Return a community class
##   -Modularity value: 0.7068849 
##   -Number of clusters: 11

Merge clusters not discriminated by DEGs

# Merge clusters which are not separable by DEGs
rd.dat <- t(norm.dat[iter.result$markers,])
merge.result <- merge_cl(norm.dat, 
                         cl = iter.result$cl, 
                         rd.dat = rd.dat,
                         de.param = de.param)

cat(length(unique(merge.result$cl))," Clusters\n")

## 18  Clusters

Plot clusters on UMAP embedding

Radial.Glia.data@ident <- as.factor(merge.result$cl)
Radial.Glia.data@meta.data$Scrattchicat.ident <- as.character(Radial.Glia.data@ident)

colors <- c("#68b041", "#e3c148", "#b7d174", "#83c3b8",  "#3e69ac", "#e46b6b","#009fda",
            "#c773a7","#ec756d", "#7293c8", "#b79f0b", "#3ca73f","#31b6bd",
            "#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d",
            "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")

DimPlot(Radial.Glia.data,
        reduction.use = "umap",
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 4,
        no.legend = F,
        no.axes = T,
        cols.use = colors)

plot <- FeaturePlot(object = Radial.Glia.data,
            features.plot =  c("Barhl2", "Otx2", "Dlk1", "Sfrp2",
                               "Foxg1", "Emx1", "Gsx2", "Olig2",
                               "Wnt3a", "Rspo3", "Fgf8", "Shh"
                               ),
            cols.use = c("grey90", brewer.pal(9,"YlGnBu")),
            reduction.use = "umap",
            no.legend = T,
            overlay = F,
            dark.theme = F,
            do.return =T,
            no.axes = T)

for (i in 1:length(plot)){
  plot[[i]]$data <- plot[[i]]$data[order(plot[[i]]$data$gene),]
}

cowplot::plot_grid(plotlist = plot[1:12], ncol = 4)

rm(list = ls()[!ls() %in% c("Radial.Glia.data")])
gc()

##             used   (Mb) gc trigger    (Mb)   max used    (Mb)
## Ncells   3247236  173.5    6057934   323.6    6057934   323.6
## Vcells 900987039 6874.0 2146939545 16379.9 2991256468 22821.5

Compute maturation and spatial scores

We use the gene signatures from Fig6 E and G to compute the maturation or spatial scores

Dorso-Ventral scores

Sgenes <- read.table("./Progenitors/Spatial.gene.clusters.csv", sep = ";", header = T, row.names = 1)

# VP enriched genes at E12
genes.list <- Sgenes %>%
              filter(Gene.Clusters %in% c("Clust.1", "Clust.2")) %>%
              pull(Gene) %>%
              as.character()

genes.list <- list(genes.list[genes.list %in% rownames(Radial.Glia.data@data)])

genes.list

## [[1]]
##   [1] "Slc39a10"      "Dner"          "Plekha6"       "Pdk1"         
##   [5] "Tmx4"          "Zcchc12"       "Arx"           "Xist"         
##   [9] "Riiad1"        "Hs3st1"        "Mthfd2l"       "Cald1"        
##  [13] "Nap1l5"        "Peg3"          "St5"           "Bnip3"        
##  [17] "Ostm1"         "Ascl1"         "Gm2694"        "Dleu7"        
##  [21] "Anxa2"         "Spsb4"         "Cyfip2"        "Gpx3"         
##  [25] "Sez6"          "Ptch1"         "Mctp1"         "Atad2b"       
##  [29] "Sybu"          "Pmm1"          "Wnt7b"         "Six3"         
##  [33] "Esco1"         "Tpgs2"         "Celf4"         "Vldlr"        
##  [37] "Spats2l"       "Klf7"          "Acadl"         "Ackr3"        
##  [41] "D1Ertd622e"    "Pam"           "Nckap5"        "Btg2"         
##  [45] "B3galt2"       "Xpr1"          "Gm5532"        "Rasal2"       
##  [49] "Atp1b1"        "Dusp10"        "Kcnk2"         "Nenf"         
##  [53] "Nacc2"         "Notch1"        "Pbx3"          "Gsn"          
##  [57] "Lypd6"         "Tnfaip6"       "Ermn"          "Gca"          
##  [61] "Fign"          "Dlx2"          "Itga6"         "Fam171b"      
##  [65] "Lrp4"          "Tspan18"       "Sema6d"        "Slc24a5"      
##  [69] "Myef2"         "Acss1"         "Stmn3"         "Fgf13"        
##  [73] "Pou3f4"        "Usp51"         "Tnik"          "Qrfpr"        
##  [77] "Anxa5"         "Fat4"          "Hspa4l"        "Slc7a11"      
##  [81] "Lhfp"          "Arhgef26"      "Ptx3"          "Fstl5"        
##  [85] "Ctso"          "Sfrp2"         "Nes"           "Igsf3"        
##  [89] "Rhoc"          "Ank2"          "Cxxc4"         "Ddit4l"       
##  [93] "Lmo4"          "Car8"          "Rgs3"          "Cntln"        
##  [97] "Jun"           "Mycl"          "Kif1b"         "Akap9"        
## [101] "Cacna2d1"      "Napepld"       "Fam126a"       "Actr3b"       
## [105] "3110082J24Rik" "Fam193a"       "Klf3"          "Chic2"        
## [109] "Gsx2"          "Epha5"         "Shroom3"       "Mn1"          
## [113] "Gatsl2"        "St7"           "Zfp800"        "Cpa4"         
## [117] "Mest"          "Ptn"           "Ndnf"          "Tcf7l1"       
## [121] "Lrrtm1"        "Prokr1"        "Wnt7a"         "Frmd4b"       
## [125] "Gpr19"         "Bbc3"          "Ttc9b"         "Zfp60"        
## [129] "Gm26604"       "Dbx1"          "Nr2f2"         "Gm2115"       
## [133] "Prss23"        "Lmo1"          "Dkk3"          "Fgfr2"        
## [137] "Fam53b"        "Adgra1"        "Dusp8"         "Cited2"       
## [141] "Map7"          "Ahi1"          "Marcks"        "Hsf2"         
## [145] "Smpdl3a"       "Jmjd1c"        "Prmt2"         "Nuak1"        
## [149] "Timp3"         "Cep290"        "Phlda1"        "Sox1"         
## [153] "Csmd1"         "Sfrp1"         "Zfp703"        "Nrg1"         
## [157] "Rwdd4a"        "Palld"         "Pgpep1"        "Neto2"        
## [161] "Tppp3"         "Abhd4"         "Ajuba"         "Dpysl2"       
## [165] "Pcdh9"         "Mycbp2"        "Sox21"         "4931406C07Rik"
## [169] "Bmper"         "Ntm"           "Ubash3b"       "Il18"         
## [173] "Herc1"         "Filip1"        "Rbp1"          "Cspg5"        
## [177] "Stac"          "Meis1"         "Sqstm1"        "Sparc"        
## [181] "Hist3h2ba"     "Myo18a"        "Spag9"         "Rdm1"         
## [185] "Tubb2a"        "Cxcl14"        "Ntrk2"         "Fam172a"      
## [189] "Nr2f1"         "Fam84a"        "Id2"           "Meg3"         
## [193] "Itgb8"         "Rai14"         "Sema5a"        "Nell2"        
## [197] "Tmem106c"      "Itgb5"         "Epha3"         "Robo2"        
## [201] "Jam2"          "Eva1c"         "Sox8"          "Ptchd4"       
## [205] "Dlgap1"        "Tmem178"       "Zeb1"          "Cxxc5"        
## [209] "Cdo1"          "Sema6a"        "Megf10"        "Setbp1"       
## [213] "Cnih2"         "Efemp2"        "Gng3"          "Gnaq"         
## [217] "Rorb"          "Fgfbp3"        "Nkx2-3"

Radial.Glia.data <- AddModuleScore(Radial.Glia.data,
                                  genes.list = genes.list,
                                  genes.pool = NULL,
                                  n.bin = 5,
                                  seed.use = 1,
                                  ctrl.size = length(genes.list),
                                  use.k = FALSE,
                                  enrich.name = "E12VP_signature",
                                  random.seed = 1)

# DP enriched genes at E12
genes.list <- Sgenes %>%
              filter(Gene.Clusters %in% c("Clust.4")) %>%
              pull(Gene) %>%
              as.character()

genes.list <- list(genes.list[genes.list %in% rownames(Radial.Glia.data@data)])

genes.list

## [[1]]
##  [1] "Igfbp2"   "Lrrfip1"  "Tgfb2"    "Fjx1"     "Lmo2"     "Mpped2"  
##  [7] "Gpc4"     "Gria2"    "Igfbpl1"  "Gm11266"  "Ubc"      "Aacs"    
## [13] "Cav1"     "Serpinh1" "Ifitm2"   "Dhcr7"    "Rspo3"    "Cd24a"   
## [19] "Hsp90b1"  "Rmst"     "Cd63"     "Calr"     "Mt1"      "Fezf2"   
## [25] "Gnl3"     "Dct"      "Cdon"     "Gsta4"    "Mat2b"    "Slc16a3" 
## [31] "Hist1h3g" "Hist1h1d" "Ctsl"     "Mycn"     "Sp8"      "Alcam"   
## [37] "Gm26917"  "Gm42418"  "Eif1a"    "Wnt8b"

Radial.Glia.data <- AddModuleScore(Radial.Glia.data,
                                  genes.list = genes.list,
                                  genes.pool = NULL,
                                  n.bin = 5,
                                  seed.use = 1,
                                  ctrl.size = length(genes.list),
                                  use.k = FALSE,
                                  enrich.name = "E12DP_signature",
                                  random.seed = 1)

Maturation scores

Tgenes <- read.table("./Progenitors/Temporal.gene.clusters.csv", sep = ";", header = T, row.names = 1)

# VP enriched genes at E12
genes.list <- Tgenes %>%
              filter(Gene.Clusters %in% c("Clust.2")) %>%
              pull(Gene) %>%
              as.character()

genes.list <- list(genes.list[genes.list %in% rownames(Radial.Glia.data@data)])

genes.list

## [[1]]
##  [1] "Nop58"    "Set"      "Hspa5"    "Lhx2"     "Flrt3"    "Wfdc2"   
##  [7] "Bex1"     "Wls"      "Bach2"    "Slc2a1"   "Mcm7"     "Kdelr2"  
## [13] "Tmem176b" "Emx1"     "Lrrn1"    "Socs2"    "Isyna1"   "Syce2"   
## [19] "Mt2"      "Chst2"    "Shmt1"    "Aldoc"    "Gas1"     "H2-K1"   
## [25] "Tgif1"    "Msh6"     "Dmrt3"    "Emx2"

Radial.Glia.data <- AddModuleScore(Radial.Glia.data,
                                  genes.list = genes.list,
                                  genes.pool = NULL,
                                  n.bin = 5,
                                  seed.use = 1,
                                  ctrl.size = length(genes.list),
                                  use.k = FALSE,
                                  enrich.name = "Young_signature",
                                  random.seed = 1)

# DP enriched genes at E12
genes.list <- Tgenes %>%
              filter(Gene.Clusters %in% c("Clust.1")) %>%
              pull(Gene) %>%
              as.character()

genes.list <- list(genes.list[genes.list %in% rownames(Radial.Glia.data@data)])

genes.list

## [[1]]
##  [1] "Gm29260"   "Pantr1"    "Zdbf2"     "Pea15a"    "Zeb2"      "Pax6"     
##  [7] "Insm1"     "Tfap2c"    "Cdh4"      "Nkain4"    "Gpm6b"     "Veph1"    
## [13] "Lxn"       "Trim2"     "Cyr61"     "Ddah1"     "Tox"       "Pou3f2"   
## [19] "Pdpn"      "Limch1"    "Ptprz1"    "Gng12"     "Prkcb"     "Ypel3"    
## [25] "Gpm6a"     "Ncan"      "Tox3"      "Mt3"       "Necab2"    "Rgcc"     
## [31] "Pcdh8"     "Ednrb"     "Elavl3"    "Scg3"      "Myo6"      "Hmgn3"    
## [37] "Eomes"     "Eif1b"     "Gm11627"   "Sox9"      "Hist1h2ap" "Atxn1"    
## [43] "Gadd45g"   "Cplx2"     "Gm17750"   "Nrcam"     "Rtn1"      "Acot1"    
## [49] "Slc1a3"    "Myo10"     "Ctnnd2"    "Snhg18"    "Ncald"     "Kif21a"   
## [55] "Lima1"     "Gap43"     "Zbtb20"    "Ccdc80"    "Cadm2"     "Lix1"     
## [61] "Vit"       "Sncaip"    "Sall3"     "Asrgl1"    "Plce1"

Radial.Glia.data <- AddModuleScore(Radial.Glia.data,
                                  genes.list = genes.list,
                                  genes.pool = NULL,
                                  n.bin = 5,
                                  seed.use = 1,
                                  ctrl.size = length(genes.list),
                                  use.k = FALSE,
                                  enrich.name = "Mature_signature",
                                  random.seed = 1)

Plots

color.pal <- c("#440154FF", "#443A83FF", "#31688EFF", "#21908CFF", "#35B779FF", "#8FD744FF", "#FDE725FF", "#F98C0AFF", "#BB3754FF", "#56106EFF", "#000004FF")

FeaturePlot(object = Radial.Glia.data,
            features.plot = c("E12VP_signature1", "E12DP_signature1", "Young_signature1", "Mature_signature1"),
            cols.use =   color.pal,
            reduction.use = "umap",
            no.axes = T,
            no.legend = T,
            overlay = F,
            dark.theme = F
            )

Plot spatial and temporal marker genes

plot <- FeaturePlot(object = Radial.Glia.data,
                    features.plot =  c("Zcchc12", "Qrfpr", "Sfrp2",
                                       "Mpped2", "Zbtb20", "Gm29260",
                                       "Tfap2c", "Lrrn1", "Nrg1"),
                    cols.use = c("grey90", brewer.pal(9,"YlGnBu")),
                    reduction.use = "umap",
                    no.legend = T,
                    overlay = F,
                    dark.theme = F,
                    do.return =T,
                    no.axes = T)

for (i in 1:length(plot)){
  plot[[i]]$data <- plot[[i]]$data[order(plot[[i]]$data$gene),]
}

cowplot::plot_grid(plotlist = plot[1:9], ncol = 3)

Session Info

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

## [1] "03 mai, 2021, 11,23"

#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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] limma_3.42.0          Rphenograph_0.99.1    igraph_1.2.5         
##  [4] matrixStats_0.55.0    RColorBrewer_1.1-2    patchwork_0.0.1      
##  [7] ggExtra_0.9           dplyr_0.8.3           reticulate_1.18      
## [10] scrattch.hicat_0.0.16 Seurat_2.3.4          Matrix_1.2-17        
## [13] cowplot_1.0.0         ggplot2_3.2.1         loomR_0.2.1.9000     
## [16] hdf5r_1.3.2.9000      R6_2.4.1             
## 
## loaded via a namespace (and not attached):
##   [1] Rtsne_0.15          colorspace_1.4-1    class_7.3-17       
##   [4] modeltools_0.2-22   ggridges_0.5.1      mclust_5.4.5       
##   [7] htmlTable_1.13.2    base64enc_0.1-3     rstudioapi_0.11    
##  [10] proxy_0.4-23        farver_2.0.1        npsurv_0.4-0       
##  [13] flexmix_2.3-15      bit64_4.0.2         codetools_0.2-16   
##  [16] splines_3.6.3       R.methodsS3_1.7.1   lsei_1.2-0         
##  [19] robustbase_0.93-5   knitr_1.26          zeallot_0.1.0      
##  [22] jsonlite_1.7.0      Formula_1.2-3       ica_1.0-2          
##  [25] cluster_2.1.0       kernlab_0.9-29      png_0.1-7          
##  [28] R.oo_1.23.0         shiny_1.4.0         compiler_3.6.3     
##  [31] httr_1.4.1          backports_1.1.5     fastmap_1.0.1      
##  [34] assertthat_0.2.1    lazyeval_0.2.2      later_1.0.0        
##  [37] lars_1.2            acepack_1.4.1       htmltools_0.5.0    
##  [40] tools_3.6.3         gtable_0.3.0        glue_1.4.1         
##  [43] reshape2_1.4.3      RANN_2.6.1          rappdirs_0.3.1     
##  [46] Rcpp_1.0.5          vctrs_0.2.0         gdata_2.18.0       
##  [49] ape_5.3             nlme_3.1-141        iterators_1.0.12   
##  [52] fpc_2.2-3           gbRd_0.4-11         lmtest_0.9-37      
##  [55] xfun_0.18           stringr_1.4.0       mime_0.7           
##  [58] miniUI_0.1.1.1      lifecycle_0.1.0     irlba_2.3.3        
##  [61] gtools_3.8.1        DEoptimR_1.0-8      MASS_7.3-53        
##  [64] zoo_1.8-6           scales_1.1.0        promises_1.1.0     
##  [67] doSNOW_1.0.18       parallel_3.6.3      yaml_2.2.1         
##  [70] pbapply_1.4-2       gridExtra_2.3       rpart_4.1-15       
##  [73] segmented_1.0-0     latticeExtra_0.6-28 stringi_1.4.6      
##  [76] foreach_1.4.7       checkmate_1.9.4     caTools_1.17.1.2   
##  [79] bibtex_0.4.2        Rdpack_0.11-0       SDMTools_1.1-221.1 
##  [82] rlang_0.4.7         pkgconfig_2.0.3     dtw_1.21-3         
##  [85] prabclus_2.3-1      bitops_1.0-6        evaluate_0.14      
##  [88] lattice_0.20-41     ROCR_1.0-7          purrr_0.3.3        
##  [91] labeling_0.3        htmlwidgets_1.5.1   bit_4.0.4          
##  [94] tidyselect_0.2.5    plyr_1.8.4          magrittr_1.5       
##  [97] snow_0.4-3          gplots_3.0.1.1      Hmisc_4.3-0        
## [100] pillar_1.4.2        foreign_0.8-72      withr_2.1.2        
## [103] fitdistrplus_1.0-14 mixtools_1.1.0      survival_2.44-1.1  
## [106] nnet_7.3-14         tsne_0.1-3          tibble_2.1.3       
## [109] crayon_1.4.0        KernSmooth_2.23-15  rmarkdown_2.5      
## [112] grid_3.6.3          data.table_1.12.6   metap_1.1          
## [115] digest_0.6.25       diptest_0.75-7      xtable_1.8-4       
## [118] httpuv_1.5.2        tidyr_1.0.0         R.utils_2.9.0      
## [121] stats4_3.6.3        munsell_0.5.0

---

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