--- title: "scEC_ana" author: "zhouruihan" date: "2025-10-27" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE) ``` # 1.Seu ```{r} library(Seurat) ``` 1106 个肝内皮细胞转录组进行分析，鉴定出 8 个可区分的簇，证实了先前报道的分区。基于普遍接受的标记进行注释，例如门静脉内皮细胞的Ly6a、Cd34 ，窦状隙内皮细胞的Bmp2、Mrc1、Gpihbp1、Aqp1 ，以及中央静脉内皮细胞的Rspo3、Wnt9b、Wnt2 ```{r} count <- fread("../00raw_data/GSE297062_endo_counts.txt") meta <- fread("../00raw_data/GSE297062_endo_meta_data_bam.txt") table(meta$p2_multilevel_cluster) ``` ```{r} count <- count %>% column_to_rownames(var = "V1") mat <- as.matrix(count) mat <- Matrix(mat, sparse = TRUE) seu <- CreateSeuratObject( counts = mat, project = "Mouse_liver_endo", min.cells = 3, # 可按需要调整 min.features = 200 # 可按需要调整 ) ``` ```{r} seu <- NormalizeData(seu) seu <- FindVariableFeatures(seu, selection.method = "vst", nfeatures = 2000, verbose = FALSE) seu <- ScaleData(seu) seu <- RunPCA(seu) ``` ```{r} n.pcs <- 20 res.used <- 0.3 seu <- FindNeighbors(seu, dims = 1:n.pcs, verbose = F) seu <- FindClusters(seu, resolution = res.used, verbose = F) seu <- RunUMAP(seu, dims = 1:n.pcs, reduction = "pca", verbose = F) seu <- RunTSNE(seu, dims = 1:n.pcs, reduction = "pca", verbose = F) ``` 文章中的e2是Ly6a和Cd34阳性的群，和我们的一致，直接提e2出来分析 ```{r} e2 <- meta[meta$p2_multilevel_cluster == 2,]$cell_ID ``` ```{r} Cd34_pos <- subset(seu, subset = Cd34 > 0) Cd34_Ly6a_pos <- subset(Cd34_pos, subset = Ly6a > 0) # 130个 ``` ```{r} seu$lab <- "neg_Cd34_Ly6a" seu@meta.data[rownames(seu@meta.data) %in% e2,]$lab <- "pos_Cd34_Ly6a" table(seu$lab) ``` ```{r} seu$lab <- factor(seu$lab,levels = rev(c("pos_Cd34_Ly6a","neg_Cd34_Ly6a"))) pdf("../02results/01Dotplot_Cd34_Ly6a.pdf",width = 4.5,height = 3.5) DotPlot(seu,features = c("Cd34","Ly6a"),group.by = "lab",cols = c("lightgrey","red"),dot.scale = 7)+theme(axis.text.x = element_text(angle = 45 , hjust = 1)) dev.off() ``` ```{r} myobj = seu Idents(myobj) <- "lab" cd34ly6DE_gene <- lapply(unique(myobj$lab),function(i){ t <- FindMarkers(myobj, ident.1 = i, assay = 'RNA',slot = 'counts', only.pos = FALSE, logfc.threshold = 0.25) t$entrez <- mapIds(org.Mm.eg.db, keys = rownames(t), column = "ENTREZID", keytype = "SYMBOL", multiVals = "first") t <- t[!is.na(t$entrez),] t$cell <- i return(t) }) names(cd34ly6DE_gene) <- unique(myobj$lab) cd34ly6DE_gene <- do.call(rbind, cd34ly6DE_gene) cd34ly6DE_gene$gene <- mapIds(org.Mm.eg.db, keys = cd34ly6DE_gene$entrez, column = "SYMBOL", keytype = "ENTREZID", multiVals = "first") cd34ly6DE_gene1 <- cd34ly6DE_gene[cd34ly6DE_gene$p_val < 0.05 & cd34ly6DE_gene$avg_log2FC > 0.5, ] write.csv(cd34ly6DE_gene1,"../01data/cd34ly6DE_gene1_log2FC0.8.csv",row.names = FALSE) write.csv(cd34ly6DE_gene,"../01data/cd34ly6DE_gene.csv",row.names = FALSE) ``` ```{r} # 初始化一个列表来存储富集结果 ego_results1 <- list() # 对每个细胞类型的差异基因进行GO富集分析 unique_cells <- unique(cd34ly6DE_gene1$cell) for (cell in unique_cells) { gene_list <- cd34ly6DE_gene1$entrez[cd34ly6DE_gene1$cell == cell] ego_results1[[cell]] <- enrichGO( gene = gene_list, OrgDb = org.Mm.eg.db, keyType = 'ENTREZID', ont = "BP", pAdjustMethod = "BH", pvalueCutoff = 0.05, qvalueCutoff = 0.05, readable = TRUE ) } ``` ```{r} ego_df1 <- do.call(rbind, lapply(names(ego_results1), function(cell) { if (nrow(ego_results1[[cell]]) > 0) { res <- as.data.frame(ego_results1[[cell]]) res$CellType <- cell return(res) } })) write.csv(ego_df1,"../01data/seu_cd34ly6_GO_log2FC0.8.csv",row.names = FALSE) ``` ```{r} ego_df1 <- read.csv("../01data/03Enrich/04Cd34_Ly6a/ctrseu_cd34ly6_GO_log2FC0.5.csv") ``` ```{r} saveRDS(seu,"../01data/Cd34_Ly6a_seu.rds") ``` ```{r} pos_ego_df <- subset(ego_df1,CellType == "pos_Cd34_Ly6a") hema_pwy <- data.frame(Type = "hema",pwy = c( "stem cell differentiation", "stem cell development", "stem cell proliferation" )) hema_pwy <- hema_pwy %>% left_join(select(pos_ego_df, Description, p.adjust, geneID), by = c("pwy" = "Description")) neo_pwy <- data.frame(Type = "neo",pwy = c( "nerve development", "peripheral nervous system development", "positive regulation of nervous system development", "central nervous system neuron axonogenesis", "neuron projection guidance", "neural crest cell development" )) neo_pwy <- neo_pwy %>% left_join(select(pos_ego_df, Description, p.adjust, geneID), by = c("pwy" = "Description")) tube_pwy <- data.frame(Type = "tube",pwy = c( "epithelial tube morphogenesis", "regulation of tube diameter", "regulation of tube size", "branching morphogenesis of an epithelial tube" )) tube_pwy <- tube_pwy %>% left_join(select(pos_ego_df, Description, p.adjust, geneID), by = c("pwy" = "Description")) ``` ```{r} liver_pwy <- rbind(hema_pwy,neo_pwy,tube_pwy) liver_pwy$`-log10(padj)` <- -log10(liver_pwy$p.adjust) write.csv(liver_pwy,"../01data/Cd34_Ly6a_pwy.csv") ``` ```{r} library(dplyr) library(ggplot2) library(patchwork, lib.loc = "/home/computingserver/anaconda3/envs/R4.3.1/lib/R/library") type_order <- liver_pwy %>% group_by(Type) %>% summarise(max_log10p = max("-log10(padj)")) %>% arrange(max_log10p) %>% pull(Type) type_order <- c("tube","neo","hema") liver_pwy <- liver_pwy %>% mutate(Type = factor(Type, levels = type_order)) %>% group_by(Type) %>% arrange(p.adjust, .by_group = TRUE) %>% ungroup() liver_pwy$pwy <- factor(liver_pwy$pwy, levels = rev(liver_pwy$pwy)) # 因为ggplot是y轴从下往上画的 # liver_pwy$Type <- factor(liver_pwy$Type, levels = rev(type_order)) # 因为ggplot是y轴从下往上画的 # --- 颜色设置 --- type_colors <- c( tube = "#BC70D8", neo = "#F08800", hema = "#0FB20A" ) p_left <- ggplot(liver_pwy, aes(x = 1, y = pwy, fill = Type)) + geom_tile() + scale_fill_manual(values = type_colors) + theme_void() + theme(legend.position = "none", axis.text.y = element_text(size = 13,hjust = 1)) p_right <- ggplot(liver_pwy, aes(x = `-log10(padj)`, y = pwy, fill = Type)) + geom_bar(stat = "identity") + scale_fill_manual(values = type_colors) + labs(x = expression(-log[10]~adjusted~p), y = "") + theme_linedraw(base_size = 18) + theme( legend.position = "bottom", plot.margin = margin(0,0,0,0,unit = "cm"), axis.title.y = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank(), panel.border = element_rect(linewidth = 1.5), panel.grid = element_blank(), ) combined_plot <- p_left + p_right + plot_layout(widths = c(0.05, 1)) ggsave("../02results/ranked_enrichment_barplot.pdf", combined_plot, width = 7, height = 7) ``` ```{r} neo <- c("Nrg1","Ednrb","Efnb1","Adgrg6","Flna","Wnt5A","Lmna","Lnc2","Igf1","Synpo") blo <- c("Ly6a","Mecom","Tcf15","Gata6","DGCR8","Ace") bile <- c("Hes1","Abcg2","Jcad","Dll4","Lgals1") genes <- c(neo,blo,bile) ``` ```{r} cd34ly6DE_gene$Regulate <- as.factor(ifelse(cd34ly6DE_gene$p_val < 0.05 & abs(cd34ly6DE_gene$avg_log2FC) >= 0.5,ifelse(cd34ly6DE_gene$avg_log2FC >= 0.5,'Up','Down'),'Stable')) table(cd34ly6DE_gene$Regulate) # Down Stable Up # 291 10952 291 ``` ```{r} cd34ly6DE_gene$Regulate <- factor(cd34ly6DE_gene$Regulate,levels = c("Up","Stable","Down")) cd34ly6DE_gene$logP <- -log10(cd34ly6DE_gene$p_val) cd34ly6DE_gene$logP[cd34ly6DE_gene$logP == "Inf"] <- max(cd34ly6DE_gene$logP[!cd34ly6DE_gene$logP == "Inf"])+10 # cd34ly6DE_gene <- cd34ly6DE_gene[!cd34ly6DE_gene$logP >= 70,] Up_nums <- nrow(cd34ly6DE_gene[cd34ly6DE_gene$Regulate == "Up",]) Down_nums <- nrow(cd34ly6DE_gene[cd34ly6DE_gene$Regulate == "Down",]) Stable_nums <- nrow(cd34ly6DE_gene[cd34ly6DE_gene$Regulate == "Stable",]) ``` ```{r} p <- ggplot(data = cd34ly6DE_gene, aes(x = avg_log2FC, y = logP, colour =Regulate)) + geom_point(shape=19, size=1.5,alpha=0.7, ) + scale_color_manual(values=c("Down"="#4389CA", "Stable" = "grey", "Up" = "#C73434"))+ geom_vline(xintercept=c(-0.5,0.5),lty=4,col="black",lwd=0.8) + geom_hline(yintercept = -log10(0.05),lty=4,col="black",lwd=0.8) + theme_linedraw()+ theme(panel.grid = element_blank(), text = element_text(size = 12,colour = "black"), axis.text = element_text(size = 15,colour = "black"), axis.title = element_text(size = 15,colour = "black"), legend.text = element_text(size = 15,colour = "black"), legend.title = element_text(size = 15,colour = "black"), panel.border = element_rect(linewidth = 1,colour = "black"), aspect.ratio = 1 )+ labs(y = "-log10(pval)" # subtitle = paste("Total genes: Up =", Up_nums , "Down =", Down_nums, "Stable =", Stable_nums) )+ scale_x_continuous(breaks = c(-10, -5, 0, 5, 10)) # annotate("text", x = -1.4, y = max(cd34ly6DE_gene$logP) + 3, label = "x= -0.5", size = 5, colour = "black") + # annotate("text", x = 1.4, y = max(cd34ly6DE_gene$logP) + 3, label = "x = 0.5", size = 5, colour = "black") # p ``` ```{r} for_label <- cd34ly6DE_gene %>% dplyr::filter(gene %in% ves & Regulate != "Down") p1 <- p + geom_point(data = for_label,size = 1.5, shape = 1, color = "black") + # stroke = 0.5, ggrepel::geom_label_repel( data = for_label, aes(label = gene), label.size = NA, fill = NA, color="black", # box.padding = 0.3, # 增加框和点之间的距离 # point.padding = 1, # 增加点与标签之间的填充距离 max.overlaps = 10, # 增加允许的重叠数量 nudge_y = 1.5, # 可以用来在 y 方向稍微移动标签 nudge_x = 0.8 # 可以用来在 x 方向稍微移动标签 )+ theme(plot.title = element_text(hjust = 0.5,size = 18))+ labs(title = "Volcano") pdf("../02results/05Volcano/02update_20250617/04ves_volcano.pdf",height = 4,width = 5) p1 dev.off() ``` ```{r} for_label <- cd34ly6DE_gene %>% dplyr::filter(gene %in% genes & Regulate != "Down") %>% mutate(category = case_when( gene %in% neo ~ "Neuro", gene %in% blo ~ "Blood", gene %in% bile ~ "Bile", )) my_colors <- c( "Neuro" = "#F08800", "Blood" = "#0FB20A", "Bile" = "#BC70D8", "Vessel" = "#54CCB6", "Down"="#4389CA", "Stable" = "grey", "Up" = "#C73434" ) p1 <- p + geom_point(data = for_label,aes(color = category),size = 1.5, shape = 1) + # stroke = 0.5, ggrepel::geom_label_repel( data = for_label, aes(label = gene,color = category), label.size = NA, fill = NA, max.overlaps = 40, size = 3, show.legend = FALSE ) + scale_color_manual(values = my_colors) + # theme_minimal(base_size = 13) + theme(plot.title = element_text(hjust = 0.5, size = 16)) + labs(title = "Volcano Plot") pdf("../02results/all_gene_volcano.pdf", height = 4, width = 5) print(p1) dev.off() ```