require(decoupleR)
require(Seurat)
require(ggsci)
require(ggplot2)
source("script/20240910.R")
net <- get_collectri(organism='human', split_complexes=FALSE)
obj <- readRDS("result/obj.pbmc3k.expression.rds")
tfs <- unique(net$source)
targets <- unique(net$target)
genes <- select.gene(obj, tfs = tfs, targets = targets, min.cells = 0)
tfs <- genes$tfs
write.table(as.matrix(obj@assays$RNA$data[tfs,]), "response/logcount.pmbc.tsv", sep = "\t")


################################################################################
#                           Analysis (CSN)
################################################################################
csn.degree <- read.table("response/csndm.tsv")
rownames(csn.degree) <- tfs
colnames(csn.degree) <- colnames(obj)
obj.csn.degree<- CreateSeuratObject(csn.degree)
LayerData(obj.csn.degree, "data") <- obj.csn.degree@assays$RNA$counts
VariableFeatures(obj.csn.degree) <- rownames(obj.csn.degree)
obj.csn.degree <- ScaleData(obj.csn.degree)
obj.csn.degree <- RunPCA(obj.csn.degree)
obj.csn.degree <- FindNeighbors(obj.csn.degree, dims = 1:10)
obj.csn.degree <- FindClusters(obj.csn.degree)
obj.csn.degree <- RunUMAP(obj.csn.degree, dims = 1:10)
obj.csn.degree$celltype <- Idents(obj)

#saveRDS(obj.csn.degree, "response/obj.csn.degree.rds")
obj.csn.degree <- readRDS("response/obj.csn.degree.rds")

DimPlot(obj.csn.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle("oCSN")

################################################################################
#                           Analysis (locCSN)
################################################################################
loccsn.degree <- read.table("response/loc.csndm.tsv") %>% t
rownames(loccsn.degree) <- tfs
colnames(loccsn.degree) <- colnames(obj)
obj.loccsn.degree<- CreateSeuratObject(loccsn.degree)
LayerData(obj.loccsn.degree, "data") <- obj.loccsn.degree@assays$RNA$counts
VariableFeatures(obj.loccsn.degree) <- rownames(obj.loccsn.degree)
obj.loccsn.degree <- ScaleData(obj.loccsn.degree)
obj.loccsn.degree <- RunPCA(obj.loccsn.degree)
obj.loccsn.degree <- FindNeighbors(obj.loccsn.degree, dims = 1:10)
obj.loccsn.degree <- FindClusters(obj.loccsn.degree)
obj.loccsn.degree <- RunUMAP(obj.loccsn.degree, dims = 1:10)
obj.loccsn.degree$celltype <- Idents(obj)

#saveRDS(obj.loccsn.degree, "response/obj.loccsn.degree.rds")
obj.loccsn.degree <- readRDS("response/obj.loccsn.degree.rds")

DimPlot(obj.loccsn.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle("locCSN")


################################################################################
#                           Analysis (NNET)
################################################################################
obj <- prepare.seurat(obj, genes = genes$genes)
obj <- prepare.graph(obj)
obj <- prepare.reg(obj, predictors = tfs,
                   responses = tfs)
obj <- run.nn.reg(obj, tfs, return.smooth = F)

nnet.degree <- apply(2*obj@misc$mod$effect^2, 1, colSums)
obj.nnet.degree<- CreateSeuratObject(t(nnet.degree))
LayerData(obj.nnet.degree, "data") <- obj.nnet.degree@assays$RNA$counts
VariableFeatures(obj.nnet.degree) <- rownames(obj.nnet.degree)
obj.nnet.degree <- ScaleData(obj.nnet.degree)
obj.nnet.degree <- RunPCA(obj.nnet.degree)
obj.nnet.degree <- FindNeighbors(obj.nnet.degree, dims = 1:10)
obj.nnet.degree <- FindClusters(obj.nnet.degree)
obj.nnet.degree <- RunUMAP(obj.nnet.degree, dims = 1:10)
obj.nnet.degree$celltype <- Idents(obj)

#saveRDS(obj.nnet.degree, "response/obj.nnet.degree.rds")
obj.nnet.degree <- readRDS("response/obj.nnet.degree.rds")

DimPlot(obj.nnet.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle("NNet")

################################################################################
#                           Plot
################################################################################
obj <- readRDS("result/obj.pbmc3k.expression.rds")
obj.csn.degree <- readRDS("response/obj.csn.degree.rds")
obj.loccsn.degree <- readRDS("response/obj.loccsn.degree.rds")
obj.nnet.degree <- readRDS("response/obj.nnet.degree.rds")

sil <- cluster::silhouette(as.numeric(obj$seurat_clusters),
                           dist(obj.csn.degree@reductions$pca@cell.embeddings[,1:10]))[,3]  %>% median
p1 <- DimPlot(obj.csn.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle(paste0("oCSN"), paste0("MSI=", round(sil,2), "; Time: 77.2 Sec"))

sil <- cluster::silhouette(as.numeric(obj$seurat_clusters),
                           dist(obj.loccsn.degree@reductions$pca@cell.embeddings[,1:10]))[,3]  %>% median
p2 <- DimPlot(obj.loccsn.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle(paste0("LocCSN"), paste0("MSI=", round(sil,2), "; Time: 4883.3 Sec"))

sil <- cluster::silhouette(as.numeric(obj$seurat_clusters),
                           dist(obj.nnet.degree@reductions$pca@cell.embeddings[,1:10]))[,3] %>% median
p3 <- DimPlot(obj.nnet.degree, group.by = "celltype", pt.size=0.2)   +
  theme_classic()+
  theme(text = element_text(size = 13, face = "bold"),
        axis.line = element_blank(),
        legend.position = "none",
        axis.text = element_blank(),
        axis.ticks = element_blank(),
  ) + xlab("")+ylab("")+
  scale_color_bmj()+
  ggtitle(paste0("NNet "), paste0("MSI=", round(sil,2), "; Time: 584.8 Sec"))

p1 +p2 + p3

genes$genes %>% length
