library("optparse")

option_list <- list(
	make_option(c("-i", "--indir"), help="indir"),
	make_option(c("-o", "--outdir"), help="outdir", default="."),
	make_option(c("--he_png"), help="he png file"),
	make_option(c("--res"), help="resolution, default %default",type="double",default=0.5),
	make_option(c("--MinCell"), help="MinCell, default %default",type="double",default=5),
	make_option(c("--MinFeatures"), help="MinFeatures, default %default",type="double",default=100),
	make_option(c("--point_size"), help="point_size, default %default",type="double",default=1),
	make_option(c("--ref"), help="1 for BlueprintEncodeData;2 for DatabaseImmuneCellExpressionData;3 for NovershternHematopoieticData;4 for MonacoImmuneData;5 for HumanPrimaryCellAtlasData;6 for ImmGenData;7 for MouseRNAseqData,default %default",type="double",default=1)
)
opt <- parse_args(OptionParser(option_list=option_list))
mode(opt$res)
#######################################
# 首先将需要用到的包load进工作环境
suppressMessages({
  library(Seurat)
  library(dplyr)
  library(tibble)
  library(ggplot2)
  library(tiff)
  library(broom)
  library(grid)
  library(gridExtra)
  library(magick)
  library(purrr)
  library(cowplot)
  library(figpatch)
  library(pdftools)
  library(cluster)
  library(ggpubr)
  library(SingleR)
  library(celldex)
})

# 记录运行时间
StartTime=proc.time()


##绘图函数
#输出ggplot对象
SavePlot <- function(od, filename, data, width = 8, height = 8){
  file.png <- paste(filename, "png", sep = ".")
  file.pdf <- paste(filename, "pdf", sep = ".")
  ggsave(filename = file.path(od, file.png), plot = data, width = width, height = height)
  ggsave(filename = file.path(od, file.pdf), plot = data, width = width, height = height)
}



# 参数初始化
#选择不同的resolution值可以获得不同的cluster数目，值越大cluster数目越多，默认值是0.5
res=opt$res
FilePath=opt$indir
he_file=opt$he_png
#过滤参数：细胞数5，基因数100
MinCell=opt$MinCell
MinFeatures=opt$MinFeatures

#输出目录
outdir=opt$outdir
if(!dir.exists(outdir)){dir.create(path=outdir)}


# 数据导入，给定输入文件路径，路径下需包括matrix.mtx,genes.tsv(或 features.tsv),和barcodes.tsv文件
expr=Read10X(FilePath,cell.column = 1)

#创建Seurat对象，数据集中测到的少于200个基因的位点（min.features = 200）和少于5个位点覆盖的基因（min.cells = 3）被过滤掉
object=CreateSeuratObject(counts = expr,assay = "Spatial",min.cells=MinCell,min.features=MinFeatures)
#标准化，鉴定出表达量高变的基因
#object <- SCTransform(object, assay = "Spatial", verbose = FALSE,return.only.var.genes = FALSE)
object <- NormalizeData(object, normalization.method = "LogNormalize", scale.factor = 10000)
object <- FindVariableFeatures(object, selection.method = "vst", nfeatures = 2000)
object <- ScaleData(object)
object1 <- RunPCA(object=object,pc.genes = VariableFeatures(object))
object1 <- FindNeighbors(object=object1, reduction = "pca",dims = 1:30,verbose = F)
object1 <- FindClusters(object=object1, verbose = F,resolution = res)
object1 <- RunUMAP(object=object1, reduction = "pca",dims = 1:30)

#整理聚类结果
#Cluster=object1@meta.data %>%
#  select(.,seurat_clusters) %>%
#  rownames_to_column(.,var = "Barcode") %>%
#  arrange(.,seurat_clusters)

#输出聚类结果
#write.csv(Cluster,file.path(outdir,"cluster.csv"),quote = F,row.names = F)

# 绘制UMAP聚类图和TSNE聚类图
umapPlot=UMAPPlot(object1, label = F)
SavePlot(outdir,"UMAP",umapPlot)

# 计录聚类时间
EndTime=proc.time()
RunTime <- EndTime - StartTime
print(paste0('ProcessingTime: ',RunTime[3][[1]],'seconds'))

#HumanPrimaryCellAtlasData、BlueprintEncodeData、DatabaseImmuneCellExpressionData、MonacoImmuneData、NovershternHematopoieticData
#ImmGenData、MouseRNAseqData
#mouseRNA <- HumanPrimaryCellAtlasData()
if(opt$ref==1){
	mouseRNA <- BlueprintEncodeData()
}else if(opt$ref==2){
	mouseRNA <- DatabaseImmuneCellExpressionData()
}else if(opt$ref==3){
	mouseRNA <- NovershternHematopoieticData()
}else if(opt$ref==4){
	mouseRNA <- MonacoImmuneData()
}else if(opt$ref==5){
        mouseRNA <- HumanPrimaryCellAtlasData()
}else if(opt$ref==6){
        mouseRNA <- ImmGenData()
}else if(opt$ref==7){
        mouseRNA <- MouseRNAseqData()
}

sce_for_SingleR<-GetAssayData(object1,slot="data")
clusters=object1@meta.data$seurat_clusters
pred.mouseRNA <- SingleR(test=sce_for_SingleR,ref=mouseRNA,labels=mouseRNA$label.fine,clusters=clusters,assay.type.test="logcounts",assay.type.ref ="logcounts")
cellType=data.frame(clusterID=levels(object1@meta.data$seurat_clusters),celltype=pred.mouseRNA$labels)
object1@meta.data$singleR=cellType[match(clusters,cellType$clusterID),"celltype"]
singleR_umapPlot=UMAPPlot(object1, label = F, group.by = c("singleR"))
SavePlot(outdir,"singleR_UMAP",singleR_umapPlot)

Cluster=object1@meta.data %>%
  select(.,c(seurat_clusters,singleR)) %>%
  rownames_to_column(.,var = "Barcode") %>%
  arrange(.,seurat_clusters)
write.csv(Cluster,file.path(outdir,"cluster.csv"),quote = F,row.names = F)

he_fig=png::readPNG(he_file)
width=ncol(he_fig)
hetght=nrow(he_fig)
cal_zoom_rate = function(width, height){
	std_width = 1000
    	std_height = std_width / (46 * 31) * (46 * 36 * sqrt(3) / 2.0)
    	if(std_width / std_height > width / height){
        	scale = width / std_width
	}
    	else{
        	scale = height / std_height
	}
    	return(scale)
}

zoom_scale = cal_zoom_rate(ncol(he_fig), nrow(he_fig))
#zoom_scale=1
bc_pos_file = gzfile(paste(FilePath, "barcodes_pos.tsv.gz", sep  = "/"),'rt')
bc_pos = read.table(bc_pos_file, header = FALSE, sep = '\t', quote = '')
names(bc_pos)=c("Barcode","x","y")
bc_pos=merge(bc_pos,Cluster,by="Barcode",all=FALSE)
bc_pos$x=bc_pos$x * zoom_scale
bc_pos$y=bc_pos$y * zoom_scale

point_size=opt$point_size
col =  c("#F56867","#FEB915","#C798EE","#59BE86","#7495D3","#D1D1D1","#6D1A9C","#15821E","#3A84E6","#997273","#787878","#DB4C6C","#9E7A7A","#554236","#AF5F3C","#93796C","#F9BD3F","#DAB370","#877F6C","#268785","#ed1299", "#09f9f5", "#246b93", "#cc8e12", "#d561dd", "#c93f00", "#ddd53e","#4aef7b", "#e86502", "#9ed84e", "#39ba30", "#6ad157", "#8249aa", "#99db27", "#e07233", "#ff523f","#ce2523", "#f7aa5d", "#cebb10", "#03827f", "#931635", "#373bbf", "#a1ce4c", "#ef3bb6", "#d66551","#1a918f", "#ff66fc", "#2927c4", "#7149af" ,"#57e559" ,"#8e3af4" ,"#f9a270" ,"#22547f", "#db5e92","#edd05e", "#6f25e8", "#0dbc21", "#280f7a", "#6373ed", "#5b910f" ,"#7b34c1" ,"#0cf29a" ,"#d80fc1")

bc_pos %>%
    group_by(singleR) %>%
    do(model = pam(.[c('x', 'y')], 1, stand = T)) %>%
    ungroup() %>% group_by(singleR) %>%
    do(map_df(.$model, broom::tidy)) %>% ### 整理模型数据
    ungroup() %>% select(singleR, x,y) %>%
    data.frame() %>%
    dplyr::rename(x.center = x, y.center = y, singleR = singleR) -> label.data
singleR_clstr_plot<-ggplot(data=bc_pos,aes(x ,dim(he_fig)[1]-y)) +
    background_image(he_fig)+
    geom_point(aes(colour = singleR),size=point_size) +
    scale_color_manual(values = col)+
    theme_bw() +
    theme(plot.title = element_text(face = 2,size = 50,hjust = 0.5)) +
    #geom_label(data = label.data, aes(label = singleR,x = x.center,y = dim(he_fig)[1]-y.center))+
    theme(axis.ticks = element_blank(), axis.text.y = element_blank(),panel.border = element_blank(), axis.text.x = element_blank()) +
    xlab('')+ylab('')+
    coord_cartesian(xlim = c(0, dim(he_fig)[2]), ylim = c(0, dim(he_fig)[1]),
                        expand = FALSE)+
    guides(colour = guide_legend(override.aes = list(size=3.5)))+
    theme(legend.position="right")

ggsave(filename = paste(outdir,"singleR_cluster.pdf",sep="/"), plot = singleR_clstr_plot, width = 10, height = 8)
ggsave(filename = paste(outdir,"singleR_cluster.png",sep="/"), plot = singleR_clstr_plot, width = 10, height = 8)
#pdf(paste(outdir,"dimplot_singleR_umap.pdf",sep="/"),width=7,height=5)
singleR_umap <- DimPlot(object1, reduction = "umap", label = FALSE, pt.size = 0.5,group.by=c('singleR'))
singleR_cluster_umap<-plot_grid(singleR_umap,singleR_clstr_plot)
ggsave(singleR_cluster_umap, filename = paste(outdir,"singleR_cluster_umap.png",sep="/"),height=8,width=16)
ggsave(singleR_cluster_umap, filename = paste(outdir,"singleR_cluster_umap.pdf",sep="/"),height=8,width=16)



