#' @description: Endothelium cells

library(Seurat)

library(harmony)

library(clustree)

library(ggpubr)

library(dplyr)

library(patchwork)

library(ComplexHeatmap)

library(circlize)

library(vegan)

set.seed(101)

library(openxlsx)

library(future)

plan("multiprocess", workers = 10) 

options(future.globals.maxSize = 50000 * 1024^2) # set 50G RAM

setwd(dir = "/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA")

source(file = "/home/longzhilin/Analysis_Code/Visualization/colorPalettes.R")

source(file = "/home/longzhilin/Analysis_Code/code/ratio.plot.R")

scRNA.data <- readRDS("data.merge.pro.rds")

###########################1.differential gene#################################

Idents(scRNA.data) <- scRNA.data$cellType_low

idents <- as.character(levels(scRNA.data))

Endothelium.markers <- FindMarkers(scRNA.data, 

                                    ident.1 = "Endothelium (VCAM1+)",

                                    ident.2 = "Endothelium (VCAM1-)",

                                    logfc.threshold = 0, 

                                    min.pct = 0.1, 

                                    test.use = "MAST", 

                                    latent.vars = "orig.ident")

Endothelium.markers <- arrange(Endothelium.markers, desc(avg_log2FC))

Endothelium.markers$gene <- rownames(Endothelium.markers)

Endothelium.DEGs <- Endothelium.markers %>% filter(abs(avg_log2FC) >= 0.25 & p_val_adj < 0.05)

write.xlsx(Endothelium.markers, file = "2.Cluster/AnnotateCellType/Endothelium.all.DEGs.xlsx", rowNames = F, overwrite = T)

saveRDS(Endothelium.markers, file = "2.Cluster/AnnotateCellType/Endothelium.all.DEGs.rds")

#### volcano map

source("/home/longzhilin/Analysis_Code/Visualization/Plot.EnhancedVolcano.R")

pdf("2.Cluster/AnnotateCellType/Endothelium.EnhancedVolcano.pdf")

Plot.EnhancedVolcano(Endothelium.DEGs, x = "avg_log2FC", y = "p_val_adj", selected.showType = "Both", select.num = 15, title = "Endothelium (VCAM1+) VS Endothelium (VCAM1-)")

dev.off()

#### pathway enrichment

source("/home/longzhilin/Analysis_Code/PathwayEnrichment/clusterProfiler.enricher.R")

up.DEGs <- filter(Endothelium.DEGs, avg_log2FC > 1 & p_val_adj < 0.05)

down.DEGs <- filter(Endothelium.DEGs, avg_log2FC< -1 & p_val_adj < 0.05)

DEGs <- list(up = up.DEGs$gene, down = down.DEGs$gene)

pdf("2.Cluster/AnnotateCellType/Endothelium.program.pathway.enrichment.pdf")

res <- lapply(names(DEGs), function(x){

    y <- DEGs[[x]]

    res <- cluterProfiler.enricher(gene = y, geneType = "SYMBOL", db.type = "MsigDB", saveDir = paste0(getwd(),"/2.Cluster/AnnotateCellType"),

                            title = x, qvalueCutoff = 0.05, pvalueCutoff = 0.05)

    # gene ratio

    res <- res$em.res.genesymbol@result %>% filter(p.adjust<0.05) #fdr adjust

    pathway.geneNum <- unlist(strsplit(res$BgRatio, "/"))[seq(1, nrow(res),by=2)]

    gene.ratio <- as.numeric(res$Count)/as.numeric(pathway.geneNum)

    res$gene.ratio <- gene.ratio

    return(res)

})

dev.off()

####cluster profiler enrichment result display

library(ggthemes)

#top10

enrich1 <- read.csv("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/up_enricherResult.csv", header = T)

enrich2 <- read.csv("/data/active_data/lzl/RenalTumor-20200713/DataAnalysis-20210803/scRNA/2.Cluster/AnnotateCellType/down_enricherResult.csv", header = T)

enrich.list <- list(up = enrich1, down = enrich2)

enrich <- list(up = enrich1[c(1,3,8,9,16,17,18,21,22,23,25,26,27,28,42,49,50,72), c("ID")],

               down = enrich2[c(2:5,8:10,11,12,19,21,22), c("ID")])

enrich.path <- Reduce(function(x,y) c(x,y), enrich)

enrich.path <- unique(enrich.path)

enrich.path <- enrich.path[-c(3,6,26,28,29,4,13,11,5,20,23,27,30)]

enrich.list.select <- lapply(enrich.list, function(x){

  idx <- which(x$ID %in% enrich.path)

  return(x[idx, c("ID", "Count", "pvalue", "p.adjust")])

})

enrich.res <- Reduce(function(x,y) rbind(x,y), enrich.list.select)

a <- sapply(enrich.list.select, function(x){nrow(x)})

enrich.res$Type <- c(rep("Endothelium (VCAM1+) ", as.numeric(a[1])), rep("Endothelium (VCAM1-) ", as.numeric(a[2])))

enrich.res$p.adjust <- -log10(enrich.res$p.adjust)

# process gene name

enrich.res$ID <- gsub("_", " ", enrich.res$ID)

# factor ID

idx1 <- na.omit(match(enrich$up, enrich.path))

idx2 <- na.omit(match(enrich$down, enrich.path))

enrich.path <- enrich.path[c(idx1, idx2)]

enrich.res$ID <- factor(enrich.res$ID, levels = gsub("_", " ", enrich.path))

#source(file = "/home/longzhilin/Analysis_Code/Visualization/ggplot.bubble.plot.R")

pdf("2.Cluster/AnnotateCellType/Endothelium.enrichment.bubble.pdf", width = 5, height = 5)

p1 <- ggplot(enrich.res, 

            aes(x = Type, 

                y = ID, 

                size = Count, 

                fill = p.adjust))

p2 <- p1 + guides(color=FALSE) + geom_point(shape = 21, alpha = 0.7) + theme_few() + scale_size_continuous(breaks = c(3,5,7,9,11), labels = c(3,5,7,9,11)) + scale_fill_gradientn(colours = c("#fed71b", "#fcc41c", "#f27620", "#e92925"), breaks = c(1.69897, 4, 6, 8), labels = c(0.02, 0.0001, 0.000001, 0.00000001))

p2 <- p2 + theme(axis.text.x = element_text(angle = 45,vjust = 1,hjust = 1, size = 8), axis.text.y = element_text(size = 8)) + xlab("") + ylab("")

print(p2)

dev.off()