Loo et all markers

library(COTAN)
library(ggplot2)
library(data.table)
library(Matrix)
library(factoextra)
library(ggrepel)
library(dendextend)

1 Gene clustering using COTAN

input_dir = "Data/"
layers = list("L1"=c("Reln","Lhx5"), "L2/3"=c("Satb2","Cux1"), "L4"=c("Rorb","Sox5") , "L5/6"=c("Bcl11b","Fezf2") , "Prog"=c("Vim","Hes1"))
objE17 = readRDS(file = "Data/E17_cortex_cl2.cotan.RDS")

g.space = get.gene.coexpression.space(objE17, 
                                      n.genes.for.marker = 25,  primary.markers = unlist(layers))

[1] "calculating gene coexpression space: output tanh of reduced coex matrix"
     L11      L12    L2/31    L2/32      L41      L42    L5/61    L5/62    Prog1    Prog2 
  "Reln"   "Lhx5"  "Satb2"   "Cux1"   "Rorb"   "Sox5" "Bcl11b"  "Fezf2"    "Vim"   "Hes1" 
[1] "Get p-values on a set of genes on columns genome wide on rows"
[1] "Using function S"
[1] "function to generate S "
[1] "Secondary markers:181"
[1] "function to generate S "
[1] "Columns (V set) number: 181 Rows (U set) number: 1236"

g.space = as.data.frame(as.matrix(g.space))

coex.pca.genes <- prcomp(t(g.space),
                 center = TRUE,
                 scale. = F) 

fviz_eig(coex.pca.genes, addlabels=TRUE,ncp = 10)

#fviz_eig(coex.pca.genes, choice = "eigenvalue", addlabels=TRUE)

1.1 Hierarchical clustering

controls =list("genes related to L5/6"=c("Foxp2","Tbr1"), "genes related to L2/3"=c("Mef2c"), "genes related to Prog"=c("Nes","Sox2") , "genes related to L1"=c() , "genes related to L4"=c())
# clustering usign Ward method
hc.norm = hclust(dist(g.space), method = "ward.D2")

# and cut the tree into 5 clusters (for example)
cut = cutree(hc.norm, k = 7, order_clusters_as_data = F)

# It crates the tree
dend <- as.dendrogram(hc.norm)

# I can use a dataframe from the pca to store some data regarding the clustering
pca_1 = as.data.frame(coex.pca.genes$rotation[,1:5])
pca_1 = pca_1[order.dendrogram(dend),]


mycolours <- c("genes related to L5/6" = "#3C5488FF","genes related to L2/3"="#F39B7FFF","genes related to Prog"="#4DBBD5FF","genes related to L1"="#E64B35FF","genes related to L4" = "#91D1C2FF", "not marked"="#B09C85FF")

# save the cluster number into the dataframe
pca_1$hclust = cut

dend =branches_color(dend,k=7,col=c("#4DBBD5FF","#91D1C2FF","#F39B7FFF","#E64B35FF","#3C5488FF","#91D1C2FF","#B09C85FF" ),groupLabels = T)
#dend =color_labels(dend,k=5,labels = rownames(pca_1),col=pca_1$colors)



dend %>%
set("labels", ifelse(labels(dend) %in% rownames(pca_1)[rownames(pca_1) %in% c(unlist(layers),unlist(controls))], labels(dend), "")) %>%
  set("branches_k_color", value = c("gray80","#4DBBD5FF","#91D1C2FF" ,"gray80","#F39B7FFF","#E64B35FF","#3C5488FF"), k = 7) %>%
plot(horiz=F, axes=T,ylim = c(0,100))

1.2 Import of markers from Loo et al. paper

Markers_Loo = read.csv("Data/Markers_Loo.csv")

genes_detected = pca_1[rownames(pca_1) %in% unlist(Markers_Loo)[!unlist(Markers_Loo) == ""],]

Gene present in the 10% of most differentially expressed genes by COTAN

dim(genes_detected)[1]

[1] 33

of total genes detected as markers by Loo et al

length(unlist(Markers_Loo)[!unlist(Markers_Loo) == ""])

[1] 48

Number of genes detected:

sum(unlist(Markers_Loo)[!unlist(Markers_Loo) == ""] %in% rownames(objE17@coex))

[1] 43

Removed becouse not detected

pure.markers = unlist(Markers_Loo)
pure.markers = pure.markers[!unlist(Markers_Loo) == ""]
pure.markers[!pure.markers %in% rownames(objE17@coex)]

    L.I2    PROG4   PROG10   PROG19   PROG22 
  "Gdf5" "Cdc25c" "Gas2l3"  "Rspo1"  "Wnt8b" 

Markers_Loo = as.list(Markers_Loo)
for (i in names(Markers_Loo)) {
    Markers_Loo[[i]] = Markers_Loo[[i]][! Markers_Loo[[i]] == ""]
}
Markers_Loo

$L.I
[1] "Ebf3"  "Gdf5"  "Lhx1"  "Lhx5"  "Ndnf"  "Reln"  "Samd3" "Trp73"

$L.II.IV
[1] "Satb2"         "3110047P20Rik" "9130024F11Rik" "Dok5"          "Inhba"         "Pou3f1"       

$L.V.VI
 [1] "Bcl11b" "Crym"   "Fezf2"  "Hs3st4" "Mc4r"   "Nfe2l3" "Nxph3"  "Plxna4" "Rwdd3"  "Sla"    "Sybu"   "Tbr1"  

$PROG
 [1] "Aldoc"    "Arhgef39" "Aspm"     "Cdc25c"   "Cdkn3"    "Cyr61"    "Dkk3"     "Ednrb"    "Gas1"    
[10] "Gas2l3"   "Hes1"     "Hes5"     "Htra1"    "Nde1"     "Nek2"     "Pax6"     "Pkmyt1"   "Plk1"    
[19] "Rspo1"    "Tcf19"    "Tk1"      "Wnt8b"   

Primary markers also used by Loo et al.

unlist(layers)[unlist(layers) %in% unique(unlist(Markers_Loo))]

     L11      L12    L2/31    L5/61    L5/62    Prog2 
  "Reln"   "Lhx5"  "Satb2" "Bcl11b"  "Fezf2"   "Hes1" 

Table

tableMarkersCOTAN = as.data.frame(matrix(nrow = 6,ncol = 4))
colnames(tableMarkersCOTAN)=c("Loo.L.I","Loo.L.II.IV","Loo.L.V.VI","Loo.PROG")
rownames(tableMarkersCOTAN)=c("COTAN.L.I","COTAN.L.II.III","COTAN.L.IV","COTAN.L.V.VI","COTAN.PROG","COTAN.Not Grouped")
tableMarkersCOTAN_no_sec = tableMarkersCOTAN

layers.cluster = c("Reln","Satb2","Rorb","Bcl11b","Vim")
names(layers.cluster) = unique(cut[unlist(layers)])
layers.cluster

       3        6        2        5        1 
  "Reln"  "Satb2"   "Rorb" "Bcl11b"    "Vim" 

1.3 Comparition between Loo et al. markers and COTAN markers

groups = list("L.I"=names(which(layers.cluster == "Reln"))
              ,"L.II.III"=names(which(layers.cluster == "Satb2")),
              "L.IV"=names(which(layers.cluster == "Rorb")),
              "L.V.VI"=names(which(layers.cluster == "Bcl11b")),
              "PROG"=names(which(layers.cluster == "Vim")))

for(layer1 in c("L.I","L.II.III", "L.IV","L.V.VI","PROG")){
    for(layer2 in c("L.I","L.II.IV","L.V.VI","PROG")){
    tableMarkersCOTAN[paste0("COTAN.",layer1),paste0("Loo.",layer2)] =
        sum(rownames(pca_1[pca_1$hclust %in% groups[[layer1]],]) %in% Markers_Loo[[layer2]])
    tableMarkersCOTAN[paste0("COTAN.","Not Grouped"),paste0("Loo.",layer2)] =
        sum(rownames(pca_1[!pca_1$hclust %in% unlist(groups),]) %in% Markers_Loo[[layer2]])
    }
}

tableMarkersCOTAN

Without secondary markers

pca_2 = pca_1[! rownames(pca_1) %in% colnames(g.space),]


for(layer1 in c("L.I","L.II.III", "L.IV","L.V.VI","PROG")){
    for(layer2 in c("L.I","L.II.IV","L.V.VI","PROG")){
    tableMarkersCOTAN_no_sec[paste0("COTAN.",layer1),paste0("Loo.",layer2)] =
        sum(rownames(pca_2[pca_2$hclust %in% groups[[layer1]],]) %in% Markers_Loo[[layer2]])
    tableMarkersCOTAN_no_sec[paste0("COTAN.","Not Grouped"),paste0("Loo.",layer2)] =
        sum(rownames(pca_2[!pca_2$hclust %in% unlist(groups),]) %in% Markers_Loo[[layer2]])
    }
}

tableMarkersCOTAN_no_sec

specific.genes.table = data.frame("genes"=c(), "COTAN"=c(),"Loo."=c())
tt1 = c()
tt2 = c()
for(layer1 in c("L.I","L.II.III", "L.IV","L.V.VI","PROG")){
    for(layer2 in c("L.I","L.II.IV","L.V.VI","PROG")){
    tt1 = data.frame("genes"= rownames(pca_1[pca_1$hclust %in% groups[[layer1]],])[rownames(pca_1[pca_1$hclust 
                          %in% groups[[layer1]],]) %in% Markers_Loo[[layer2]]])
    if (dim(tt1)[1] > 0) {
                   tt1 = cbind(tt1,  "COTAN"=layer1, "Loo."=layer2)
        }
    tt2 = data.frame("genes"= 
        rownames(pca_1[!pca_1$hclust %in% unlist(groups),])[rownames(pca_1[!pca_1$hclust %in%                   unlist(groups),]) %in% Markers_Loo[[layer2]]])
    if (dim(tt2)[1] > 0) {
        tt2 = cbind(tt2, "COTAN"= "Not Grouped", "Loo."=layer2)
    }
    specific.genes.table = rbind(specific.genes.table,tt1,tt2)
    }
}

specific.genes.table[!(duplicated(specific.genes.table)) , ]

NA

dend %>%
set("labels", ifelse(labels(dend) %in% rownames(pca_1)[rownames(pca_1) %in% specific.genes.table$genes], labels(dend), "")) %>%
  set("branches_k_color", value = c("gray80","#4DBBD5FF","#91D1C2FF" ,"gray80","#F39B7FFF","#E64B35FF","#3C5488FF"), k = 7) %>%
    plot(horiz=T, axes=T,xlim = c(0,80),cex = 10, dLeaf = -7)

2 WGCNA

library(WGCNA)
library(cluster)
library(data.table)
library(Matrix)
library(Seurat)
library(utils)
library(dplyr)
library(patchwork)
library(graphics)
options(stringsAsFactors = FALSE)
data_dir = "Data/"
library(gplots)
myheatcol = colorpanel(250,'red',"orange",'lemonchiffon')

2.1 Test with the 2000 most varied genes

This seems the best option for the analysis with WGCNA.

data = as.data.frame(fread(paste(input_dir,"E175_only_cortical_cells.txt.gz", sep = "/"),sep = "\t"))
data = as.data.frame(data)
rownames(data) = data$V1
data = data[,2:ncol(data)]
data[1:10,1:10]

E17 <- CreateSeuratObject(counts = data, project = "Cortex E17.5", min.cells = 3, min.features = 200)

Feature names cannot have underscores ('_'), replacing with dashes ('-')

E17[["percent.mt"]] <- PercentageFeatureSet(E17, pattern = "^mt-")
VlnPlot(E17, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

hist(E17$nFeature_RNA/E17$nCount_RNA, breaks = 100)

plot1 <- FeatureScatter(E17, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(E17, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
plot1 + plot2

E17 <- NormalizeData(E17, normalization.method = "LogNormalize", scale.factor = 10000)

Performing log-normalization
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

E17 <- FindVariableFeatures(E17, selection.method = "vst", nfeatures = 2000)

Calculating gene variances
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|
Calculating feature variances of standardized and clipped values
0%   10   20   30   40   50   60   70   80   90   100%
[----|----|----|----|----|----|----|----|----|----|
**************************************************|

# Identify the 10 most highly variable genes
top10 <- head(VariableFeatures(E17), 10)

# plot variable features with and without labels
plot1 <- VariableFeaturePlot(E17)
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)

When using repel, set xnudge and ynudge to 0 for optimal results

plot2

all.genes <- rownames(E17)
E17 <- ScaleData(E17, features = all.genes)

Centering and scaling data matrix

  |                                                                                                            
  |                                                                                                      |   0%
  |                                                                                                            
  |========                                                                                              |   8%
  |                                                                                                            
  |================                                                                                      |  15%
  |                                                                                                            
  |========================                                                                              |  23%
  |                                                                                                            
  |===============================                                                                       |  31%
  |                                                                                                            
  |=======================================                                                               |  38%
  |                                                                                                            
  |===============================================                                                       |  46%
  |                                                                                                            
  |=======================================================                                               |  54%
  |                                                                                                            
  |===============================================================                                       |  62%
  |                                                                                                            
  |=======================================================================                               |  69%
  |                                                                                                            
  |==============================================================================                        |  77%
  |                                                                                                            
  |======================================================================================                |  85%
  |                                                                                                            
  |==============================================================================================        |  92%
  |                                                                                                            
  |======================================================================================================| 100%

E17 <- RunPCA(E17, features = VariableFeatures(object = E17))

PC_ 1 
Positive:  Fabp7, Aldoc, Mfge8, Dbi, Ednrb, Vim, Slc1a3, Mt3, Apoe, Ttyh1 
       Tnc, Sox2, Atp1a2, Ddah1, Hes5, Sparc, Mlc1, Ppap2b, Rgcc, Bcan 
       Ndrg2, Qk, Lxn, Id3, Phgdh, Slc9a3r1, Nr2e1, Aldh1l1, Gpx8, Mt1 
Negative:  Tubb3, Stmn2, Neurod6, Stmn4, Map1b, Stmn1, Myt1l, Mef2c, Thra, 4930506M07Rik 
       Bcl11a, Gap43, Bhlhe22, Syt4, Cntn2, Nell2, Hs6st2, 9130024F11Rik, Olfm1, Satb2 
       Akap9, Ptprd, Rbfox1, Clmp, Ina, Enc1, Camk2b, Dync1i1, Dab1, Atp2b1 
PC_ 2 
Positive:  Sstr2, Mdk, Meis2, Pou3f2, Eomes, Zbtb20, Unc5d, Sema3c, Fos, Tead2 
       Palmd, Mfap4, Nhlh1, Ulk4, H1f0, Uaca, Neurog2, Neurod1, Ezr, Ier2 
       Nrn1, Baz2b, Pdzrn3, Btg2, Egr1, Mfap2, Loxl1, H2afv, Hbp1, Nnat 
Negative:  Gap43, Sybu, Dync1i1, Meg3, Mef2c, Map1b, Fezf2, Camk2b, Ina, Stmn2 
       Cdh13, Thra, Nin, Rac3, Igfbp3, Ssbp2, Neto2, Cd200, Hmgcs1, Tuba1b 
       Syt1, Slc6a15, Mapre2, Plk2, Rprm, Atp1b1, Cadm2, Arpp21, Kitl, Ntrk2 
PC_ 3 
Positive:  Meg3, Smpdl3a, Slc9a3r1, Slc15a2, Timp3, Tmem47, Ndrg2, Apoe, Ttyh1, Fmo1 
       Mlc1, Scrg1, Islr2, Malat1, Gstm1, Gja1, Ndnf, Aldh1l1, Mt3, Sparc 
       Serpinh1, Paqr7, Asrgl1, Sepp1, S100a1, Atp1b1, Ctsl, Cpe, S100a16, Lhx5 
Negative:  Birc5, Top2a, Cenpm, Pbk, Tpx2, Cenpe, Mki67, Cdca8, Gmnn, Cks2 
       Ccnb1, Ccnb2, Spc24, Hmgb2, Cenpf, Tk1, Hmmr, Prc1, Kif11, Ccna2 
       2810417H13Rik, C330027C09Rik, Cdca2, Ect2, Nusap1, Cenpa, Uhrf1, Plk1, Spc25, Knstrn 
PC_ 4 
Positive:  Lhx5, Nhlh2, Snhg11, Reln, 1500016L03Rik, Trp73, Cacna2d2, Ndnf, Car10, Lhx1 
       Islr2, Pcp4, Meg3, RP24-351J24.2, Rcan2, Pnoc, Mab21l1, Zic1, E330013P04Rik, Emx2 
       Malat1, Ebf3, Nr2f2, Zcchc12, Zbtb20, Celf4, Tmem163, Ache, Calb2, Unc5b 
Negative:  Ptn, Satb2, 9130024F11Rik, Neurod6, Mef2c, Dab1, Limch1, Hs6st2, Abracl, Dok5 
       Gucy1a3, Nell2, Ptprz1, Syt4, Ttc28, Clmp, Macrod2, Fam19a2, Smpdl3a, Ndrg1 
       Gstm1, 4930506M07Rik, Paqr7, Aldh1l1, Myt1l, Hmgcs1, Slc15a2, Pdzrn4, Slc9a3r1, Aldoc 
PC_ 5 
Positive:  Fam210b, Sfrp1, Pax6, Enkur, Tubb3, Tuba1b, Mcm3, Veph1, Stmn1, Eif1b 
       Map1b, Hopx, Abracl, Cdk2ap2, Tfap2c, Rps27l, 2810025M15Rik, Slc14a2, Prdx1, Hells 
       Gap43, Sept11, Egln3, Gm1840, Ezr, Cpne2, 9130024F11Rik, Nes, Efnb2, Cux1 
Negative:  Serpine2, Id1, Olig1, Sparcl1, Igfbp3, Fam212b, Ccnb2, Ppic, Gng12, Ccnb1 
       Bcan, Cenpe, Pbk, Id3, Rasl11a, Plk1, Aqp4, Aspm, Hmmr, Slc6a1 
       Slc4a4, Malat1, Myo6, Timp3, Meg3, Cdk1, Prrx1, Npy, B2m, Cspg4 

DimPlot(E17, reduction = "pca")

all.genes <- rownames(E17)
E17 <- ScaleData(E17, features = all.genes)

Centering and scaling data matrix

  |                                                                                                            
  |                                                                                                      |   0%
  |                                                                                                            
  |========                                                                                              |   8%
  |                                                                                                            
  |================                                                                                      |  15%
  |                                                                                                            
  |========================                                                                              |  23%
  |                                                                                                            
  |===============================                                                                       |  31%
  |                                                                                                            
  |=======================================                                                               |  38%
  |                                                                                                            
  |===============================================                                                       |  46%
  |                                                                                                            
  |=======================================================                                               |  54%
  |                                                                                                            
  |===============================================================                                       |  62%
  |                                                                                                            
  |=======================================================================                               |  69%
  |                                                                                                            
  |==============================================================================                        |  77%
  |                                                                                                            
  |======================================================================================                |  85%
  |                                                                                                            
  |==============================================================================================        |  92%
  |                                                                                                            
  |======================================================================================================| 100%

seurat.data = as.matrix(E17[["RNA"]]@data)

Var.genes = VariableFeatures(E17)[1:2000]

datExpr0 = t(seurat.data[rownames(seurat.data) %in% Var.genes,])
gsg = goodSamplesGenes(datExpr0, verbose = 3)

 Flagging genes and samples with too many missing values...
  ..step 1

gsg$allOK

[1] TRUE

if (!gsg$allOK){
    # Optionally, print the gene and sample names that were removed:
    if (sum(!gsg$goodGenes)>0)
        printFlush(paste("Removing genes:", paste(names(datExpr0)[!gsg$goodGenes], collapse = ", ")));
    if (sum(!gsg$goodSamples)>0)
        printFlush(paste("Removing samples:", paste(rownames(datExpr0)[!gsg$goodSamples], collapse = ", ")));
    # Remove the offending genes and samples from the data:
    datExpr0 = datExpr0[gsg$goodSamples, gsg$goodGenes]
}



sampleTree = hclust(dist(datExpr0), method = "average");
# Plot the sample tree: Open a graphic output window of size 12 by 9 inches
# The user should change the dimensions if the window is too large or too small.
sizeGrWindow(12,9)

par(cex = 0.6);
par(mar = c(0,4,2,0))
plot(sampleTree, main = "Sample clustering to detect outliers", sub="", xlab="", cex.lab = 1.5,
     cex.axis = 1.5, cex.main = 2)

# Plot a line to show the cut
#abline(h = 400, col = "red")

No outliner detected

# Automatic network construction and module detection
# Choose a set of soft-thresholding powers
powers = c(c(1:10), seq(from = 10, to=25, by=2))
# Call the network topology analysis function
sft = pickSoftThreshold(datExpr0, powerVector = powers, verbose = 5)

pickSoftThreshold: will use block size 2000.
 pickSoftThreshold: calculating connectivity for given powers...
   ..working on genes 1 through 2000 of 2000

Plot the results:

sizeGrWindow(9, 5)
par(mfrow = c(1,2));
cex1 = 0.9;
# Scale-free topology fit index as a function of the soft-thresholding power
plot(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     xlab="Soft Threshold (power)",ylab="Scale Free Topology Model Fit,signed R^2",type="n",
     main = paste("Scale independence"));
text(sft$fitIndices[,1], -sign(sft$fitIndices[,3])*sft$fitIndices[,2],
     labels=powers,cex=cex1,col="red");

# this line corresponds to using an R^2 cut-off of h
abline(h=0.90,col="red")
# Mean connectivity as a function of the soft-thresholding power
plot(sft$fitIndices[,1], sft$fitIndices[,5],
     xlab="Soft Threshold (power)",ylab="Mean Connectivity", type="n",
     main = paste("Mean connectivity"))

text(sft$fitIndices[,1], sft$fitIndices[,5], labels=powers, cex=cex1,col="red")

Tested with 5, 3 and 2 and 4. The best seems 2

net = blockwiseModules(datExpr0, power = 2, maxBlockSize = 20000,
                       TOMType = "signed", minModuleSize = 30,
                       reassignThreshold = 0, mergeCutHeight = 0.25,
                       numericLabels = TRUE, pamRespectsDendro = FALSE,
                       saveTOMs = TRUE,
                       saveTOMFileBase = "E17.5",
                       verbose = 3)

 Calculating module eigengenes block-wise from all genes
   Flagging genes and samples with too many missing values...
    ..step 1
 ..Working on block 1 .
    TOM calculation: adjacency..
    ..will not use multithreading.
     Fraction of slow calculations: 0.000000
    ..connectivity..
    ..matrix multiplication (system BLAS)..
    ..normalization..
    ..done.
   ..saving TOM for block 1 into file E17.5-block.1.RData
 ....clustering..
 ....detecting modules..
 ....calculating module eigengenes..
 ....checking kME in modules..
     ..removing 552 genes from module 1 because their KME is too low.
     ..removing 154 genes from module 2 because their KME is too low.
     ..removing 86 genes from module 3 because their KME is too low.
     ..removing 42 genes from module 4 because their KME is too low.
     ..removing 20 genes from module 5 because their KME is too low.
 ..merging modules that are too close..
     mergeCloseModules: Merging modules whose distance is less than 0.25
       Calculating new MEs...

# open a graphics window
sizeGrWindow(12, 9)
# Convert labels to colors for plotting
mergedColors = labels2colors(net$colors)
# Plot the dendrogram and the module colors underneath
plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)

2.2 Comparition between Loo et al. markers and WGCNA markers

plotNetworkHeatmap(datExpr0, plotGenes = unlist(Markers_Loo),
networkType="signed", useTOM=TRUE,
power=2, main="D. TOM in an signed network")

Warning: Not all gene names were recognized. Only the following genes were recognized. 
    Ebf3,     Lhx1,     Lhx5,     Ndnf,     Reln,     Samd3,     Trp73,     Satb2,     9130024F11Rik,     Dok5,     Inhba,     Bcl11b,     Crym,     Fezf2,     Hs3st4,     Nfe2l3,     Nxph3,     Rwdd3,     Sla,     Sybu,     Aldoc,     Arhgef39,     Aspm,     Cdc25c,     Cyr61,     Ednrb,     Gas1,     Hes1,     Hes5,     Htra1,     Nek2,     Pax6,     Plk1,     Tcf19,     Tk1
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

plotNetworkHeatmap(datExpr0, plotGenes = unlist(Markers_Loo),

There were 15 warnings (use warnings() to see them)

networkType="unsigned", useTOM=TRUE,
power=2, main="D. TOM in an signed network")

Warning: Not all gene names were recognized. Only the following genes were recognized. 
    Ebf3,     Lhx1,     Lhx5,     Ndnf,     Reln,     Samd3,     Trp73,     Satb2,     9130024F11Rik,     Dok5,     Inhba,     Bcl11b,     Crym,     Fezf2,     Hs3st4,     Nfe2l3,     Nxph3,     Rwdd3,     Sla,     Sybu,     Aldoc,     Arhgef39,     Aspm,     Cdc25c,     Cyr61,     Ednrb,     Gas1,     Hes1,     Hes5,     Htra1,     Nek2,     Pax6,     Plk1,     Tcf19,     Tk1
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

plotNetworkHeatmap(datExpr0, plotGenes = unique(unlist(layers)),
networkType="signed", useTOM=TRUE,
power=2, main="D. TOM in an signed network")

..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

plotNetworkHeatmap(datExpr0, plotGenes = unique(unlist(layers)),
networkType="unsigned", useTOM=TRUE,
power=2, main="D. TOM in an unsigned network")

..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

# Calculate topological overlap anew: this could be done more efficiently by saving the TOM
# calculated during module detection, but let us do it again here.
dissTOM = 1-TOMsimilarityFromExpr(datExpr0, power = 2);

TOM calculation: adjacency..
..will not use multithreading.
 Fraction of slow calculations: 0.000000
..connectivity..
..matrix multiplication (system BLAS)..
..normalization..
..done.

# Transform dissTOM with a power to make moderately strong connections more visible in the heatmap
plotTOM = dissTOM^7;
# Set diagonal to NA for a nicer plot
diag(plotTOM) = NA;
rownames(dissTOM)=colnames(datExpr0)
colnames(dissTOM)=colnames(datExpr0)

selectTOM = dissTOM[rownames(dissTOM) %in% unlist(Markers_Loo),colnames(dissTOM) %in% unlist(Markers_Loo)];
# There’s no simple way of restricting a clustering tree to a subset of genes, so we must re-cluster.
selectTree = hclust(as.dist(selectTOM), method = "average")
moduleColors = mergedColors
names(moduleColors) = rownames(dissTOM)
selectColors = moduleColors[names(moduleColors) %in% unlist(Markers_Loo)]
# Open a graphical window
sizeGrWindow(9,9)
# Taking the dissimilarity to a power, say 10, makes the plot more informative by effectively changing
# the color palette; setting the diagonal to NA also improves the clarity of the plot
plotDiss = selectTOM^7;
diag(plotDiss) = NA;
TOMplot(plotDiss, selectTree, selectColors, main = "Network heatmap plot, selected genes",col=myheatcol)

plotDendroAndColors(selectTree,selectColors,
                     "Module colors",
                     dendroLabels = NULL, hang = 0.03,
                     addGuide = TRUE, guideHang = 0.05)

net$colors[names(net$colors) %in% unlist(Markers_Loo)]

9130024F11Rik         Aldoc      Arhgef39          Aspm        Bcl11b        Cdc25c          Crym         Cyr61 
            0             1             0             2             0             2             0             1 
         Dok5          Ebf3         Ednrb         Fezf2          Gas1          Hes1          Hes5        Hs3st4 
            0             4             1             0             1             1             1             0 
        Htra1         Inhba          Lhx1          Lhx5          Ndnf          Nek2        Nfe2l3         Nxph3 
            1             0             4             4             4             2             0             0 
         Pax6          Plk1          Reln         Rwdd3         Samd3         Satb2           Sla          Sybu 
            1             2             4             0             0             0             0             0 
        Tcf19           Tk1         Trp73 
            0             2             4 

There were 15 warnings (use warnings() to see them)

WGCNA using the 2000 genes most varied by Seurat, detects the following number of markers.

length(net$colors[names(net$colors) %in% unlist(Markers_Loo)])

[1] 35

tableMarkersWGCNA = as.data.frame(matrix(nrow = 5,ncol = 4))
colnames(tableMarkersWGCNA)=c("Loo.L.I","Loo.L.II.IV","Loo.L.V.VI","Loo.PROG")
rownames(tableMarkersWGCNA)=c("WGCNA.L.I","WGCNA.Not Grouped","WGCNA.PROG","WGCNA.unknown1", "WGCNA.unknown2")

net$colors[unique(unlist(layers))]

  Reln   Lhx5  Satb2   Cux1   Rorb   Sox5 Bcl11b  Fezf2    Vim   Hes1 
     4      4      0      0      0      0      0      0      1      1 

# Attention! The next list need to be updated by hand!
groups = list("L.I"=4,"unknown1"=3,"PROG"=1, "unknown2"=2 , "Not Grouped"=0 )

for(layer1 in c("L.I","unknown1","PROG","unknown2","Not Grouped")){
    for(layer2 in c("L.I","L.II.IV", "L.V.VI","PROG")){
    tableMarkersWGCNA[paste0("WGCNA.",layer1),paste0("Loo.",layer2)] =
        sum(names(net$colors[net$colors %in% groups[[layer1]]]) %in% Markers_Loo[[layer2]])
    #tableMarkersWGCNA[paste0("WGCNA.","Not Grouped"),paste0("Loo.",layer2)] =
     #   sum(names(net$colors[!net$colors %in% unlist(groups)]) %in% Markers_Loo[[layer2]])
    }
}

tableMarkersWGCNA

dend %>%
set("labels", ifelse(labels(dend) %in% rownames(pca_1)[rownames(pca_1) %in% unlist(Markers_Loo)], labels(dend), "")) %>%
  set("branches_k_color", value = c("gray80","#4DBBD5FF","#91D1C2FF" ,"gray80","#F39B7FFF","#E64B35FF","#3C5488FF"), k = 7) %>%
    set("labels_cex"=0.1) %>%
    set("branches_lwd", 0.5) %>%
plot_horiz.dendrogram(horiz=T, axes=T,xlim = c(0,80), dLeaf = -15,text_pos = 3)

sessionInfo()

R version 4.0.4 (2021-02-15)
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/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.2.20.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C             LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] gplots_3.1.1          patchwork_1.1.1       dplyr_1.0.4           SeuratObject_4.0.0   
 [5] Seurat_4.0.0          cluster_2.1.1         WGCNA_1.70-3          fastcluster_1.1.25   
 [9] dynamicTreeCut_1.63-1 dendextend_1.14.0     ggrepel_0.9.1         factoextra_1.0.7     
[13] Matrix_1.3-2          data.table_1.13.6     ggplot2_3.3.3         COTAN_0.1.0          

loaded via a namespace (and not attached):
  [1] R.utils_2.10.1        reticulate_1.18       tidyselect_1.1.0      RSQLite_2.2.3        
  [5] AnnotationDbi_1.52.0  htmlwidgets_1.5.3     grid_4.0.4            Rtsne_0.15           
  [9] munsell_0.5.0         preprocessCore_1.52.1 codetools_0.2-18      ica_1.0-2            
 [13] future_1.21.0         miniUI_0.1.1.1        withr_2.4.1           colorspace_2.0-0     
 [17] Biobase_2.50.0        filelock_1.0.2        knitr_1.31            rstudioapi_0.13      
 [21] stats4_4.0.4          ROCR_1.0-11           ggsignif_0.6.1        tensor_1.5           
 [25] listenv_0.8.0         labeling_0.4.2        polyclip_1.10-0       bit64_4.0.5          
 [29] farver_2.0.3          basilisk_1.2.1        parallelly_1.23.0     vctrs_0.3.6          
 [33] generics_0.1.0        xfun_0.20             R6_2.5.0              doParallel_1.0.16    
 [37] clue_0.3-58           bitops_1.0-6          spatstat.utils_2.0-0  cachem_1.0.3         
 [41] gridGraphics_0.5-1    assertthat_0.2.1      promises_1.1.1        scales_1.1.1         
 [45] nnet_7.3-15           gtable_0.3.0          Cairo_1.5-12.2        globals_0.14.0       
 [49] goftest_1.2-2         rlang_0.4.10          GlobalOptions_0.1.2   splines_4.0.4        
 [53] rstatix_0.7.0         lazyeval_0.2.2        impute_1.64.0         checkmate_2.0.0      
 [57] broom_0.7.5           BiocManager_1.30.10   yaml_2.2.1            reshape2_1.4.4       
 [61] abind_1.4-5           backports_1.2.1       httpuv_1.5.5          Hmisc_4.5-0          
 [65] tools_4.0.4           ggplotify_0.0.5       ellipsis_0.3.1        jquerylib_0.1.3      
 [69] RColorBrewer_1.1-2    BiocGenerics_0.36.0   ggridges_0.5.3        Rcpp_1.0.6           
 [73] plyr_1.8.6            base64enc_0.1-3       purrr_0.3.4           basilisk.utils_1.2.2 
 [77] ggpubr_0.4.0          rpart_4.1-15          deldir_0.2-10         pbapply_1.4-3        
 [81] GetoptLong_1.0.5      viridis_0.5.1         cowplot_1.1.1         S4Vectors_0.28.1     
 [85] zoo_1.8-8             haven_2.3.1           magrittr_2.0.1        scattermore_0.7      
 [89] openxlsx_4.2.3        circlize_0.4.12       lmtest_0.9-38         RANN_2.6.1           
 [93] fitdistrplus_1.1-3    matrixStats_0.58.0    hms_1.0.0             mime_0.9             
 [97] evaluate_0.14         xtable_1.8-4          jpeg_0.1-8.1          rio_0.5.16           
[101] readxl_1.3.1          IRanges_2.24.1        gridExtra_2.3         shape_1.4.5          
[105] compiler_4.0.4        tibble_3.0.6          KernSmooth_2.23-18    crayon_1.4.0         
[109] R.oo_1.24.0           htmltools_0.5.1.1     mgcv_1.8-33           later_1.1.0.1        
[113] Formula_1.2-4         tidyr_1.1.2           DBI_1.1.1             ComplexHeatmap_2.6.2 
[117] MASS_7.3-53.1         rappdirs_0.3.3        car_3.0-10            R.methodsS3_1.8.1    
[121] parallel_4.0.4        igraph_1.2.6          forcats_0.5.1         pkgconfig_2.0.3      
[125] rvcheck_0.1.8         foreign_0.8-81        plotly_4.9.3          foreach_1.5.1        
[129] bslib_0.2.4           stringr_1.4.0         digest_0.6.27         sctransform_0.3.2    
[133] RcppAnnoy_0.0.18      spatstat.data_2.0-0   rmarkdown_2.7         cellranger_1.1.0     
[137] leiden_0.3.7          htmlTable_2.1.0       uwot_0.1.10           curl_4.3             
[141] gtools_3.8.2          shiny_1.6.0           rjson_0.2.20          lifecycle_0.2.0      
[145] nlme_3.1-152          jsonlite_1.7.2        carData_3.0-4         viridisLite_0.3.0    
[149] pillar_1.4.7          lattice_0.20-41       fastmap_1.1.0         httr_1.4.2           
[153] survival_3.2-7        GO.db_3.12.1          glue_1.4.2            zip_2.1.1            
[157] spatstat_1.64-1       png_0.1-7             iterators_1.0.13      bit_4.0.4            
[161] sass_0.3.1            stringi_1.5.3         blob_1.2.1            caTools_1.18.1       
[165] latticeExtra_0.6-29   memoise_2.0.0         irlba_2.3.3           future.apply_1.7.0