require(Seurat)
require(magrittr)
require(decoupleR)
require(reshape2)
require(tidyverse)
require(tidyr)
require(broom)
require(limma)
require(ggsci)
require(fossil)
require(clusterProfiler)
require(msigdbr)
require(wordcloud)
require(org.Hs.eg.db)
require(dplyr)
require(plyr)
require(tidyr)
require(pryr)
require(ggplot2)
source("script/20240910.R")

find.npc <- function(sd){
  k <- length(sd)
  s <- abs(diff(sd))
  mu <- mean(s[(0.8*k):k-1])
  sigma <- sd(s[(0.8*k):k-1])
  sk <- mu + 6*sigma
  npcs <- max(which(s > sk))+1
  npcs
}
################################# Load data and Seurat preprocessing #################################
counts <- read.delim("data/hsc/GSM3305366_LinNegCD34PosCD164Pos.raw_counts.tsv.gz", header = T, row.names = "Barcode")[,-1]
LinNegCD34PosCD164Pos <- CreateSeuratObject(counts = counts, project = "LinNegCD34PosCD164Pos")
counts <- read.delim("data/hsc/GSM3305367_LinNegCD34NegCD164high.raw_counts.tsv.gz", header = T, row.names = "Barcode")[,-1]
LinNegCD34NegCD164high <- CreateSeuratObject(counts = counts, project = "LinNegCD34NegCD164high")
counts <- read.delim("data/hsc/GSM3305368_LinNegCD34lowCD164high.raw_counts.tsv.gz", header = T, row.names = "Barcode")[,-1]
LinNegCD34lowCD164high <- CreateSeuratObject(counts = counts, project = "LinNegCD34lowCD164high")
counts <- read.delim("data/hsc/GSM3305369_LinNegCD34NegCD164low.raw_counts.tsv.gz", header = T, row.names = "Barcode")[,-1]
LinNegCD34NegCD164low <- CreateSeuratObject(counts = counts, project = "LinNegCD34NegCD164low")
obj <- merge(LinNegCD34PosCD164Pos, y = c(LinNegCD34NegCD164high, LinNegCD34lowCD164high,
                                          LinNegCD34NegCD164low),
             add.cell.ids = c("LinNegCD34PosCD164Pos", "LinNegCD34NegCD164high", "LinNegCD34lowCD164high",
                              "LinNegCD34NegCD164low"), project = "LinNeg")

# Seurat pipeline
obj <- NormalizeData(obj)
obj <- FindVariableFeatures(obj)
obj <- ScaleData(obj)
obj <- RunPCA(obj, npcs = 100)
obj <- FindNeighbors(obj, dims = 1:10)
obj <- FindClusters(obj, resolution = 0.4)
obj <- RunUMAP(obj, dims = 1:10)

# Annotation
new.cluster.ids <- c("P", "ErP", "Ery", "MEP", "GMP", "MP",
                     "Ery", "Lym", "Neu", "B/E/Mast", "Ery")
names(new.cluster.ids) <- levels(obj)
obj <- RenameIdents(obj, new.cluster.ids)
DimPlot(obj, group.by = "orig.ident")
DimPlot(obj, label = T)

# Save Seurat object
# saveRDS(obj, "result/obj.LinNeg.rds")

################################# Run NNet #################################
obj <- readRDS("result/obj.LinNeg.rds")
obj$celltype <- Idents(obj)
obj <- JoinLayers(obj)
rt.ppr <- get.ppr()
genes <- select.gene(obj)
obj <- prepare.seurat(obj, genes = genes$genes)

mem_change({
  obj <- prepare.graph(obj)
  # obj <- select.cell(obj, n = 500)
  cells <- readLines("result/LinNeg.cells")
  obj@misc$setting$cells <- as.numeric(cells)
  names(obj@misc$setting$cells) <- colnames(obj)[as.numeric(cells)]

  t1 <- system.time({
    obj <- prepare.reg(obj, predictors = genes$targets,
                       responses = genes$tfs)
  })
  responses <- genes$tfs

  t2 <- system.time({
    obj <- run.nn.reg(obj, responses = responses, return.smooth = T, return.prune = T)
  })
  #obj <- set.defaults(obj, defaults = list(cutoff = 0.99))

  t3 <- system.time({
    obj <- build.meta.network(obj, return.p.val = F)
    cells <- obj@misc$setting$cells
  })
})
t1
t2
t3
# Save cells used for the analysis
#writeLines(obj@misc$setting$cells %>% as.character, "result/LinNeg.cells")

################################# Figure A1 Dim plot #################################
p <- DimPlot(obj, pt.size = 0.8)  + theme(text = element_text(face = "bold", size = 13),
                                          legend.position = "none",
                                          axis.text = element_blank(), axis.ticks = element_blank()
) + xlab("UMAP1")+ylab("UMAP2")+
  scale_color_bmj()+
  ggtitle("Pellin et al., 2019", "Early hematopoiesis atlas")

LabelClusters(p, id = "ident", fontface = "bold", color = "black",size = 4, box = 1.5, alpha=0.8)

################################# Potentially supplementary #################################
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2), size = 1, col = "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,], aes(umap_1, umap_2), col = "black", size = 1)+
  theme_classic()+
  scale_color_gradientn(colors = gg.color.spec(11), limits = c(0,1), )+ theme(text = element_text(face = "bold"),
                                                                              legend.position = "none",
                                                                              axis.line = element_blank(),
                                                                              axis.text = element_blank(), axis.ticks = element_blank()
  )+ xlab("")+ylab("")

################################# Supplementary Figure A nPCA loading heatmap #################################
obj.sub <- CreateSeuratObject(t(obj@misc$mod$meta.network$npca.loadings),
                              meta.data = cbind(obj@meta.data[cells,], obj@misc$mod$meta.network$npca.loadings))
Idents(obj.sub) <- obj.sub$celltype
LayerData(obj.sub, "data") <- obj.sub@assays$RNA$counts
LayerData(obj.sub, "scale.data") <- obj.sub@assays$RNA$counts
npcs <- obj@misc$mod$meta.network$npca.loadings[order(obj.sub$celltype),]
row.anno <- data.frame(Celltype = obj.sub$celltype[rownames(npcs)],  row.names = rownames(npcs))
colnames(npcs) <-  paste("MN",1:20, sep = "")
row.color <- ggsci::pal_bmj()(9)
names(row.color) <- levels(row.anno$Celltype)
pheatmap::pheatmap(t(npcs), color = colorRampPalette(RColorBrewer::brewer.pal(9,"BuPu"))(100),
                   cluster_rows = F, cluster_cols = F, show_colnames = F, annotation_col = row.anno,
                   annotation_colors = list(Celltype = row.color), fontsize = 11
)

################################# Figure C, D, Supplementary Figure B Meta-network and nPCA loading dimplot #################################
# Select central genes to show on heatmaps representing adjacency matrices of meta-networks
central.genes <- select.central.genes(obj, n.net = 2, k = 1, n.per.component = 10)
row.anno <- data.frame(Module = central.genes$response.module, row.names = central.genes$central.responses)
col.anno <- data.frame(Module = central.genes$predictor.module, row.names = central.genes$central.predictors)
anno.colors <- list(Module = rainbow(2))
names(anno.colors$Module) <- paste("M",1:2, sep = "")

################################# Figure C Net1 #################################
mat <- get.network(obj, i = 1, assay = "meta.network")
new_row_labels <- paste0("italic(", central.genes$central.responses, ")")
new_row_labels <- as.expression(sapply(new_row_labels, str2lang))
new_col_labels <- paste0("italic(", central.genes$central.predictors, ")")
new_col_labels <- as.expression(sapply(new_col_labels, str2lang))

# Heatmap
pheatmap::pheatmap(mat[central.genes$central.responses,central.genes$central.predictors],
                   cluster_rows = F, cluster_cols = F, color = RColorBrewer::brewer.pal(9,"Greys"), show_colnames = T,
                   show_rownames = T, annotation_row = row.anno, annotation_col = col.anno, annotation_colors =anno.colors,
                   gaps_col = table(col.anno$Module) %>% cumsum, gaps_row = table(row.anno$Module) %>% cumsum, fontsize = 11,
                   main = "Meta-network 1 adjacency matrix",
                   labels_row = new_row_labels,
                   labels_col =new_col_labels )

# Loading dimplot
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2),size=1, col= "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj@misc$mod$meta.network$npca.loadings[,1]),size=1)+
  scale_color_gradientn(colors = RColorBrewer::brewer.pal(9,"BuPu"))+
  theme_classic()+ theme(text = element_text(face = "bold", size=13),
                         legend.position = "none",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank()
  )+ xlab("")+ylab("") + ggtitle("Meta-network 1")

################################# Figure C Net2 #################################
mat <- get.network(obj, i = 2, assay = "meta.network")
new_row_labels <- paste0("italic(", central.genes$central.responses, ")")
new_row_labels <- as.expression(sapply(new_row_labels, str2lang))
new_col_labels <- paste0("italic(", central.genes$central.predictors, ")")
new_col_labels <- as.expression(sapply(new_col_labels, str2lang))

# Heatmap
pheatmap::pheatmap(mat[central.genes$central.responses,central.genes$central.predictors],
                   cluster_rows = F, cluster_cols = F, color = RColorBrewer::brewer.pal(9,"Greys"),
                   show_colnames = T, show_rownames = T,
                   annotation_row = row.anno, annotation_col = col.anno, annotation_colors =anno.colors,
                   gaps_col = table(col.anno$Module) %>% cumsum, gaps_row = table(row.anno$Module) %>% cumsum, fontsize = 11,
                   main = "Meta-network 2 adjacency matrix",
                   labels_row = new_row_labels,
                   labels_col =new_col_labels)

# Loading dimplot
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2),size=1, col= "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj@misc$mod$meta.network$npca.loadings[,2]),size=1)+
  scale_color_gradientn(colors = RColorBrewer::brewer.pal(9,"BuPu"))+
  theme_classic()+ theme(text = element_text(face = "bold", size = 13),
                         legend.position = "none",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank()
  )+ xlab("")+ylab("") + ggtitle("Meta-network 2")

# Legend
ggplot() + geom_col(aes(x = c(1,2), y = c(1,2), fill = c("M1","M2"), col = c(1,12))) +
  theme_classic()+ theme(text = element_text(face = "bold", size = 13),
                         legend.position = "bottom",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank()
  ) + scale_fill_manual("Central gene of",values = anno.colors$Module) +
  scale_color_gradientn("Weight", colours = RColorBrewer::brewer.pal(9,"Greys"))

################################# Supplementary Figure B1 Net3 #################################
central.genes <- select.central.genes(obj, n.net = 4, k = 1, n.per.component = 10)
row.anno <- data.frame(Module = central.genes$response.module, row.names = central.genes$central.responses)
col.anno <- data.frame(Module = central.genes$predictor.module, row.names = central.genes$central.predictors)
anno.colors <- list(Module = rainbow(4))
names(anno.colors$Module) <- paste("M",1:4, sep = "")
row.anno <- subset(row.anno, Module %in% c("M3","M4"))
col.anno <- subset(col.anno, Module %in% c("M3","M4"))
new_row_labels <- paste0("italic(", rownames(row.anno), ")")
new_row_labels <- as.expression(sapply(new_row_labels, str2lang))
new_col_labels <- paste0("italic(", rownames(col.anno), ")")
new_col_labels <- as.expression(sapply(new_col_labels, str2lang))

mat <- get.network(obj, i = 3, assay = "meta.network")
# Heatmap
pheatmap::pheatmap(mat[rownames(row.anno),rownames(col.anno)],
                   cluster_rows = F, cluster_cols = F, color = RColorBrewer::brewer.pal(9,"Greys"), show_colnames = T,
                   show_rownames = T, annotation_row = row.anno, annotation_col = col.anno, annotation_colors =anno.colors,
                   gaps_col = table(col.anno$Module) %>% cumsum, gaps_row = table(row.anno$Module) %>% cumsum, fontsize = 11,
                   main = "Meta-network 3 adjacency matrix",
                   labels_row = new_row_labels,
                   labels_col = new_col_labels)

# Loading dimplot
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2),size=1, col= "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj@misc$mod$meta.network$npca.loadings[,3]),size=1)+
  scale_color_gradientn(colors = RColorBrewer::brewer.pal(9,"BuPu"))+
  theme_classic()+ theme(text = element_text(face = "bold", size = 13),
                         legend.position = "none",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank()
  )+ xlab("")+ylab("") + ggtitle("Meta-network 3")

################################# Supplementary Figure B2 Net4 #################################
mat <- get.network(obj, i = 4, assay = "meta.network")
# Heatmap
pheatmap::pheatmap(mat[rownames(row.anno),rownames(col.anno)],
                   cluster_rows = F, cluster_cols = F, color = RColorBrewer::brewer.pal(9,"Greys"), show_colnames = T,
                   show_rownames = T, annotation_row = row.anno, annotation_col = col.anno, annotation_colors =anno.colors,
                   gaps_col = table(col.anno$Module) %>% cumsum, gaps_row = table(row.anno$Module) %>% cumsum, fontsize = 11,
                   main = "Meta-network 4 adjacency matrix")

# Loading dimplot
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2),size=1, col= "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = obj@misc$mod$meta.network$npca.loadings[,4]),size=1)+
  scale_color_gradientn(colors = RColorBrewer::brewer.pal(9,"BuPu"))+
  theme_classic()+ theme(text = element_text(face = "bold", size = 13),
                         legend.position = "none",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank()
  )+ xlab("")+ylab("") + ggtitle("Meta-network 4")

################################# Supplementary Figure C Pseudotime analysis #################################
dm <- destiny::DiffusionMap(obj@reductions$pca@cell.embeddings[,1:10], n_pcs = NA)
dpt <- destiny::DPT(dm)
time <- dpt$DPT2666[cells]
DimPlot(obj, cells.highlight = colnames(obj)[2666])

################################# Supplementary Figure C1: GATA degree vs DPT #################################
mat <- get.network(obj, responses = "GATA1")
mat1 <- get.network(obj, responses = "GATA2")

p1 <- ggplot() +
  geom_path(aes(x = time[order(time)], y = colSums(mat[,order(time)]) %>% log, col = obj$celltype[cells][order(time)]), alpha = 0.5) +
  geom_point(aes(x = time[order(time)], y = colSums(mat[,order(time)]) %>% log, col = obj$celltype[cells][order(time)])) +
  scale_color_bmj() +
  theme_classic()+
  theme(text = element_text(face = "bold", size = 13), plot.title = element_text(hjust = 0.5),
        axis.title.y = element_text(face = "bold.italic")
  ) + xlab("") + ylab("GATA1") +
  guides(color=guide_legend(title="Celltype")) + ggtitle("Log degree")

p2 <- ggplot() +
  geom_path(aes(x = time[order(time)], y = colSums(mat1[,order(time)]) %>% log, col = obj$celltype[cells][order(time)]), alpha = 0.5) +
  geom_point(aes(x = time[order(time)], y = colSums(mat1[,order(time)]) %>% log, col = obj$celltype[cells][order(time)])) +
  scale_color_bmj() +
  theme_classic()+
  theme(text = element_text(face = "bold", size = 13),
        axis.title.y = element_text(face = "bold.italic")
  ) + xlab("") + ylab("GATA2") +
  guides(color=guide_legend(title="Celltype"))

p3 <- ggplot() +
  geom_path(aes(x = time[order(time)], y = obj@assays$RNA$data["GATA1", cells][order(time)], col = obj$celltype[cells][order(time)]), alpha = 0.5) +
  geom_point(aes(x = time[order(time)], y = obj@assays$RNA$data["GATA1", cells][order(time)], col = obj$celltype[cells][order(time)])) +
  scale_color_bmj() +
  theme_classic()+
  theme(text = element_text(face = "bold", size = 13), plot.title = element_text(hjust = 0.5),
        axis.title.y = element_text(face = "bold.italic")
  ) + xlab("") + ylab("GATA1") +
  guides(color=guide_legend(title="Celltype")) + ggtitle("Log expression")

p4 <- ggplot()  +
  geom_path(aes(x = time[order(time)], y = obj@assays$RNA$data["GATA2", cells][order(time)], col = obj$celltype[cells][order(time)]), alpha = 0.5) +
  geom_point(aes(x = time[order(time)], y = obj@assays$RNA$data["GATA2", cells][order(time)], col = obj$celltype[cells][order(time)])) +
  scale_color_bmj() +
  theme_classic()+
  theme(text = element_text(face = "bold", size = 13),
        axis.title.y = element_text(face = "bold.italic")
  ) + xlab("Pseudotime") + ylab("GATA2") +
  guides(color=guide_legend(title="Celltype"))

ggarrange(p1,p2,p3,p4, ncol = 1, align = "hv", common.legend = T,legend = "bottom" )

################################# Supplementary Figure D Meta-TF and nPCA loading dimplot #################################
meta.response <- build.meta.response(obj, scale = T, return.p.val = F)
npcs <- meta.response$npca.loadings
colnames(npcs) <-  paste("MT",1:20, sep = "")
top.tfs <- apply(npcs, 2, function(x) rownames(npcs)[rank(-x)<=50]) %>% as.vector()
breaks <- rep(paste("MT",1:20, sep = ""), each = 50)
breaks <- breaks[!duplicated(top.tfs)]
top.tfs <- top.tfs[!duplicated(top.tfs)]
breaks <- factor(breaks, levels = paste("MT",1:20, sep = ""))
col.anno <- data.frame(module = breaks, row.names = top.tfs)
col.colors <- list(module = rainbow(20))
names(col.colors$module) <- paste("MT",1:20, sep = "")
pheatmap::pheatmap(t(npcs[top.tfs,]), colorRampPalette(RColorBrewer::brewer.pal(9,"Reds"))(100),
                   cluster_rows = F, cluster_cols = F, show_colnames = F,show_rownames = T, fontsize = 13 , annotation_col = col.anno,
                   annotation_colors = col.colors, gaps_col = table(breaks) %>% cumsum()
)

################################# Meta-network degree #################################
tf.degree <- apply(meta.response$meta.network, 3, rowSums)

obj.sub <- CreateSeuratObject(tf.degree,
                              meta.data = obj@meta.data[cells,])
Idents(obj.sub) <- obj.sub$celltype
tf.degree <- tf.degree[,order(obj.sub$celltype)]
row.anno <- data.frame(Celltype = obj.sub$celltype[colnames(tf.degree)],  row.names = colnames(tf.degree))
rownames(tf.degree) <-  paste("MT",1:20, sep = "")
row.color <- ggsci::pal_bmj()(9)
names(row.color) <- levels(row.anno$Celltype)
pheatmap::pheatmap(tf.degree, color = colorRampPalette(RColorBrewer::brewer.pal(9,"Reds"))(100),
                   cluster_rows = F, cluster_cols = F, show_colnames = F, annotation_col = row.anno,
                   annotation_colors = list(Celltype = row.color), fontsize = 11
)



################################# Figure D Meta-TF 8 degree dimplot #################################
# Meta-TF 7 (sometime 8) is the meta-TF that show the highest degree during monopoesis lineage.
tf.degree <- apply(meta.response$meta.network, 1, colSums)
ggplot() + geom_point(data = obj@reductions$umap@cell.embeddings, aes(umap_1, umap_2),size=1, col= "lightgrey") +
  geom_point(data = obj@reductions$umap@cell.embeddings[cells,],
             aes(umap_1, umap_2, col = tf.degree[,7]),size=1)+
  scale_color_gradientn(colors = RColorBrewer::brewer.pal(9,"Reds"))+
  ggtitle("Meta-TF 7 degree", "Meta-TF central in MP lineage")+
  theme_classic()+ theme(text = element_text(face = "bold", size = 13),
                         legend.position = "none",
                         axis.line = element_blank(),
                         axis.text = element_blank(), axis.ticks = element_blank() , title = element_text()
  )+ xlab("")+ylab("")

################################# Figure E MSigDB enrichment analysis #################################
# Save Meta-TF anlaysis nPCA loadings
# saveRDS(glist, file = "result/LinNeg.nmf.loading.rds")
# glist <- readRDS("result/LinNeg.nmf.loading.rds")

# Over-representation analysis
glist <- sort(npcs[,7], decreasing = T)
em <- enrichGO(gene          = names(glist)[1:50],
               universe = names(glist),
               OrgDb         = org.Hs.eg.db,
               ont           = "BP",
               pAdjustMethod = "fdr",
               readable      = TRUE,
               keyType = "SYMBOL")

# Calculate Pi-value
em@result <- em@result %>% mutate(GeneRatio1 = sapply(GeneRatio, function(x) eval(parse(text = x)))) %>%
  mutate(BgRatio1 = sapply(BgRatio, function(x) eval(parse(text = x)))) %>%
  mutate(fc = GeneRatio1/BgRatio1) %>%
  mutate(pi = -log10(p.adjust) * fc)
em <- simplify(em, by = "pi", select_fun =  max)
em@result$Description <- toupper(em@result$Description)

# Simplify terms
showCategory <- em@result %>% filter(!duplicated(geneID)) %>% arrange(desc(pi)) %>% dplyr::select(Description)
cnetplot(em, categorySize="pi", foldChange=glist, showCategory = showCategory[1:5,], face = "bold", cex_label_gene = 1.2, cex_label_category = 1.2)+
  theme(text = element_text(size=15, face = "bold"))+
  guides(size = "none", color = "none")

# Title
ggplot() + ggtitle("GO over-representation on \n top 50 weighted TFs of Meta-TF 7") +
  theme(text = element_text(size=20, face = "bold"), plot.title = element_text(hjust = 0.5))

################################# Supplementary Figure F Meta-TF 8 module identification #################################
# Top 500 targets most connected to Meta-TF 8
top.targets <- apply(meta.response$meta.network[7,,], 1, sum)
top.targets <- names(sort(top.targets) %>% tail(n = 500))

# Calculate these targets' correlation in co-expression with meta-TF 8 and visualise correlation heatmap
cor.mat <- cor(t(meta.response$meta.network[7,top.targets,]))
hclu.mod <- hclust(dist(cor.mat))
tree.cluster <- cutree(hclu.mod, k = 4)
tree.cluster <- data.frame(Cluster = factor(tree.cluster))
anno.col <- ggsci::pal_lancet()(4)
anno.col[4] <- "orange"
anno.col <- anno.col[c(3,2,1,4)]
names(anno.col) <- 1:4
pheatmap::pheatmap(cor.mat, breaks = seq(-1,1, length.out = 100), fontsize = 11,
                   colorRampPalette(gg.color.spec(11))(100), show_rownames = F, show_colnames = F,
                   annotation_row = tree.cluster,
                   annotation_col = tree.cluster, annotation_colors = list(Cluster = anno.col))

################################# Figure G Trace plot #################################
celltype <- obj$celltype[cells]
lineage <- as.character(celltype)
lineage[celltype%in%c("GMP","Mono/DC")] <- "MP"
lineage[celltype%in%c("ErP","MEP", "Ery")] <- "Ery"
plot.data <- data.frame(time, t( meta.response$meta.network[7,top.targets,]), celltype = celltype, lineage = lineage)
plot.data <- plot.data[order(time),]
plot.data <- reshape2::melt(plot.data, id.vars = c("time", "celltype", "lineage"))
plot.data$Cluster <- 1
plot.data$Cluster[plot.data$variable %in% rownames(tree.cluster)[tree.cluster==2]] <- 2
plot.data$Cluster[plot.data$variable %in% rownames(tree.cluster)[tree.cluster==3]] <- 3
plot.data$Cluster[plot.data$variable %in% rownames(tree.cluster)[tree.cluster==4]] <- 4
plot.data$Cluster <- as.factor(plot.data$Cluster)
subset(plot.data, lineage %in% c("P", "MP")) %>% ggplot() +
  geom_point(aes(time, sqrt(value), col = Cluster, shape = celltype), alpha = 0.5) +
  geom_path(aes(time, sqrt(value), col = Cluster, group = variable), alpha = 0.1) +
  theme_classic()+
  scale_shape_manual("Celltype", values = c(1, 16, 8, 3))+
  scale_color_manual(values = anno.col)+ theme(text = element_text(face = "bold", size = 13),
                                               plot.title = element_text(hjust = 0.5))+
  xlab("Pseudotime (along MP lineage)")+ylab("Coexpression") +
  ggtitle("Coexpression with Meta-TF 7")+
  guides(colour = guide_legend(order = 1),
         shape = guide_legend(order = 2))

################################# Figure G Check regulon composition #################################
top.targets[tree.cluster == 1] %>% cat(., sep = ",")

# Check regulon 4 hub genes
module.act <- meta.response$meta.network[7,,]
v <- RSpectra::svds(module.act[top.targets[tree.cluster == 4],], k = 1)$v
module.membership1 <- cor(t(module.act), v)[,1]

module.membership1 <- module.membership1[top.targets[tree.cluster == 4]]
module.membership1[abs(module.membership1) > 0.9] %>%
  sort(.,decreasing = F)  %>% names %>% cat(., sep = ", ")

# Check regulon 1 hub genes
v <- RSpectra::svds(module.act[top.targets[tree.cluster == 1],], k = 1)$v
module.membership2 <- cor(t(module.act), v)[,1]
module.genes <- c(top.targets[tree.cluster == 4] , top.targets[tree.cluster == 1])

module.membership2 <- module.membership2[top.targets[tree.cluster == 1]]
module.membership2[abs(module.membership2) > 0.9] %>%
  sort(.,decreasing = F)  %>% names %>% cat(., sep = ", ")

top.targets.no.tfs <- top.targets[!top.targets %in% genes$tfs]

################################# Discarded Figure G Trace plot of regulon mean co-expression #################################
celltype <- obj$celltype[cells]
lineage <- as.character(celltype)
lineage[celltype%in%c("GMP","Mono")] <- "MP"
lineage[celltype%in%c("ErP","MEP", "Ery")] <- "Ery"
plot.data.mean <- data.frame(cluster = factor(tree.cluster[,1]), meta.response$meta.network[7,top.targets,])
plot.data.mean <- plot.data.mean %>% group_by(cluster) %>%
  summarise_each(funs(mean)) %>% dplyr::select(-cluster) %>%
  t %>% data.frame %>% mutate(time = time, celltype = celltype, lineage = lineage) %>%
  reshape2::melt(id.vars = c("time", "celltype", "lineage")) %>%
  mutate(Cluster = gsub("X","", variable)) %>% arrange(time)

subset(plot.data.mean, lineage %in% c("HSC", "Mye")) %>% ggplot() +
  geom_point(aes(time, value, col = Cluster, shape = celltype), alpha = 1) +
  geom_smooth(aes(time, value, col = Cluster), alpha = 0.,  span  = 0.2) +
  theme_classic()+
  scale_shape_manual("Celltype", values = c(1, 16, 8, 3))+
  scale_color_manual(values = anno.col)+ theme(text = element_text(face = "bold", size = 15),
                                               legend.position = "bottom",
                                               , legend.box="vertical"
  )+ xlab("Pseudotime")+ylab("Co-expression") +
  ggtitle("Co-expression with Meta-TF 7 withikn mononuclear phagocyte lineage")+
  guides(colour = guide_legend(order = 1),
         shape = guide_legend(order = 2)) +
  ylim(c(0,NA))

