gc()
library(optparse)
library(Seurat)
library(SingleR)
library(scRNAseq)
library(RColorBrewer)
library(viridis)
library(gridExtra)
library(ggplot2)

options(future.globals.maxSize= 1000000000000)

option_list=list(
  make_option("--p",action="store",default="huPr_PdPb",type='character',help="Project Name"),
  make_option("--s",action="store",default="hu",type='character',help="Species"),
  make_option("--o",action="store",default="pr",type='character',help="Organ")
)
opt=parse_args(OptionParser(option_list=option_list))
rm(option_list)


setwd("../")

source("./r.scripts/sc-TissueMapper_functions.R")
source("./r.scripts/sc-TissueMapper_process.R")

if (opt$s == "hu"){
  load("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/Travaglini_Nature2020/assets/Analysis/r objects/droplet_normal_lung_blood_seurat_ntiss10x.P1.anno.20191002.RC4.Robj")
  lin.P1 <- UpdateSeuratObject(ntiss10x.P1.anno)
  lin.P1$label.main <- lin.P1$magnetic.selection
  lin.P1$label.fine <- lin.P1$free_annotation
  lin.P1@meta.data <- lin.P1@meta.data[,c("label.main","label.fine")]
  lin.P1 <- as.SingleCellExperiment(lin.P1,assay="RNA")
  rm(ntiss10x.P1.anno)
  load("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/Travaglini_Nature2020/assets/Analysis/r objects/droplet_normal_lung_seurat_ntiss10x.P2.anno.20191002.RC4.Robj")
  lin.P2 <- UpdateSeuratObject(ntiss10x.P2.anno)
  lin.P2$label.main <- lin.P2$magnetic.selection
  lin.P2$label.fine <- lin.P2$free_annotation
  lin.P2@meta.data <- lin.P2@meta.data[,c("label.main","label.fine")]
  lin.P2 <- as.SingleCellExperiment(lin.P2,assay="RNA")
  rm(ntiss10x.P2.anno)
  load("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/Travaglini_Nature2020/assets/Analysis/r objects/droplet_normal_lung_blood_seurat_ntiss10x.P3.anno.20191002.RC4.Robj")
  lin.P3 <- UpdateSeuratObject(ntiss10x.P3.anno)
  lin.P3$label.main <- lin.P3$magnetic.selection
  lin.P3$label.fine <- lin.P3$free_annotation
  lin.P3@meta.data <- lin.P3@meta.data[,c("label.main","label.fine")]
  lin.P3 <- as.SingleCellExperiment(lin.P3,assay="RNA")
  rm(ntiss10x.P3.anno)
  
  pop <- readRDS("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/huPr_Pd_all.rds")
  try(
    if (as.numeric(substring(pop@version,1,1))<3){
      pop <- UpdateSeuratObject(pop)
    }
  )
  Idents(pop) <- "Lineage"
  pop.epi <- subset(pop,idents="Epi")
  Idents(pop.epi) <- "Population"
  pop.epi <- subset(pop.epi,idents=c("BE","LE","Hillock","Club"))
  pop.st <- subset(pop,idents=c("FMSt","Endo"))
  pop.epi$label.main <- pop.epi$Lineage
  pop.st$label.main <- pop.st$Lineage
  pop.epi$label.fine <- pop.epi$Population
  pop.st$label.fine <- pop.st$Population
  pop.epi@meta.data <- pop.epi@meta.data[,c("label.main","label.fine")]
  pop.st@meta.data <- pop.st@meta.data[,c("label.main","label.fine")]
  pop.epi <- as.SingleCellExperiment(pop.epi,assay="RNA")
  pop.st <- as.SingleCellExperiment(pop.st,assay="RNA")
  rm(pop)
} else if (opt$s == "mu"){
  load("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/Travaglini_Nature2020/assets/Analysis/r objects/mouse_droplet_TMS_tiss10x.mouse.anno.ncbi.20200222.RC4.Robj")
  lin <- UpdateSeuratObject(tiss10x.mouse.anno)
  lin$label.main <- lin$free_annotation
  lin$label.fine <- lin$free_annotation
  lin@meta.data <- lin@meta.data[,c("label.main","label.fine")]
  lin <- as.SingleCellExperiment(lin,assay="RNA")
  rm(tiss10x.mouse.anno)
}

sc10x <- readRDS(paste0("./analysis/",opt$p,"_raw.rds"))

sc10x.se <- as.SingleCellExperiment(sc10x)
if (opt$s == "hu"){
  common <- Reduce(intersect, list(rownames(sc10x.se),rownames(lin.P1),rownames(lin.P2),rownames(lin.P3)))
  lin.P1 <- lin.P1[common,]
  lin.P2 <- lin.P1[common,]
  lin.P3 <- lin.P3[common,]
} else if (opt$s == "mu"){
  common <- Reduce(intersect, list(rownames(sc10x.se),rownames(lin)))
  lin <- lin[common,]
}
sc10x.se <- sc10x.se[common,]
rm(common)

if (opt$s == "hu"){
  ref <- list(lung.P1=lin.P1,lung.P2=lin.P2,lung.P3=lin.P3)
  labels <- list(lung.P1=lin.P1$label.fine,lung.P2=lin.P2$label.fine,lung.P3=lin.P3$label.fine)
  rm(lin.P1)
  rm(lin.P2)
  rm(lin.P3)
} else if (opt$s == "mu"){
  ref <- lin
  labels <- lin$label.fine
  rm(lin)
}
  
#singler <- SingleR(sc10x.se,ref=ref,method="single",labels=labels,de.method="wilcox",BPPARAM=BiocParallel::MulticoreParam(workers=20))
#labs <- singler$labels
singler <- SingleR(sc10x.se,ref=ref,clusters=sc10x.se$integrated_snn_res.1,labels=labels,de.method="wilcox",BPPARAM=BiocParallel::MulticoreParam(workers=20))
labs <- singler$labels[match(sc10x.se$integrated_snn_res.1,singler@rownames)]
sc10x$pop <- labs
labs.raw <- labs
labs[labs %in% c("Differentiating Basal","Proliferating Basal","Basal","Mesothelial","Alveolar Epithelial Type 2","Club","Alveolar Epithelial Type 1","Ciliated","Signaling Alveolar Epithelial Type 2","Proximal Basal","Neuroendocrine","Mucous","Ionocyte","Serous","Proximal Ciliated","Goblet")] <- "Epithelia"
labs[labs %in% c("Adventitial Fibroblast","Alveolar Fibroblast","Lipofibroblast")] <- "Fibroblast"
labs[labs %in% c("Airway Smooth Muscle","Pericyte","Myofibroblast","Fibromyocyte","Vascular Smooth Muscle")] <- "Smooth Muscle"
labs[labs %in% c("Capillary Aerocyte","Capillary","Bronchial Vessel 2","Vein","Artery","Lymphatic","Bronchial Vessel 1","Lympatic")] <- "Endothelia"
labs[labs %in% c("CD4+ Memory/Effector T","CD4+ Naive T","CD8+ Memory/Effector T","CD8+ Naive T","Natural Killer","B","Plasmacytoid Dendritic","Plasma","Proliferating NK/T","Natural Killer T","CD4+ T","CD8+ T","Proliferating T","Regulatory T","Ly6g5b+ T","Proliferating NK","Alox5+ Lymphocytes","Zbtb32+ B")] <- "Lymphoid"
labs[labs %in% c("Classical Monocyte","Nonclassical Monocyte","Myeloid Dendritic Type 2","IGSF21+ Dendritic","EREG+ Dendritic","Myeloid Dendritic Type 1","TREM2+ Dendritic","Macrophage","Proliferating Macrophage","Proliferating Classical Monocyte","Intermediate Monocyte","Ccr7+ Dendritic","Proliferating Dendritic","Alveolar Macrophage","Interstitial Macrophage","Proliferating Alveolar Macrophage")] <- "Myeloid"
labs[labs %in% c("Basophil/Mast 1","Basophil/Mast 2","Basophil","Neutrophil")] <- "Granulocyte"
sc10x$lin <- labs
sc10x$lung <- labs.raw
#DimPlot(sc10x,group.by="lin",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

Idents(sc10x) <- "lin"
sc10x.epi <- subset(sc10x, idents="Epithelia")
sc10x.fmst <- subset(sc10x, idents=unique(Idents(sc10x))[unique(Idents(sc10x)) %in% c("Fibroblast","Smooth Muscle")])
sc10x.fib <- subset(sc10x, idents="Fibroblast")
sc10x.sm <- subset(sc10x, idents="Smooth Muscle")
sc10x.leu <- subset(sc10x, idents=unique(Idents(sc10x))[unique(Idents(sc10x)) %in% c("Lymphoid", "Myeloid", "Granulocyte")])

res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))

if (ncol(sc10x.epi) > 1000) {
  pc.calc.epi <- 1000
} else if (ncol(sc10x.epi) > 500) {
  pc.calc.epi <- 500
} else if (ncol(sc10x.epi) > 100) {
  pc.calc.epi <- ncol(sc10x.epi)-1
}
results <- scPC(sc10x.epi,pc=pc.calc.epi,hpc=0.9,file="epi",print="umap")
sc10x.epi <- results[[1]]
pc.use.epi <- results[[2]]
rm(results)
sc10x.epi <- scCluster(sc10x.epi,res=res,red="pca",dim=pc.use.epi,print="umap",folder="epi")
#DimPlot(sc10x.epi,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

if (ncol(sc10x.fmst) > 1000) {
  pc.calc.fmst <- 1000
} else if (ncol(sc10x.fmst) > 500) {
  pc.calc.fmst <- 500
} else if (ncol(sc10x.fmst) > 100) {
  pc.calc.fmst <- ncol(sc10x.fmst)-1
}
results <- scPC(sc10x.fmst,pc=pc.calc.fmst,hpc=0.9,file="fmst",print="umap")
sc10x.fmst <- results[[1]]
pc.use.fmst <- results[[2]]
rm(results)
sc10x.fmst <- scCluster(sc10x.fmst,res=res,red="pca",dim=pc.use.fmst,print="umap",folder="fmst")
#DimPlot(sc10x.fmst,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

if (ncol(sc10x.fib) > 1000) {
  pc.calc.fib <- 1000
} else if (ncol(sc10x.fib) > 500) {
  pc.calc.fib <- 500
} else if (ncol(sc10x.fib) > 100) {
  pc.calc.fib <- ncol(sc10x.fib)-1
}
results <- scPC(sc10x.fib,pc=pc.calc.fib,hpc=0.9,file="all",print="umap")
sc10x.fib <- results[[1]]
pc.use.fib <- results[[2]]
rm(results)
sc10x.fib <- scCluster(sc10x.fib,res=res,red="pca",dim=pc.use.fib,print="umap",folder="fib")
#DimPlot(sc10x.fib,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

if (ncol(sc10x.sm) > 1000) {
  pc.calc.sm <- 1000
} else if (ncol(sc10x.sm) > 500) {
  pc.calc.sm <- 500
} else if (ncol(sc10x.sm) > 100) {
  pc.calc.sm <- ncol(sc10x.sm)-1
}
results <- scPC(sc10x.sm,pc=pc.calc.sm,hpc=0.9,file="sm",print="umap")
sc10x.sm <- results[[1]]
pc.use.sm <- results[[2]]
rm(results)
sc10x.sm <- scCluster(sc10x.sm,res=res,red="pca",dim=pc.use.sm,print="umap",folder="sm")
#DimPlot(sc10x.sm,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

if (ncol(sc10x.leu) > 1000) {
  pc.calc.leu <- 1000
} else if (ncol(sc10x.leu) > 500) {
  pc.calc.leu <- 500
} else if (ncol(sc10x.leu) > 100) {
  pc.calc.leu <- ncol(sc10x.leu)-1
}
results <- scPC(sc10x.leu,pc=pc.calc.leu,hpc=0.9,file="leu",print="umap")
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="umap",folder="leu")
#DimPlot(sc10x.leu,group.by="integrated_snn_res.0.1",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

if (opt$o == "pr" && opt$s == "hu") {
  sc10x.se.epi <- as.SingleCellExperiment(sc10x.epi,assay="RNA")
  common <- intersect(rownames(sc10x.se.epi),rownames(pop.epi))
  sc10x.se.epi <- sc10x.se.epi[common,]
  pop.epi <- pop.epi[common,]
  rm(common)
  singler.epi <- SingleR(sc10x.se.epi,ref=pop.epi,method="single",labels=pop.epi$label.fine,de.method="wilcox",BPPARAM=BiocParallel::MulticoreParam(workers=20))
  sc10x.epi$pop <- singler.epi$labels
  Idents(sc10x.epi) <- "pop"
  sc10x[["pop"]][row.names(sc10x.epi[["pop"]])] <- sc10x.epi[["pop"]][,"pop"]
  #DimPlot(sc10x.epi,group.by="pop",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
} else if (opt$o == "pr" && opt$s == "mu") {
  sc10x.epi$pop <- sc10x.epi$lin
} else {
  sc10x.epi$pop <- sc10x.epi$lin
}

sc10x$pop <- sc10x$lin

sc10x$pop[names(sc10x.epi$pop)] <- sc10x.epi$pop
sc10x$pop[names(sc10x.leu$pop)] <- sc10x.leu$pop
#DimPlot(sc10x,group.by="pop",reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")

Idents(sc10x) <- "pop"
Idents(sc10x.epi) <- "pop"
Idents(sc10x.fmst) <- "lin"
Idents(sc10x.fib) <- "lin"
Idents(sc10x.sm) <- "lin"
Idents(sc10x.leu) <- "pop"

if (!dir.exists("./analysis/vis/singler/")){
  dir.create("./analysis/vis/singler/")
}

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("epi","fmst","fib","sm","leu")){
  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()
}

DefaultAssay(object=sc10x) <- "SCT"
DefaultAssay(object=sc10x.epi) <- "SCT"
DefaultAssay(object=sc10x.fmst) <- "SCT"
DefaultAssay(object=sc10x.fib) <- "SCT"
DefaultAssay(object=sc10x.sm) <- "SCT"
DefaultAssay(object=sc10x.leu) <- "SCT"
for (i in res){
  sc10x[[paste0("resolution_",i)]] <- sc10x[[paste0("integrated_snn_res.",i)]]
  sc10x.epi[[paste0("resolution_",i)]] <- sc10x.epi[[paste0("integrated_snn_res.",i)]]
  sc10x.fmst[[paste0("resolution_",i)]] <- sc10x.fmst[[paste0("integrated_snn_res.",i)]]
  sc10x.fib[[paste0("resolution_",i)]] <- sc10x.fib[[paste0("integrated_snn_res.",i)]]
  sc10x.sm[[paste0("resolution_",i)]] <- sc10x.sm[[paste0("integrated_snn_res.",i)]]
  sc10x.leu[[paste0("resolution_",i)]] <- sc10x.leu[[paste0("integrated_snn_res.",i)]]
}
sc10x@meta.data <- sc10x@meta.data[,c("lin","pop","lung","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.epi@meta.data <- sc10x.epi@meta.data[,c("lin","pop","lung","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.fmst@meta.data <- sc10x.fmst@meta.data[,c("lin","lung","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.fib@meta.data <- sc10x.fib@meta.data[,c("lin","lung","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.sm@meta.data <- sc10x.sm@meta.data[,c("lin","lung","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
sc10x.leu@meta.data <- sc10x.leu@meta.data[,c("lin","lung","pop","samples","resolution_0.1","resolution_0.2","resolution_0.3","resolution_0.4","resolution_0.5","resolution_0.75","resolution_1","resolution_2","resolution_3","resolution_4","resolution_5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]
saveRDS(sc10x,paste0("./analysis/",opt$p,"_id_all.rds"))
saveRDS(sc10x.epi,paste0("./analysis/",opt$p,"_id_epi.rds"))
saveRDS(sc10x.fmst,paste0("./analysis/",opt$p,"_id_fmst.rds"))
saveRDS(sc10x.fib,paste0("./analysis/",opt$p,"_id_fib.rds"))
saveRDS(sc10x.sm,paste0("./analysis/",opt$p,"_id_sm.rds"))
saveRDS(sc10x.leu,paste0("./analysis/",opt$p,"_id_leu.rds"))
library(sceasy)
library(reticulate)
use_condaenv('sceasy')
convertFormat(sc10x,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_all.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.epi,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_epi.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.fmst,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_fmst.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.fib,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_fib.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.sm,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_sm.h5ad"),assay="SCT",main_layer="scale.data")
convertFormat(sc10x.leu,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",opt$p,"_id_leu.h5ad"),assay="SCT",main_layer="scale.data")
saveRDS(sc10x,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_all.rds"))
saveRDS(sc10x.epi,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_epi.rds"))
saveRDS(sc10x.fmst,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_fmst.rds"))
saveRDS(sc10x.fib,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_fib.rds"))
saveRDS(sc10x.sm,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_sm.rds"))
saveRDS(sc10x.leu,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",opt$p,"_id_leu.rds"))

save(list=ls(pattern="^singler"),file='./analysis/singler_objects.RData')
rm(list=ls(pattern="^sc10x"))
rm(list=ls(pattern="^singler"))
save.image(file="./analysis/sc10x.id.RData")