Switch to side-by-side view

--- a
+++ b/RNA-seq/AnalysisPipeline/Endothelial.R
@@ -0,0 +1,101 @@
+#' @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()