seurat

rm(list = ls())
#安装相关的包
# install.packages("BiocManager") 
# BiocManager::install(c('Seurat','dplyr','ggplot2','patchwork',
#                      'magrittr','gtools','Matrix',tidyverse))
if(!require(Seurat))install.packages("Seurat")
if(!require(dplyr))install.packages("dplyr")
if(!require(ggplot2))install.packages("ggplot2")
if(!require(magrittr))install.packages("magrittr")
if(!require(gtools))install.packages("gtools")
if(!require(stringr))install.packages("stringr")
if(!require(Matrix))install.packages("panel.border = element_rect(color = "black", fill = NA),  # 设置外框为黑色线条
Matrix")
if(!require(tidyverse))install.packages("tidyverse")
if(!require(patchwork))install.packages("patchwork")

setwd("D:/scrna-seq/Ppscdata/")
load('WT.rds')
##========================1.多样本整合创建Seurat对象========================
dir = c('WT_1/filtered_feature_bc_matrix', 
        'WT_2/filtered_feature_bc_matrix',
        'WT_3/filtered_feature_bc_matrix')
names(dir) = c('WT1', 'WT2', 'WT3')
counts <- Read10X(data.dir = dir)
PPdata = CreateSeuratObject(counts, project = "zsz", min.cells=3, min.features = 100)
dim(PPdata)
table([email protected]$seurat_clusters)

##========================2.数据质控与标准化================================
#计算细胞中线粒体基因比例
PPdata[["percent.mt"]]<- PercentageFeatureSet(PPdata, pattern = "^M")
##   add   "^M” “mt”小鼠样本  “^MT”人源样本
PPdata[["percent.C"]]<- PercentageFeatureSet(PPdata, pattern = "^C")
## View(PPdata[["percent.C"]])
col.num <- length(levels([email protected]$orig.ident))
violin <- VlnPlot(PPdata,
                  features = c("nFeature_RNA", "percent.mt","percent.C"), 
                  pt.size = 0.1, #不需要显示点，可以设置pt.size = 0
                  ncol = 3)
                  
ggsave("QC/vlnplot-before-qc_all.pdf", plot = violin, width = 15, height = 6) 
ggsave("QC/vlnplot-before-qc_all.png", plot = violin, width = 15, height = 6) 

plot1 <- FeatureScatter(PPdata, feature1 = "nCount_RNA", feature2 = "percent.mt")
plot2 <- FeatureScatter(PPdata, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
pearplot <- CombinePlots(plots = list(plot1, plot2), nrow=1, legend="none") 
ggsave("QC/pearplot-before-qc_all.pdf", plot = pearplot, width = 12, height = 5) 
ggsave("QC/pearplot-before-qc_all.png", plot = pearplot, width = 12, height = 5)
##设置质控标准
PPdata<-subset(PPdata,subset=nFeature_RNA>500& nFeature_RNA<5000 &percent.mt<0.5&percent.C<10)
##add  只筛选掉细胞，不会筛选掉gene
VlnPlot(PPdata,features = "percent.mt")
PPdata<-subset(PPdata,subset=nFeature_RNA>100& nFeature_RNA<5000 &percent.mt<5&percent.C<10)
PPdata<-subset(PPdata,subset=nFeature_RNA>200& nFeature_RNA<2500 &percent.mt<5&percent.C<10)
col.num <- length(levels([email protected]))
## 绘制质量控制后的图
violin <-VlnPlot(PPdata,
                 features = c("nFeature_RNA","percent.mt","percent.C"), 
                 group.by = "orig.ident",
                 pt.size = 0.1, 
                 ncol = 3)
ggsave("QC/vlnplot-after-qc-all.pdf", plot = violin, width = 15, height = 6) 
ggsave("QC/vlnplot-after-qc-all.png", plot = violin, width = 15, height = 6)
## 基因表达量标准化
## 它的作用是让测序数据量不同的细胞的基因表达量具有可比性。计算公式如下：
## 标准化后基因表达量 = log1p（10000*基因counts/细胞总counts）
##细胞标准化
PPdata <- NormalizeData(PPdata, normalization.method = "LogNormalize", scale.factor = 10000)
##========================3.数据降维与聚类==================================
## 寻找高变基因
PPdata <- FindVariableFeatures(PPdata, selection.method = "vst", nfeatures = 2000)
top10 <- head(VariableFeatures(PPdata), 10) 
plot1 <- VariableFeaturePlot(PPdata) 
plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE, size=2.5) 
plot <- CombinePlots(plots = list(plot1, plot2),legend="bottom")
## 横坐标是某基因在所有细胞中的平均表达值，纵坐标是此基因的方差。
## 红色的点是被选中的高变基因，黑色的点是未被选中的基因，变异程度最高的10个基因在如图中标注了基因名称。
ggsave("cluster/VariableFeatures_all.pdf", plot = plot, width = 8, height = 6) 
ggsave("cluster/VariableFeatures_all.png", plot = plot, width = 8, height = 6)

## 在进行PCA降维之前还有要对数据进行中心化（必选）和细胞周期回归分析（可选），
## 它们可以使用ScaleData()函数一步完成。
##每个细胞里基因的标准化
scale.genes <-  rownames(PPdata)
PPdata <- ScaleData(PPdata, features = scale.genes)## 看周期相关基因对细胞聚类的影响
## 如果需要消除细胞周期的影响
#PPdata <- ScaleData(PPdata, vars.to.regress = c("S.Score", "G2M.Score"), features = rownames(PPdata))

## scRNA对象中原始表达矩阵经过标准化和中心化之后，已经产生了三套基因表达数据，可以通过以下命令获得：
# 原始表达矩阵
GetAssayData(PPdata,slot="counts",assay="RNA")  
#标准化之后的表达矩阵              
GetAssayData(scRNA,slot="data",assay="RNA")
#中心化之后的表达矩阵 
GetAssayData(scRNA,slot="scale.data",assay="RNA")

## cc.genes$s.genes=c("Pp3c8-22990","Pp3c1-32870","Pp3c12-18160","Pp3c4-11630","Pp3c16-6040",
##                    "Pp3c15-19470","Pp3c25-2980","Pp3c1-34820","Pp3c19-18510",
##                    "Pp3c23-390","Pp3c22-3180","Pp3c6-14810","Pp3c5-8430","Pp3c14-15380",
##                    "Pp3c24-19790","Pp3c15-21840","Pp3c10-16970","Pp3c3-20850",
##                    "Pp3c3-1164","Pp3c13-1670","Pp3c17-19520","Pp3c16-21800","Pp3c5-14880",
##                   "Pp3c26-14310","Pp3c8-14970","Pp3c24-5960","Pp3c14-680","Pp3c10-2020",
##                   "Pp3c7-4800","Pp3c18-8230","Pp3c27-7450","Pp3c24-12860","Pp3c22-850",
##                   "Pp3c11-22000","Pp3c7-7920","Pp3c25-15290","Pp3c25-4870","Pp3c5-12400",
##                   "Pp3c23-10930","Pp3c6-23160","Pp3c25-1550","Pp3c5-4530",
##                    "Pp3c23-19980", "Pp3c13-15990","Pp3c2-14050","Pp3c8-18140","Pp3c5-22720",
##                   "Pp3c2-1760","Pp3c1-21120","Pp3c1-14663","Pp3c10-650")
##cc.genes$g2m.genes=c("Pp3c16-3910","Pp3c8-2160","Pp3c11-10380","Pp3c7-18430","Pp3c27-6070",
##                     "Pp3c14-12680","Pp3c1-5080","Pp3c1-31370","Pp3c6-12520","Pp3c11-11860",
##                     "Pp3c23-8670","Pp3c8-8040","Pp3c17-11160","Pp3c2-24850","Pp3c19-15920",
##                     "Pp3c20-13190","Pp3c24-16930","Pp3c1-25950","Pp3c24-4900","Pp3c21-13710",
##                     "Pp3c8-17230","Pp3c22-17190","Pp3c1-11980","Pp3c5-10270","Pp3c23-4540",
##                     "Pp3c23-8650","Pp3c23-3390","Pp3c11-11580","Pp3c7-8870", "Pp3c12-6220", 
##                     "Pp3c18-16050","Pp3c20-12130","Pp3c24-5040","Pp3c3-17710","Pp3c4-17450",
##                     "Pp3c26-10570","Pp3c13-14170","Pp3c8-6330","Pp3c21-1630","Pp3c24-8550",
##                     "Pp3c21-3130","Pp3c2-10200","Pp3c21-3140","Pp3c6-320","Pp3c16-20900", 
##                     "Pp3c7-9140","Pp3c11-16930","Pp3c15-23870","Pp3c6-21870","Pp3c3-10820",
##                     "Pp3c22-20430","Pp3c11-22070","Pp3c19-17530","Pp3c25-10280","Pp3c8-15890",
##                     "Pp3c5-10610","Pp3c27-2950","Pp3c2-15020","Pp3c1-19220","Pp3c14-5110",
##                     "Pp3c21-8410","Pp3c9-5880","Pp3c26-14480", "Pp3c7-18190","Pp3c16-8790",
##                     "Pp3c24-2350","Pp3c11-22000","Pp3c25-15290","Pp3c5-12400","Pp3c23-10930","Pp3c7-7920")
## 查看我们选择的高变基因中有哪些细胞周期相关基因
## CaseMatch(c(cc.genes$s.genes,cc.genes$g2m.genes),VariableFeatures(PPdata))

## 细胞周期评分
## g2m_genes = CaseMatch(search = cc.genes$g2m.genes, match = rownames(PPdata))
## s_genes = CaseMatch(search = cc.genes$s.genes, match = rownames(PPdata))
## PPdata <- CellCycleScoring(object=PPdata,  g2m.features=g2m_genes,  s.features=s_genes)

## 以上代码运行之后会在[email protected]中添加S.Score、G2M.Score和Phase三列有关细胞周期的信息。
## [email protected]

## 查看细胞周期基因对细胞聚类的影响
## PPdata <- RunPCA(PPdata, features = c(s_genes, g2m_genes))
## p <- DimPlot(PPdata, reduction = "pca", group.by = "Phase")
## ggsave("cluster/cellcycle_pca_all.png", p, width = 8, height = 6)
## 影响不大，继续分析
##凸显出我们想要的细胞
DimPlot(integrated, label = T, group.by = "Treat", 
        cells.highlight = WhichCells(integrated, 
                                     idents = c("group1_untreated", "group1_treated")), 
        cols.highlight = c("darkblue", "darkred"), cols = "grey")

## PCA降维并提取主成分
PPdata <- RunPCA(PPdata, features = VariableFeatures(PPdata)) 
##查看分组信息
print(pbmc[["pca"]], dims = 1:5, nfeatures = 5)
VizDimLoadings(PPdata, dims = 1:2, reduction = "pca")
plot1 <- DimPlot(PPdata, reduction = "pca", group.by="orig.ident") 
plot2 <- ElbowPlot(PPdata, ndims=30, reduction="pca") 
plotc <- plot1+plot2
ggsave("cluster/pca-all.pdf", plot = plotc, width = 8, height = 4) 
ggsave("cluster/pca-all.png", plot = plotc, width = 8, height = 4)
## 后续分析要根据Elbow选择提取的pc轴数量，一般选择斜率平滑的点之前的所有pc轴，此图我的建议是选择前20个pc轴。
#查看决定pc值的top10基因在500个细胞中的热图，查看pc1-pc9轴
pdf("cluster/PC_heatmap-all.pdf",height = 9,width = 18)
DimHeatmap(PPdata, dims = 1:20, nfeatures=10, cells = 500, balanced = TRUE)
dev.off()
png("cluster/PC_heatmap-all.png")
DimHeatmap(PPdata, dims = 1:9, nfeatures=10, cells = 500, balanced = TRUE)
dev.off()
##add    dim值
PPdata <- JackStraw(PPdata,num.replicate = 100)
PPdata <- ScoreJackStraw(PPdata,dims = 1:20)
JackStrawPlot(PPdata,dims = 1:20)
ElbowPlot(PPdata)
## 细胞聚类
## 此步利用 细胞-PC值 矩阵计算细胞之间的距离，
## 然后利用距离矩阵来聚类。其中有两个参数需要人工选择，
## 第一个是FindNeighbors()函数中的dims参数，需要指定哪些pc轴用于分析，选择依据是之前介绍的cluster/pca.png文件中的右图。
## 第二个是FindClusters()函数中的resolution参数，需要指定0.1-0.9之间的一个数值，用于决定clusters的相对数量，数值越大cluters越多。
PPdata <- FindNeighbors(PPdata, dims = 1:20) 
PPdata <- FindClusters(PPdata, resolution = 0.5)
##add   umap
PPdata <- RunUMAP(PPdata,dims = 1:20)
DimPlot(PPdata,reduction = "umap")
DimPlot(PPdata,reduction = "tsne")
##add find maker
cluster5.maker <- FindMarkers(PPdata,ident.1 = 5,ident.2 = c(0,3),min.pct = 0.25)
PPdata.makers <- FindAllMarkers(PPdata,only.pos = TRUE,min.pct = 0.25,logfc.threshold = 0.25)
##add pct.1:在当前cluster细胞中检测到该基因表达的细胞比例
##add pct.2:在其它cluster细胞中检测到该基因表达的细胞比例
##avg_logFC:两组间平均logFC,Seurat v4默认log2。正值表示特征在第一组中表达得更高
PPdata.makers %>% group_by(cluster) %>% top_n(n=10,wt=avg_logFC)

metadata <- [email protected]
cell_cluster <- data.frame(cell_ID=rownames(metadata), cluster_ID=metadata$seurat_clusters)
write.csv(cell_cluster,'cluster/cell_cluster_all.csv',row.names = F)

## 非线性降维
## tsne
PPdata <- RunTSNE(PPdata, dims = 1:20)
embed_tsne <- Embeddings(PPdata, 'tsne')
write.csv(embed_tsne,'cluster/embed_tsne_all.csv')
## 根据处理显示
plot = DimPlot(PPdata, reduction = "tsne" ,group.by = "orig.ident",pt.size = 0.5,label.size = 4)
ggsave("cluster/tSNE_group_all.pdf", plot = plot, width = 8, height = 7)
ggsave("cluster/tSNE_group_all.png", plot = plot, width = 8, height = 7)
## 根据聚类显示
plot1 = DimPlot(PPdata, reduction = "tsne" ,label = "T", pt.size = 0.5,label.size = 4)
ggsave("cluster/tSNE_cluster_all.pdf", plot = plot1, width = 8, height = 7)
ggsave("cluster/tSNE_cluster_all.png", plot = plot1, width = 8, height = 7)
## UMAP
PPdata <- RunUMAP(PPdata,n.neighbors = 30,dims = 1:20)
embed_umap <- Embeddings(PPdata, 'umap')
write.csv(embed_umap,'cluster/embed_umap_all.csv')
## 根据处理显示
plot = DimPlot(PPdata, reduction = "umap",group.by = "orig.ident",pt.size = 0.5,label.size = 4) 
ggsave("cluster/UMAP_group_all.pdf", plot = plot, width = 8, height = 7)
ggsave("cluster/UMAP_group_all.png", plot = plot, width = 8, height = 7)
## 根据聚类显示
plot2 = DimPlot(PPdata, reduction = "umap",label = "T", pt.size = 0.5,label.size = 4) 
ggsave("cluster/UMAP_cluster_all.pdf", plot = plot2, width = 8, height = 7)
ggsave("cluster/UMAP_cluster_all.png", plot = plot2, width = 8, height = 7)
##修改细胞颜色
plot3 <- DimPlot(PPdata, label = T, pt.size = 1)+
  NoLegend()+labs(x = "UMAP1", y = "UMAP2",title = "Celltype")##去掉后显示图注信息 +
  theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"),
        axis.text.y = element_blank(),
        axis.ticks.y = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())

library(ggsci)
library(ggplot2)
#nature版本
plot4 = plot3 + scale_color_npg()
#science版本
plot5 = plot3 +scale_color_aaas()
plot_grid(plot4,plot5)
## 合并tSNE与UMAP
plotc <- plot1+plot2+ plot_layout(guides = 'collect')
ggsave("cluster/tSNE_UMAP_all.pdf", plot = plotc, width = 10, height = 5)
ggsave("cluster/tSNE_UMAP_all.png", plot = plotc, width = 10, height = 5)
load(file="WT.rda")
##提取单基因表达矩阵
expr <- FetchData(object = PPdata8, vars = stem1)
output <- colnames(PPdata8[, which(expr > 0)])
expr[expr=="0"] <- NA
expr1 <- as.data.frame(t(expr))
expr1 <- select(expr1,where(~!all(is.na(.))))
as <- intersect(colnames(expr1),cell_347)
##========================4.细胞类型鉴定=====================================
## Cluster差异基因
## 通过marker基因鉴定细胞类型目前仍是最常用的方法，
## 应用此方法的基础是找到各个cluster的显著高表达的基因，这就需要差异分析提供候选基因列表。
## Seura提供多种差异分析的方法，默认wilcox方法，MASK方法和DESeq2方法也可以留意一下。
dir.create("cell_identify")
diff.wilcox = FindAllMarkers(PPdata,min.pct = 0.25,logfc.threshold = 0.58)
all.markers = diff.wilcox %>% select(gene, everything()) %>% subset(p_val<0.05)
top20<-all.markers %>% group_by(cluster) %>% top_n(n = 20, wt = avg_logFC)
write.csv(all.markers, "cell_identify/diff_genes_wilcox_12.14.csv", row.names = F)
write.csv(top20, "cell_identify/top20_diff_genes_wilcox_12.14.csv", row.names = F)
save(PPdata, file="WTnew.rda")
table([email protected]$seurat_clusters)
top10 = all.markers %>% group_by(cluster) %>% top_n(n = 10, wt = avg_logFC)
##rename cluster
new.cluster.id <- c("A","B","C","D","E","F","G","H","I")
names(new.cluster.id) <- levels(PPdata)
PPdata <- RenameIdents(PPdata, new.cluster.id)
write.csv(top10, "cell_identify/top10_diff_genes_wilcox_12.14.csv", row.names = F)
top10_genes <- read.csv("cell_identify/top10_diff_genes_wilcox_12.14.csv")
top10_genes = CaseMatch(search = as.vector(top10_genes$gene), match = rownames(PPdata)) 
plot1 = DoHeatmap(PPdata, features = top10_genes, group.by = "seurat_clusters", group.bar = T, size = 4)
ggsave("cell_identify/top10_markers_12.14.pdf", plot=plot1, width=8, height=6) 
ggsave("cell_identify/top10_markers_12.14.png", plot=plot1, width=8, height=6)

VlnPlot(PPdata,features = makertop10$gene[1:20],pt.size = 0,ncol = 1)
FeaturePlot(PPdata, features = c("Pp3c1-27440","Pp3c7-24260"))
  ##========================5.makergene查找0======================
pdf('cell_identify/cluster0_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c("Pp3c18-14350","Pp3c24-9690","Pp3c19-16280","Pp3c15-6160"),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c("Pp3c19-5170","Pp3c3-15000","Pp3c5-11880","Pp3c13-23660"),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c("Pp3c2-35200","Pp3c3-29550","Pp3c3-28050","Pp3c12-22320"),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c("Pp3c17-11500","Pp3c7-6560","Pp3c7-22090","Pp3c12-16880"),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c("Pp3c1-3610","Pp3c19-3750","Pp3c12-9550","Pp3c9-23190"),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist0<-c("Pp3c18-14350","Pp3c24-9690","Pp3c19-16280","Pp3c15-6160","Pp3c19-5170","Pp3c3-15000","Pp3c5-11880","Pp3c13-23660",
             "Pp3c2-35200","Pp3c3-29550","Pp3c3-28050","Pp3c12-22320","Pp3c17-11500","Pp3c7-6560","Pp3c7-22090","Pp3c12-16880",
             "Pp3c1-3610","Pp3c19-3750","Pp3c12-9550","Pp3c9-23190")
pdf('cell_identify/cluster0_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist0)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster0")
dev.off()

FeaturePlot(PPdata,features='Pp3c19-6370',cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
##========================6.cluster1=================================
pdf('cell_identify/cluster1_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c4-17500','Pp3c13-14480','Pp3c22-12650','Pp3c2-20930'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-32410','Pp3c22-4630','Pp3c1-36410','Pp3c1-23000'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c18-15840','Pp3c6-3133','Pp3c15-12690','Pp3c1-9290'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c5-18360','Pp3c5-4970','Pp3c16-23160','Pp3c19-20640'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c20-3770','Pp3c21-15570','Pp3c21-9360','Pp3c14-8120'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()

FeaturePlot(PPdata,features=c('Pp3c21-5600'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 1,label.size = 4)
genelist1<-c('Pp3c4-17500','Pp3c13-14480','Pp3c22-12650','Pp3c2-20930',
             'Pp3c2-32410','Pp3c22-4630','Pp3c1-36410','Pp3c1-23000',
             'Pp3c18-15840','Pp3c6-3133','Pp3c15-12690','Pp3c1-9290',
             'Pp3c5-18360','Pp3c5-4970','Pp3c16-23160','Pp3c19-20640',
             'Pp3c20-3770','Pp3c21-15570','Pp3c21-9360','Pp3c14-8120')
pdf('cell_identify/cluster1_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist1)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster1")
dev.off()

##========================7.cluster2===================================
pdf('cell_identify/cluster2_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c6-6610','Pp3c4-21680','Pp3c6-6560','Pp3c14-22890'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c9-2920','Pp3c11-4360','Pp3c20-10830','Pp3c9-1620'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c15-5450','Pp3c6-6545','Pp3c10-25030','Pp3c5-22560'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-23840','Pp3c2-36840','Pp3c7-23740','Pp3c12-8440'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c8-3600','Pp3c1-19610','Pp3c10-1100','Pp3c8-15410'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist2<-c('Pp3c6-6610','Pp3c4-21680','Pp3c6-6560','Pp3c14-22890',
             'Pp3c9-2920','Pp3c11-4360','Pp3c20-10830','Pp3c9-1620',
             'Pp3c15-5450','Pp3c6-6545','Pp3c10-25030','Pp3c5-22560',
             'Pp3c13-23840','Pp3c2-36840','Pp3c7-23740','Pp3c12-8440',
             'Pp3c8-3600','Pp3c1-19610','Pp3c10-1100','Pp3c8-15410')
pdf('cell_identify/cluster2_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist2)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster2")
dev.off()
##========================8.cluster3===================================
pdf('cell_identify/cluster3_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c16-26080','Pp3c3-10140','Pp3c26-5100','Pp3c17-5320'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c1-23970','Pp3c25-6580','Pp3c19-20160','Pp3c23-4180'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c12-6960','Pp3c1-36360','Pp3c3-18320','Pp3c16-21900'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c5-4560','Pp3c14-11840','Pp3c9-4170','Pp3c4-25130'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c12-7050','Pp3c12-12300','Pp3c21-6920','Pp3c12-5390'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist3<-c('Pp3c16-26080','Pp3c3-10140','Pp3c26-5100','Pp3c17-5320',
             'Pp3c1-23970','Pp3c25-6580','Pp3c19-20160','Pp3c23-4180',
             'Pp3c12-6960','Pp3c1-36360','Pp3c3-18320','Pp3c16-21900',
             'Pp3c5-4560','Pp3c14-11840','Pp3c9-4170','Pp3c4-25130',
             'Pp3c12-7050','Pp3c12-12300','Pp3c21-6920','Pp3c12-5390')
pdf('cell_identify/cluster3_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist3)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster3")
dev.off()
##========================9.cluster4===================================
pdf('cell_identify/cluster4_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c19-20840','Pp3c25-3840','Pp3c26-1520','Pp3c6-25270'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c10-4900','Pp3c5-7570','Pp3c14-22870','Pp3c2-12080'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c26-1680','Pp3c7-25080','Pp3c16-19530','Pp3c19-1760'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c14-5560','Pp3c5-960','Pp3c21-10620','Pp3c11-24130'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c4-11030','Pp3c23-16520','Pp3c1-30410','Pp3c24-18690'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist4<-c('Pp3c19-20840','Pp3c25-3840','Pp3c26-1520','Pp3c6-25270',
             'Pp3c10-4900','Pp3c5-7570','Pp3c14-22870','Pp3c2-12080',
             'Pp3c26-1680','Pp3c7-25080','Pp3c16-19530','Pp3c19-1760',
             'Pp3c14-5560','Pp3c5-960','Pp3c21-10620','Pp3c11-24130',
             'Pp3c4-11030','Pp3c23-16520','Pp3c1-30410','Pp3c24-18690')
pdf('cell_identify/cluster4_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist4)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster4")
dev.off()
##========================10.cluster5===================================
pdf('cell_identify/cluster5_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c15-23900','Pp3c24-11870','Pp3c3-26770','Pp3c16-23620'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c20-18130','Pp3c8-230','Pp3c2-20840','Pp3c5-3730'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c9-4170','Pp3c21-13190','Pp3c8-25400','Pp3c4-24420'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-3400','Pp3c21-19080','Pp3c16-7110','Pp3c23-18630'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c11-14810','Pp3c26-10790','Pp3c1-7030','Pp3c8-6770'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist5<-c('Pp3c15-23900','Pp3c24-11870','Pp3c3-26770','Pp3c16-23620',
             'Pp3c20-18130','Pp3c8-230','Pp3c2-20840','Pp3c5-3730',
             'Pp3c9-4170','Pp3c21-13190','Pp3c8-25400','Pp3c4-24420',
             'Pp3c13-3400','Pp3c21-19080','Pp3c16-7110','Pp3c23-18630',
             'Pp3c11-14810','Pp3c26-10790','Pp3c1-7030','Pp3c8-6770')
pdf('cell_identify/cluster5_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist5)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster5")
dev.off()
##========================11.cluster6===================================
pdf('cell_identify/cluster6_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c19-9000','Pp3c16-630','Pp3c14-12270','Pp3c2-30370'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-27350','Pp3c2-29800','Pp3c18-17700','Pp3c18-1550'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-9000','Pp3c1-18940','Pp3c25-13010','Pp3c19-13690'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c4-25770','Pp3c23-11380','Pp3c7-22740','Pp3c21-3730'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c16-11530','Pp3c3-6500','Pp3c17-3060','Pp3c20-1250'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist6<-c('Pp3c19-9000','Pp3c16-630','Pp3c14-12270','Pp3c2-30370',
             'Pp3c2-27350','Pp3c2-29800','Pp3c18-17700','Pp3c18-1550',
             'Pp3c13-9000','Pp3c1-18940','Pp3c25-13010','Pp3c19-13690',
             'Pp3c4-25770','Pp3c23-11380','Pp3c7-22740','Pp3c21-3730',
             'Pp3c16-11530','Pp3c3-6500','Pp3c17-3060','Pp3c20-1250')
pdf('cell_identify/cluster6_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist6)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster6")
dev.off()
##========================12.cluster7===================================
pdf('cell_identify/cluster7_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c20-15801','Pp3c6-1200','Pp3c6-3133','Pp3c14-10980'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c22-12650','Pp3c2-31241','Pp3c19-6839','Pp3c3-4440'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c10-19100','Pp3c19-23300','Pp3c21-2790','Pp3c6-26320'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c22-4640','Pp3c3-26890','Pp3c13-23990','Pp3c11-16051'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c22-1680','Pp3c9-18211','Pp3c8-23460','Pp3c14-1230'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist7<-c('Pp3c20-15801','Pp3c6-1200','Pp3c6-3133','Pp3c14-10980',
             'Pp3c22-12650','Pp3c2-31241','Pp3c19-6839','Pp3c3-4440',
             'Pp3c10-19100','Pp3c19-23300','Pp3c21-2790','Pp3c6-26320',
             'Pp3c22-4640','Pp3c3-26890','Pp3c13-23990','Pp3c11-16051',
             'Pp3c22-1680','Pp3c9-18211','Pp3c8-23460','Pp3c14-1230')
pdf('cell_identify/cluster7_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist7)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster7")
dev.off()
##========================13.cluster8===================================

pdf('cell_identify/cluster8_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c10-3020','Pp3c7-14450','Pp3c22-5610','Pp3c13-5750'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-5930','Pp3c1-30410','Pp3c15-19390','Pp3c7-12850'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c5-4200','Pp3c3-10740','Pp3c2-35930','Pp3c15-5850'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-15980','Pp3c3-18750','Pp3c9-19520','Pp3c19-21160'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c21-4230','Pp3c5-17290','Pp3c6-12510','Pp3c7-25080'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
genelist8<-c('Pp3c10-3020','Pp3c7-14450','Pp3c22-5610','Pp3c13-5750',
             'Pp3c13-5930','Pp3c1-30410','Pp3c15-19390','Pp3c7-12850',
             'Pp3c5-4200','Pp3c3-10740','Pp3c2-35930','Pp3c15-5850',
             'Pp3c13-15980','Pp3c3-18750','Pp3c9-19520','Pp3c19-21160',
             'Pp3c21-4230','Pp3c5-17290','Pp3c6-12510','Pp3c7-25080')
pdf('cell_identify/cluster8_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist8)
p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("cluster8")
dev.off()
##========================16.特异基因==================================
pdf('cell_identify/specific_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c12-12140','Pp3c6-7260','Pp3c2-3730','Pp3c21-13190'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c18-14350','Pp3c24-9690','Pp3c7-16750','Pp3c19-5170'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-30370','Pp3c21-7690','Pp3c17-12450','Pp3c3-14700'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c24-6680','Pp3c21-10030','Pp3c22-17930','Pp3c23-8600'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c26-5100','Pp3c11-5340','Pp3c10-3020','Pp3c7-14450'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c18-1550','Pp3c2-29800'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c7-14450'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 1,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c10-25070'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 1,label.size = 4)
dev.off()
genelist_s<-c('Pp3c19-16280','Pp3c15-6160','Pp3c4-17500','Pp3c13-14480',
              'Pp3c6-6160','Pp3c4-21680','Pp3c16-26080','Pp3c3-10140',
              'Pp3c19-20840','Pp3c25-3840','Pp3c15-23900','Pp3c24-11870',
              'Pp3c19-9000','Pp3c16-630','Pp3c20-15801','Pp3c6-1200',
              'Pp3c10-3020','Pp3c7-14450')
pdf('cell_identify/cluster_s-new_dotplot.pdf',height = 10,width = 12)
p<-DotPlot(PPdata,cols=c("lightyellow","darkblue"),features = genelist_s)
p1<-p+theme(axis.text.x=element_text(angle=45,hjust=1))+xlab("specific")
dev.off()
ggsave("cell_identify/cluster_s-new_dotplot.png", plot = p1, width = 12, height = 10)
##========================16.转录因子==================================
pdf('cell_identify/TF_specific_feature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c6-26890','Pp3c15-15020','Pp3c11-20170','Pp3c11-25660'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-31200','Pp3c4-16270','Pp3c7-1850','Pp3c7-7950'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c9-3300','Pp3c1-38440','Pp3c11-3960','Pp3c14-17330'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c17-11650','Pp3c7-11360','Pp3c2-6080','Pp3c24-2500'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c7-19440','Pp3c3-35530','Pp3c2-31580','Pp3c5-12180'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c7-25600','Pp3c7-25600','Pp3c7-25600','Pp3c7-5280'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-21190','Pp3c15-3180','Pp3c3-26370','Pp3c4-850'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c5-14440','Pp3c1-7740','Pp3c1-7800','Pp3c10-20000'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c10-6570','Pp3c11-14690','Pp3c11-23290','Pp3c12-25330'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-24020','Pp3c16-13260','Pp3c19-6370','Pp3c2-25760'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c20-3770','Pp3c21-16540','Pp3c23-21450','Pp3c27-5490'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c3-31490','Pp3c3-6830','Pp3c3-7020','Pp3c4-23930'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c4-25130','Pp3c4-2660','Pp3c5-820','Pp3c5-880'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c6-28290','Pp3c7-20170','Pp3c7-2300','Pp3c1-41530'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c17-23620','Pp3c19-2940','Pp3c24-7690','Pp3c7-20870'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c10-22480','Pp3c1-2540','Pp3c17-13550','Pp3c2-20930'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c15-13310','Pp3c9-18520','Pp3c12-21760','Pp3c17-18130'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c6-2730','Pp3c7-7720','Pp3c21-16050','Pp3c23-14630'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-24590','Pp3c12-21500','Pp3c23-7530','Pp3c9-15970'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c14-5130','Pp3c4-5410','Pp3c13-20650','Pp3c16-13150'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c20-18130','Pp3c24-8270','Pp3c3-12890','Pp3c3-8880'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c4-20430','Pp3c6-27690','Pp3c8-6080','Pp3c15-5200'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c3-14260','Pp3c9-5650','Pp3c18-8920','Pp3c21-16440'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-15520','Pp3c14-17020','Pp3c2-2510'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-9700','Pp3c4-26880','Pp3c7-24490'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()
pdf('cell_identify/NAC_fature.pdf',height = 10,width = 12)
FeaturePlot(PPdata,features=c('Pp3c17-11420','Pp3c9-3350','Pp3c15-17590','Pp3c2-2370'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c2-1200','Pp3c20-9770','Pp3c23-6860','Pp3c23-15540'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c8-19190','Pp3c1-35660','Pp3c22-130','Pp3c24-8270'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c26-4660','Pp3c26-13950','Pp3c4-20430','Pp3c4-29470'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c4-11490','Pp3c7-6590','Pp3c13-6470','Pp3c13-10800'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c13-20650','Pp3c2-32530','Pp3c19-1580','Pp3c3-12890'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c3-8880','Pp3c3-25610','Pp3c3-2470','Pp3c20-18130'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c23-9520','Pp3c27-7430','Pp3c27-7560','Pp3c27-7510'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c5-740','Pp3c5-630','Pp3c5-570','Pp3c8-6080'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c12-18020','Pp3c12-15190','Pp3c6-28230','Pp3c6-1310'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
FeaturePlot(PPdata,features=c('Pp3c6-27690','Pp3c16-19830','Pp3c16-23260'),cols=c("grey","yellow","red","brown")
            ,reduction = 'umap',pt.size = 0.5,label.size = 4)
dev.off()

##========================17.拟时间分析================================
dir.create("pseudotime")
  library(monocle)
id<-c("0","3")
Cell.sub3 <- subset([email protected],seurat_clusters==id)
scRNAsub <- subset(PPdata, cells=row.names(Cell.sub))
data <- as(as.matrix(scRNAsub@assays$RNA@counts), 'sparseMatrix')
pd <- new('AnnotatedDataFrame', data = [email protected])
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new('AnnotatedDataFrame', data = fData)
mycds <- newCellDataSet(data,
                        phenoData = pd,
                        featureData = fd,
                        expressionFamily = negbinomial.size())
##估计size factor和分散度
mycds <- estimateSizeFactors(mycds)
mycds <- estimateDispersions(mycds, cores=10, relative_expr = TRUE)
##过滤低表达的细胞 optional
mycds <- detectGenes(mycds,min_expr = 0.1) 
disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.1 & dispersion_empirical >= 1 * dispersion_fit)$gene_id
mycds <- setOrderingFilter(mycds, disp.genes)
p3 <- plot_ordering_genes(mycds)
#降维
mycds <- reduceDimension(mycds, max_components = 2, method = 'DDRTree')
#排序
mycds <- orderCells(mycds)
+##顺序相反 
#State轨迹分布图
plot1 <- plot_cell_trajectory(mycds, color_by = "State")
ggsave("pseudotime/State.pdf", plot = plot1, width = 6, height = 5)
ggsave("pseudotime/State.png", plot = plot1, width = 6, height = 5)
##Cluster轨迹分布图
plot2 <- plot_cell_trajectory(mycds, color_by = "seurat_clusters")
##把每一簇都画出来
plot_cell_trajectory(mycds,color_by = "seurat_clusters")+facet_wrap(~seurat_clusters, nrow = 1)
ggsave("pseudotime/Cluster.pdf", plot = plot2, width = 6, height = 5)
ggsave("pseudotime/Cluster.png", plot = plot2, width = 6, height = 5)
##Pseudotime轨迹图
plot3 <- plot_cell_trajectory(mycds, color_by = "Pseudotime")
ggsave("pseudotime/Pseudotime.pdf", plot = plot3, width = 6, height = 5)
ggsave("pseudotime/Pseudotime.png", plot = plot3, width = 6, height = 5)
##合并作图
plotc <- plot1|plot2|plot3
ggsave("pseudotime/Combination.pdf", plot = plotc, width = 10, height = 3.5)
ggsave("pseudotime/Combination.png", plot = plotc, width = 10, height = 3.5)
##基因作图
plot_genes_in_pseudotime(mycds["Pp3c21-15370",],color_by = "seurat_clusters")
plot_genes_in_pseudotime(mycds["Pp3c21-15370",],color_by = "Pseudotime")
plot_genes_in_pseudotime(mycds["Pp3c17-6500",],color_by = "State")
##单个基因
pData(mycds)$HEC2 = log2( exprs(mycds)['Pp3c17-15380',]+1)
p1=plot_cell_trajectory(mycds, color_by = "HEC2") + scale_color_gradient(na.value = "grey",low="blue",high="red")

##保存结果
write.csv(pData(mycds), "pseudotime/pseudotime.csv")

disp_table <- dispersionTable(mycds)
disp.genes <- subset(disp_table, mean_expression >= 0.5&dispersion_empirical >= 1*dispersion_fit)
disp.genes <- as.character(disp.genes$gene_id)
diff_test <- differentialGeneTest(mycds[disp.genes,], cores = 4, 
                                  fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_gene_names <- row.names(subset(diff_test, qval < 1e-04))
p2 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=3,
                             show_rownames=T, return_heatmap=T)
clusters <- cutree(p2$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
write.csv(clustering,"pseudotime/gene_cluster.csv")

##热图
p3 = plot_pseudotime_heatmap(mycds[sig_gene_names,], num_clusters=3,
                             show_rownames=F, return_heatmap=T)
ggsave("pseudotime/pseudotime_heatmap2.png", plot = p3, width = 7, height = 7)

clusters <- cutree(p2$tree_row, k = 3)
clustering <- data.frame(clusters)
clustering[,1] <- as.character(clustering[,1])
colnames(clustering) <- "Gene_Clusters"
write.csv(clustering, "pseudotime/heatmap_gene.csv")

gene<-t(as.matrix(scRNAsub@assays$RNA@counts[c('Pp3c3-6830','Pp3c4-25130',"Pp3c2-32160",
                                               'Pp3c19-6370','Pp3c9-5650','Pp3c20-18130'),
                                             colnames(scRNAsub@assays$RNA@counts)]))
colnames(gene)<-c('Pp3c3_6830','Pp3c4_25130',"Pp3c2_32160",
                  'Pp3c19_6370','Pp3c9_5650','Pp3c20_18130')
mycds@phenoData@data<-cbind(mycds@phenoData@data,gene)
mycds@phenoData@data[mycds@phenoData@data == 0] <- NA
p1<-plot_cell_trajectory(mycds, color_by = "Pp3c3_6830")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
p2<-plot_cell_trajectory(mycds, color_by = "Pp3c4_25130")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
p3<-plot_cell_trajectory(mycds, color_by = "Pp3c19_6370")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
p4<-plot_cell_trajectory(mycds, color_by = "Pp3c9_5650")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
p5<-plot_cell_trajectory(mycds, color_by = "Pp3c20_18130")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
p6<-plot_cell_trajectory(mycds, color_by = "Pp3c2_32160")+scale_color_gradient(na.value = "grey",low="yellow",high="red")
plotc<-p1|p2|p3|p4|p5|p6
ggsave("pseudotime/TF.pdf", plot = plotc, width = 30, height = 5)
ggsave("pseudotime/TF.png", plot = plotc, width = 30, height = 5)