---
title: "mouse liver"
author: "zhouruihan"
date: "2024-12-06"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


```{r}
library(Seurat)
library(SeuratData)
library(Matrix)
library(RSpectra)
library(ggplot2)
library(magrittr)
library(dplyr)
library(tidyr)
library(patchwork)
library(org.Mm.eg.db)
library(enrichplot)
library(clusterProfiler)
library(rrvgo)
library(tibble)
library(ComplexHeatmap)
library(reshape2)
library(ggrepel)
```


```{r}
setwd(dir = "/data/share/zhouruihan/01proj/14Liver/01Mouse/02ScRNA/01liver_endothelial_CMGH_PRJNA615641/")
```


文献：CMGH 《Single-Cell Transcriptomics Reveals Zone-Specific Alterations of Liver Sinusoidal Endothelial Cells in Cirrhosis》

```{r}
ctr_count <- read.csv("../00rawdata/GSM4434565_control_count.csv",row.names = 1)
```


1. 过滤掉至少需要在 3 个细胞中表达的基因
2. 过滤掉表达少于 200 个基因的细胞
3. 计算线粒体基因的比例并过滤掉高线粒体基因表达的细胞
4. 过滤掉不表达 GFP 的细胞
#1.create seu
```{r}
ctrseu <- CreateSeuratObject(
  counts = ctr_count, 
  project = "Liver Endothelial", 
  min.cells = 3,         
  min.features = 200    
)

```


```{r}
ctrseu[["pct.mt"]] <- PercentageFeatureSet(ctrseu, pattern = "^Mt-")
```

#2.QC
```{r}
VlnPlot(ctrseu, 
        fill.by = "feature", # "feature", "ident"
        features = c("nFeature_RNA","nCount_RNA","pct.mt"),
        ncol = 4, pt.size = 0.1) +
  theme(plot.title = element_text(size=12))
```

#3.fillter
```{r}
# 根据 GFP 表达筛选阳性细胞
rownames(ctrseu)[grep("GFP", rownames(ctrseu))]

# 筛选 GFP 基因的表达大于 0 的细胞
gfp_gene <- grep("GFP", rownames(ctrseu), value = TRUE)
ctrseu <- ctrseu[, which(ctrseu@assays$RNA@counts[gfp_gene, ] > 0)]
```

#4.process seu
```{r}
ctrseu <- NormalizeData(ctrseu)
ctrseu <- FindVariableFeatures(ctrseu, 
                     selection.method = "vst",
                     nfeatures = 2000,     
                     verbose = FALSE)
ctrseu <- ScaleData(ctrseu)
ctrseu <- RunPCA(ctrseu)

```

```{r}
ElbowPlot(ctrseu)
```


#5.reduction
```{r}
n.pcs <- 20
res.used <- 0.3

ctrseu <- FindNeighbors(ctrseu, dims = 1:n.pcs, verbose = F)
ctrseu <- FindClusters(ctrseu, resolution = res.used, verbose = F)

ctrseu <- RunUMAP(ctrseu, dims = 1:n.pcs, reduction = "pca", verbose = F)
ctrseu <- RunTSNE(ctrseu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```


```{r}
sort(table(ctrseu$RNA_snn_res.0.3))
```

#6.plot
##6.1.UMAP
```{r}
p1 <- DimPlot(ctrseu, reduction = 'umap',label = TRUE,raster=FALSE,repel = TRUE,pt.size =0.05,label.size = 8)+theme(plot.title = element_text(hjust = 0.5),text = element_text(size = 25),axis.text = element_text(size = 25))+labs(title = "UMAP")
pdf("../02results/01reduction/01ctr_filter_UMAP.pdf",width = 5,height = 5)
p1
dev.off()
```

##6.2.Dotplot
```{r fig.height=4, fig.width=5}
pdf("../02results/03Dotpot/fitter_Dotplot.pdf",width = 4,height = 4)
DotPlot(ctrseu,features = c("ly6a","Cd36"),cols = c("grey","red"))
dev.off()
```





##6.3.statistics
```{r fig.height=3, fig.width=6}
table <- data.frame(table(ctrseu$seurat_clusters)) %>% 
  rename(cluster = Var1)
p <- ggplot(table,aes(cluster,Freq,fill = cluster))+
        geom_col(position = "dodge")+
        theme_linedraw()+
        theme(legend.position = "none",
              panel.grid = element_blank(),
              text = element_text(size = 20,colour = "black")
          )+
  geom_text(aes(label = Freq),size = 5)
p 
```




```{r}
saveRDS(ctrseu,"../01data/01seu/ctr_seu.rds")
```


##6.4.细胞注释
```{r}
ctrseu <- readRDS("../01data/01seu/01ctr_seu.rds")
```

```{r}

n.pcs <- 20
res.used <- 0.5

ctrseu <- FindNeighbors(ctrseu, dims = 1:n.pcs, verbose = F)
ctrseu <- FindClusters(ctrseu, resolution = res.used, verbose = F)

ctrseu <- RunUMAP(ctrseu, dims = 1:n.pcs, reduction = "pca", verbose = F)
ctrseu <- RunTSNE(ctrseu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```



```{r}
table(ctrseu$RNA_snn_res.0.5)
```



##6.5.DEGs
```{r}
myobj = ctrseu
Idents(myobj) <- "RNA_snn_res.0.5"
seuDE_gene <- lapply(unique(myobj$RNA_snn_res.0.5),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(seuDE_gene) <- unique(myobj$RNA_snn_res.0.5)
seuDE_gene <- do.call(rbind, seuDE_gene)

seuDE_gene$gene <- mapIds(org.Mm.eg.db,
        keys = seuDE_gene$entrez,
        column = "SYMBOL",
        keytype = "ENTREZID",
        multiVals = "first")

seuDE_gene1 <- seuDE_gene[seuDE_gene$p_val_adj < 0.05 & seuDE_gene$avg_log2FC > 0.5, ]
write.csv(seuDE_gene1,"../01data/02DEGs/08cluster_DEGs/sign_seuDE_gene.csv",row.names = FALSE)
```



```{r}
p1 <- DimPlot(ctrseu, reduction = 'umap',label = TRUE,raster=FALSE,repel = TRUE,pt.size =0.05,label.size = 8)+theme(plot.title = element_text(hjust = 0.5),text = element_text(size = 25),axis.text = element_text(size = 25))+labs(title = "UMAP")
pdf("../02results/01Reduction/03res0.5_ctr_UMAP.pdf",width = 5,height = 5)
p1
dev.off()
```




##6.6.大类注释
```{r}
markers <- c("Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
            "Lyve1","Mrc1","Stab2","Fcgr2b",  # 1 (EC3)Midzonal 中间区域
            "Wnt9b","Lhx6","Thbd","Ramp3",   # 2 (EC5)中央静脉
            "Alb","Ttr","Apoa1",            # 3 肝细胞
            "Cd68","C1qa","C1qb",        # 4 8 巨噬
            "Wnt2","Kit","Lgals1","Plpp1", # 5 (EC4)近中央静脉
            "Colec11","Dcn",       # 6 HSC
            "Nkg7","Trbc1","Cxcr6",       # 7 T cell
            "Gja5","Gmkir1","Sdc1","Ly6a","Cd34", # 9 (EC1)门静脉
            "Spp1","Krt8","Ktr18"       # 10 胆管细胞
            )
```


```{r}
ctrseu1@meta.data$celltype <- plyr::mapvalues(
                from = c(0,
                         1,
                         2,
                         3,
                         4,8,
                         5,
                         6,
                         7,
                         9,
                         10
                
                ),
                to = c(rep("Periportal",1),
                       rep("Midzonal",1),
                       rep("Central_vein",1),
                       rep("Hepatocyte",1),
                       rep("Macrophage",2),
                       rep("Pericentral",1),
                       rep("HSC",1),
                       rep("T_cell",1),
                       rep("Portal",1),
                       rep("cholangiocyte",1)
                ),
                x = ctrseu1@meta.data$RNA_snn_res.0.5)
table(ctrseu1$celltype)
```
```{r}
# ctrseu1 <- ctrseu[,!ctrseu$RNA_snn_res.0.5 %in% c(11)]
```



```{r}
pdf("../02results/03Dotplot/tmp.pdf",width = 10,height = 4)
DotPlot(ctrseu1,features = markers,cols = c("grey","red"),group.by = "celltype")+theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust = 1))
dev.off()
```

##6.7.Periportal小类分群
```{r}
Periportal <- ctrseu1[,ctrseu1$celltype %in% c("Periportal")]
```


```{r}
n.pcs <- 20
res.used <- 0.8

Periportal <- FindNeighbors(Periportal, dims = 1:n.pcs, verbose = F)
Periportal <- FindClusters(Periportal, resolution = res.used, verbose = F)

Periportal <- RunUMAP(Periportal, dims = 1:n.pcs, reduction = "pca", verbose = F)
# Periportal <- RunTSNE(Periportal, dims = 1:n.pcs, reduction = "pca", verbose = F)

```


```{r}
table(Periportal$seurat_clusters)
```

```{r}
p1 <- DimPlot(Periportal, reduction = 'umap',label = TRUE,raster=FALSE,repel = TRUE,pt.size =0.05,label.size = 8)+theme(plot.title = element_text(hjust = 0.5),text = element_text(size = 25),axis.text = element_text(size = 25))+labs(title = "UMAP")
pdf("../02results/01Reduction/04Periportal_UMAP.pdf",width = 5,height = 5)
p1
dev.off()
```



```{r fig.height=4, fig.width=5}
pdf("../02results/03Dotplot//Periportal_Dotplot.pdf",width = 5,height = 5)
DotPlot(Periportal,features = c("Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
                                "Lyve1","Mrc1","Stab2","Fcgr2b"),  # 1 (EC3)Midzonal 中间区域),
                   cols = c("grey","red"))
dev.off()
```


```{r}
Periportal@meta.data$celltype <- plyr::mapvalues(
    from = c(1,2,5,
             3
    ),
    to = c(rep("Midzonal",3),
           rep("Periportal",1)
          ),
    x = Periportal@meta.data$RNA_snn_res.0.8)
table(Periportal$celltype)
```


```{r}
Periportal1 <- Periportal[,!Periportal$celltype %in% c(0,4)]
```


```{r}

pdf("../02results/03Dotplot/Periportal_Dotplot.pdf",width = 10,height = 4)
DotPlot(Periportal1,features = c("Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
                                "Lyve1","Mrc1","Stab2","Fcgr2b"),  # 1 (EC3)Midzonal 中间区域),
                    cols = c("grey","red"),group.by = "celltype",scale = FALSE)+theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust = 1))
dev.off()
```


```{r}
Periportal$celltype <- as.character(Periportal$celltype)
ctrseu1$celltype <- as.character(ctrseu1$celltype)
ctrseu1@meta.data$celltype[match(rownames(Periportal@meta.data), rownames(ctrseu1@meta.data))] <- Periportal$celltype
sort(table(ctrseu1$celltype))
```

```{r}
ctrseu2 <- ctrseu1[,!ctrseu1$celltype %in% c(0,4)]
```



```{r}
markers <- c(
            "Colec11","Dcn",       # 6 HSC
            "Nkg7","Trbc1","Cxcr6",       # 7 T cell
            "Cd68","C1qa","C1qb",        # 4 8 巨噬
            "Alb","Ttr","Apoa1",            # 3 肝细胞
            "Spp1","Krt8","Ktr18",       # 10 胆管细胞
            "Pecam1","Cdh5",
            "Gja5","Gmkir1","Sdc1", # 9 (EC1)门静脉
            "Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
            "Lyve1","Mrc1","Stab2",  # 1 (EC3)Midzonal 中间区域
            "Wnt2","Kit", # 5 (EC4)近中央静脉
            "Thbd","Lhx6","Wnt9b"   # 2 (EC5)中央静脉

            )
```

```{r}
ctrseu2$celltype <- gsub("Portal","EC_Portal",ctrseu2$celltype)
ctrseu2$celltype <- gsub("Periportal","EC_Periportal",ctrseu2$celltype)
ctrseu2$celltype <- gsub("Midzonal","EC_Midzonal",ctrseu2$celltype)
ctrseu2$celltype <- gsub("Pericentral","EC_Pericentral",ctrseu2$celltype)
ctrseu2$celltype <- gsub("Central_vein","EC_Central_vein",ctrseu2$celltype)
ctrseu2$celltype <- gsub("cholangiocyte","Cholangiocyte",ctrseu2$celltype)
```


```{r}
ctrseu2$celltype <- factor(ctrseu2$celltype,levels = rev(c("HSC","T_cell","Macrophage","Hepatocyte","Cholangiocyte","EC_Portal","EC_Periportal","EC_Midzonal","EC_Pericentral","EC_Central_vein")))
pdf("../02results/03Dotplot/ctrseu2.pdf",width = 10,height = 4)
DotPlot(ctrseu2,features = markers,cols = c("#63B9F2","red"),group.by = "celltype")+
        # scale_color_viridis_c()+
        theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust = 1))
dev.off()
```


```{r}
saveRDS(ctrseu2,"../01data/01seu/04ctrseu2.rds")
```





##6.8.细胞类型UMAP
```{r}
col <- c("HSC" = "#FC7C03",
         "T_cell" = "#F4B114",
         "Macrophage" = "#3DC148",
         "Hepatocyte" = "#78DBDF",
         "Cholangiocyte" = "#5598A6",
         "EC_Portal" = "#E42A2A",
         "EC_Periportal" = "#86B8CF",
         "EC_Midzonal" = "#9A86CF",
         "EC_Pericentral" = "pink",
         "EC_Central_vein" = "#BE569A")
```

```{r}
p <- DimPlot(ctrseu2, reduction = 'umap',group.by = "celltype",cols = col,label = FALSE,label.size = 15,raster=FALSE,repel = TRUE,pt.size = 0.1)+
     guides(color = guide_legend(override.aes = list(size = 6),ncol = 1))+
     theme(plot.title = element_text(color = "black",size = 20,hjust = 0.5),
           text = element_text(color = "black",size = 20),
           axis.text = element_text(color = "black",size = 25),
           axis.title  = element_text(color = "black",size = 25),
           legend.position = "right",
           legend.text = element_text(color = "black",size = 25),
           legend.title = element_text(color = "black",size = 25),
           aspect.ratio = 1,
           plot.margin = margin(t =0,b=0,r=0,l=0,unit = "cm"))+
     labs(title = "Liver EC celltype UMAP")
```


```{r}
pdf("../02results/01Reduction/02ctr_celltype_umap.pdf",width = 7,height = 7)
p
dev.off()
```


##6.9.statistics
```{r fig.height=3.5, fig.width=2}
table <- data.frame(table(ctrseu2$celltype)) %>% 
  rename(cluster = Var1)
table1 <- subset(table,cluster %in% c("EC_Central_vein","EC_Pericentral","EC_Midzonal","EC_Periportal","EC_Portal"))
p <- ggplot(table1,aes(reorder(cluster,-Freq),Freq,fill = cluster))+
        geom_col(position = "dodge",width = 0.8)+
        scale_fill_manual(values  = col)+
        theme_linedraw()+
        theme(legend.position = "none",
              panel.grid = element_blank(),
              text = element_text(size = 15,colour = "black"),
              axis.text.x = element_text(size = 15,angle = 30,hjust = 1,colour = "black"),
              plot.margin = margin(t =0.5,b=0.5,r=0.5,l=1,unit = "cm")
          )+
  geom_text(aes(label = Freq),size = 5)+
  labs(x = "")

```

```{r}
pdf("../02results/06Statistic/02EC_col_statistics.pdf",width = 4,height = 4)
p
dev.off()
```

#7.ECseu

```{r}
ECseu <- ctrseu2[,ctrseu2$celltype %in% c("EC_Central_vein","EC_Pericentral","EC_Midzonal","EC_Periportal","EC_Portal")]
```

```{r}
ECseu <- NormalizeData(ECseu)
ECseu <- FindVariableFeatures(ECseu, 
                     selection.method = "vst",
                     nfeatures = 2000,     
                     verbose = FALSE)
ECseu <- ScaleData(ECseu)
ECseu <- RunPCA(ECseu)

```

```{r}
n.pcs <- 20
res.used <- 0.3

ECseu <- FindNeighbors(ECseu, dims = 1:n.pcs, verbose = F)
ECseu <- FindClusters(ECseu, resolution = res.used, verbose = F)

ECseu <- RunUMAP(ECseu, dims = 1:n.pcs, reduction = "pca", verbose = F)
ECseu <- RunTSNE(ECseu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```


##7.1.细分过度区域
```{r}
EC_Midzonal <- ECseu[,ECseu$celltype %in% c("EC_Midzonal")]
```


```{r}
n.pcs <- 20
res.used <- 1.2

EC_Midzonal <- FindNeighbors(EC_Midzonal, dims = 1:n.pcs, verbose = F)
EC_Midzonal <- FindClusters(EC_Midzonal, resolution = res.used, verbose = F)

EC_Midzonal <- RunUMAP(EC_Midzonal, dims = 1:n.pcs, reduction = "pca", verbose = F)
# EC_Midzonal <- RunTSNE(EC_Midzonal, dims = 1:n.pcs, reduction = "pca", verbose = F)

```


```{r}
table(EC_Midzonal$seurat_clusters)
```


```{r}
p1 <- DimPlot(EC_Midzonal, reduction = 'umap',label = TRUE,raster=FALSE,repel = TRUE,pt.size =0.05,label.size = 8)+theme(plot.title = element_text(hjust = 0.5),text = element_text(size = 25),axis.text = element_text(size = 25))+labs(title = "UMAP")
pdf("../02results/01Reduction/04EC_Midzonal_UMAP.pdf",width = 5,height = 5)
p1
dev.off()
```



```{r fig.height=4, fig.width=5}
pdf("../02results/03Dotplot//EC_Midzonal_Dotplot.pdf",width = 10,height = 5)
DotPlot(EC_Midzonal,features = markers,
                   cols = c("blue","red"))
dev.off()
```


```{r}
EC_Midzonal@meta.data$celltype <- plyr::mapvalues(
    from = c(1,
             3,
             4    ),
    to = c(rep("EC_Pericentral",1),
           rep("EC_Periportal",1),
           rep("EC_Midzonal",1)
          ),
    x = EC_Midzonal@meta.data$RNA_snn_res.1.2)
table(EC_Midzonal$celltype)
```


```{r}
EC_Midzonal$celltype <- as.character(EC_Midzonal$celltype)
ECseu$celltype <- as.character(ECseu$celltype)
ECseu@meta.data$celltype[match(rownames(EC_Midzonal@meta.data), rownames(ECseu@meta.data))] <- EC_Midzonal$celltype
sort(table(ECseu$celltype))
```

```{r}
ECseu1 <- ECseu[,!ECseu$celltype %in% c(0,2,5,6,7)]
```





```{r}
p1 <- DimPlot(ECseu1, reduction = 'umap',group.by = "celltype",label = TRUE,raster=FALSE,repel = TRUE,pt.size =0.05,label.size = 8)+theme(plot.title = element_text(hjust = 0.5),text = element_text(size = 25),axis.text = element_text(size = 25))+labs(title = "UMAP")
pdf("../02results/01Reduction/06EC1_UMAP.pdf",width = 8,height = 5)
p1
dev.off()
```


```{r}
markers <- c(
            "Pecam1","Cdh5",
            "Gja5","Gmkir1","Sdc1", # 9 (EC1)门静脉
            "Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
            "Lyve1","Mrc1","Stab2",  # 1 (EC3)Midzonal 中间区域
            "Wnt2","Kit", # 5 (EC4)近中央静脉
            "Thbd","Lhx6","Wnt9b"   # 2 (EC5)中央静脉

            )
```



```{r}
ECseu1$celltype <- factor(ECseu1$celltype,levels = c("EC_Central_vein","EC_Pericentral","EC_Midzonal","EC_Periportal","EC_Portal"))
pdf("../02results/03Dotplot/ECseu.pdf",width = 10,height = 4)
DotPlot(ECseu1,features = markers,cols = c("grey","red"),group.by = "celltype")+theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust = 1))
dev.off()
```


```{r}
saveRDS(ECseu1,"../01data/01seu/04ctrseu2_EC.rds")
```



##7.2.差异分析
```{r}
myobj = ECseu1
Idents(myobj) <- "celltype"
ECDE_gene <- lapply(unique(myobj$celltype),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(ECDE_gene) <- unique(myobj$celltype)
ECDE_gene <- do.call(rbind, ECDE_gene)

ECDE_gene$gene <- mapIds(org.Mm.eg.db,
        keys = ECDE_gene$entrez,
        column = "SYMBOL",
        keytype = "ENTREZID",
        multiVals = "first")

ECDE_gene1 <- ECDE_gene[ECDE_gene$p_val_adj < 0.05 & ECDE_gene$avg_log2FC > 0.5, ]
write.csv(ECDE_gene1,"../01data/02DEGs/08cluster_DEGs/EC_sign_seuDE_gene.csv",row.names = FALSE)

```



##7.3.韦恩交集
门静脉和近门静脉区域的交集 
差异基因出来的基因太多了，这里选的是高变基因里的差异基因也就是在2000个里选的
所以最后大概全部的900多个，up的330多个，剩下的就是down的，
我们关注的细胞群体大概就是门静脉和近门静脉的交界处，所以他们之间的交集基因就比较重要，选择了两群top100的交集，大概15个，也就包括了Cd34和Ly6a
```{r}
all_de_genes <- unique(ECDE_gene1$gene)
Idents(ECseu1) <- "celltype"
avg_mat <- AverageExpression(ECseu1, assays = "RNA", slot = "scale.data")$RNA   # scale.data
mat <- avg_mat[intersect(all_de_genes, rownames(avg_mat)), ]
```


```{r}
mat1 <- as.data.frame(mat) %>% 
  rownames_to_column(var = "gene")
melt_mat <- reshape2::melt(mat1)
write.csv(melt_mat,"../01data/02DEGs/08cluster_DEGs/var_with_DEGs_sign_ECDE_gene.csv",row.names = FALSE)
```


```{r}
topper10_portal <- melt_mat %>%
  dplyr::filter(variable == "EC-Portal") %>%
  arrange(desc(value)) %>%
  slice_head(n = 100) %>%
  pull(gene)  

topper10_periportal <- melt_mat %>%
  dplyr::filter(variable == "EC-Periportal") %>%
  arrange(desc(value)) %>%
  slice_head(n = 100) %>%
  pull(gene)

inter_genes <- intersect(topper10_portal, topper10_periportal)
inter_genes
```

```{r}
# library(VennDiagram)
# library(grid)

inter_labels <- paste(inter_genes, collapse = "\n")
pdf("../02results/07DEGs/02venn.pdf", width = 6, height = 6)
grid.newpage()
venn.plot <- draw.pairwise.venn(
  area1 = 100,
  area2 = 100,
  cross.area = 50,
  category = c("EC-Portal", "EC-Periportal"),
  fill = c("#E42A2A", "#86B8CF"),
  lwd = 2,
  cat.cex = 1.5,
  cex = 0,              
  ind = TRUE
)
grid.text(inter_labels, x = 0.5, y = 0.5, gp = gpar(fontsize = 10))
dev.off()

```



```{r}
# gene_cluster_map <- ECDE_gene1 %>%
#   group_by(gene) %>%
#   slice_max(avg_log2FC, n = 1, with_ties = FALSE) %>%  # 取 log2FC 最大的记录
#   ungroup() %>%
#   select(gene, cell)  # 只保留基因名和归属 cluster
# gene_cluster_map$cell <- gsub("_","-",gene_cluster_map$cell)
# gene_cluster_map <- gene_cluster_map[gene_cluster_map$gene %in% rownames(mat), ]
# 
# mat <- mat[gene_cluster_map$gene, ]
# 
# row_split <- factor(gene_cluster_map$cell, levels = c("EC-Portal","EC-Periportal","EC-Midzonal","EC-Pericentral","EC-Central-vein"))  # 自定义排序

```



9-门静脉：Sdc1，Adgrg6，Gja5，Vwf，Ly6a，
0-近门静脉：Socs3，Ltbp4，Klf4，Apold1，Msr1
1：中间区域：Efnb1，Esm1，Pcdh17，GM26917，Insr
5：近中央静脉：Wnt2，Kit，Gas6，Gja4，Ltc4s
2：中央静脉：Cdh13，Thbd，Wnt9b，Cd9，Plvap

```{r}
# library(ComplexHeatmap)
# library(circlize)
# 
# highlight_genes <- c("Sdc1","Adgrg6","Gja5","Vwf","Ly6a","Cd34",   # 门静脉
#                      "Socs3","Ltbp4","Klf4","Apold1","Msr1",  # 近门静脉
#                      "Efnb1","Esm1","Pcdh17","GM26917","Insr", # 中央区域
#                      "Wnt2","Kit","Gas6","Gja4","Ltc4s",  # 近中央静脉
#                      "Cdh13","Thbd","Wnt9b","Cd9","Plvap"
#                      
# )  
# gene_labels <- ifelse(rownames(mat) %in% highlight_genes, rownames(mat), "")
# 
# cluster_order <- rev(c("EC-Portal","EC-Periportal","EC-Midzonal","EC-Pericentral","EC-Central_vein"))
# cluster_colors <- c(
#          "EC-Portal" = "#E42A2A",
#          "EC-Periportal" = "#86B8CF",
#          "EC-Midzonal" = "#9A86CF",
#          "EC-Pericentral" = "pink",
#          "EC-Central_vein" = "#BE569A"
# )
# 
# top_anno <- HeatmapAnnotation(
#   ClusterBar = cluster_order,
#   col = list(ClusterBar = cluster_colors),
#   show_annotation_name = FALSE,
#   simple_anno_size = unit(0.4, "cm")
# )
# 
# 
# gene_labels <- ifelse(rownames(mat) %in% highlight_genes, rownames(mat), "")
# mat <- mat[rowSums(abs(mat)) > 0, ]
# 
# highlight_genes <- c("Sdc1","Adgrg6","Gja5","Vwf","Ly6a","Cd34",   # 门静脉
#                      "Socs3","Ltbp4","Klf4","Apold1","Msr1",  # 近门静脉
#                      "Efnb1","Esm1","Pcdh17","GM26917","Insr", # 中央区域
#                      "Wnt2","Kit","Gas6","Gja4","Ltc4s",  # 近中央静脉
#                      "Cdh13","Thbd","Wnt9b","Cd9","Plvap"
#                      
# )  
# mat <- as.matrix(mat)
# Ht <- Heatmap(as.matrix(mat),
#         name = "Expression",
#         row_labels = gene_labels,
#         top_annotation = top_anno,
#         column_names_side = "top",    
#         column_names_rot = 0, 
#         row_split = row_split,
#         cluster_rows = FALSE,
#         # cluster_rows = TRUE,
#         cluster_columns = FALSE,
#         # cluster_columns = TRUE,
#         show_row_names = TRUE,
#         show_row_dend = FALSE,
#         show_column_dend = FALSE,
#         row_names_gp = gpar(fontsize = 8),
#         col = colorRamp2(c(-1, 0, 1), c("#0B57A3", "white", "#D10909")))

```



```{r}
# pdf("../02results/07DEGs_heatmap/Ht_v3.pdf",width = 6,height = 6)
# Ht
# dev.off()
```

门静脉差异基因实在太多了，只选top20
##7.4.差异表达热图
```{r}
# 选top20
top20_genes <- melt_mat %>%
  group_by(variable) %>%
  top_n(n = 20, wt = value) %>%  
  ungroup()
# 细胞类型排序：从门静脉到中央静脉
celltype_order <- c("EC-Portal", "EC-Periportal", "EC-Midzonal", "EC-Pericentral", "EC-Central-vein")
top20_genes$variable <- factor(top20_genes$variable, levels = celltype_order)
# 便于热图的基因排序，不使用自动聚类
gene_ordered <- top20_genes %>%
  arrange(variable, desc(value)) %>%
  distinct(gene, .keep_all = TRUE) %>%  # 去重，保留第一次出现（即最高值归属的 celltype）
  pull(gene)

# 基因块
# 去重并保留基因归属顺序（比如按表达最大 celltype）
gene_celltype <- top20_genes %>%
  arrange(variable, desc(value)) %>%
  distinct(gene, .keep_all = TRUE) %>%
  dplyr::filter(gene %in% gene_ordered) %>%
  mutate(variable = factor(variable, levels = celltype_order)) %>%
  arrange(variable)
gene_split <- gene_celltype$variable
names(gene_split) <- gene_celltype$gene

```


```{r}
# pdf("../02results/07DEGs_heatmap/Ht_v4.pdf",width = 12,height = 20)
# DoHeatmap(ECseu1, features = gene_list, group.by = "celltype") +
#   scale_fill_gradientn(colors = c("blue", "white", "red"))
# dev.off()
```



```{r}
# Idents(ECseu1) <- "celltype"
# avg_expr <- AverageExpression(ECseu1, group.by = "celltype", slot = "data")$RNA
# 
# genes_to_plot <- gene_list
# 
# mat <- avg_expr[genes_to_plot, ]
# mat_scaled <- t(scale(t(mat)))  # 每行Z-score标准化
# 
# # 确定要标注的基因在矩阵中的行位置
# marker_positions <- which(rownames(mat_scaled) %in% genes_to_plot)
# 
# # 创建行注释对象（右侧加线加标签）
# row_anno <- rowAnnotation(marker = anno_mark(
#   at = marker_positions,
#   labels = rownames(mat_scaled)[marker_positions],
#   side = "right",
#   labels_gp = gpar(fontsize = 9)
# ))
# 
# pdf("../02results/07DEGs_heatmap/Ht_v5.pdf",width = 12,height = 20)
# 
# Heatmap(mat_scaled,
#         name = "z-score",
#         col = colorRamp2(c(-2, 0, 2), c("blue", "white", "red")),
#         show_row_names = FALSE,
#         cluster_rows = TRUE,
#         cluster_columns = FALSE,
#         column_names_gp = gpar(fontsize = 10, fontface = "italic"),
#         column_names_rot = 45,
#         right_annotation = row_anno)
# dev.off()
```


```{r}
# library(ComplexHeatmap)
# library(circlize)
# library(grid)
Idents(ECseu1) <- "celltype"
avg_expr <- AverageExpression(ECseu1, group.by = "celltype", slot = "data")$RNA
genes_to_plot <- gene_list
mat <- avg_expr[genes_to_plot, ]
df2 <- mat  
df2_scaled <- t(scale(t(df2)))  # 每行Z-score标准化
df2_scaled <- df2_scaled[gene_ordered, ]
df2_scaled <- df2_scaled[,celltype_order]
# 用线标注的基因
select <- c(
  "Cdh13","Wnt2","Kit","Lhx6","Thbd",                      # 中央静脉
  "Gas6","Lgals1","Rab3b","Ptgs1","Tcim",                  # 近中央静脉
  "Lyve1","Dnm3os","Hgf", "Esm1", "Ltbp4",                 # 过渡区域
  "Msr1","Ntn4","Efnb1","Cd36","Ltbp4","Glul",             # 近门静脉
  "Sdc1","Adgrg6","Gja5","Vwf","Ly6a","Cavin3","Atp13a3"   # 门静脉
)

mat_split <- df2_scaled[names(gene_split), ]

# 提取在表达矩阵中的基因 & 位置
mark_genes <- select[select %in% rownames(df2_scaled)]
mark_indices <- which(rownames(df2_scaled) %in% mark_genes)

# 构建右侧标注对象
ha <- rowAnnotation(
  mark = anno_mark(
    at = mark_indices,
    labels = rownames(df2_scaled)[mark_indices],
    labels_gp = gpar(fontsize = 9)
  )
)

# 与UMAP一致
celltype_colors <- c(
 "EC-Portal" = "#E42A2A",
 "EC-Periportal" = "#86B8CF",
 "EC-Midzonal" = "#9A86CF",
 "EC-Pericentral" = "pink",
 "EC-Central-vein" = "#BE569A"
)

# # 去掉前缀 “EC_” 或 “EC–”
# celltype_labels <- gsub("^EC[_\\-]", "", celltype_order)

# 再用于 topan 注释标签
topan <- HeatmapAnnotation(
  block = anno_block(
    gp = gpar(fill = celltype_colors[celltype_order],col = "white"),
    labels = NULL
  )
)
```


```{r}
pdf("../02results/07DEGs_heatmap/Ht_v5.pdf",width = 4,height = 5)
Heatmap(
  df2_scaled,
  column_title = NULL,
  name = "z-score",
  column_split = factor(colnames(mat_split), levels = celltype_order),
  top_annotation = topan,
  right_annotation = ha,
  row_title = NULL,
  cluster_rows = FALSE,
  show_column_dend = FALSE,
  show_row_dend = FALSE,
  cluster_columns = FALSE,
  show_row_names = FALSE,
  show_column_names = TRUE,
  gap = unit(0.5, "mm"),
  column_gap = unit(0.5, "mm"),
  # column_split = column_groups,
  row_split = gene_split,        
  # cluster_rows = TRUE,           
  column_names_side = "top",
  column_names_rot = 45,
  col = colorRamp2(c(-2, 0, 2), c("#78BEDF", "white", "firebrick3")),
  row_names_gp = gpar(fontsize = 10),
  column_names_gp = gpar(fontsize = 10)
)
dev.off()
```


##7.5.重新画UMAP
因为之前从EC里又删了些细胞，所以需要重新对齐一下数据，重新画全部细胞的Dotplot和UMAP，以及统计图
```{r}
ctrseu2 <- readRDS("../01data/01seu/04ctrseu2.rds")
```

```{r}
Idents(ctrseu2) <- "celltype"
non_EC <- subset(ctrseu2, idents = setdiff(Idents(ctrseu2), c("EC_Central_vein", "EC_Pericentral", "EC_Midzonal", "EC_Periportal", "EC_Portal")))
table(non_EC$celltype)
combined_seu <- merge(non_EC, ECseu1)
Idents(combined_seu) <- combined_seu$celltype
table(combined_seu$celltype)
```


```{r}
combined_seu <- NormalizeData(combined_seu)
combined_seu <- FindVariableFeatures(combined_seu, 
                     selection.method = "vst",
                     nfeatures = 2000,     
                     verbose = FALSE)
combined_seu <- ScaleData(combined_seu)
combined_seu <- RunPCA(combined_seu)

```

```{r}
ElbowPlot(combined_seu)
```



```{r}
n.pcs <- 20
res.used <- 0.3

combined_seu <- FindNeighbors(combined_seu, dims = 1:n.pcs, verbose = F)
combined_seu <- FindClusters(combined_seu, resolution = res.used, verbose = F)

combined_seu <- RunUMAP(combined_seu, dims = 1:n.pcs, reduction = "pca", verbose = F)
combined_seu <- RunTSNE(combined_seu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```




```{r}
col <- c("HSC" = "#FC7C03",
         "T_cell" = "#F4B114",
         "Macrophage" = "#3DC148",
         "Hepatocyte" = "#78DBDF",
         "Cholangiocyte" = "#5598A6",
         "EC_Portal" = "#E42A2A",
         "EC_Periportal" = "#86B8CF",
         "EC_Midzonal" = "#9A86CF",
         "EC_Pericentral" = "pink",
         "EC_Central_vein" = "#BE569A")
```

```{r}
p <- DimPlot(combined_seu, reduction = 'umap',group.by = "celltype",cols = col,label = FALSE,label.size = 15,raster=FALSE,repel = TRUE,pt.size = 0.1)+
     guides(color = guide_legend(override.aes = list(size = 6),ncol = 1))+
     theme(plot.title = element_text(color = "black",size = 20,hjust = 0.5),
           text = element_text(color = "black",size = 20),
           axis.text = element_text(color = "black",size = 25),
           axis.title  = element_text(color = "black",size = 25),
           legend.position = "right",
           legend.text = element_text(color = "black",size = 25),
           legend.title = element_text(color = "black",size = 25),
           aspect.ratio = 1,
           plot.margin = margin(t =0,b=0,r=0,l=0,unit = "cm"))+
     labs(title = "Liver EC celltype UMAP")
```


```{r}
pdf("../02results/01Reduction/02ctr_celltype_umap.pdf",width = 7,height = 7)
p
dev.off()
```



```{r}
markers <- c(
            "Colec11","Dcn",       # 6 HSC
            "Nkg7","Trbc1","Cxcr6",       # 7 T cell
            "Cd68","C1qa","C1qb",        # 4 8 巨噬
            "Alb","Ttr","Apoa1",            # 3 肝细胞
            "Spp1","Krt8","Ktr18",       # 10 胆管细胞
            "Pecam1","Cdh5",
            "Gja5","Gmkir1","Sdc1", # 9 (EC1)门静脉
            "Ntn4","Acer2","Itga9",  # 0 (EC2)近门静脉
            "Lyve1","Mrc1","Stab2",  # 1 (EC3)Midzonal 中间区域
            "Wnt2","Kit", # 5 (EC4)近中央静脉
            "Thbd","Lhx6","Wnt9b"   # 2 (EC5)中央静脉

            )
```



```{r}
combined_seu$celltype <- factor(combined_seu$celltype,levels = rev(c("HSC","T_cell","Macrophage","Hepatocyte","Cholangiocyte","EC_Portal","EC_Periportal","EC_Midzonal","EC_Pericentral","EC_Central_vein")))
pdf("../02results/03Dotplot/ctrseu2.pdf",width = 10,height = 4)
DotPlot(combined_seu,features = markers,cols = c("#63B9F2","red"),group.by = "celltype")+
        # scale_color_viridis_c()+
        theme(axis.text.x = element_text(angle = 90,vjust = 0.5,hjust = 1))
dev.off()
```




```{r fig.height=3.5, fig.width=2}
table <- data.frame(table(combined_seu$celltype)) %>% 
  rename(cluster = Var1)
table1 <- subset(table,cluster %in% c("EC_Central_vein","EC_Pericentral","EC_Midzonal","EC_Periportal","EC_Portal"))
p <- ggplot(table1,aes(reorder(cluster,-Freq),Freq,fill = cluster))+
        geom_col(position = "dodge",width = 0.8)+
        scale_fill_manual(values  = col)+
        theme_linedraw()+
        theme(legend.position = "none",
              panel.grid = element_blank(),
              text = element_text(size = 15,colour = "black"),
              axis.text.x = element_text(size = 15,angle = 30,hjust = 1,colour = "black"),
              plot.margin = margin(t =0.5,b=0.5,r=0.5,l=1,unit = "cm")
          )+
  geom_text(aes(label = Freq),size = 5)+
  labs(x = "")

```


```{r}
pdf("../02results/06Statistic/02EC_col_statistics.pdf",width = 4,height = 4)
p
dev.off()
```


```{r}
saveRDS(combined_seu,"../01data/01seu/04ctrseu2_filter.rds")
```





前面得到了Cd34和Ly6a过后，取双阳的，只在门静脉和近门静脉两群里取Cd34+和Ly6a+的细胞，进行后续的分析
#8.Cd34和Scal-1

```{r}
sub_seu <- ECseu1[,ECseu1$celltype %in% c("EC_Periportal","EC_Portal")]
```

```{r}
# ECseu1$barcode <- rownames(ECseu1@meta.data)
Cd34_pos <- subset(ECseu1, subset = Cd34 > 0)
Cd34_Ly6a_pos <- subset(Cd34_pos, subset = Ly6a > 0)  # 130个
```



```{r}
ECseu1@meta.data <- ECseu1@meta.data %>% 
  mutate(class=case_when(
    ECseu1$barcode %in% Cd34_Ly6a_pos$barcode ~ "Cd34_Ly6a_pos",
    TRUE~"Cd34_Ly6a_neg"
  ))
table(ECseu1$class)
# Cd34_Ly6a_neg Cd34_Ly6a_pos 
#    1718            181
```


##8.1.Dot/Freature
```{r fig.height=3, fig.width=4}
ECseu1$class <- factor(ECseu1$class,levels = rev(c("Cd34_Ly6a_pos","Cd34_Ly6a_neg")))
pdf("../02results/03Dotplot/03Cd34_Ly6a.pdf",width = 4.5,height = 3.5)
DotPlot(ECseu1,features =  c("Cd34","Ly6a"),group.by = "class",cols = c("lightgrey","red"),dot.scale = 7)+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```


```{r}
pdf("../02results/02FeaturePlot/03Cd34_Ly6a.pdf",width = 9,height = 4)
FeaturePlot(ECseu1,features = c("Cd34","Ly6a"),cols = c("grey","red"),pt.size = 0.01,label = TRUE,repel = TRUE)
dev.off()

```

##8.2.umap
```{r}
col <- c("Cd34_Ly6a_neg" = "grey",
         "Cd34_Ly6a_pos" = "#D13131")
```

```{r}
p <- DimPlot(ECseu1, reduction = 'umap',group.by = "celltype",cols = col,label = FALSE,label.size = 15,raster=FALSE,repel = TRUE,pt.size = 0.3)+
     guides(color = guide_legend(override.aes = list(size = 6),ncol = 1))+
     theme(plot.title = element_text(color = "black",size = 20,hjust = 0.5),
           text = element_text(color = "black",size = 20),
           axis.text = element_text(color = "black",size = 25),
           axis.title  = element_text(color = "black",size = 25),
           legend.position = "right",
           legend.text = element_text(color = "black",size = 25),
           legend.title = element_text(color = "black",size = 25),
           aspect.ratio = 1,
           plot.margin = margin(t =0,b=0,r=0,l=0,unit = "cm"))+
     labs(title = "Liver EC celltype UMAP")
```


```{r}
pdf("../02results/01Reduction/02ctr_class_umap.pdf",width = 6,height = 6)
p
dev.off()

pdf("../02results/01Reduction/02EC_celltype_umap.pdf",width = 8,height = 6)
p
dev.off()
```




##8.3.statistics
```{r fig.height=3.5, fig.width=2}
table <- data.frame(table(ECseu1$class)) %>% 
  rename(cluster = Var1)
p <- ggplot(table,aes(cluster,Freq,fill = cluster))+
        geom_col(position = "dodge",width = 0.8)+
        scale_fill_manual(values  = col)+
        theme_linedraw()+
        theme(legend.position = "none",
              panel.grid = element_blank(),
              text = element_text(size = 15,colour = "black"),
              axis.text.x = element_text(size = 15,angle = 30,hjust = 1,colour = "black"),
              plot.margin = margin(t =0.5,b=0.5,r=0.5,l=1,unit = "cm")
          )+
  geom_text(aes(label = Freq),size = 5)
p 
```

```{r}
pdf("../02results/06Statistic/01col_statistics.pdf",width = 2.5,height = 4)
p
dev.off()
```


##8.4.DEGs
```{r}
myobj = ECseu1
Idents(myobj) <- "class"
ly6aDE_gene <- lapply(unique(myobj$class),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(ly6aDE_gene) <- unique(myobj$class)
ly6aDE_gene <- do.call(rbind, ly6aDE_gene)

ly6aDE_gene$gene <- mapIds(org.Mm.eg.db,
        keys = ly6aDE_gene$entrez,
        column = "SYMBOL",
        keytype = "ENTREZID",
        multiVals = "first")

ly6aDE_gene1 <- ly6aDE_gene[ly6aDE_gene$p_val < 0.05 & ly6aDE_gene$avg_log2FC > 0.3, ]
write.csv(ly6aDE_gene1,"../01data/02DEGs/06Ly6a/ly6aDE_gene1_log2FC0.3.csv",row.names = FALSE)
write.csv(ly6aDE_gene,"../01data/02DEGs/06Ly6a/ly6aDE_gene.csv",row.names = FALSE)
```




##8.5.enrich
```{r}
# 初始化一个列表来存储富集结果
ego_results <- list()

# 对每个细胞类型的差异基因进行GO富集分析
unique_cells <- unique(ly6aDE_gene1$cell)
for (cell in unique_cells) {
    gene_list <- ly6aDE_gene1$entrez[ly6aDE_gene1$cell == cell]
    ego_results[[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_df <- do.call(rbind, lapply(names(ego_results), function(cell) {
    if (nrow(ego_results[[cell]]) > 0) {
        res <- as.data.frame(ego_results[[cell]])
        res$CellType <- cell
        return(res)
    }
}))

write.csv(ego_df,"../01data/03Enrich/05Ly6a/ctrseu_ly6a_GO_log2FC0.3.csv",row.names = FALSE)
```


```{r}
pos_ego <- subset(ego_df,CellType == "Cd34_Ly6a_pos")
```



```{r}
saveRDS(ECseu1,"../01data/01seu/04ctrseu2_EC.rds")
```


```{r}
pos_ego_df <- subset(ego_df,CellType == "Cd34_Ly6a_pos")
hema_pwy <- data.frame(Type = "hema",pwy = c(
                "stem cell development",
                "stem cell differentiation",
                "regulation of hematopoietic stem cell proliferation",
                "regulation of stem cell proliferation",
                "hematopoietic stem cell homeostasis",
                "mesenchymal stem cell differentiation"
                ))
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(
                "neuron projection extension",
                "neuron projection guidance",
                "regulation of neuron projection regeneration",
                "positive regulation of neuron projection development",
                "central nervous system neuron differentiation",
                "neuron migration"
                ))

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(
                "hepaticobiliary system development",
                "epithelial tube morphogenesis",
                "branching morphogenesis of an epithelial tube",
                "regulation of tube diameter",
                "regulation of tube size",
                "epithelial tube formation"
                ))
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(pagj)` <- -log10(liver_pwy$p.adjust)
write.csv(liver_pwy,"../01data/03Enrich/04Cd34_Ly6a/plot_pwy.csv")
```


```{r}
my_colors <- c(
  "Neuro"  = "#BF00FF",  
  "Blood"  = "#00C8FF",  
  "Bile"   = "#309427",  
  "Vessel" = "#ff7f0e",   
  "Down"="#4389CA", 
  "Stable" = "grey",
  "Up" = "#C73434"
)
```


##8.6.pwy plot
```{r}
text_size = 25
c = "hema"
pwy_dat <- subset(liver_pwy,Type == c)
p <- ggplot(pwy_dat, aes(x = `-log10(padj)`, y = reorder(pwy,`-log10(padj)`), fill = `-log10(padj)`)) +
        geom_bar(stat = "identity", show.legend = FALSE) +  
        # geom_text(aes(label = Count), hjust = 0.5, size = 15) +  
        labs(x = "-log10(padj)", y = "", title = c," pathway") +
        scale_fill_gradient(low = "#E6D2A3", high = "#0FB20A") +
  # low = "#F6E01E","#B8E0F1"  "#8FC0CA" high = "#2BA866","#C82890"  "#0700D3"
        theme_bw() +
        theme(
          axis.text.y = element_text(color = "black",size = text_size+5),  
          axis.text.x = element_text(color = "black",size = text_size), 
          axis.title.y = element_text(size = text_size),
          axis.title.x = element_text(size = text_size+5),
          plot.title = element_text(hjust = 0.5, size = text_size + 5),
          panel.border = element_rect(size = 3),
          panel.grid = element_blank(),
          plot.margin = margin(1,1,1,1,unit = "cm")
        ) +
        geom_vline(xintercept = 1.30103, color = "black", linetype = "dashed", linewidth = 2) +
        scale_y_discrete(labels = function(i) stringr::str_wrap(i, width = 50))  # 换行

pdf(paste0("../02results/04Ego/04Cd34_Ly6a/",c,"_pwy.pdf"),width = 11,height = 6)
p
dev.off()


```




汇总的通路
```{r}
library(dplyr)
library(ggplot2)
library(patchwork)

type_order <- liver_pwy %>%
  group_by(Type) %>%
  summarise(max_log10p = max("-log10(padj)")) %>%
  arrange(max_log10p) %>%
  pull(Type)

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 = "none",
    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/04Ego/04Cd34_Ly6a/ranked_enrichment_barplot.pdf", combined_plot, width = 7, height = 6)  

```



##8.7.Volcano
```{r}
neo <- c("Epha4","Nrg1","Ednrb","Efnb1","Adgrg6","Flna","Wnt5A","Unc5b")
blo <- c("Ly6a","Mecom","Tcf15","Bap1","Gata6","Pdgfra","DGCR8","Ace")
bile <- c("Hes1","Abcg2","Jcad","Dll4","Heg1")
ves <- c("Gja5","Ednrb","Vwf","Lrg1","Jam2","Ier2")

# lym <- c("Vegfc","Sox18","Itga9")
# dy <- c("Cavin3")
# rebir <- c("Efnb1","IL6st","Id1","Ahnak")
# met <- c("Nrg1","Tsc22d1","Igfbp5","Insr","Timp3","Cd9","Plvap","Acer2","Foxp1")
# fib <- c("Sdc1","Ltbp4","Ece1","Mfge8","Vim","Bmpr2","Fbln2")
# immu <- c("Msr1","Clec1a","Fabp5","Cebpd","Id3")

genes <- c(neo,blo,bile,ves)

```



```{r}
ECseu1$class <- factor(ECseu1$class,levels = c("Cd34_Ly6a_pos","Cd34_Ly6a_neg"))
DotPlot(ECseu1,features = genes,group.by = "class",scale = FALSE)
```


```{r}
ly6aDE_gene$Regulate <- as.factor(ifelse(ly6aDE_gene$p_val < 0.05 & abs(ly6aDE_gene$avg_log2FC) >= 0.4,ifelse(ly6aDE_gene$avg_log2FC >= 0.4,'Up','Down'),'Stable'))
table(ly6aDE_gene$Regulate)
# Down Stable     Up 
#  291   10952    291
```



```{r}
ly6aDE_gene$Regulate <- factor(ly6aDE_gene$Regulate,levels = c("Up","Stable","Down"))
ly6aDE_gene$logP <- -log10(ly6aDE_gene$p_val)
ly6aDE_gene$logP[ly6aDE_gene$logP == "Inf"] <- max(ly6aDE_gene$logP[!ly6aDE_gene$logP == "Inf"])+10
# ly6aDE_gene <- ly6aDE_gene[!ly6aDE_gene$logP >= 70,]
Up_nums <- nrow(ly6aDE_gene[ly6aDE_gene$Regulate == "Up",])
Down_nums <- nrow(ly6aDE_gene[ly6aDE_gene$Regulate == "Down",])
Stable_nums <- nrow(ly6aDE_gene[ly6aDE_gene$Regulate == "Stable",])
```


```{r}
p <- ggplot(data = ly6aDE_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 <- ly6aDE_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 <- ly6aDE_gene %>%
  dplyr::filter(gene %in% genes & Regulate != "Down") %>%
  mutate(category = case_when(
    gene %in% neo  ~ "Neuro",
    gene %in% blo  ~ "Blood",
    gene %in% bile ~ "Bile",
    gene %in% ves  ~ "Vessel"
  ))

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 = 20,
    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/01control/05Volcano/02update_20250617/05all_gene_volcano.pdf", height = 4, width = 5)
print(p1)
dev.off()

```


#9.四氯化碳处理的组

```{r}
cir_count <- read.csv("../00rawdata/GSM4434566_cirrhosis_count.csv",row.names = 1)
```


1. 过滤掉至少需要在 3 个细胞中表达的基因
2. 过滤掉表达少于 200 个基因的细胞
3. 计算线粒体基因的比例并过滤掉高线粒体基因表达的细胞
4. 过滤掉不表达 GFP 的细胞
##9.1.create seu
```{r}
cirseu <- CreateSeuratObject(
  counts = cir_count, 
  project = "Liver Endothelial", 
  min.cells = 3,         
  min.features = 200    
)

```


```{r}
cirseu[["pct.mt"]] <- PercentageFeatureSet(cirseu, pattern = "^Mt-")
```

##9.2.QC
```{r}
pdf("../02results/02cirrhosis/01QC.pdf")
VlnPlot(cirseu, 
        fill.by = "feature", # "feature", "ident"
        features = c("nFeature_RNA","nCount_RNA","pct.mt"),
        ncol = 4, pt.size = 0.1) +
  theme(plot.title = element_text(size=12))
dev.off()
```

##9.3.filter
```{r}
# 根据 GFP 表达筛选阳性细胞
rownames(cirseu)[grep("GFP", rownames(cirseu))]

# 筛选 GFP 基因的表达大于 0 的细胞
gfp_gene <- grep("GFP", rownames(cirseu), value = TRUE)
cirseu <- cirseu[, which(GetAssayData(cirseu, slot = "counts")[gfp_gene, ] > 0)]
```


```{r}
cirseu <- NormalizeData(cirseu)
cirseu <- FindVariableFeatures(cirseu, 
                     selection.method = "vst",
                     nfeatures = 2000,     
                     verbose = FALSE)
cirseu <- ScaleData(cirseu)
cirseu <- RunPCA(cirseu)

```

```{r}

n.pcs <- 20
res.used <- 0.5

cirseu <- FindNeighbors(cirseu, dims = 1:n.pcs, verbose = F)
cirseu <- FindClusters(cirseu, resolution = res.used, verbose = F)

cirseu <- RunUMAP(cirseu, dims = 1:n.pcs, reduction = "pca", verbose = F)
cirseu <- RunTSNE(cirseu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```


##9.4.内皮marker
有一部分GFP标注的，不是内皮细胞，所以还是marker标注一下
```{r fig.height=3, fig.width=4}
pdf("../02results/02cirrhosis/02EC_Dotplot.pdf",width = 4.5,height = 3.5)
DotPlot(cirseu,features =  c("Pecam1","Cdh5"),cols = c("lightgrey","red"),dot.scale = 7)+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```


```{r fig.height=3, fig.width=4}
pdf("../02results/02cirrhosis/02umap.pdf",width = 4.5,height = 3.5)
DimPlot(cirseu,reduction = "umap",label = TRUE)+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```


```{r}
cirseu@meta.data$celltype <- plyr::mapvalues(
    from = c(0,1,2,3,4,5,6,9,11,12,
             7,8,10,13
             
    ),
    to = c(rep("ECs",10),
           rep("non_ECs",4)
          ),
    x = cirseu@meta.data$RNA_snn_res.0.5)
table(cirseu$celltype)
```


```{r fig.height=3, fig.width=4}
cirseu$celltype <- factor(cirseu$celltype,levels = rev(c("ECs","non_ECs")))
pdf("../02results/02cirrhosis/02EC_Dotplot.pdf",width = 4.5,height = 4)
DotPlot(cirseu,features =  c("Pecam1","Cdh5"),group.by = "celltype",cols = c("lightgrey","red"),dot.scale = 7,scale = TRUE)+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```

```{r}
col <- c("ECs" = "#058563",
         "non_ECs" = "#C2B9A0"
        )
```

```{r}
p <- DimPlot(cirseu, reduction = 'umap',group.by = "celltype",cols = col,label = FALSE,label.size = 15,raster=FALSE,repel = TRUE,pt.size = 0.1)+
     guides(color = guide_legend(override.aes = list(size = 6),ncol = 1))+
     theme(plot.title = element_text(color = "black",size = 20,hjust = 0.5),
           text = element_text(color = "black",size = 20),
           axis.text = element_text(color = "black",size = 25),
           axis.title  = element_text(color = "black",size = 25),
           legend.position = "right",
           legend.text = element_text(color = "black",size = 25),
           legend.title = element_text(color = "black",size = 25),
           aspect.ratio = 1,
           plot.margin = margin(t =0,b=0,r=0,l=0,unit = "cm"))+
     labs(title = "Liver cirEC celltype UMAP")
```


```{r}
pdf("../02results/02cirrhosis/03cir_celltype_umap.pdf",width = 7,height = 7)
p
dev.off()
```

```{r}
saveRDS(cirseu,"../01data/02cirrhosis/01cirrhosis_seu.rds")
```


#10.四氯化碳处理的EC
```{r}
ECseu <- cirseu[,cirseu$celltype == "ECs"]
```


```{r}

n.pcs <- 20
res.used <- 0.5

ECseu <- FindNeighbors(ECseu, dims = 1:n.pcs, verbose = F)
ECseu <- FindClusters(ECseu, resolution = res.used, verbose = F)

ECseu <- RunUMAP(ECseu, dims = 1:n.pcs, reduction = "pca", verbose = F)
ECseu <- RunTSNE(ECseu, dims = 1:n.pcs, reduction = "pca", verbose = F)

```



```{r fig.height=3, fig.width=4}
pdf("../02results/02cirrhosis/02cir_ECumap.pdf",width = 4.5,height = 3.5)
DimPlot(ECseu,reduction = "umap",label = TRUE,group.by = )+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```



```{r}
ECseu$barcode <- rownames(ECseu@meta.data)
cir_Cd34_pos <- subset(ECseu, subset = Cd34 > 0)
cir_Cd34_Ly6a_pos <- subset(cir_Cd34_pos, subset = Ly6a > 0)  # 130个
```



```{r}
ECseu@meta.data <- ECseu@meta.data %>% 
  mutate(class=case_when(
    ECseu$barcode %in% cir_Cd34_Ly6a_pos$barcode ~ "Cd34_Ly6a_pos",
    TRUE~"Cd34_Ly6a_neg"
  ))
table(ECseu$class)
# Cd34_Ly6a_neg Cd34_Ly6a_pos 
#    1792           3246 
```


```{r}
saveRDS(ECseu,"../01data/02cirrhosis/01cirrhosis_ECseu.rds")
```


##11.1.Dot/Freature
```{r fig.height=3, fig.width=4}
ECseu$class <- factor(ECseu$class,levels = rev(c("Cd34_Ly6a_pos","Cd34_Ly6a_neg")))
pdf("../02results/02cirrhosis/04Dotplot_Cd34_Ly6a.pdf",width = 4.5,height = 3.5)
DotPlot(ECseu,features =  c("Cd34","Ly6a"),group.by = "class",cols = c("lightgrey","red"),dot.scale = 7)+theme(axis.text.x = element_text(angle = 45 , hjust = 1))
dev.off()
```


```{r}
pdf("../02results/02cirrhosis/05FeaturePlot_03Cd34_Ly6a.pdf",width = 9,height = 4)
FeaturePlot(ECseu,features = c("Cd34","Ly6a"),cols = c("grey","red"),pt.size = 0.01,label = TRUE,repel = TRUE)
dev.off()

```

##11.2.umap
```{r}
col <- c("Cd34_Ly6a_neg" = "grey",
         "Cd34_Ly6a_pos" = "#D13131")
```

```{r}
p <- DimPlot(ECseu, reduction = 'umap',group.by = "class",cols = col,label = FALSE,label.size = 15,raster=FALSE,repel = TRUE,pt.size = 0.1)+
     guides(color = guide_legend(override.aes = list(size = 6),ncol = 1))+
     theme(plot.title = element_text(color = "black",size = 20,hjust = 0.5),
           text = element_text(color = "black",size = 20),
           axis.text = element_text(color = "black",size = 25),
           axis.title  = element_text(color = "black",size = 25),
           legend.position = "right",
           legend.text = element_text(color = "black",size = 25),
           legend.title = element_text(color = "black",size = 25),
           aspect.ratio = 1,
           plot.margin = margin(t =0,b=0,r=0,l=0,unit = "cm"))+
     labs(title = "Liver cirEC celltype UMAP")
```


```{r}
pdf("../02results/02cirrhosis/06cir_class_umap.pdf",width = 8,height = 8)
p
dev.off()

```



```{r}
cir <- data.frame(freq = table(ECseu$class)[2],lab = "cirrhosis")
ctr <- data.frame(freq = table(ECseu1$class)[2],lab = "control")
ECtable <- rbind(ctr,cir)
ECtable$lab <- factor(ECtable$lab,levels = c("control","cirrhosis"))
```

```{r}
p <- ggplot(ECtable,aes(lab,freq,fill = lab))+
        geom_col(position = "dodge",width = 0.8)+
        scale_fill_manual(values  = c("cirrhosis" = "#C061CF",
                                      "control" = "#ECB40C"))+
        theme_linedraw()+
        theme(legend.position = "none",
              panel.grid = element_blank(),
              text = element_text(size = 15,colour = "black"),
              axis.text.x = element_text(size = 15,angle = 30,hjust = 1,colour = "black"),
              plot.margin = margin(t =0.5,b=0.5,r=0.5,l=1,unit = "cm")
          )+
  geom_text(aes(label = freq),size = 5)+
  labs(x = "")

```

```{r}
pdf("../02results/02cirrhosis/07compare_statistics_v2.pdf",width = 2.5,height = 4)
p
dev.off()
```


##11.3.两个组Cd34-Cd36阳性的差异基因
```{r}
ctr1 <- ECseu1[,ECseu1$class == "Cd34_Ly6a_pos"]
cir1 <- ECseu[,ECseu$class == "Cd34_Ly6a_pos"]
```

```{r}
ctr1$group <- "control"
cir1$group <- "cirrhosis"
merged <- merge(ctr1, y = cir1, add.cell.ids = c("ctr", "cir"), project = "Merged_posEC")

merged <- NormalizeData(merged)
merged <- FindVariableFeatures(merged)
merged <- ScaleData(merged)
merged <- RunPCA(merged)

```


```{r}
Idents(merged) <- merged$group
merged <- JoinLayers(merged)
markers <- FindMarkers(merged, ident.1 = "cirrhosis", ident.2 = "control", 
                       logfc.threshold = 0.25, min.pct = 0.1)
markers$gene <- rownames(markers)
sign_markers <- markers[markers$p_val < 0.05 & markers$avg_log2FC > 0.5,]
write.csv(markers,"../01data/02cirrhosis/04all_DEGs.csv",row.names = FALSE)
write.csv(sign_markers,"../01data/02cirrhosis/02sign_DEGs.csv",row.names = FALSE)
```




##11.4Volcano
```{r}
neo <- c("Nes","Msx2","Ednrb","Sox4","Cnp","Dbn1","Lmna","Igf1","Synpo")   # 神经
bile <- c("Hgf","Hoxa2","Tubb3","Ncald","Spp1","Krt8","Krt18","Lgals1")        # 胆管
genes <- c(neo,bile)
```


```{r}
markers$Regulate <- as.factor(ifelse(markers$p_val < 0.05 & abs(markers$avg_log2FC) >= 0.5,ifelse(markers$avg_log2FC >= 0.5,'Up','Down'),'Stable'))
table(markers$Regulate)
# Down Stable     Up 
#  687   3074    721
```



```{r}
markers$Regulate <- factor(markers$Regulate,levels = c("Up","Stable","Down"))
markers$logP <- -log10(markers$p_val)
markers$logP[markers$logP == "Inf"] <- max(markers$logP[!markers$logP == "Inf"])+10
markers <- markers[!markers$logP >= 70,]
markers <- markers[!markers$avg_log2FC < -8,]
Up_nums <- nrow(markers[markers$Regulate == "Up",])
Down_nums <- nrow(markers[markers$Regulate == "Down",])
Stable_nums <- nrow(markers[markers$Regulate == "Stable",])
```


```{r}
p <- ggplot(data = markers, 
            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( -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 <- markers %>%
  dplyr::filter(gene %in% genes & Regulate != "Down") %>%
  mutate(category = case_when(
    gene %in% neo  ~ "Neuro",
    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 = 20,
    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/02cirrhosis/08DEGs_volcano.pdf", height = 4, width = 5)
print(p1)
dev.off()

```


```{r}
for_label <- markers %>% 
  dplyr::filter(gene %in% neo & 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 = 20,   # 增加允许的重叠数量
        nudge_y = 1, # 可以用来在 y 方向稍微移动标签
        nudge_x = 0.8 # 可以用来在 x 方向稍微移动标签
        
      )+
      theme(plot.title = element_text(hjust = 0.5,size = 18))+
      labs(title = "Volcano")
pdf("../02results/02cirrhosis/09neo_volcano.pdf",height = 4,width = 5)
p1
dev.off()
```