Investigation of the apical progenitors diversity

Load libraries and QCFiltered dataset

# Load libraries
library(Seurat)
library(monocle)
library(princurve)
library(cluster)
library(parallel)
library(scrattch.hicat)
library(scrattch.vis)
library(ggplot2)
library(ggExtra)
library(dplyr)
library(reshape)
library(gridExtra)
library(RColorBrewer)
library(wesanderson)
library(viridis)

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

# Load the full annotated dataset
Allcells.data <- readRDS("./QC.filtered.cells.RDS")

Extract Apical progenitors and compute Pseudo-DorsoVentral score

Extract apical progenitors and filter the gene expression matrix

We perform K-means clustering on the 2 cell state scores :

• Apical progenitors AP
• Basal progenitors BP

set.seed(100)
# K-means clustering based on AP, BP scores across cells
cl <- kmeans(cbind(Allcells.data@meta.data$AP_signature1, Allcells.data@meta.data$BP_signature1), 3)
Allcells.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)

col.pal <- wes_palette("GrandBudapest1", 3, 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") ; rm(p1)

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 Apical progenitor as the cluster with the highest mean Apical progenitor signature

# Find clusters wiht the highest mean APscore
MeanKclust.APscore <- aggregate(AP_signature1 ~ kmeanClust, Allcells.data@meta.data, mean)
APclust <- MeanKclust.APscore %>% filter(AP_signature1 == max(AP_signature1)) %>% pull(kmeanClust)

# Extract apical progenitors 
barcodes <- Allcells.data@meta.data %>% filter(kmeanClust == APclust) %>% pull(Barcodes)
AP.data <-  SubsetData(Allcells.data, cells.use = barcodes , subset.raw = T,  do.clean = F)

# Further filter the 3 outlier cells based on spring coordinates
cells <- rownames(AP.data@dr$spring@cell.embeddings[AP.data@dr$spring@cell.embeddings[,2] > 250,])
AP.data <-  SubsetData(AP.data, cells.use = cells , subset.raw = T,  do.clean = F)

DimPlot(AP.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)

Filter gene counts matrix

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

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

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

AP.data <- ScaleData(object = AP.data, vars.to.regress = c("CC.Difference","percent.mito", "nUMI"), display.progress = F)

Import Spring dimensionality reduction

Cell cycle associated genes were excluded for PCA dimensionality reduction and Spring plot was generating using these parameters:

Number of cells: 1648
Number of genes that passed filter: 857
Min expressing cells (gene filtering): 3
Min number of UMIs (gene filtering): 3
Gene variability %ile (gene filtering): 90
Number of principal components: 7
Number of nearest neighbors: 20
Number of force layout iterations: 500

# Import Spring coordinates calculated without cell cycle genes
Coordinates <- read.table("./Progenitors/E12.AP.Coordinates.txt", sep=",", header = F)[,c(2,3)]
rownames(Coordinates) <- rownames(AP.data@meta.data)
AP.data <- SetDimReduction(AP.data,
                           reduction.type = "spring.AP",
                           slot = "cell.embeddings",
                           new.data = as.matrix(Coordinates))
AP.data@dr$spring.AP@key <- "spring.AP"
colnames(AP.data@dr$spring.AP@cell.embeddings) <- paste0(GetDimReduction(object = AP.data, reduction.type = "spring.AP",slot = "key"), c(1,2))

Fit a principal curve over the AP in the Spring space

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

# Fit the principal curve
fit <- principal_curve(as.matrix(data[,1:2]),
                       smoother='lowess',
                       trace=TRUE,
                       f = 0.7,
                       stretch=0,
                       plot_iterations = F)

## Starting curve---distance^2: 2022449778
## Iteration 1---distance^2: 745172
## Iteration 2---distance^2: 696270.5
## Iteration 3---distance^2: 677509.1
## Iteration 4---distance^2: 668251.7
## Iteration 5---distance^2: 663393.6
## Iteration 6---distance^2: 660928.6
## Iteration 7---distance^2: 659586.4
## Iteration 8---distance^2: 658917.8
## Iteration 9---distance^2: 658638.8

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

data$Phase <- as.character(AP.data@meta.data$Phase)
data$DorsoVentral.Score <- DorsoVentral.Score

# Direction of the maturation score
if (cor(data$DorsoVentral.Score, AP.data@data['Zbtb20', ]) > 0) { data$DorsoVentral.Score <- -(data$DorsoVentral.Score - max(data$DorsoVentral.Score))}

AP.data@meta.data$DorsoVentral.Score <- data$DorsoVentral.Score

# Plot Cell onto PC1 and PC2 with principal curve
ggplot(data, aes(springAP.1, springAP.2)) + 
  geom_point(aes(color=Phase), size=2, shape=16, values=col.pal) +
  geom_line(data=pc.line, color='red', size=0.77)

# Plot pseudotime color gradient on the cell cycle filtered Spring embbeding
ggplot(data, aes(springAP.1, springAP.2)) +
  geom_point(aes(color=DorsoVentral.Score), size=2, shape=16) +
  scale_color_viridis(direction = -1) +
  geom_line(data=pc.line, color='red', size=0.77)

Manuscript Fig. 5A

# Plot pseudotime color gradient on the Spring embbeding calculated from full dataset
ggplot(data, aes(spring1, spring2)) + 
  geom_point(aes(color=DorsoVentral.Score), size=2, shape=16) +
  scale_color_viridis(direction = -1) 

Find differentially expressed genes along the pseudo DV axis

Initialize a monocle object

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

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

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

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

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

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

Test each gene trend over pseudo-DV score

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

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

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

Smooth significative gene expression pseudo-DV axis

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

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

Cluster cells by similar expression profiles on the pseudo-DV axis

set.seed(100)
# Cluster cells using the Partitioning Around Medoids algorithm
Cells.Clust <- pam(as.dist((1 - cor(Smooth.curve.matrix ,method = "spearman"))), k=7)
Domaines.Clust <- data.frame(Domaines = paste0("Clust.",Cells.Clust$clustering))

Cluster genes by similar profiles over pseudo-DV axis

Perform hierarchical clustering

# Build genes dendrogram
dst <- dist(t(scale(t(Smooth.curve.matrix))), method = "manhattan")
hc <- hclust(dst, method = "ward.D")

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

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

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

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

write.table(Gene.dynamique, "./Progenitors/Gene.dynamique.csv", sep = ";", quote = F)

Plot heatmap

anno.colors <- list(Domaines = c(Clust.1="#83c3b8", Clust.2="#009fda", Clust.3="#3e69ac", Clust.4="#e46b6b", Clust.5="#e3c148", Clust.6="#b7d174", Clust.7="#68b041"),
                    Gene.Clusters = c(Clust.1 ="#ec756d" , Clust.2="#ebcb2e", Clust.3="#b79f0b", Clust.4="#7293c8", Clust.5="#3ca73f",
                                      Clust.6="#9ec22f", Clust.7="#cc3a1b", Clust.8="#4cabdc", Clust.9="#cc8778", Clust.10="#4784a2"))

pheatmap::pheatmap(Smooth.curve.matrix,
                   cluster_rows= as.hclust(dend),
                   scale = "row",
                   cluster_cols = F,
                   gaps_col = cumsum(as.numeric(table(Domaines.Clust$Domaines))),
                   annotation_row = Gene.dynamique %>% dplyr::select(Gene.Clusters),
                   annotation_col = Domaines.Clust,
                   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 = "Genes expression along Dorso-Ventral axis")

Manuscript Fig. 5C

For visualization purposes we manually rotated the two last branches

Plot gene clusters trends

source("./functions/GenesTrendPlots.R")
Clusters.trend(AP.data,
               Which.cluster = 1:10,
               clust.list = clusters,
               group.by = "global",
               span = 1,
               Smooth.method = "auto",
               Use.scale.data = T)

Manuscript Fig. 5C

Assign domain identity

We assign domains identity based on clusters’ transcriptional profile by setting boundaries over pseudo-DV score

# Set the boundaries over pseudotime score
new_data$cluster <- Domaines.Clust$Domaines
Infered.Domain.boundary <- aggregate(DorsoVentral.Score ~ cluster, new_data, max) %>% pull(DorsoVentral.Score)

# Assign identities based on the position of the cell on the pseudo-dv axis
Domaine.Ident <- sapply(AP.data@meta.data$DorsoVentral.Score,
                        function(x){ if(x<Infered.Domain.boundary[1]){ x = "Sub.Pallium.1"
                        } else if(x> Infered.Domain.boundary[1] & x< Infered.Domain.boundary[2]){ x ="Sub.Pallium.2"
                        } else if(x> Infered.Domain.boundary[2] & x< Infered.Domain.boundary[3]){ x = "Sub.Pallium.3" 
                        } else if(x> Infered.Domain.boundary[3] & x< Infered.Domain.boundary[4]){ x = "Ventral.Pallium"
                        } else if(x> Infered.Domain.boundary[4] & x< Infered.Domain.boundary[5]){ x = "lateral.Pallium.1"
                        } else if(x> Infered.Domain.boundary[6]){ x = "Dorsal.Pallium"
                        } else x="lateral.Pallium.2"})

# Transfert the identity to the Seurat object
AP.data@meta.data$Domaine <- Domaine.Ident
AP.data <- SetAllIdent(AP.data, id = "Domaine")

Plot.Genes.trend(AP.data,
                 genes = c("Gsx2", "Dbx1", "Gm29260", "Tfap2c", "Emx1", "Lrrn1"),
                 Use.scale.data = F)

Manuscript Fig. 5D

DimPlot(AP.data,
         group.by = "Domaine",
         reduction.use = "spring",
         dim.1 = 1,
         dim.2 = 2,
         do.label=F,
         label.size = 4,
         no.legend = F,
         cols.use = tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B"))
         )

Plot Figure 5B

# Load full dataset
Allcells.data <- readRDS("./QC.filtered.cells.RDS")

# Transfer the identities
Rename.Clust <-  function(Clustdata, RawQCdata) {
  unClustered.cells <- RawQCdata@meta.data$Barcodes
  RawQCdata <- SetIdent(RawQCdata, cells.use = unClustered.cells, ident.use = "All.Unclustered.Cells")
  
  for(i in unique(Clustdata@meta.data$Domaine)){
  New.ident <- i
  Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$Domaine == i))
  print(paste0("Cluster_",i,": ",length(Barcodes), " Cells"))
  Barcodes <- Barcodes[Barcodes %in% rownames(RawQCdata@meta.data)]
  RawQCdata <- SetIdent(RawQCdata, cells.use = Barcodes ,ident.use = paste0("AP.",i))
  }
  return(RawQCdata)
}

Allcells.data <- Rename.Clust(Clustdata = AP.data, RawQCdata = Allcells.data)

## [1] "Cluster_Sub.Pallium.2: 363 Cells"
## [1] "Cluster_Ventral.Pallium: 312 Cells"
## [1] "Cluster_Sub.Pallium.3: 197 Cells"
## [1] "Cluster_lateral.Pallium.1: 219 Cells"
## [1] "Cluster_Dorsal.Pallium: 201 Cells"
## [1] "Cluster_lateral.Pallium.2: 113 Cells"
## [1] "Cluster_Sub.Pallium.1: 241 Cells"

colors <-  c("#969696",tolower(c("#68B041", "#E3C148", "#B7D174", "#83C3B8", "#009FDA", "#3E69AC", "#E46B6B")))

DimPlot(Allcells.data,
        reduction.use = "spring", 
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        label.size = 2,
        no.legend = T,
        cols.use = colors)

Manuscript Fig. 5B

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

Plot the representative gene of the Figure S6

# Load custom plotting functions
source("./functions/GenesTrendPlots.R")

Plot.Genes.trend(AP.data,
                 genes = c("Dlx1", "Ano1", "Dlx2", "Olig2"),
                 Use.scale.data = F)

Manuscript Fig.S6A

Plot.Genes.trend(AP.data,
                 genes = c("Otx2", "Zbtb20", "Six3", "Ascl1"),
                 Use.scale.data = F)

Manuscript Fig.S6B

Plot.Genes.trend(AP.data,
                 genes = c("Meis1", "Rorb", "Epha3"),
                 Use.scale.data = F)

Manuscript Fig.S6C

Plot.Genes.trend(AP.data,
                 genes = c("Sfrp2", "Eva1c", "Sema5a", "Fat4"),
                 Use.scale.data = F)

Manuscript Fig.S6D

Plot.Genes.trend(AP.data,
                 genes = c("Dmrta2","Lypd6", "Pax6", "Dmrta1"),
                 Use.scale.data = F)

Manuscript Fig.S6E

Plot.Genes.trend(AP.data,
                 genes = c("Emx1", "Arx", "Lmo3"),
                 Use.scale.data = F)

Manuscript Fig.S6F

Plot.Genes.trend(AP.data,
                 genes = c("Fezf2", "Emx2", "Sp8"),
                 Use.scale.data = F)

Manuscript Fig.S6G

Transfert identities on the full dataset (Manuscript Fig. 2A)

# load full dataset
Allcells.data <- readRDS("./Clustered.cells.RDS")

#Transfer the identities
Rename.Clust <-  function(Clustdata, RawQCdata) {
  for(i in unique(Clustdata@meta.data$Domaine)){
    New.ident <- i
    Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$Domaine == i))
    print(paste0("Cluster_",i,": ",length(Barcodes), " Cells"))
    Barcodes <- Barcodes[Barcodes %in% rownames(RawQCdata@meta.data)]
    RawQCdata <- SetIdent(RawQCdata, cells.use = Barcodes ,ident.use = paste0("AP.",i))
  }
  return(RawQCdata)
}

Allcells.data <- Rename.Clust(Clustdata = AP.data, RawQCdata = Allcells.data)

## [1] "Cluster_Sub.Pallium.2: 363 Cells"
## [1] "Cluster_Ventral.Pallium: 312 Cells"
## [1] "Cluster_Sub.Pallium.3: 197 Cells"
## [1] "Cluster_lateral.Pallium.1: 219 Cells"
## [1] "Cluster_Dorsal.Pallium: 201 Cells"
## [1] "Cluster_lateral.Pallium.2: 113 Cells"
## [1] "Cluster_Sub.Pallium.1: 241 Cells"

colors2 <-  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 = T,
        cols.use = colors2)

Manuscript Fig. 2A

Allcells.data <- StashIdent(Allcells.data, save.name = "Cluster.ident")
saveRDS(Allcells.data, "./Clustered.cells.RDS")
saveRDS(AP.data, "./AP.data.RDS")

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

Perform iterative graph-based clustering to find if discrete pallial sub-domains exist

Extract pallial progenitors

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

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

Perform iterative clustering implemented in scrattch.hicat package

For more detail on the scrattch.hicat please refer to the package page.

Prepare the dataset for clustering with scrattch.hicat

Independent gene filtering

We decide 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 = "\t", header = F)[,1])

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

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(AP.data@raw.data)[rownames(AP.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("AP.data", "Data.matrix")])

Select features used to remove PCs which would correlate with unwanted sources of variation

gene.counts <- log2(colSums(as.matrix(Data.matrix$norm.dat) > 0))
nUMI <- log2(colSums(Data.matrix$raw.dat))
perctMito <- AP.data@meta.data$percent.mito
perctRibo <- AP.data@meta.data$percent.ribo
PseudoDVscore <- AP.data@meta.data$DorsoVentral.Score

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

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

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

rm(gene.counts, nUMI, perctMito, perctRibo,PseudoDVscore)

Iterative clustering

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

Run the iterative clustering

The default iter_clust function use in this version of the scrattch.hicat package does not allow to set the k.param argument. We modified this function to allow this argument to be set to other values.

# Perform the iterative clustering 
iter.result <- iter_clust(norm.dat, 
                          counts = raw.dat,
                          dim.method = "pca",
                          max.dim = 15,
                          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 ~ 0.022 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.047 s
##   Build undirected graph from the weighted links...DONE ~ 0.052 s
##   Run louvain clustering on the graph ...DONE ~ 0.031 s
##   Return a community class
##   -Modularity value: 0.6822028 
##   -Number of clusters: 12[1] "test-iter_clust.1"
##   Finding nearest neighbors...DONE ~ 0.016 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.043 s
##   Build undirected graph from the weighted links...DONE ~ 0.046 s
##   Run louvain clustering on the graph ...DONE ~ 0.021 s
##   Return a community class
##   -Modularity value: 0.6939065 
##   -Number of clusters: 12[1] "test-iter_clust.2"
##   Finding nearest neighbors...DONE ~ 0.001 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.008 s
##   Build undirected graph from the weighted links...DONE ~ 0.008 s
##   Run louvain clustering on the graph ...DONE ~ 0.003 s
##   Return a community class
##   -Modularity value: 0.7585217 
##   -Number of clusters: 7

Merge clusters not seperable by DEGs

# Merge clusters which are not seperable 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")

## 2  Clusters

cat(length(unique(merge.result$markers)), " DE genes")

## 10  DE genes

Import cell labels into the Seurat obj

AP.data@ident <- as.factor(merge.result$cl)
AP.data@meta.data$Pallial.subDomaine <-as.character(AP.data@ident)

colors <-  c("#969696", "#e7823a")

DimPlot(AP.data,
        reduction.use = "spring",
        dim.1 = 1,
        dim.2 = 2,
        do.label=T,
        pt.size = 2,
        label.size = 4,
        no.legend = F,
        cols.use = colors)

Find sub-domain marker genes

# Find genes which best distinguish the subdomains
Allmarkers <- FindAllMarkers(object = AP.data,
                             test.use = "roc",
                             min.pct = 0.3,
                             logfc.threshold = 0.6,
                             print.bar = T,
                             only.pos = T)

Allmarkers$q.diff <- abs(Allmarkers$pct.1 - Allmarkers$pct.2)/max(Allmarkers$pct.1, Allmarkers$pct.2)
Allmarkers$Specificity.index <- Allmarkers$avg_logFC/ (1-Allmarkers$q.diff)

Selected.markers <- Allmarkers %>%
  filter(Specificity.index > 1) %>% pull(gene)

# Prepare annotation for plotting
colors <-  c("#969696", "#e7823a")
colorsident <- cbind(ident = unique(as.character(AP.data@ident)),
                     colors = colors,
                     id = unique(as.character(AP.data@ident)))

# Create annotation data.frame
anno.df <- as.data.frame(cbind(
  sample_name = row.names(AP.data@meta.data),
  primary_type_id = colorsident[match(as.character(AP.data@ident), colorsident[,1]),3],
  primary_type_label = as.character(AP.data@ident),
  primary_type_color = colorsident[match(as.character(AP.data@ident), colorsident[,1]),2]
))


# Plot cells barplot
data <- cbind(sample_name = colnames(AP.data@data),
              as.data.frame(t(as.matrix(AP.data@data))))

sample_bar_plot(data, 
                anno.df, 
                genes = Selected.markers,
                group_order = c("1","3"),
                grouping = "primary_type",
                log_scale = FALSE,
                font_size = 7,
                label_height = 10,
                label_type = "angle",
                bg_color ="#f7f7f7")

Session Info

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

## [1] "30 novembre, 2020, 10,52"

#Packages used
sessionInfo()

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

---

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