hu-mu comparison and DIY for muPrUr

r.scripts/diy_muPrUr.R

+library(Seurat)
+library(SingleR)
+library(ggplot2)
+library(cowplot)
+library(gridExtra)
+library(RColorBrewer)
+library(viridis)
+
+options(bitmapType="cairo")
+
+setwd("../")
+
+load("./analysis/sc10x.id.rda")
+load("./analysis/singler_objects.RData")
+
+if (!dir.exists(paste0("./analysis/vis/diy"))){
+  dir.create(paste0("./analysis/vis/diy"))
+}
+
+
+sc10x$Region <- sc10x$samples
+Idents(sc10x) <- "Region"
+sc10x <- RenameIdents(sc10x,"musAd001_PrF"="Prostate")
+sc10x <- RenameIdents(sc10x,"musAd002_PrF"="Prostate")
+sc10x <- RenameIdents(sc10x,"musAd003_PrF_St"="Prostate")
+sc10x <- RenameIdents(sc10x,"musAd007_PrFcol_GEX"="Prostate")
+sc10x <- RenameIdents(sc10x,"musAd002_UrF"="Urethra")
+sc10x <- RenameIdents(sc10x,"musAd004n5_UrF_GEX"="Prostate")
+sc10x$Region <- Idents(sc10x)
+sc10x$Region <- factor(sc10x$Region,levels=c("Prostate","Urethra"))
+
+sc10x$`Cell Lineage` <- sc10x$lin
+Idents(sc10x) <- "Cell Lineage"
+sc10x$`Cell Lineage` <- Idents(sc10x)
+sc10x$`Cell Lineage` <- factor(sc10x$`Cell Lineage`,levels=c("Epi","FMSt","Endo","Leu"))
+Idents(sc10x) <- "Cell Lineage"
+
+postscript(paste0("./analysis/vis/diy/UMAP_lin.eps"))
+DimPlot(sc10x)+scale_color_viridis(discrete=TRUE,option="magma")+theme_cowplot()
+dev.off()
+postscript(paste0("./analysis/vis/diy/UMAP_lin.split.eps"))
+DimPlot(sc10x,split.by="Region")+scale_color_viridis(discrete=TRUE,option="magma")+theme_cowplot()
+dev.off()
+
+labs.lin <- singler.lin$labels
+labs.lin <- replace(labs.lin,labs.lin=="Epithelial cells","Epi")
+labs.lin <- replace(labs.lin,labs.lin %in% c("Stromal cells","Fibroblasts"),"FMSt")
+labs.lin <- replace(labs.lin,labs.lin == "Endothelial cells","Endo")
+labs.lin <- replace(labs.lin,labs.lin %in% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia"),"Leu")
+labs.lin <- data.frame(Cell.Lineage=labs.lin,row.names=rownames(singler.lin))
+labs.lin$Cell.Lineage <- factor(labs.lin$Cell.Lineage,levels=c("Epi","FMSt","Endo","Leu"))
+labs.lin.col =list(
+  Cell.Lineage=c(Epi=magma(4)[1],FMSt=magma(4)[2],Endo=magma(4)[3],Leu=magma(4)[4])
+)
+singler.lin@listData$scores <- singler.lin@listData$scores[,c(14,7,5,9,10,16,12,3,11,2,1,4,8,13,15,6)]
+plot <- plotScoreHeatmap(singler.lin,show.labels=FALSE,show.pruned=FALSE,annotation_col=labs.lin,cluster_rows=FALSE,annotation_colors=labs.lin.col,cells.order= as.integer(c(1,2,3,4,5,16,14,8,9,10,15,13,12,7,11,6)))
+postscript(paste0("./analysis/vis/diy/Corr_lin.eps"))
+print(plot)
+dev.off()
+
+postscript(paste0("./analysis/vis/diy/Dot_lin.anchor.eps"))
+DotPlot(sc10x,features=rev(c("Epcam","Cdh1","Dcn","Myl9","Pecam1","Vwf","Ptprc")),dot.scale=10,cols=rev(heat.colors(2)))
+dev.off()
+
+sc10x.epi$Region <- sc10x.epi$samples
+Idents(sc10x.epi) <- "Region"
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd001_PrF"="Prostate")
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd002_PrF"="Prostate")
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd003_PrF_St"="Prostate")
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd007_PrFcol_GEX"="Prostate")
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd002_UrF"="Urethra")
+sc10x.epi <- RenameIdents(sc10x.epi,"musAd004n5_UrF_GEX"="Prostate")
+sc10x.epi$Region <- Idents(sc10x.epi)
+sc10x.epi$Region <- factor(sc10x.epi$Region,levels=c("Prostate","Urethra"))
+
+
+Idents(sc10x.epi) <- "integrated_snn_res.0.1"
+cells.ap <- WhichCells(sc10x.epi,idents="0")
+cells.vp <- WhichCells(sc10x.epi,idents="3")
+cells.ed <- WhichCells(sc10x.epi,idents="4")
+cells.dlp <- WhichCells(sc10x.epi,idents="5")
+
+sc10x.epi$`Cell Type` <- sc10x.epi$hu_pops
+Idents(sc10x.epi) <- "Cell Type"
+sc10x.epi <- SetIdent(sc10x.epi,cell=cells.ap,value="APr LE")
+sc10x.epi <- SetIdent(sc10x.epi,cell=cells.vp,value="VPr LE")
+sc10x.epi <- SetIdent(sc10x.epi,cell=cells.ed,value="ED")
+sc10x.epi <- SetIdent(sc10x.epi,cell=cells.dlp,value="DLPr LE")
+sc10x.epi$`Cell Type` <- Idents(sc10x.epi)
+sc10x.epi$`Cell Type` <- factor(sc10x.epi$`Cell Type`,levels=c("BE","Ur","VPr LE","DLPr LE","APr LE","ED"))
+Idents(sc10x.epi) <- "Cell Type"
+
+postscript(paste0("./analysis/vis/diy/UMAP_epi.eps"))
+DimPlot(sc10x.epi)+scale_color_viridis(discrete=TRUE)+theme_cowplot()
+dev.off()
+postscript(paste0("./analysis/vis/diy/UMAP_epi.split.eps"))
+DimPlot(sc10x.epi,split.by="Region")+scale_color_viridis(discrete=TRUE)+theme_cowplot()
+dev.off()
+
+labs.epi <- data.frame(ID=singler.epi$labels,`Cell Type`=c("APr LE","BE","Ur","VPr LE","ED","DLPr LE"),row.names=rownames(singler.epi))
+labs.epi$ID <- factor(labs.epi$ID,levels=c("BE","Ur","LE"))
+labs.epi.col =list(
+  ID=c(BE="brown4",Ur="brown3",LE="brown2"),
+  Cell.Type=c(BE=viridis(6)[1],Ur=viridis(6)[2],`VPr LE`=viridis(6)[3],`DLPr LE`=viridis(6)[4],`APr LE`=viridis(6)[5],ED=viridis(6)[6])
+)
+labs.epi$Cell.Type <- factor(labs.epi$Cell.Type,levels=c("BE","Ur","VPr LE","DLPr LE","APr LE","ED"))
+singler.epi@listData$scores <- singler.epi@listData$scores[,c(1,3,2)]
+plot <- plotScoreHeatmap(singler.epi,show.labels=FALSE,show.pruned=FALSE,annotation_col=labs.epi,cluster_rows=FALSE,annotation_colors=labs.epi.col,cells.order= as.integer(c(5,1,2,3,6,4)))
+postscript(paste0("./analysis/vis/diy/Corr_epi.eps"))
+print(plot)
+dev.off()
+
+postscript(paste0("./analysis/vis/diy/Dot_epi.anchor.eps"))
+DotPlot(sc10x.epi,features=rev(c("Krt5","Krt4","Sbp","Msmb","Tgm4","Svs2")),dot.scale=10,cols=rev(heat.colors(2)))
+dev.off()
+
+deg.lin <- FindAllMarkers(sc10x,assay="SCT",slot="data",logfc.threshold=0,test.use="MAST",only.pos=TRUE)
+deg.epi <- FindAllMarkers(sc10x.epi,assay="SCT",slot="data",logfc.threshold=0,test.use="MAST",only.pos=TRUE)
+
+write.table(deg.lin,file="./analysis/vis/diy/DEG_lin.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=TRUE)
+write.table(deg.epi,file="./analysis/vis/diy/DEG_epi.csv",sep=",",quote=FALSE,row.names=FALSE,col.names=TRUE)
+
+Idents(sc10x) <- "integrated_snn_res.0.5"
+postscript(paste0("./analysis/vis/diy/UMAP_lin.res0.5.eps"))
+DimPlot(sc10x)+theme_cowplot()
+dev.off()
+Idents(sc10x.epi) <- "integrated_snn_res.0.1"
+postscript(paste0("./analysis/vis/diy/UMAP_epi.res0.1.eps"))
+DimPlot(sc10x.epi)+theme_cowplot()
+dev.off()
\ No newline at end of file

r.scripts/huPrmuPr.R

@@ -19,9 +19,8 @@ library(ggplot2)
 options(bitmapType="cairo")
 
 option_list=list(
-  make_option("--p",action="store",default="muPr",type='character',help="Project Name"),
-  make_option("--s",action="store",default="mu",type='character',help="Species"),
-  make_option("--o",action="store",default="pr",type='character',help="Organ")
+  make_option("--p",action="store",default="muPrUr",type='character',help="Project Name"),
+  make_option("--r",action="store",default="epi",type='character',help="Reference Subset")
 )
 opt=parse_args(OptionParser(option_list=option_list))
 rm(option_list)
@@ -32,173 +31,68 @@ setwd("../")
 source("./r.scripts/sc-TissueMapper_functions.R")
 source("./r.scripts/sc-TissueMapper_process.R")
 
-if (opt$s == "hu"){
-  hpca.se <- HumanPrimaryCellAtlasData()
-  lin.se <- hpca.se
-  leu.se <- hpca.se
-} else if (opt$s == "mu"){
-  igd.se <- ImmGenData()
-  #lin.se <- igd.se
-  #leu.se <- igd.se
-  lin.se <- subset(igd.se,select=(as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) >= 920648 & as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) <= 2112477))
-  lin.se <- subset(lin.se,select=!(as.numeric(gsub("_","",substr(colnames(lin.se),4,10))) %in% 1398469:1398471))
-  leu.se <- lin.se
-  
-  convert <- read.delim("/work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/RUN/muPr/genesets/Ensemble.mus-hum.txt")
-}
+igd.se <- ImmGenData()
+lin.se <- subset(igd.se,select=(as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) >= 920648 & as.numeric(gsub("_","",substr(colnames(igd.se),4,10))) <= 2112477))
+lin.se <- subset(lin.se,select=!(as.numeric(gsub("_","",substr(colnames(lin.se),4,10))) %in% 1398469:1398471))
+leu.se <- lin.se
+rm(igd.se)
 
-if (opt$o == "pr"){
-  load("./genesets/scDWShuPr.rda")
-  
-  try(
-    if (as.numeric(substring(sc10x@version,1,1))<3){
-      sc10x <- UpdateSeuratObject(sc10x)
-    }
-  )
-  Idents(sc10x) <- "pops"
-  sc10x <- subset(sc10x,idents=c("NE","Endo"),invert=TRUE)
-  Idents(sc10x) <- "lin"
-  sc10x <- subset(sc10x,idents=c("Leu"),invert=TRUE)
-  scDWSpr.se <- as.SingleCellExperiment(sc10x)
-  rm(sc10x)
+convert <- read.delim("/work/urology/ghenry/RNA-Seq/SingleCell/PIPELINE/RUN/muPr/genesets/Ensemble.mus-hum.txt")
+
+if (opt$r == "epi"){
+  load("./genesets/sc10x.epi.id.rda")
+  sc10x <- sc10x.epi
   rm(sc10x.epi)
-  rm(sc10x.fmst)
-  rm(sc10x.leu)
+  Idents(sc10x) <- "pops"
+  sc10x <- subset(sc10x,idents=c("NE"),invert=TRUE)
+  sc10x <- RenameIdents(sc10x,"Hillock"="Ur")
+  sc10x <- RenameIdents(sc10x,"Club"="Ur")
+  sc10x$pops <- Idents(sc10x)
+  sc10x$pops <- factor(sc10x$pops,levels=c("BE","LE","Ur"))
+  Idents(sc10x) <- "pops"
+} else if (opt$r == "fmst"){
+  load("./genesets/sc10x.fmst.id.rda")
+  sc10x <- sc10x.st
   rm(sc10x.st)
-  
-  if (opt$s == "mu"){
-    scDWSpr.se <- scDWSpr.se[rownames(scDWSpr.se) %in% convert$Human.gene.name,]
-    rownames(scDWSpr.se) <- convert$Gene.name[match(rownames(scDWSpr.se),convert$Human.gene.name)]
-  }
+  Idents(sc10x) <- "pops"
 }
 
+try(
+  if (as.numeric(substring(sc10x@version,1,1))<3){
+    sc10x <- UpdateSeuratObject(sc10x)
+  }
+)
+Idents(sc10x) <- "pops"
+
+scDWSpr.se <- as.SingleCellExperiment(sc10x,assay="RNA")
+rm(sc10x)
+scDWSpr.se <- scDWSpr.se[rownames(scDWSpr.se) %in% convert$Human.gene.name,]
+rownames(scDWSpr.se) <- convert$Gene.name[match(rownames(scDWSpr.se),convert$Human.gene.name)]
+
 load("./analysis/sc10x.raw.rda")
 
-sc10x.se <- as.SingleCellExperiment(sc10x)
+sc10x.se <- as.SingleCellExperiment(sc10x,assay="RNA")
 common <- intersect(rownames(sc10x.se),rownames(lin.se))
 lin.se <- lin.se[common,]
 sc10x.se <- sc10x.se[common,]
 rm(common)
 
-#singler.lin <- SingleR(sc10x.se,ref=lin.se,method="cluster",clusters=sc10x.se$integrated_snn_res.0.1,labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
-#sc10x$lin <- singler.lin$labels[match(sc10x.se$integrated_snn_res.0.1,singler.lin@rownames)]
-singler.lin <- SingleR(sc10x.se,ref=lin.se,method="single",labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
-sc10x$lin <- singler.lin$labels
-if (opt$s == "hu"){
-  sc10x$lin[sc10x$lin %in% c("DC","B_cell","Neutrophil","T_cells","Monocyte","Macrophage","NK_cell","Neutrophils","CMP","GMP","MEP","Myelocyte","Pre-B_cell_CD34-","Pro-B_cell_CD34+","Pro-Myelocyte","HSC_-G-CSF","HSC_CD34+")] <- "Leu"
-  sc10x$lin[sc10x$lin %in% c("Endothelial_cells","Erythroblast","Platelets")] <- "Endo"
-  sc10x$lin[sc10x$lin %in% c("Smooth_muscle_cells","Fibroblasts","Chondrocytes","Osteoblasts","MSC","Tissue_stem_cells")] <- "FMSt"
-  sc10x$lin[sc10x$lin %in% c("Epithelial_cells","Keratinocytes","Neuroepithelial_cell")] <- "Epi"
-} else if (opt$s == "mu"){
-  sc10x$lin[sc10x$lin %in% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")] <- "Leu"
-  sc10x$lin[sc10x$lin %in% c("Endothelial cells")] <- "Endo"
-  sc10x$lin[sc10x$lin %in% c("Stromal cells","Fibroblasts")] <- "FMSt"
-  sc10x$lin[sc10x$lin %in% c("Epithelial cells")] <- "Epi"
-  #sc10x$leu <- sc10x$lin
-  #sc10x$leu[sc10x$leu != "Leu"] <- "non-Leu"
-}
-DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
+singler.lin <- SingleR(sc10x.se,ref=lin.se,method="cluster",clusters=sc10x.se$integrated_snn_res.0.5,labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
+sc10x$lin <- singler.lin$labels[match(sc10x.se$integrated_snn_res.0.5,singler.lin@rownames)]
+#singler.lin <- SingleR(sc10x.se,ref=lin.se,method="single",labels=lin.se$label.main,BPPARAM=MulticoreParam(workers=10))
+#sc10x$lin <- singler.lin$labels
+sc10x$lin[sc10x$lin %in% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")] <- "Leu"
+sc10x$lin[sc10x$lin %in% c("Endothelial cells")] <- "Endo"
+sc10x$lin[sc10x$lin %in% c("Stromal cells","Fibroblasts")] <- "FMSt"
+sc10x$lin[sc10x$lin %in% c("Epithelial cells")] <- "Epi"
+#DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
 
-#Idents(sc10x) <- "leu"
-#sc10x.leu <- subset(sc10x,idents="Leu")
-#sc10x.nonleu <- subset(sc10x,idents="non-Leu")
 Idents(sc10x) <- "lin"
-sc10x.EpiFMSt <- subset(sc10x,idents=c("Leu","Endo"),invert=TRUE)
+sc10x.epi <- subset(sc10x,idents=c("Epi"))
+sc10x.fmst <- subset(sc10x,idents=c("FMSt"))
 
 res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
 
-try({
-  results <- scPC(sc10x.EpiFMSt,pc=100,hpc=0.9,file="EpiFMSt",print="2")
-  sc10x.EpiFMSt <- results[[1]]
-  pc.use.EpiFMSt <- results[[2]]
-  rm(results)
-  sc10x.EpiFMSt <- scCluster(sc10x.EpiFMSt,res=res,red="pca",dim=pc.use.EpiFMSt,print="2",folder="EpiFMSt")
-})
-
-sc10x.se.EpiFMSt <- as.SingleCellExperiment(sc10x.EpiFMSt)
-common <- intersect(rownames(sc10x.se.EpiFMSt),rownames(scDWSpr.se))
-scDWSpr.se <- scDWSpr.se[common,]
-sc10x.se.EpiFMSt <- sc10x.se.EpiFMSt[common,]
-rm(common)
-
-singler.EpiFMSt <- SingleR(sc10x.se.EpiFMSt,ref=scDWSpr.se,method="cluster",clusters=sc10x.se.EpiFMSt$integrated_snn_res.1,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
-sc10x.EpiFMSt$hu_pops <- singler.EpiFMSt$first.labels[match(sc10x.se.EpiFMSt$integrated_snn_res.1,singler.EpiFMSt@rownames)]
-#singler.EpiFMSt <- SingleR(sc10x.se.EpiFMSt,ref=scDWSpr.se,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
-#sc10x.EpiFMSt$hu_pops <- singler.EpiFMSt$labels
-Idents(sc10x.EpiFMSt) <- "hu_pops"
-#sc10x.nonleu$scDWSpr[sc10x.nonleu$scDWSpr == "OE2"] <- "Hillock"
-#sc10x.nonleu$scDWSpr[sc10x.nonleu$scDWSpr == "OE1"] <- "Club"
-sc10x.EpiFMSt$hu_pops <- factor(sc10x.EpiFMSt$hu_pops)
-DimPlot(sc10x.EpiFMSt,group.by="hu_pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-
-
-
-
-Idents(sc10x.EpiFMSt) <- "hu_pops"
-sc10x.epi <- subset(sc10x.EpiFMSt,idents=levels(sc10x.EpiFMSt)[levels(sc10x.EpiFMSt) %in% c("BE","LE","Club","Hillock")])
-sc10x.fmst <- subset(sc10x.EpiFMSt,idents=levels(sc10x.EpiFMSt)[levels(sc10x.EpiFMSt) %in% c("Fib","SM")])
-#sc10x.st <- subset(sc10x.EpiFMSt,idents=setdiff(levels(sc10x.EpiFMSt),c("BE","LE","Club","Hillock")))
-
-if (opt$o == "pr" && opt$s == "hu"){
-  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi)
-  scDWSpr.se.epi <- scDWSpr.se[,scDWSpr.se$Lin=="Epi",]
-  common <- intersect(rownames(sc10x.se.epi),rownames(scDWSpr.se.epi))
-  scDWSpr.se.epi <- scDWSpr.se.epi[common,]
-  sc10x.se.epi <- sc10x.se.epi[common,]
-  rm(common)
-  
-  sc10x.se.fmst <- as.SingleCellExperiment(sc10x.fmst)
-  scDWSpr.se.fmst <- scDWSpr.se[,scDWSpr.se$Merge_Epi.dws_St.go %in% c("Fib","SM"),]
-  common <- intersect(rownames(sc10x.se.fmst),rownames(scDWSpr.se.fmst))
-  scDWSpr.se.fmst <- scDWSpr.se.fmst[common,]
-  sc10x.se.fmst <- sc10x.se.fmst[common,]
-  rm(common)
-  
-  singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se.epi,labels=scDWSpr.se.epi$ident,BPPARAM=MulticoreParam(workers=10))
-  sc10x.epi$scDWSpr <- singler.epi$labels
-  Idents(sc10x.epi) <- "scDWSpr"
-  sc10x.epi$scDWSpr[sc10x.epi$scDWSpr == "OE2"] <- "Hillock"
-  sc10x.epi$scDWSpr[sc10x.epi$scDWSpr == "OE1"] <- "Club"
-  
-  singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se.fmst,labels=scDWSpr.se.fmst$ident,BPPARAM=MulticoreParam(workers=10))
-  sc10x.fmst$scDWSpr <- singler.fmst$labels
-  Idents(sc10x.fmst) <- "scDWSpr"
-}
-
-#sc10x.se.leu <- as.SingleCellExperiment(sc10x.leu)
-#leu.se <- leu.se[,leu.se$label.main %in% c("DC","B_cell","Neutrophil","T_cells","Monocyte","Macrophage","NK_cell","Neutrophils","CMP","GMP","MEP","Myelocyte","Pre-B_cell_CD34-","Pro-B_cell_CD34+","Pro-Myelocyte","Macrophages","Monocytes","B cells","Eosinophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia")]
-#common <- intersect(rownames(sc10x.se.leu),rownames(leu.se))
-#leu.se <- leu.se[common,]
-#sc10x.se.leu <- sc10x.se.leu[common,]
-#rm(common)
-
-#singler.leu <- SingleR(sc10x.se.leu,ref=leu.se,labels=leu.se$label.main,BPPARAM=MulticoreParam(workers=10))
-#sc10x.leu$leu <- singler.leu$labels
-#Idents(sc10x.leu) <- "leu"
-#DimPlot(sc10x.leu,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-#sc10x$pops <- sc10x$lin
-#sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
-#sc10x.epi$pops <- sc10x.epi$lin
-#sc10x.fmst$pops <- sc10x.fmst$lin
-#sc10x.st$pops <- sc10x.st$lin
-#sc10x.leu$pops <- sc10x.leu$leu
-if (opt$o == "pr" && opt$s == "hu"){
-  sc10x$pops[names(sc10x.epi$scDWSpr)] <- sc10x.epi$scDWSpr
-  sc10x.epi$pops <- sc10x.epi$scDWSpr
-  sc10x$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-  sc10x.fmst$pops <- sc10x.fmst$scDWSpr
-  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-  sc10x$pops[names(sc10x.leu$leu)] <- sc10x.leu$leu
-  sc10x.st$pops[names(sc10x.fmst$scDWSpr)] <- sc10x.fmst$scDWSpr
-}
-#DimPlot(sc10x,group.by="pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
-#sc10x$leu <- "non-Leu"
-#sc10x$leu[names(sc10x.leu$leu)] <- sc10x.leu$leu
-#DimPlot(sc10x,group.by="leu",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-
 try({
   results <- scPC(sc10x.epi,pc=100,hpc=0.9,file="epi",print="2")
   sc10x.epi <- results[[1]]
@@ -206,7 +100,6 @@ try({
   rm(results)
   sc10x.epi <- scCluster(sc10x.epi,res=res,red="pca",dim=pc.use.epi,print="2",folder="epi")
 })
-
 try({
   results <- scPC(sc10x.fmst,pc=100,hpc=0.9,file="fmst",print="2")
   sc10x.fmst <- results[[1]]
@@ -215,26 +108,41 @@ try({
   sc10x.fmst <- scCluster(sc10x.fmst,res=res,red="pca",dim=pc.use.fmst,print="2",folder="fmst")
 })
 
-#try({
-#  results <- scPC(sc10x.st,pc=100,hpc=0.9,file="st",print="2")
-#  sc10x.st <- results[[1]]
-#  pc.use.st <- results[[2]]
-#  rm(results)
-#  sc10x.st <- scCluster(sc10x.st,res=res,red="pca",dim=pc.use.st,print="2",folder="st")
-#})
-
-#try({
-#  results <- scPC(sc10x.leu,pc=100,hpc=0.9,file="leu",print="2")
-#  sc10x.leu <- results[[1]]
-#  pc.use.leu <- results[[2]]
-#  rm(results)
-#  sc10x.leu <- scCluster(sc10x.leu,res=res,red="pca",dim=pc.use.leu,print="2",folder="leu")
-#})
+if (opt$r == "epi"){
+  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi,assay="RNA")
+  common <- intersect(rownames(sc10x.se.epi),rownames(scDWSpr.se))
+  scDWSpr.se <- scDWSpr.se[common,]
+  sc10x.se.epi <- sc10x.se.epi[common,]
+  rm(common)
+  
+  singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se,method="cluster",clusters=sc10x.se.epi$integrated_snn_res.0.1,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
+  sc10x.epi$hu_pops <- singler.epi$labels[match(sc10x.se.epi$integrated_snn_res.0.1,singler.epi@rownames)]
+  #singler.epi <- SingleR(sc10x.se.epi,ref=scDWSpr.se,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
+  #sc10x.epi$hu_pops <- singler.epi$labels
+  Idents(sc10x.epi) <- "hu_pops"
+  sc10x.epi$hu_pops <- factor(sc10x.epi$hu_pops)
+  #DimPlot(sc10x.epi,group.by="hu_pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
+  sc10x.fmst$hu_pops <- sc10x.fmst$lin
+} else if (opt$r == "fmst"){
+  sc10x.se.fmst <- as.SingleCellExperiment(sc10x.fmst,assay="RNA")
+  common <- intersect(rownames(sc10x.se.fmst),rownames(scDWSpr.se))
+  scDWSpr.fmse <- scDWSpr.se[common,]
+  sc10x.se.fmst <- sc10x.se.fmst[common,]
+  rm(common)
+  
+  singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se,method="cluster",clusters=sc10x.se.fmst$integrated_snn_res.0.1,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
+  sc10x.fmst$hu_pops <- singler.fmst$labels[match(sc10x.se.fmst$integrated_snn_res.0.1,singler.fmst@rownames)]
+  #singler.fmst <- SingleR(sc10x.se.fmst,ref=scDWSpr.se,labels=scDWSpr.se$pops,BPPARAM=MulticoreParam(workers=10))
+  #sc10x.fmst$hu_pops <- singler.fmst$labels
+  Idents(sc10x.fmst) <- "hu_pops"
+  sc10x.fmst$hu_pops <- factor(sc10x.fmst$hu_pops)
+  #DimPlot(sc10x.fmst,group.by="hu_pops",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
+  sc10x.epi$hu_pops <- sc10x.epi$lin
+}
 
 Idents(sc10x) <- "lin"
-sc10x$pops <- Idents(sc10x)
-levels(sc10x$pops) <- c(levels(sc10x$pops),levels(sc10x.EpiFMSt$hu_pops))
-sc10x$pops[names(sc10x.EpiFMSt$hu_pops)] <- sc10x.EpiFMSt$hu_pops
+Idents(sc10x.epi) <- "hu_pops"
+Idents(sc10x.fmst) <- "hu_pops"
 
 if (!dir.exists(paste0("./analysis/vis/singler"))){
   dir.create(paste0("./analysis/vis/singler"))
@@ -245,39 +153,35 @@ postscript(paste0("./analysis/vis/singler/UMAP_all.lin.eps"))
 print(plot)
 dev.off()
 
-Idents(sc10x) <- "pops"
-Idents(sc10x.EpiFMSt) <- "hu_pops"
-Idents(sc10x.epi) <- "hu_pops"
-Idents(sc10x.fmst) <- "hu_pops"
-#Idents(sc10x.st) <- "pops"
-#Idents(sc10x.leu) <- "leu"
-
-plot <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
-postscript(paste0("./analysis/vis/singler/UMAP_all.eps"))
-print(plot)
-dev.off()
-
-for (i in c("EpiFMSt","epi","fmst")){
+for (i in c("epi","fmst")){
   plot <- DimPlot(get(paste0("sc10x.",i)),reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
   postscript(paste0("./analysis/vis/singler/UMAP_",i,".eps"))
   print(plot)
   dev.off()
 }
 
-plot <- plotScoreHeatmap(singler.lin,show.labels=TRUE,annotation_col=data.frame(ID=sc10x$lin,row.names=rownames(singler.lin)))
+labs <- singler.lin$labels
+labs <- replace(labs,labs=="Epithelial cells","Epi")
+labs <- replace(labs,labs %in% c("Stromal cells","Fibroblasts"),"FMSt")
+labs <- replace(labs,labs == "Endothelial cells","Endo")
+labs <- replace(labs,labs %in% c("Macrophages","Monocytes","B cells","DC","Eosinophils","Neutrophils","T cells","ILC","NK cells","Basophils","Mast cells","Tgd","NKT","B cells, pro","Microglia"),"Leu")
+plot <- plotScoreHeatmap(singler.lin,show.labels=FALSE,annotation_col=data.frame(ID=labs,row.names=rownames(singler.lin)))
 postscript(paste0("./analysis/vis/singler/ScoreHeatmap_lin.norm.eps"))
 print(plot)
 dev.off()
+rm(labs)
 
-plot <- plotScoreHeatmap(singler.EpiFMSt,show.labels=TRUE,show.pruned=TRUE,annotation_col=data.frame(first.label=singler.EpiFMSt$first.labels,row.names=rownames(singler.EpiFMSt)))
-postscript(paste0("./analysis/vis/singler/ScoreHeatmap_EpiFMSt.norm.eps"))
-print(plot)
-dev.off()
-
-plot <- plotScoreHeatmap(singler.EpiFMSt,show.labels=TRUE,show.pruned=TRUE,annotation_col=data.frame(first.label=singler.EpiFMSt$first.labels,row.names=rownames(singler.EpiFMSt)),normalize=FALSE)
-postscript(paste0("./analysis/vis/singler/ScoreHeatmap_EpiFMSt.eps"))
-print(plot)
-dev.off()
+if (opt$r == "epi"){
+  plot <- plotScoreHeatmap(singler.epi,show.labels=FALSE,show.pruned=FALSE,annotation_col=data.frame(ID=singler.epi$labels,row.names=rownames(singler.epi)))
+  postscript(paste0("./analysis/vis/singler/ScoreHeatmap_epi.norm.eps"))
+  print(plot)
+  dev.off()
+} else if (opt$r == "fmst"){
+plot <- plotScoreHeatmap(singler.fmst,show.labels=FALSE,show.pruned=FALSE,annotation_col=data.frame(ID=singler.fmst$labels,row.names=rownames(singler.fmst)))
+  postscript(paste0("./analysis/vis/singler/ScoreHeatmap_fmst.norm.eps"))
+  print(plot)
+  dev.off()
+}
 
 rm(list=ls(pattern="^sc10x.se"))