R Seurat 单细胞处理pipline 代码

阅读原文时间：2021年04月21日阅读：1

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")
}
}

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

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

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

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

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"))

手机扫一扫

移动阅读更方便

你可能感兴趣的文章