R Seurat 单细胞处理pipline 代码

提示：扫一扫查出行
【扫一扫了解最新限行尾号】
复制提示
 options(stringsAsFactors = F )

 rm(list = ls())

 library(Seurat)

 library(dplyr)

 library(ggplot2)

 library(Hmisc)

 library(pheatmap)

 #读入数据

 #合并gene去batch

 expr_1 <- readRDS("C:/Gu_lab/PA/result/pipline_results/P1_normal/expr.RDS")

 expr_1 <- RenameCells(expr_1,add.cell.id = "P1",for.merge = T )

 expr_8 <- readRDS("C:/Gu_lab/PA/result/pipline_results/P8_normal/expr.RDS")

 expr_8 <- RenameCells(expr_8,add.cell.id = "P8",for.merge = T )

 expr_all_P1_P8_nobatch <- FindIntegrationAnchors(c(expr_1,expr_8))

 expr_all_P1_P8_nobatch <- IntegrateData(anchorset = expr_all_P1_P8_nobatch, dims = :)

 expr_all_P1_P8_nobatch <- FindVariableFeatures(expr_all_P1_P8_nobatch)

 all.genes <- rownames(expr_all_P1_P8_nobatch)

 expr_all_P1_P8_nobatch <- ScaleData(expr_all_P1_P8_nobatch,features = all.genes)

 expr_all_P1_P8_nobatch <- RunPCA(expr_all_P1_P8_nobatch,features = VariableFeatures(object = expr_all_P1_P8_nobatch))

 ElbowPlot(expr_all_P1_P8_nobatch)

 expr_all_P1_P8_nobatch <- FindNeighbors(expr_all_P1_P8_nobatch,dims = :)

 expr_all_P1_P8_nobatch <- FindClusters(expr_all_P1_P8_nobatch,resolution = 0.5)

 expr_all_P1_P8_nobatch <- RunUMAP(expr_all_P1_P8_nobatch, dims = :)

 DimPlot(expr_all_P1_P8_nobatch,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 cell_soruce <-rep(c('P1','P8'), times = c(,))

 newident_1 <- factor(cell_soruce)

 expr_all_P1_P8_nobatch_save <- expr_all_P1_P8_nobatch

 Idents(expr_all_P1_P8_nobatch) <- newident_1

 #expr_nobatch_seurat@active.ident <- newident_1

 DimPlot(expr_all_P1_P8_nobatch,reduction = "umap",cols =c("#F0E442", "#999999", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73","#FF0000","#333333") )

 expr_all_P1_P8_nobatch <- expr_all_P1_P8_nobatch_save

 outdir <- "C:/Gu_lab/PA/result/Pred_result/P1_P8_normal/"

 saveRDS(expr_all_P1_P8_nobatch_save, file =paste0(outdir,"expr_all_P1_P8.RDS"))

 #根据Marker去掉免疫细胞

 pdf(file = paste0(outdir,"Tcell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"),pt.size = )#T-cell Markers

 dev.off()

 #FeaturePlot(expr_all_P1_P8_nobatch,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"))

 pdf(file = paste0(outdir,"Bcell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("PTPRC","CD79A","CD19","CD20","CD45"),pt.size = )#T-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"NK_cell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("PTPRC","NKG7","CD56"),pt.size = )#T-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Macrophage-cell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("CD14","CD163", "CD68", "CSF1R","CD33","HLA-DR","AIF1","FCER1G", "FCGR3A", "TYROBP"),pt.size = )#Macrophage-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Fibroblasts-cell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("ACTA2","FAP", "PDPN","COL1A2","DCN", "COL3A1", "COL6A1","S100A4","COL1A1","THY1"),pt.size = )#Fibroblasts-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Endothelial-cell_feature.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch,features = c("PLVAP","PECAM1","VWF","ENG","MCAM","CD146"),pt.size = )#Endothelial-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Somatotrope_feature.pdf"),width = ,height = )

 Somatotrope_feature =  c("Pappa2","Ceacam10","Tcerg1l","Cabp2","Wnt10a","RP24-294M17.1","Mmp7","Ghrhr","Gh1")

 Somatotrope_feature = toupper(Somatotrope_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Somatotrope_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Somatotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Lactotrope_feature.pdf"),width = ,height = )

 Lactotrope_feature =  c("Prl","Agtr1a","Syndig1","Angpt1","Edil3","6030419C18Rik","Hepacam2","Cilp","Akr1c14")

 Lactotrope_feature = toupper(Lactotrope_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Lactotrope_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Lactotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Corticotrope_feature.pdf"),width = ,height = )

 Corticotrope_feature =  c("Gpc5","Tnnt1","Tnni3","Gm15543","Cplx3","Egr4","Cdh8","Galnt9","Pomc")

 Corticotrope_feature = toupper(Corticotrope_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Corticotrope_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Corticotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Melanotrope_feature.pdf"),width = ,height = )

 Melanotrope_feature = c("Oacyl","Esm1","Pkib","Gulo","Megf11","Rbfox3","Scn2a1","Pax7","Pcsk2")

 Melanotrope_feature = toupper(Melanotrope_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Melanotrope_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Melanotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Thyrotrope_feature.pdf"),width = ,height = )

 Thyrotrope_feature = c("Tshb","Trhr","Dio2")

 Thyrotrope_feature = toupper(Thyrotrope_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Thyrotrope_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Thyrotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Stem_cell_feature.pdf"),width = ,height = )

 Stem_cell_feature =  c("Cyp2f2","Lcn2","Mia","Cpxm2","Rbpms","Gm266","Cdh26","Aqp3","Aqp4","CD8A")

 Stem_cell_feature = toupper(Stem_cell_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Stem_cell_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Stem_cell_feature)

 dev.off()

 pdf(file = paste0(outdir,"Proliferatring_Pou1f1_feature.pdf"),width = ,height = )

 Proliferatring_Pou1f1_feature = c("POU1f1","Pbk","Spc25","Ube2c","Hmmr","Ckap2l","Spc25","Cdca3","Birc5","Ccnb1")

 Proliferatring_Pou1f1_feature = toupper(Proliferatring_Pou1f1_feature)

 VlnPlot(expr_all_P1_P8_nobatch, features =Proliferatring_Pou1f1_feature,pt.size = )

 #FeaturePlot(expr_all_P1_P8_nobatch,features = Proliferatring_Pou1f1_feature)

 dev.off()

 #激素的Gene

 pdf(file = paste0(outdir,"VlnPlot.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch, features =c("TSHB","TBX19","PCSK2","FHSB","GH1", "PRL", "POMC", "CGA", "LHB", "POU1F1","MIA"),pt.size = )

 dev.off()

 #Cell_Type

 expr_all_P1_P8_nobatch <- readRDS(paste0(outdir,"expr_all_P1_P8.RDS"))

 expr <- RenameIdents(expr_all_P1_P8_nobatch,'' = 'T_NK','' = 'Fibroblasts','' = 'Endothelial','' = 'Macrophage','' = 'Somatotrope','' = 'Lactotrope','' = 'Corticotrope','' = 'Gonadotrope','' = 'Stem cell')

 DimPlot(expr,cols = c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333"))

 ##Fig2 DEgene heatmap

 expr_all_P1_P8_nobatch.markers <- FindAllMarkers(expr_all_P1_P8_nobatch)

 top10 <-  expr_all_P1_P8_nobatch.markers %>% group_by(cluster) %>% top_n(n = , wt = avg_logFC)

 DoHeatmap(expr_all_P1_P8_nobatch, features = top10$gene) + NoLegend()

 top5 <-  expr_all_P1_P8_nobatch.markers %>% group_by(cluster) %>% top_n(n = , wt = avg_logFC)

 DoHeatmap(expr_all_P1_P8_nobatch, features = top5$gene) + NoLegend()

 pdf(file = paste0(outdir,"vlnplot_new.pdf"),width = ,height = )

 VlnPlot(expr_all_P1_P8_nobatch, features = top5$gene, pt.size = )

 #for(i in :length(top5$gene)){

   #print(top5$gene[i])

   #VlnPlot(expr_all_P1_P8_nobatch, features = top5$gene, pt.size = )

   #FeaturePlot(expr_all_P1_P8_nobatch,features = top5$gene[i])

 #}

 dev.off()

 #mice cell type, reclus

 mice_expr_init <- readRDS(file = "C:/Gu_lab/PA/data/mouse/mice_exp.rds")

 DimPlot(mice_expr_init,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 #mice_corr_P1&P8

 mice_expr<- readRDS(file = "C:/Gu_lab/PA/data/mouse/mice_exp_final.rds")

 DimPlot(mice_expr,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 #top10 <-  expr_all_P1_P8_nobatch.markers %>% group_by(cluster) %>% top_n(n = , wt = avg_logFC)

 #

 top10 <-  expr_all_P1_P8_nobatch.markers %>% group_by(cluster) %>% top_n(n = , wt = p_val)

 top5 <-  expr_all_P1_P8_nobatch.markers %>% group_by(cluster) %>% top_n(n = , wt = p_val)

 #定义corr矩阵,row_mice,col_human

 featureGene_mice <- mice_expr@assays$RNA@var.features

 featureGene_human <- toupper(featureGene_mice)

 corr_matrix = matrix(nrow = , ncol = )

 mice_ident <- mice_expr@active.ident

 expr_all_P1_P8_nobatch_ident <- expr_all_P1_P8_nobatch@active.ident

 expr_ident <- expr@active.ident

 scale_expr_matrix_mice <- mice_expr@assays$RNA@scale.data

 scale_expr_matrix_human <- expr@assays$integrated@scale.data

 tm_rowname <- rownames(scale_expr_matrix_mice)

 tm_rowname_human<- rownames(scale_expr_matrix_human)

 for( i in :){

   featrue_x = array(rep(,length(featureGene_mice)),dim = length(featureGene_mice))

   x_name = levels(mice_ident)[i]

   expr_mice <- scale_expr_matrix_mice[,which(mice_ident==x_name)]

   for(k in :length(featureGene_mice)){

     id = which(tm_rowname == featureGene_mice[k])

     featrue_x[k] = sum(expr_mice[id,])

   }

   featrue_x = featrue_x/length(colnames(scale_expr_matrix_mice))

   for(j in :){

     featrue_y = array(rep(,length(featureGene_mice)),dim = length(featureGene_mice))

     y_name = levels(expr_ident)[j]

     expr_human <- scale_expr_matrix_human[,which(expr_ident==y_name)]

     for(k in :length(featureGene_human)){

       if(featureGene_human[k] %in% tm_rowname_human){

         id = which(tm_rowname_human == featureGene_human[k])

         featrue_y[k] = sum(expr_human[id,])

       }

     }

     featrue_y = featrue_y/length(colnames(scale_expr_matrix_human))

     print(featrue_y[:])

     corr_matrix[i,j] = cor(featrue_x,featrue_y,method = "spearman")

   }

 }

 rownames(corr_matrix) <- levels(mice_ident)

 colnames(corr_matrix) <- levels(expr_ident)

 #pheatmap(corr_matrix,clustering_distance_rows = "correlation",display_numbers = T,number_format = "%.2f")

 pheatmap(corr_matrix,clustering_distance_rows = "correlation",cluster_rows = T,cluster_cols =T,display_numbers = T,number_format = "%.2f")

 #########################################################################################################

 #老鼠的数据处理见 mouse_data_init.R

 #use human var gene

 corr_matrix = matrix(nrow = , ncol = )

 mice_ident <- mice_expr@active.ident

 expr_ident <- expr@active.ident

 #scale_expr_matrix_mice <- mice_expr@assays$RNA@scale.data

 #scale_expr_matrix_human <- expr@assays$integrated@scale.data

 tm_rowname <- rownames(scale_expr_matrix_mice)

 tm_rowname_human<- rownames(scale_expr_matrix_human)

 featureGene_human <- expr@assays$integrated@var.features

 featureGene_mice <- tolower(featureGene_human)

 featureGene_mice <- capitalize(featureGene_mice)

 for( i in :){

   featrue_x = array(rep(,length(featureGene_mice)),dim = length(featureGene_mice))

   x_name = levels(mice_ident)[i]

   expr_mice <- scale_expr_matrix_mice[,which(mice_ident==x_name)]

   for(k in :length(featureGene_mice)){

     if(featureGene_mice[k] %in% tm_rowname){

       id = which(tm_rowname == featureGene_mice[k])

       featrue_x[k] = sum(expr_mice[id,])

     }

   }

   featrue_x = featrue_x/length(colnames(scale_expr_matrix_mice))

   for(j in :){

     featrue_y = array(rep(,length(featureGene_mice)),dim = length(featureGene_mice))

     y_name = levels(expr_ident)[j]

     expr_human <- scale_expr_matrix_human[,which(expr_ident==y_name)]

     for(k in :length(featureGene_human)){

       if(featureGene_human[k] %in% tm_rowname_human){

         id = which(tm_rowname_human == featureGene_human[k])

         featrue_y[k] = sum(expr_human[id,])

       }

     }

     featrue_y = featrue_y/length(colnames(scale_expr_matrix_human))

     print(featrue_y[:])

     corr_matrix[i,j] = cor(featrue_x,featrue_y,method = "spearman")

   }

 }

 rownames(corr_matrix) <- levels(mice_ident)

 colnames(corr_matrix) <- levels(expr_ident)

 #pheatmap(corr_matrix,clustering_distance_rows = "correlation",display_numbers = T,number_format = "%.2f")

 pheatmap(corr_matrix,clustering_distance_rows = "correlation",cluster_rows = T,cluster_cols =T,display_numbers = T,number_format = "%.2f")

 #########################################################################################################

 expr_all_P1_P8_nobatch_filter <- SubsetData(expr_all_P1_P8_nobatch, ident.remove = '')

 saveRDS(expr_all_P1_P8_nobatch_filter, file =paste0(outdir,"expr_all_P1_P8_filter.RDS"))

 expr_all_P1_P8_nobatch_filter <- SubsetData(expr_all_P1_P8_nobatch_filter,ident.remove = '')

 pdf(file = paste0(outdir,"Dimplot_filter_cell.pdf"),width = ,height = )

 DimPlot(expr_all_P1_P8_nobatch_filter,cols = c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333"))

 dev.off()

 Gene_train <- read.table("C:/Gu_lab/PA/result/mouse_scaleFeatures_rf/Train_Gene_capital_1.txt",sep = "\n",stringsAsFactors = F)

 Gene_train <- Gene_train[,]

 Gene_train

 scale_data <- expr_all_P1_P8_nobatch_filter@assays$integrated@scale.data

 rownames(scale_data)

 var_gene_exp <- scale_data[which(rownames(scale_data) %in% Gene_train),]

 write.table(var_gene_exp,paste0(outdir,"exp_train_gene_scale_cells.txt"),row.names = T,col.names = T, quote = F)

 #

 #python 预测

 #

 expr <- readRDS("C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/expr_all_P1_P8_filter.RDS")

 Gene_list <- read.table("C:/Gu_lab/PA/result/mouse_scaleFeatures_rf/Train_Gene_capital_1.txt")

 #expr<- RunPCA(expr,features = Gene_list$V1)

 #DimPlot(expr, reduction = "pca")

 #ElbowPlot(expr)

 #expr <- FindNeighbors(expr, dims = :)

 #expr <- FindClusters(expr, resolution = 0.5)

 #expr <- RunUMAP(expr,dims = :)

 #Plot 预测结果

 outdir <- "C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/Replot/"

 #expr <- readRDS("C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/expr_filter_cells.rds")

 Pred_data <- read.table('C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/Pred_results_scale_cell_resample.txt')

 newident <- Pred_data['V2'][[]]

 expr[["CellType"]] <- Pred_data['V2'][[]]

 newident <- factor(newident)

 Idents(expr) <- newident

 expr[["PredScore"]] <- Pred_data['V3'][[]]

 #expr_filter_unknown<- SubsetData(expr,ident.use = c("Lactotrope","Thyrotrope","Somatotrope"))

 expr_filter_unknown <- subset(expr,idents = c("Lactotrope","Somatotrope"))

 expr_filter_unknown <- RunPCA(expr_filter_unknown,features = Gene_list$V1[:])

 DimPlot(expr_filter_unknown, reduction = "pca")

 ElbowPlot(expr_filter_unknown)

 expr_filter_unknown <- RunUMAP(expr_filter_unknown,dims = :,n.neighbors = )

 pdf(file = paste0(outdir,"Dimplot_Pred_results.pdf"),width = ,height = )

 #去掉数目特别少的类

 DimPlot(expr_filter_unknown,pt.size = ,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 dev.off()

 #DimPlot(expr,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 Alapha = expr_filter_unknown@meta.data$PredScore

 umap_pos <- expr_filter_unknown@reductions$umap@cell.embeddings

 umap_pos_frame <- data.frame(umap_pos)

 #(tsne_pos_frame, aes(x = tSNE_1,y = tSNE_2,color = newident)) + geom_point(size =,shape = ,alpha = Alapha)

 #输出文件的位置

 dirout <- "C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/Replot/"

 #不同类别细胞分开画

 iden_cell = expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Lactotrope") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 #gg

 pdf(file = paste0(dirout,"Lactotrope.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell")+scale_colour_manual(values=c( "red","grey"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell =expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Melanotrope") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Melanotrope.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Melanotrope_like cell")+scale_colour_manual(values=c( "red","grey"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell = expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Corticotrope") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Corticotrope.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Corticotrope_like cell")+scale_colour_manual(values=c( "red","grey"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell = expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Proliferating_Pou1f1") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Proliferating_Pou1f1.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Proliferating_Pou1f1_like cell")+scale_colour_manual(values=c( "grey","red"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell = expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Somatotrope") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Somatotrope.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Somatotrope_like cell")+scale_colour_manual(values=c( "grey","red"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell = Pred_data['V2'][[]]

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Stem_cell") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Stem_cell.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type )) + geom_point(size =,shape = ,alpha =Alapha)+labs(title = "Stem_cell_like cell")+scale_colour_manual(values=c( "grey","red"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #不同类别细胞分开画

 iden_cell = expr_filter_unknown@meta.data$CellType

 for(i in :length(iden_cell)){

   if(iden_cell[i]!="Thyrotrope") {

     iden_cell[i] = "Others"

   }

 }

 Cell_type = iden_cell

 #gg <- ggplot(data = tsne_pos_frame)+geom_point(mapping = aes(x = tSNE_1,y = tSNE_2,color = CC),size =,shape = ,alpha = Alapha)+labs(title = "Lactotrope_like cell") + scale_color_hue(name = "Cell_type",labels = c("Lactotrope","Others"))

 pdf(file = paste0(dirout,"Thyrotrope.pdf"),width = ,height = )

 ggplot(umap_pos_frame, aes(x = UMAP_1,y = UMAP_2,color = Cell_type)) + geom_point(size =,shape = ,alpha = Alapha)+labs(title = "Thyrotrope_like cell")+scale_colour_manual(values=c( "grey","red"))+theme(plot.title = element_text(hjust = 0.5))

 dev.off()

 #看一下两坨细胞来源是否分别是两个样本

 cell_source <-colnames(expr_filter_unknown)

 cell_type <- strsplit(cell_source,split = '_')

 ct <-rep(c('P1'), times = length(cell_source))

 for( i in :length(cell_source)){

   if(cell_type[[i]][]=='P4'){

     ct[i]='P4'

   }

 }

 newident <- factor(ct)

 Idents(expr_filter_unknown) <- newident

 DimPlot(expr_filter_unknown,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 #统计一下，两中细胞来源的细胞的细胞类型比例

 sum1_S =

 sum1_L =

 sum4_S =

 sum4_L = 

 for(i in :length(cell_source)){

   if(cell_type[[i]][]=='P1' && Pred_data$V2[i]=='Somatotrope'){sum1_S = sum1_S + }

   else if(cell_type[[i]][]=='P1' && Pred_data$V2[i]=='Lactotrope'){sum1_L = sum1_L + }

   else if(cell_type[[i]][]=='P4' && Pred_data$V2[i]=='Lactotrope'){sum4_L = sum4_L + }

   else{sum4_S = sum4_S + }

 }

 #泌乳型的分的并不好，可能是去除了batch的原因，因为MNN的假设是不公的细胞类型相同

 #但因为P1和P4是不同的细胞类型，所以可以不用做MNN，可以用简单的回归，把表达量拉到同一尺度即可

 expr_1 <- readRDS("C:/Gu_lab/PA/result/pipline_results/P1_tumor/expr.RDS")

 expr_1 <- RenameCells(expr_1,add.cell.id = "P1_",for.merge = T )

 expr_4 <- readRDS("C:/Gu_lab/PA/result/pipline_results/P4_tumor/expr.RDS")

 expr_4 <- RenameCells(expr_4,add.cell.id = "P4_",for.merge = T )

 expr_1_4 <- merge(expr_1,expr_4)

 expr_1_4 <- FindVariableFeatures(expr_1_4)

 all.genes <- rownames(expr_1_4)

 expr_1_4 <- ScaleData(expr_1_4,features = all.genes)

 expr_1_4 <- RunPCA(expr_1_4,features = VariableFeatures(object = expr_1_4))

 ElbowPlot(expr_1_4)

 expr_1_4 <- FindNeighbors(expr_1_4,dims = :)

 expr_1_4 <- FindClusters(expr_1_4,resolution = 0.5)

 expr_1_4 <- RunUMAP(expr_1_4, dims = :)

 DimPlot(expr_1_4,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 cell_soruce <-rep(c('P','N'), times = c(,))

 newident_1 <- factor(cell_soruce)

 expr_1_4_save <- expr_1_4

 Idents(expr_1_4) <- newident_1

 #expr_nobatch_seurat@active.ident <- newident_1

 DimPlot(expr_1_4,reduction = "umap",cols =c("#F0E442", "#999999", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73","#FF0000","#333333") )

 expr_1_4 <- expr_1_4_save

 outdir <- "C:/Gu_lab/PA/result/Pred_result/P1_P8_Tumor/"

 saveRDS(expr_1_4, file =paste0(outdir,"expr_all_P1_P8_new.RDS"))

 #feature plot

 outdir <- "C:/Gu_lab/PA/result/Pred_result/P1_P8_Tumor/with_batch/"

 pdf(file = paste0(outdir,"Somatotrope_feature.pdf"))

 Somatotrope_feature =  c("Pappa2","Ceacam10","Tcerg1l","Cabp2","Wnt10a","RP24-294M17.1","Mmp7","Ghrhr","Gh")

 Somatotrope_feature = toupper(Somatotrope_feature)

 VlnPlot(expr_1_4, features =Somatotrope_feature)

 FeaturePlot(expr_1_4,features = Somatotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Lactotrope_feature.pdf"))

 Lactotrope_feature =  c("Prl","Agtr1a","Syndig1","Angpt1","Edil3","6030419C18Rik","Hepacam2","Cilp","Akr1c14")

 Lactotrope_feature = toupper(Lactotrope_feature)

 VlnPlot(expr_1_4, features =Lactotrope_feature)

 FeaturePlot(expr_1_4,features = Lactotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Corticotrope_feature.pdf"))

 Corticotrope_feature =  c("Gpc5","Tnnt1","Tnni3","Gm15543","Cplx3","Egr4","Cdh8","Galnt9","Pomc")

 Corticotrope_feature = toupper(Corticotrope_feature)

 VlnPlot(expr_1_4, features =Corticotrope_feature)

 FeaturePlot(expr_1_4,features = Corticotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Melanotrope_feature.pdf"))

 Melanotrope_feature = c("Oacyl","Esm1","Pkib","Gulo","Megf11","Rbfox3","Scn2a1","Pax7","Pcsk2")

 Melanotrope_feature = toupper(Melanotrope_feature)

 VlnPlot(expr_1_4, features =Melanotrope_feature)

 FeaturePlot(expr_1_4,features = Melanotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Thyrotrope_feature.pdf"))

 Thyrotrope_feature = c("Tshb","Trhr","Dio2")

 Thyrotrope_feature = toupper(Thyrotrope_feature)

 VlnPlot(expr_1_4, features =Thyrotrope_feature)

 FeaturePlot(expr_1_4,features = Thyrotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Stem_cell_feature.pdf"))

 Stem_cell_feature =  c("Cyp2f2","Lcn2","Mia","Cpxm2","Rbpms","Gm266","Cdh26","Aqp3","Aqp4")

 Stem_cell_feature = toupper(Stem_cell_feature)

 VlnPlot(expr_1_4, features =Stem_cell_feature)

 FeaturePlot(expr_1_4,features = Stem_cell_feature)

 dev.off()

 pdf(file = paste0(outdir,"Proliferatring_Pou1f1_feature.pdf"))

 Proliferatring_Pou1f1_feature = c("POU1f1","Pbk","Spc25","Ube2c","Hmmr","Ckap2l","Spc25","Cdca3","Birc5","Ccnb1")

 Proliferatring_Pou1f1_feature = toupper(Proliferatring_Pou1f1_feature)

 VlnPlot(expr_1_4, features =Proliferatring_Pou1f1_feature)

 FeaturePlot(expr_1_4,features = Proliferatring_Pou1f1_feature)

 dev.off()

 pdf(file = paste0(outdir,"T_B_NK_feature.pdf"))

 VlnPlot(expr_1_4,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"))#T-cell Markers

 FeaturePlot(expr_1_4,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"))

 VlnPlot(expr_1_4,features = c("PTPRC","CD79A","CD19","CD20","CD45"))#B-cell Markers

 FeaturePlot(expr_1_4,features = c("PTPRC","CD79A","CD19","CD20","CD45"))

 VlnPlot(expr_1_4,features = c("PTPRC","NKG7","CD56"))#NK-cell Markers

 FeaturePlot(expr_1_4,features = c("PTPRC","NKG7","CD56"))

 dev.off()

 pdf(file = paste0(outdir,"Macrophage-cell_feature.pdf"))

 VlnPlot(expr_1_4,features = c("CD14","CD163", "CD68", "CSF1R","CD33","HLA-DR","AIF1","FCER1G", "FCGR3A", "TYROBP"))#Macrophage-cell Markers

 FeaturePlot(expr_1_4,features = c("CD14","CD163", "CD68", "CSF1R","CD33","HLA-DR","AIF1","FCER1G", "FCGR3A", "TYROBP"))

 dev.off()

 pdf(file = paste0(outdir,"Fibroblasts-cell_feature.pdf"))

 VlnPlot(expr_1_4,features = c("ACTA2","FAP", "PDPN","COL1A2","DCN", "COL3A1", "COL6A1","S100A4","COL1A1","THY1"))#Fibroblasts-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Endothelial-cell_feature.pdf"))

 VlnPlot(expr_1_4,features = c("PLVAP","PECAM1","VWF","ENG","MCAM","CD146"))#Endothelial-cell Markers

 FeaturePlot(expr_1_4,features = c("PLVAP","PECAM1","VWF","ENG","MCAM","CD146"))

 dev.off()

 VlnPlot(expr_1_4,features = c("GH1","PRL","POU1F1"))

 #下面是为了去掉非垂体细胞：

 expr_1_4 <- SubsetData(expr_1_4,ident.use = c( '','','','','','','',''))

 pdf(file = paste0(outdir,"Dimplot_filter_cell.pdf"),width = ,height = )

 DimPlot(expr_1_4,cols = c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333"))

 dev.off()

 saveRDS(expr_1_4,file = paste0(outdir,"expr_1_4_filter.RDS"))

 Gene_train <- read.table("C:/Gu_lab/PA/result/mouse_scaleFeatures_rf/Train_Gene_capital.txt",sep = "\n",stringsAsFactors = F)

 Gene_train <- Gene_train[,]

 Gene_train

 scale_data <- expr_1_4@assays$RNA@scale.data

 rownames(scale_data)

 var_gene_exp <- scale_data[which(rownames(scale_data) %in% Gene_train),]

 write.table(var_gene_exp,"C:/Gu_lab/PA/result/Pred_result/P1_P8_Tumor/exp_train_gene_scale_cells.txt",row.names = T,col.names = T, quote = F)

 #去掉垂体细胞以后看一下C('GH1','PRL','POU1F1','TSHB')四个基因的分布

 pdf(file = paste0(outdir,"GH1_PRL_POU1F1_TSHB.pdf"),width = ,height = )

 VlnPlot(expr_1_4,features = c('GH1','PRL','POU1F1','TSHB'))#Endothelial-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"GH1_PRL_POU1F1_TSHB_FeaturePlot.pdf"),width = ,height = )

 FeaturePlot(expr_1_4,features = c('GH1','PRL','POU1F1','TSHB'))

 dev.off()

 #Plot 预测结果

 expr <- readRDS(file = paste0(outdir,"expr_1_4_filter.RDS"))

 Pred_data <- read.table('C:/Gu_lab/PA/result/Pred_result/P1_P8_tumor/Pred_results_scale_cell_resample.txt')

 newident <- Pred_data['V2'][[]]

 newident <- factor(newident)

 Idents(expr) <- newident

 pdf(file = paste0(outdir,"Dimplot_Pred_results.pdf"),width = ,height = )

 #去掉数目特别少的类

 DimPlot(expr,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 dev.off()

 #DimPlot(expr,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 Alapha = Pred_data['V3'][[]]

 tsne_pos <- expr@reductions$tsne@cell.embeddings

 tsne_pos_frame <- data.frame(tsne_pos)

 col = c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442",  "#CC79A7", "#0072B2", "#D55E00","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333")

 cid <- Pred_data['V4'][[]]

 CC = col[cid+]

 #去掉unknown

 expr_filter_unknown <- SubsetData(expr,ident.remove = "Unknown")

 DimPlot(expr_filter_unknown,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 expr_filter_unknown <- SubsetData(expr,ident.use = c("Lactotrope","Somatotrope"))

 DimPlot(expr_filter_unknown,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 expr_filter_unknown <- FindVariableFeatures(expr_filter_unknown)

 all.genes <- rownames(expr_filter_unknown)

 expr_filter_unknown <- ScaleData(expr_filter_unknown,features = all.genes)

 expr_filter_unknown <- RunPCA(expr_filter_unknown,features = VariableFeatures(object = expr_filter_unknown))

 ElbowPlot(expr_filter_unknown)

 #expr_filter_unknown <- FindNeighbors(expr_filter_unknown,dims = :)

 #expr_filter_unknown <- FindClusters(expr_filter_unknown,resolution = 0.5)

 expr_filter_unknown <- RunUMAP(expr_filter_unknown, dims = :)

 DimPlot(expr_filter_unknown,reduction = "umap",cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#000099","#660066","#333333") )

 expr_filter_unknown_save <- expr_filter_unknown

 #看一下细胞来源

 #看一下两坨细胞来源是否分别是两个样本

 cell_source <-colnames(expr_filter_unknown)

 cell_type <- strsplit(cell_source,split = '_')

 ct <-rep(c('P1'), times = length(cell_source))

 for( i in :length(cell_source)){

   if(cell_type[[i]][]=='P4'){

     ct[i]='P4'

   }

 }

 newident <- factor(ct)

 Idents(expr_filter_unknown) <- newident

 DimPlot(expr_filter_unknown,cols =c("#999999","#FF0099", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7","#990000","#9900cc","#66FF66","#663300","#0000FF","#CC0033","#FF0000","#333333") )

 #统计一下，两中细胞来源的细胞的细胞类型比例

 sum1_S =

 sum1_L =

 sum4_S =

 sum4_L = 

 predresult <- expr_filter_unknown@active.ident

 for(i in :length(cell_type)){

   if(cell_type[[i]][]=='P1' && predresult[i]=='Somatotrope'){sum1_S = sum1_S + }

   else if(cell_type[[i]][]=='P1' && predresult[i]=='Lactotrope'){sum1_L = sum1_L + }

   else if(cell_type[[i]][]=='P4' && predresult[i]=='Lactotrope'){sum4_L = sum4_L + }

   else{sum4_S = sum4_S + }

 }

 #看一下细胞的marker表达情况

 VlnPlot(expr_filter_unknown,features = c("GH1","PRL","TSHB","POU1F1"))#Endothelial-cell Markers

 FeaturePlot(expr_filter_unknown,features = c("GH1","PRL","TSHB","POU1F1"))

 pdf(file = paste0(outdir,"Somatotrope_feature.pdf"))

 Somatotrope_feature =  c("Pappa2","Ceacam10","Tcerg1l","Cabp2","Wnt10a","RP24-294M17.1","Mmp7","Ghrhr","Gh")

 Somatotrope_feature = toupper(Somatotrope_feature)

 VlnPlot(expr, features =Somatotrope_feature)

 FeaturePlot(expr,features = Somatotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Lactotrope_feature.pdf"))

 Lactotrope_feature =  c("Prl","Agtr1a","Syndig1","Angpt1","Edil3","6030419C18Rik","Hepacam2","Cilp","Akr1c14")

 Lactotrope_feature = toupper(Lactotrope_feature)

 VlnPlot(expr, features =Lactotrope_feature)

 FeaturePlot(expr,features = Lactotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Corticotrope_feature.pdf"))

 Corticotrope_feature =  c("Gpc5","Tnnt1","Tnni3","Gm15543","Cplx3","Egr4","Cdh8","Galnt9","Pomc")

 Corticotrope_feature = toupper(Corticotrope_feature)

 VlnPlot(expr, features =Corticotrope_feature)

 FeaturePlot(expr,features = Corticotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Melanotrope_feature.pdf"))

 Melanotrope_feature = c("Oacyl","Esm1","Pkib","Gulo","Megf11","Rbfox3","Scn2a1","Pax7","Pcsk2")

 Melanotrope_feature = toupper(Melanotrope_feature)

 VlnPlot(expr, features =Melanotrope_feature)

 FeaturePlot(expr,features = Melanotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Thyrotrope_feature.pdf"))

 Thyrotrope_feature = c("Tshb","Trhr","Dio2")

 Thyrotrope_feature = toupper(Thyrotrope_feature)

 VlnPlot(expr, features =Thyrotrope_feature)

 FeaturePlot(expr,features = Thyrotrope_feature)

 dev.off()

 pdf(file = paste0(outdir,"Stem_cell_feature.pdf"))

 Stem_cell_feature =  c("Cyp2f2","Lcn2","Mia","Cpxm2","Rbpms","Gm266","Cdh26","Aqp3","Aqp4")

 Stem_cell_feature = toupper(Stem_cell_feature)

 VlnPlot(expr, features =Stem_cell_feature)

 FeaturePlot(expr,features = Stem_cell_feature)

 dev.off()

 pdf(file = paste0(outdir,"Proliferatring_Pou1f1_feature.pdf"))

 Proliferatring_Pou1f1_feature = c("POU1f1","Pbk","Spc25","Ube2c","Hmmr","Ckap2l","Spc25","Cdca3","Birc5","Ccnb1")

 Proliferatring_Pou1f1_feature = toupper(Proliferatring_Pou1f1_feature)

 VlnPlot(expr, features =Proliferatring_Pou1f1_feature)

 FeaturePlot(expr,features = Proliferatring_Pou1f1_feature)

 dev.off()

 pdf(file = paste0(outdir,"T_B_NK_feature.pdf"))

 VlnPlot(expr,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"))#T-cell Markers

 FeaturePlot(expr,features = c("CD2","CD3D","CD3E","CD3G","CD4","CD8","CD45"))

 VlnPlot(expr,features = c("PTPRC","CD79A","CD19","CD20","CD45"))#B-cell Markers

 FeaturePlot(expr,features = c("PTPRC","CD79A","CD19","CD20","CD45"))

 VlnPlot(expr,features = c("PTPRC","NKG7","CD56"))#NK-cell Markers

 FeaturePlot(expr,features = c("PTPRC","NKG7","CD56"))

 dev.off()

 pdf(file = paste0(outdir,"Macrophage-cell_feature.pdf"))

 VlnPlot(expr,features = c("CD14","CD163", "CD68", "CSF1R","CD33","HLA-DR","AIF1","FCER1G", "FCGR3A", "TYROBP"))#Macrophage-cell Markers

 FeaturePlot(expr,features = c("CD14","CD163", "CD68", "CSF1R","CD33","HLA-DR","AIF1","FCER1G", "FCGR3A", "TYROBP"))

 dev.off()

 pdf(file = paste0(outdir,"Fibroblasts-cell_feature.pdf"))

 VlnPlot(expr,features = c("ACTA2","FAP", "PDPN","COL1A2","DCN", "COL3A1", "COL6A1","S100A4","COL1A1","THY1"))#Fibroblasts-cell Markers

 dev.off()

 pdf(file = paste0(outdir,"Endothelial-cell_feature.pdf"))

 VlnPlot(expr,features = c("PLVAP","PECAM1","VWF","ENG","MCAM","CD146"))#Endothelial-cell Markers

 FeaturePlot(expr,features = c("PLVAP","PECAM1","VWF","ENG","MCAM","CD146"))

 dev.off()

 VlnPlot(expr,features = c("GH1","CHGA","SLPI"))

 VlnPlot(expr,features = c("GH1","CHGA","SLPI"))

 VlnPlot(expr,features = c("GH"))
正文

相关阅读

病毒程序算法C语言，病毒程序算法c语言代码

html5标记绘图数学图形，html小图标代码

js抽奖代码抽奖统计，javascript抽奖小程序

html5好玩的代码，H5代码

js中为什么要用到闭包，JS闭包的示例代码

html5密码强度，html密码框代码怎么写

京东搜索框代码html5，京东搜索栏

html5代码下拉箭头显示出来，html5下拉菜单选项代码

目录[+]
