We extract the most mature pallial cells in two steps :
# Cell state scores use to select
FeaturePlot(object = QCFiltered.data,
features.plot = c("AP_signature1", "SP_signature1", "Pal_signature1"),
cols.use = rev(brewer.pal(10,"Spectral")),
reduction.use = "spring",
no.legend = T,
overlay = F,
dark.theme = F
)
We perform Kmeans clustering on the 3 cell state scores :
AP
SP
Pal
# K-means clustering based on AP, SP and Pal signature scores
set.seed(100)
cl <- kmeans(cbind(QCFiltered.data@meta.data$AP_signature1,
QCFiltered.data@meta.data$SP_signature1,
QCFiltered.data@meta.data$Pal_signature1), 3)
QCFiltered.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
col.pal <- wes_palette("GrandBudapest1", 3, type = "discrete")
p1 <- ggplot(QCFiltered.data@meta.data, aes(x=Pal_signature1, y=SP_signature1, colour = kmeanClust)) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey")
DimPlot(QCFiltered.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 = T)
We then extract the glutamatergic neuron branch as the cluster with the highest mean Pallial neurons
signature
# Find the k-means cluster with the highest mean Pal_signature1 score
MeanKclust.Palscore <- aggregate(Pal_signature1 ~ kmeanClust, QCFiltered.data@meta.data, mean)
Palclust <- MeanKclust.Palscore %>% filter(Pal_signature1 == max(Pal_signature1)) %>% pull(kmeanClust)
# Extract barcodes and filter the seurat object
Glut.cells <- QCFiltered.data@meta.data %>% filter(kmeanClust == Palclust) %>% pull(Barcodes)
Glut.QCFiltered.data <- SubsetData(QCFiltered.data, cells.use = Glut.cells , subset.raw = T, do.clean = F)
# K-means clustering based on EN and LN signature scores
cl <- kmeans(cbind(Glut.QCFiltered.data@meta.data$EN_signature1, Glut.QCFiltered.data@meta.data$LN_signature1), 2)
Glut.QCFiltered.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
# Find the k-means cluster with the highest mean Late_signature1 score
MeanKclust.LNscore <- aggregate(LN_signature1 ~ kmeanClust, Glut.QCFiltered.data@meta.data, mean)
LN.clust <- MeanKclust.LNscore %>% filter(LN_signature1 == max(LN_signature1)) %>% pull(kmeanClust)
# Extract barcodes
LN.cells <- Glut.QCFiltered.data@meta.data %>% filter(kmeanClust == LN.clust) %>% pull(Barcodes)
col.pal <- wes_palette("GrandBudapest1", 2, type = "discrete")
p1 <- ggplot(Glut.QCFiltered.data@meta.data, aes(x=LN_signature1, y=EN_signature1, colour = kmeanClust)) +
annotate(geom = "text", label = paste0(length(LN.cells), " Late Neurons"), x = 1, y =1.5) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey") ; rm(p1)
DimPlot(Glut.QCFiltered.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 = T)
We then extract the glutamatergic neuron with the highest mean Late neurons
signature
# Filter the seurat object
Glut.QCFiltered.data <- SubsetData(Glut.QCFiltered.data, cells.use = LN.cells , subset.raw = T, do.clean = F)
For more detail on the scrattch.hicat please refer to the package page.
We decide to exclude Cell cycle, ribosomal and mitochondrial associated genes, as well as Xist for the clustering step.
# Exclude cell cycle associated genes
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")
# Exclude genes detected in less than 3 cells
num.cells <- Matrix::rowSums(Glut.QCFiltered.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", s.genes, g2m.genes)
genes.use <- genes.use[!genes.use %in% GenesToRemove]
dgeMatrix_count <- as.matrix(Glut.QCFiltered.data@raw.data)[rownames(Glut.QCFiltered.data@raw.data) %in% genes.use,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1)
gene.counts <- log2(colSums(as.matrix(Data.matrix$norm.dat) > 0))
nUMI <- log2(colSums(Data.matrix$raw.dat))
perctMito <- Glut.QCFiltered.data@meta.data$percent.mito
perctRibo <- Glut.QCFiltered.data@meta.data$percent.ribo
rm.eigen <- as.matrix(cbind(gene.counts,
nUMI,
perctMito,
perctRibo))
row.names(rm.eigen) <- names(gene.counts)
colnames(rm.eigen) <- c("log2nGenes",
"log2nUMI",
"perctMito",
"perctRibo")
rm(gene.counts, nUMI, perctMito, perctRibo)
# 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 = 30,
min.cells = 10)
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 = 8,
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.006 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.011 s
## Build undirected graph from the weighted links...DONE ~ 0.024 s
## Run louvain clustering on the graph ...DONE ~ 0.012 s
## Return a community class
## -Modularity value: 0.8698811
## -Number of clusters: 21[1] "test-iter_clust.1"
## [1] "test-iter_clust.2"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.3210654
## -Number of clusters: 2[1] "test-iter_clust.3"
## [1] "test-iter_clust.4"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.002 s
## Build undirected graph from the weighted links...DONE ~ 0.004 s
## Run louvain clustering on the graph ...DONE ~ 0.002 s
## Return a community class
## -Modularity value: 0.7275928
## -Number of clusters: 9[1] "test-iter_clust.4.1"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.5114614
## -Number of clusters: 3[1] "test-iter_clust.4.2"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.002 s
## Build undirected graph from the weighted links...DONE ~ 0.004 s
## Run louvain clustering on the graph ...DONE ~ 0.002 s
## Return a community class
## -Modularity value: 0.7128587
## -Number of clusters: 7[1] "test-iter_clust.5"
## [1] "test-iter_clust.6"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.5893365
## -Number of clusters: 4[1] "test-iter_clust.7"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.4688921
## -Number of clusters: 3[1] "test-iter_clust.8"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.6325457
## -Number of clusters: 5[1] "test-iter_clust.9"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.002 s
## Build undirected graph from the weighted links...DONE ~ 0.004 s
## Run louvain clustering on the graph ...DONE ~ 0.002 s
## Return a community class
## -Modularity value: 0.7003194
## -Number of clusters: 6[1] "test-iter_clust.10"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.626505
## -Number of clusters: 4[1] "test-iter_clust.11"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.003 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.7033592
## -Number of clusters: 7[1] "test-iter_clust.12"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.5403066
## -Number of clusters: 3
# 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")
## 13 Clusters
## 481 DE genes
Glut.QCFiltered.data@ident <- as.factor(merge.result$cl)
Glut.QCFiltered.data@meta.data$LN.ident <-as.character(Glut.QCFiltered.data@ident)
colors <- c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#cc8778" , "#d14c8d", "#4cabdc", "#5ab793", "#e7823a","#e6bb9b", "#046c9a", "#4784a2" , "#4990c9")
DimPlot(Glut.QCFiltered.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 4,
no.legend = T,
cols.use = colors)
# 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$LN.ident)){
New.ident <- i
Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$LN.ident == 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("LN.Glut.",i))
}
return(RawQCdata)
}
Allcells.data <- Rename.Clust(Clustdata = Glut.QCFiltered.data, RawQCdata = Allcells.data)
## [1] "Cluster_23: 109 Cells"
## [1] "Cluster_18: 86 Cells"
## [1] "Cluster_24: 45 Cells"
## [1] "Cluster_26: 26 Cells"
## [1] "Cluster_19: 48 Cells"
## [1] "Cluster_14: 16 Cells"
## [1] "Cluster_1: 33 Cells"
## [1] "Cluster_22: 37 Cells"
## [1] "Cluster_25: 66 Cells"
## [1] "Cluster_21: 24 Cells"
## [1] "Cluster_20: 31 Cells"
## [1] "Cluster_16: 24 Cells"
## [1] "Cluster_13: 16 Cells"
Glut.QCFiltered.data@meta.data %>% ggplot(aes(x= LN.ident, y= LN_signature1, colour= LN.ident)) +
scale_color_manual(values=colors[-1]) +
geom_boxplot() +
geom_hline(yintercept = 1)
# Extract clusters with median LN_score > 1
Median_LN_score <- aggregate(LN_signature1 ~ LN.ident, Glut.QCFiltered.data@meta.data, median)
mature_LN <- Median_LN_score %>% filter(LN_signature1 > 1) %>% pull(LN.ident)
Glut.QCFiltered.data <- SubsetData(Glut.QCFiltered.data, ident.use = mature_LN , subset.raw = T, do.clean = F)
colors <- c("#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#d14c8d", "#4cabdc", "#5ab793", "#e7823a", "#046c9a", "#4990c9") #colors[-c(1,5,10,12)]
DimPlot(Glut.QCFiltered.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 5,
no.legend = T,
cols.use = colors)
# Prepare annotation for hicat pipeline
colorsident <- cbind(ident = unique(as.character(Glut.QCFiltered.data@ident)),
colors = colors,
id = unique(as.character(Glut.QCFiltered.data@ident)))
# Create annotation data.frame
anno.df <- as.data.frame(cbind(
sample_name = row.names(Glut.QCFiltered.data@meta.data),
primary_type_id = colorsident[match(as.character(Glut.QCFiltered.data@ident), colorsident[,1]),3],
primary_type_label = as.character(Glut.QCFiltered.data@ident),
primary_type_color = colorsident[match(as.character(Glut.QCFiltered.data@ident), colorsident[,1]),2]
))
# Make a data.frame of unique cluster id, type, color, and broad type
cl.df <- anno.df %>%
select(primary_type_id,
primary_type_label,
primary_type_color) %>%
unique()
colnames(cl.df)[1:3] <- c("cluster_id", "cluster_label", "cluster_color")
# Sort by cluster_id
cl.df <- arrange(cl.df, cluster_id)
row.names(cl.df) <- cl.df$cluster_id
cl.fact <- setNames(factor(anno.df$primary_type_id), anno.df$sample_name)
# Filter genes
num.cells <- Matrix::rowSums(Glut.QCFiltered.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")
genes.use <- genes.use[!genes.use %in% GenesToRemove] ; rm(GenesToRemove, num.cells)
Glut.QCFiltered.data@raw.data <- Glut.QCFiltered.data@raw.data[genes.use, ]
Glut.QCFiltered.data <- NormalizeData(object = Glut.QCFiltered.data,
normalization.method = "LogNormalize",
scale.factor = round(median(Glut.QCFiltered.data@meta.data$nUMI)),
display.progress = F)
# Find all var genes
Glut.QCFiltered.data <- FindVariableGenes(object = Glut.QCFiltered.data,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.02,
x.high.cutoff = 3,
y.cutoff = 1,
do.plot = F, display.progress = F)
dgeMatrix_count <- as.matrix(Glut.QCFiltered.data@raw.data)[rownames(Glut.QCFiltered.data@raw.data) %in% Glut.QCFiltered.data@var.genes,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1) ; rm(dgeMatrix_cpm)
Data.matrix <- list(raw.dat=dgeMatrix_count, norm.dat=norm.dat) ; attach(Data.matrix)
rm(dgeMatrix_count,norm.dat)
Scrattch.hicat perform hierarchical clustering on a cluster correlation matrix based on median expression values for the top 50 most DEGs between every pair of clusters. It estimates branch confidence level using a bootstrap approach implemented by the pvclust package
# Build the dendrogram
dend.result <- build_dend(cl.med[,levels(cl.fact)],
l.color= setNames(as.character(cl.df$cluster_color), row.names(cl.df)),
nboot = 100)
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
# Attach cluster labels to the leaves of the tree
dend.labeled <- dend.result$dend
labels(dend.labeled) <- cl.df[labels(dend.labeled), "cluster_label"]
plot(dend.labeled)
Because some cells display intermediate transcriptomic signatures between two closely related clusters we used only clusters’ signatures on their core cells
# Select marker genes
de.param <- de_param(padj.th = 0.01,
lfc.th = 0.6,
low.th = 0.5,
q1.th = 0.25,
q2.th = NULL,
q.diff.th = 0.6,
de.score.th = 20,
min.cells = 5)
display.result <- display_cl(Glut.QCFiltered.data@ident,
norm.dat,
plot=F,
de.param = de.param,
n.markers = 50)
select.markers <- select_markers(norm.dat,
cl.fact,
de.genes = display.result$de.genes,
n.markers = 50)
marker.genes <- select.markers$markers
# Run the centroide classifier
map.results <- map_sampling(train.dat = norm.dat,
train.cl = Glut.QCFiltered.data@ident,
test.dat = norm.dat,
markers = marker.genes,
markers.perc = 0.8,
iter = 100)
##
Running iteration 1 of 100.
Running iteration 2 of 100.
Running iteration 3 of 100.
Running iteration 4 of 100.
Running iteration 5 of 100.
Running iteration 6 of 100.
Running iteration 7 of 100.
Running iteration 8 of 100.
Running iteration 9 of 100.
Running iteration 10 of 100.
Running iteration 11 of 100.
Running iteration 12 of 100.
Running iteration 13 of 100.
Running iteration 14 of 100.
Running iteration 15 of 100.
Running iteration 16 of 100.
Running iteration 17 of 100.
Running iteration 18 of 100.
Running iteration 19 of 100.
Running iteration 20 of 100.
Running iteration 21 of 100.
Running iteration 22 of 100.
Running iteration 23 of 100.
Running iteration 24 of 100.
Running iteration 25 of 100.
Running iteration 26 of 100.
Running iteration 27 of 100.
Running iteration 28 of 100.
Running iteration 29 of 100.
Running iteration 30 of 100.
Running iteration 31 of 100.
Running iteration 32 of 100.
Running iteration 33 of 100.
Running iteration 34 of 100.
Running iteration 35 of 100.
Running iteration 36 of 100.
Running iteration 37 of 100.
Running iteration 38 of 100.
Running iteration 39 of 100.
Running iteration 40 of 100.
Running iteration 41 of 100.
Running iteration 42 of 100.
Running iteration 43 of 100.
Running iteration 44 of 100.
Running iteration 45 of 100.
Running iteration 46 of 100.
Running iteration 47 of 100.
Running iteration 48 of 100.
Running iteration 49 of 100.
Running iteration 50 of 100.
Running iteration 51 of 100.
Running iteration 52 of 100.
Running iteration 53 of 100.
Running iteration 54 of 100.
Running iteration 55 of 100.
Running iteration 56 of 100.
Running iteration 57 of 100.
Running iteration 58 of 100.
Running iteration 59 of 100.
Running iteration 60 of 100.
Running iteration 61 of 100.
Running iteration 62 of 100.
Running iteration 63 of 100.
Running iteration 64 of 100.
Running iteration 65 of 100.
Running iteration 66 of 100.
Running iteration 67 of 100.
Running iteration 68 of 100.
Running iteration 69 of 100.
Running iteration 70 of 100.
Running iteration 71 of 100.
Running iteration 72 of 100.
Running iteration 73 of 100.
Running iteration 74 of 100.
Running iteration 75 of 100.
Running iteration 76 of 100.
Running iteration 77 of 100.
Running iteration 78 of 100.
Running iteration 79 of 100.
Running iteration 80 of 100.
Running iteration 81 of 100.
Running iteration 82 of 100.
Running iteration 83 of 100.
Running iteration 84 of 100.
Running iteration 85 of 100.
Running iteration 86 of 100.
Running iteration 87 of 100.
Running iteration 88 of 100.
Running iteration 89 of 100.
Running iteration 90 of 100.
Running iteration 91 of 100.
Running iteration 92 of 100.
Running iteration 93 of 100.
Running iteration 94 of 100.
Running iteration 95 of 100.
Running iteration 96 of 100.
Running iteration 97 of 100.
Running iteration 98 of 100.
Running iteration 99 of 100.
Running iteration 100 of 100.
Glut.QCFiltered.data@meta.data$Core_cells <- ifelse(map.results$map.df$prob > 0.95, "Core.Cells", "Intermediate")
DimPlot(Glut.QCFiltered.data,
group.by = "Core_cells",
reduction.use = "spring",
cols.use = wes_palette("Royal1", 2, type = "discrete"),
dim.1 = 1,
dim.2 = 2,
do.label=F,
label.size = 2,
no.legend = F)
# Remove intermediate cells
Core.cells <- Glut.QCFiltered.data@meta.data %>% dplyr::filter(Core_cells == "Core.Cells") %>% pull(Barcodes)
Glut.QCFiltered.data <- SubsetData(Glut.QCFiltered.data, cells.use = Core.cells , subset.raw = T, do.clean = F)
anno.df <- anno.df %>% filter(sample_name %in% Core.cells)
DimPlot(Glut.QCFiltered.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 3,
no.legend = T,
cols.use = colors)
data <- cbind(sample_name = colnames(Glut.QCFiltered.data@data),
as.data.frame(t(as.matrix(Glut.QCFiltered.data@data))))
Selected.markers <- c("Slc17a6","Tbr1", "Gap43", "Foxg1","Mn1","Reln",
"Lhx5","Calb2","Samd3", "Dmrta2", "Nr2f2", "Ebf2",
"Gm27199", "Chl1","Unc5c", "Grm1", "Insm2","Trp73",
"Ebf3", "Cdkn1a", "Lhx1", "Lhx1os","Zfp503", "Cacna2d2",
"Serpine2", "Rspo3", "Grm2", "Cxcr4","Lamp5","Neurod2",
"Neurod6", "Tcf4","Lrfn5", "Fezf1", "Fezf2", "Pcp4",
"Nfia", "Nfib", "Nfix", "Bhlhe22","Ppp1r14c", "Scg2",
"Meis2","Pbx3", "Mab21l1", "Barhl2", "Runx1t1","Mef2c",
"Nxph1", "Meis1", "Tshz2","Lhx9", "Nr4a2", "Cck", "Foxp2",
"Pou3f2", "Sox6","Ebf1","Nrip3", "Zic1", "Etv1", "Tfap2e", "Pax6")
sample_bar_plot(data,
anno.df,
genes = Selected.markers,
group_order = levels(cl.fact),
grouping = "primary_type",
log_scale = FALSE,
font_size = 7,
label_height = 10,
label_type = "angle",
bg_color ="#f7f7f7")
Load QC filtered dataset
As for the glutamaterigic neurons, we extract the most mature sub-pallial cells in two steps :
FeaturePlot(object = QCFiltered.data,
features.plot = c("AP_signature1", "Pal_signature1", "SP_signature1"),
cols.use = rev(brewer.pal(10,"Spectral")),
reduction.use = "spring",
no.legend = T,
overlay = F,
dark.theme = F
)
We perform K-means clustering on the 3 cell state scores :
AP
SP
Pal
# K-means clustering based on AP, SP and Pal signature scores
set.seed(100)
cl <- kmeans(x = cbind(QCFiltered.data@meta.data$AP_signature1,
QCFiltered.data@meta.data$SP_signature1,
QCFiltered.data@meta.data$Pal_signature1),
centers = 3)
QCFiltered.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
col.pal <- wes_palette("GrandBudapest1", 3, type = "discrete")
p1 <- ggplot(QCFiltered.data@meta.data, aes(x=Pal_signature1, y=SP_signature1, colour = kmeanClust)) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey")
DimPlot(QCFiltered.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 = T)
We then extract the GABAergic neuron branche as beeing the K-means cluster with the highest mean Sub-Pallial neurons
signature
# Find the k-means clusters with the highest mean SP_signature1 score
MeanKclust.SPscore <- aggregate(SP_signature1 ~ kmeanClust, QCFiltered.data@meta.data, mean)
SPclust <- MeanKclust.SPscore %>% filter(SP_signature1 == max(SP_signature1)) %>% pull(kmeanClust)
# Extract barcodes and filter the seurat object
SP.cells <- QCFiltered.data@meta.data %>% filter(kmeanClust == SPclust) %>% pull(Barcodes)
SP.QCFiltered.data <- SubsetData(QCFiltered.data, cells.use = SP.cells , subset.raw = T, do.clean = F)
# Calculate SP early neurons signature score based on selected marker genes
LNgenes <- c("Dlx6os1", "Nudt4", "Abracl", "Arl4d", "Tmem123", "Ccdc109b", "Hmgn2", "E130006D01Rik", "Cdca7")
genes.list <- list(LNgenes)
enrich.name <- "EN_signature"
SP.QCFiltered.data <- AddModuleScore(SP.QCFiltered.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)
# K-means clustering based on EN and LN signature scores
cl <- kmeans(cbind(SP.QCFiltered.data@meta.data$EN_signature1, SP.QCFiltered.data@meta.data$LN_signature1), 3)
SP.QCFiltered.data@meta.data$kmeanClust <- paste0("Clust.",cl$cluster)
# Find the k-means clusters with the highest mean Late_signature1 score
MeanKclust.LNscore <- aggregate(LN_signature1 ~ kmeanClust, SP.QCFiltered.data@meta.data, mean)
LN.clust <- MeanKclust.LNscore %>% filter(LN_signature1 == max(LN_signature1)) %>% pull(kmeanClust)
# Extract barcodes
LN.cells <- SP.QCFiltered.data@meta.data %>% filter(kmeanClust == LN.clust) %>% pull(Barcodes)
col.pal <- wes_palette("GrandBudapest1", 3, type = "discrete")
p1 <- ggplot(SP.QCFiltered.data@meta.data, aes(x=LN_signature1, y=EN_signature1, colour = kmeanClust)) +
annotate(geom = "text", label = paste0(length(LN.cells), " Late Neurons"), x = 1, y =1.5) +
scale_color_manual(values=col.pal) +
geom_point() +
theme(legend.position="none")
ggMarginal(p1, type = "histogram", fill="lightgrey") ; rm(p1)
DimPlot(SP.QCFiltered.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 = T)
We perform the clustering procedure as for the glutamatergic neurons
We decide to exclude Cell cycle, ribosomal and mitochondrial associated genes, as well as Xist for the clustering step.
# Exclude cell cycle associated genes
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")
# Exclude genes detected in less than 3 cells
num.cells <- Matrix::rowSums(SP.QCFiltered.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), s.genes, g2m.genes, "Xist")
genes.use <- genes.use[!genes.use %in% GenesToRemove]
dgeMatrix_count <- as.matrix(SP.QCFiltered.data@raw.data)[rownames(SP.QCFiltered.data@raw.data) %in% genes.use,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1)
gene.counts <- log2(colSums(as.matrix(Data.matrix$norm.dat) > 0))
nUMI <- log2(colSums(Data.matrix$raw.dat))
perctMito <- SP.QCFiltered.data@meta.data$percent.mito
perctRibo <- SP.QCFiltered.data@meta.data$percent.ribo
rm.eigen <- as.matrix(cbind(gene.counts,
nUMI,
perctMito,
perctRibo))
row.names(rm.eigen) <- names(gene.counts)
colnames(rm.eigen) <- c("log2nGenes",
"log2nUMI",
"perctMito",
"perctRibo")
rm(gene.counts, nUMI, perctMito, perctRibo)
# 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 = 30,
min.cells = 10)
# 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,
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.002 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.008 s
## Build undirected graph from the weighted links...DONE ~ 0.013 s
## Run louvain clustering on the graph ...DONE ~ 0.006 s
## Return a community class
## -Modularity value: 0.7894448
## -Number of clusters: 11[1] "test-iter_clust.1"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.4405265
## -Number of clusters: 3[1] "test-iter_clust.2"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.002 s
## Build undirected graph from the weighted links...DONE ~ 0.003 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.6305993
## -Number of clusters: 4[1] "test-iter_clust.3"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0 s
## Return a community class
## -Modularity value: 0.4473667
## -Number of clusters: 2[1] "test-iter_clust.4"
## Finding nearest neighbors...DONE ~ 0.001 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.001 s
## Build undirected graph from the weighted links...DONE ~ 0.002 s
## Run louvain clustering on the graph ...DONE ~ 0.001 s
## Return a community class
## -Modularity value: 0.2015328
## -Number of clusters: 2[1] "test-iter_clust.5"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.003 s
## Build undirected graph from the weighted links...DONE ~ 0.007 s
## Run louvain clustering on the graph ...DONE ~ 0.003 s
## Return a community class
## -Modularity value: 0.6771333
## -Number of clusters: 7[1] "test-iter_clust.6"
## Finding nearest neighbors...DONE ~ 0 s
## Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 0.002 s
## Build undirected graph from the weighted links...DONE ~ 0.004 s
## Run louvain clustering on the graph ...DONE ~ 0.002 s
## Return a community class
## -Modularity value: 0.6456297
## -Number of clusters: 4
# 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")
## 6 Clusters
## 164 DE genes
SP.QCFiltered.data@ident <- as.factor(merge.result$cl)
SP.QCFiltered.data@meta.data$LN.ident <-as.character(SP.QCFiltered.data@ident)
# Load full dataset
Allcells.data <- readRDS("./Clustered.cells.RDS")
# Transfer the identities
Rename.Clust <- function(Clustdata, RawQCdata) {
for(i in unique(Clustdata@meta.data$LN.ident)){
New.ident <- i
Barcodes <- rownames(subset(Clustdata@meta.data, Clustdata@meta.data$LN.ident == 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("LN.GABA.",i))
}
return(RawQCdata)
}
Allcells.data <- Rename.Clust(Clustdata = SP.QCFiltered.data, RawQCdata = Allcells.data)
## [1] "Cluster_10: 113 Cells"
## [1] "Cluster_1: 26 Cells"
## [1] "Cluster_7: 45 Cells"
## [1] "Cluster_11: 54 Cells"
## [1] "Cluster_9: 16 Cells"
## [1] "Cluster_8: 23 Cells"
colors2 <- c("#969696",
"#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)
# Extract late neurons clusters
glut.clusters <- paste0(c("LN.Glut."),c("13", "1", "14", "22", "16", "19", "24", "26", "20", "21"))
gaba.clusters <- paste0( c("LN.GABA."),c("1", "7", "8", "9", "10", "11"))
clusters <- c(gaba.clusters, glut.clusters)
All.LN.data <- SubsetData(Allcells.data, ident.use = clusters, subset.raw = T, do.clean = F)
rm(glut.clusters, gaba.cluster, clusters)
colors <- c("#ec756d", "#c773a7", "#7293c8", "#b79f0b", "#3ca73f","#31b6bd",
"#ebcb2e", "#9ec22f", "#a9961b", "#cc3a1b", "#d14c8d", "#4cabdc", "#5ab793", "#e7823a", "#046c9a", "#4990c9")
DimPlot(All.LN.data,
reduction.use = "spring",
dim.1 = 1,
dim.2 = 2,
do.label=T,
label.size = 2,
no.legend = T,
cols.use = colors)
# Sub-pallial neurons
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.GABA.", c(8,11))), ident.use = "dLGE")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.GABA.", c(7,9,10))), ident.use = "Striatal_IN")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.GABA.", c(1))), ident.use = "Cortical_IN")
# Pallial neurons
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.Glut.", c(1,13,14))), ident.use = "CR")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.Glut.", c(22))), ident.use = "Fezf1")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.Glut.", c(16,19))), ident.use = "Pcp4")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.Glut.", c(24))), ident.use = "Nr4a2")
All.LN.data <- SetIdent(All.LN.data, cells.use = WhichCells(All.LN.data, ident = paste0("LN.Glut.", c(26,20,21))), ident.use = "Foxp2")
colors <- c("#ea6569", "#7694c8", "#e6a0c4", "#e2d203", "#096b9a", "#49adc7", "#e58606", "#b31021")
# Prepare annotation for hicat pipeline
colorsident <- cbind(ident = unique(as.character(All.LN.data@ident)),
colors = colors,
id = unique(as.character(All.LN.data@ident)))
# Create annotation data.frame
anno.df <- as.data.frame(cbind(
sample_name = row.names(All.LN.data@meta.data),
primary_type_id = colorsident[match(as.character(All.LN.data@ident), colorsident[,1]),3],
primary_type_label = as.character(All.LN.data@ident),
primary_type_color = colorsident[match(as.character(All.LN.data@ident), colorsident[,1]),2]
))
# Make a data.frame of unique cluster id, type, color, and broad type
cl.df <- anno.df %>%
select(primary_type_id, primary_type_label, primary_type_color) %>%
unique()
colnames(cl.df)[1:3] <- c("cluster_id", "cluster_label", "cluster_color")
# Sort by cluster_id
cl.df <- arrange(cl.df, cluster_id)
row.names(cl.df) <- cl.df$cluster_id
cl.fact <- setNames(factor(anno.df$primary_type_id), anno.df$sample_name)
# Filter genes
num.cells <- Matrix::rowSums(All.LN.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")
genes.use <- genes.use[!genes.use %in% GenesToRemove] ; rm(GenesToRemove, num.cells)
All.LN.data@raw.data <- All.LN.data@raw.data[genes.use, ]
All.LN.data <- NormalizeData(object = All.LN.data,
normalization.method = "LogNormalize",
scale.factor = round(median(All.LN.data@meta.data$nUMI)),
display.progress = F)
# Find all var genes
All.LN.data <- FindVariableGenes(object = All.LN.data,
mean.function = ExpMean,
dispersion.function = LogVMR,
x.low.cutoff = 0.02,
x.high.cutoff = 3,
y.cutoff = 1,
do.plot = F, display.progress = F)
dgeMatrix_count <- as.matrix(All.LN.data@raw.data)[rownames(All.LN.data@raw.data) %in% All.LN.data@var.genes,]
dgeMatrix_cpm <- cpm(dgeMatrix_count)
norm.dat <- log2(dgeMatrix_cpm + 1) ; rm(dgeMatrix_cpm)
Data.matrix <- list(raw.dat=dgeMatrix_count, norm.dat=norm.dat) ; attach(Data.matrix)
rm(dgeMatrix_count,norm.dat)
Scrattch.hicat perform hierarchical clustering on a cluster correlation matrix based on median expression values for the top 50 most DEGs between every pair of clusters. It estimates branch confidence level using a bootstrap approach implemented by the pvclust package
# Build the dendrogram
dend.result <- build_dend(cl.med[,levels(cl.fact)],
l.color= setNames(as.character(cl.df$cluster_color), row.names(cl.df)),
nboot = 100)
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
# Attach cluster labels to the leaves of the tree
dend.labeled <- dend.result$dend
labels(dend.labeled) <- cl.df[labels(dend.labeled), "cluster_label"]
# Rotate dendrogramme leafs
NewOrder <- c("Cortical_IN","Striatal_IN","dLGE","CR","Fezf1","Pcp4","Nr4a2","Foxp2")
l.rank <- setNames(1:nrow(cl.df), NewOrder) #set the cluster order
# Color of the leaf nodes.
l.color <- setNames(as.character(cl.df$cluster_color), NewOrder)
# Re-build the dendrogram
cl.med <- get_cl_medians(norm.dat, cl.fact)
dend.result <- build_dend(cl.med[,levels(cl.fact)],
l.rank,
l.color=l.color,
nboot = 100)
## Bootstrap (r = 0.5)... Done.
## Bootstrap (r = 0.6)... Done.
## Bootstrap (r = 0.7)... Done.
## Bootstrap (r = 0.8)... Done.
## Bootstrap (r = 0.9)... Done.
## Bootstrap (r = 1.0)... Done.
## Bootstrap (r = 1.1)... Done.
## Bootstrap (r = 1.2)... Done.
## Bootstrap (r = 1.3)... Done.
## Bootstrap (r = 1.4)... Done.
# Import dendrogram into Seurat object
All.LN.data@cluster.tree[[1]] <- ape::as.phylo(dend.labeled)
labels(All.LN.data@cluster.tree[[1]]) <- cl.df[as.numeric(labels(All.LN.data@cluster.tree[[1]])), "cluster_label"]
All.LN.data@ident <- factor(All.LN.data@ident, levels = c("Cortical_IN","Striatal_IN","dLGE","CR","Fezf1","Pcp4","Nr4a2","Foxp2"))
p1 <- ggdendrogram(dendro_data(as.hclust(dend.labeled)), labels = F, rotate = T) + scale_y_reverse()
p2 <- Cluster.dotplot(All.LN.data,
Marker.genes = rev(c("Foxg1","Tbr1", "Slc17a6",
"Reln", "Lhx5", "Neurod6",
"Pcp4", "Fezf1", "Pbx3",
"Foxp2", "Nr4a2", "Dlx5", "Gad2",
"Six3", "Sp8", "Zfp503", "Isl1", "Lhx6", "Sst")),
min.expression = 0.7, percent.mi=0.15, maxdot.size = 5)
plot_grid(plotlist = list(p1,p2), ncol=2, align='h', rel_widths = c(0.2, 1.5))
## [1] "30 novembre, 2020, 10,34"
## 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] pvclust_2.2-0 limma_3.42.0 Rphenograph_0.99.1
## [4] igraph_1.2.5 wesanderson_0.3.6 ggExtra_0.9
## [7] tidyr_1.0.0 ggdendro_0.1-20 scrattch.vis_0.0.210
## [10] purrr_0.3.3 ggbeeswarm_0.6.0 RColorBrewer_1.1-2
## [13] Seurat_2.3.4 cowplot_1.0.0 ggplot2_3.2.1
## [16] Matrix_1.2-17 matrixStats_0.55.0 dplyr_0.8.3
## [19] dendextend_1.12.0 scrattch.hicat_0.0.16
##
## 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 splines_3.6.3
## [7] digest_0.6.25 foreach_1.4.7 htmltools_0.5.0
## [10] viridis_0.5.1 lars_1.2 gdata_2.18.0
## [13] magrittr_1.5 checkmate_1.9.4 cluster_2.1.0
## [16] mixtools_1.1.0 ROCR_1.0-7 R.utils_2.9.0
## [19] colorspace_1.4-1 xfun_0.18 crayon_1.3.4
## [22] jsonlite_1.7.0 zeallot_0.1.0 survival_2.44-1.1
## [25] zoo_1.8-6 iterators_1.0.12 ape_5.3
## [28] glue_1.4.1 gtable_0.3.0 kernlab_0.9-29
## [31] prabclus_2.3-1 DEoptimR_1.0-8 scales_1.1.0
## [34] bibtex_0.4.2 miniUI_0.1.1.1 Rcpp_1.0.5
## [37] metap_1.1 dtw_1.21-3 viridisLite_0.3.0
## [40] xtable_1.8-4 htmlTable_1.13.2 reticulate_1.13
## [43] foreign_0.8-72 bit_4.0.4 proxy_0.4-23
## [46] mclust_5.4.5 SDMTools_1.1-221.1 Formula_1.2-3
## [49] stats4_3.6.3 tsne_0.1-3 htmlwidgets_1.5.1
## [52] httr_1.4.1 gplots_3.0.1.1 fpc_2.2-3
## [55] ellipsis_0.3.0 acepack_1.4.1 modeltools_0.2-22
## [58] ica_1.0-2 farver_2.0.1 pkgconfig_2.0.3
## [61] R.methodsS3_1.7.1 flexmix_2.3-15 nnet_7.3-14
## [64] tidyselect_0.2.5 labeling_0.3 rlang_0.4.7
## [67] reshape2_1.4.3 later_1.0.0 munsell_0.5.0
## [70] tools_3.6.3 ggridges_0.5.1 evaluate_0.14
## [73] stringr_1.4.0 fastmap_1.0.1 yaml_2.2.1
## [76] npsurv_0.4-0 knitr_1.26 bit64_4.0.2
## [79] fitdistrplus_1.0-14 robustbase_0.93-5 caTools_1.17.1.2
## [82] RANN_2.6.1 pbapply_1.4-2 nlme_3.1-141
## [85] mime_0.7 R.oo_1.23.0 hdf5r_1.3.2.9000
## [88] compiler_3.6.3 rstudioapi_0.11 beeswarm_0.2.3
## [91] png_0.1-7 lsei_1.2-0 tibble_2.1.3
## [94] stringi_1.4.6 highr_0.8 lattice_0.20-41
## [97] vctrs_0.2.0 pillar_1.4.2 lifecycle_0.1.0
## [100] Rdpack_0.11-0 lmtest_0.9-37 data.table_1.12.6
## [103] bitops_1.0-6 irlba_2.3.3 gbRd_0.4-11
## [106] httpuv_1.5.2 R6_2.4.1 latticeExtra_0.6-28
## [109] promises_1.1.0 KernSmooth_2.23-15 gridExtra_2.3
## [112] vipor_0.4.5 codetools_0.2-16 MASS_7.3-53
## [115] gtools_3.8.1 assertthat_0.2.1 withr_2.1.2
## [118] diptest_0.75-7 parallel_3.6.3 doSNOW_1.0.18
## [121] grid_3.6.3 rpart_4.1-15 class_7.3-17
## [124] rmarkdown_2.5 segmented_1.0-0 Rtsne_0.15
## [127] shiny_1.4.0 base64enc_0.1-3
Institute of Psychiatry and Neuroscience of Paris, INSERM U1266, 75014, Paris, France, matthieu.moreau@inserm.fr↩