gc()
library(methods)
library(optparse)
library(Seurat)
library(readr)
library(readxl)
library(fBasics)
library(pastecs)
library(qusage)
library(RColorBrewer)
library(monocle)
library(dplyr)
library(viridis)
library(gridExtra)

options(bitmapType="cairo")

setwd("../")

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

project.name="PdPbPc"

scFolders()

sc10x <- scLoad(p=project.name)

results <- scQC(sc10x,lg=500,hg=2500,hm=0.1,sub=FALSE)
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
counts.gene.raw <- results[[3]]
counts.cell.filtered <- results[[4]]
counts.gene.filtered <- results[[5]]
rm(results)

results <- scCellCycle(sc10x,sub=FALSE)
sc10x <- results[[1]]
genes.s <- results[[2]]
genes.g2m <- results[[3]]
rm(results)

sc10x.l <- list()
sc10x.l[["Pd"]] <- subset(sc10x,subset= Donor=="Donor")
sc10x.l[["Pb"]] <- subset(sc10x,subset= BPH=="BPH")
sc10x.l[["Pc"]] <- subset(sc10x,subset= Cancer=="Cancer")

sc10x <- scCCA(sc10x.l)
rm(sc10x.l)

#gc()
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","S.Score","G2M.Score"),do.par=TRUE,num.cores=45,verbose=FALSE)
#sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),do.par=TRUE,num.cores=45,verbose=FALSE)
#gc()

results <- scPC(sc10x,pc=50,hpc=0.9,file="pre.stress",print="2",cca=TRUE)
sc10x <- results[[1]]
pc.use.prestress <- results[[2]]
rm(results)

sc10x <- scCluster(sc10x,res=0.5,red="pca",dim=pc.use.prestress,print="2",folder="pre.stress")

genes.stress <- read_delim("./genesets/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
genes.stress <- genes.stress[1]
colnames(genes.stress) <- "scDWS.Stress"

results <- scScore(sc10x,score="Stress",geneset=as.list(genes.stress),cut.pt=0.9,anchor=c("EGR1","FOS","JUN"))
sc10x.preStress <- results[[1]]
sc10x <- results[[2]]
rm(results)

#gc()
sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito","S.Score","G2M.Score"),do.par=TRUE,num.cores=45,verbose=FALSE)
#sc10x <- ScaleData(object=sc10x,vars.to.regress=c("nFeature_RNA","percent.mito"),do.par=TRUE,num.cores=45,verbose=FALSE)
#gc()

results <- scPC(sc10x,pc=50,hpc=0.9,file="post.stress",print="2",cca=TRUE)
sc10x <- results[[1]]
pc.use.poststress <- results[[2]]
rm(results)

res <- c(seq(0.1,0.5,0.1),0.75,seq(1,5,1))
sc10x <- scCluster(sc10x,res=res,red="pca",dim=pc.use.poststress,print="2",folder="ALL")

# gene.set1 <- read_delim("./genesets/DEG_Epi_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
# gene.set1 <- as.list(gene.set1)
# names(gene.set1) <- "Epi"
# gene.set <- c(gene.set1)
# gene.set1 <- read_delim("./genesets/DEG_FMSt_5FC.txt","\t",escape_double=FALSE,trim_ws=TRUE,col_names=FALSE)
# gene.set1 <- as.list(gene.set1)
# names(gene.set1) <- "St"
# gene.set <- c(gene.set,gene.set1)

gene.set1 <- read_delim("./genesets/genes.deg.BE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "BE"
#gene.set <- c(gene.set,gene.set1)
gene.set <- c(gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.LE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "LE"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE1.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Club"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.OE2.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Hillock"
gene.set <- c(gene.set,gene.set1)

gene.set1 <- read_delim("./genesets/genes.deg.NE.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "NE"
gene.set <- c(gene.set,gene.set1)

gene.set1 <- read_delim("./genesets/genes.deg.Endo.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Endo"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.SM.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "SM"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.Fib.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Fib"
gene.set <- c(gene.set,gene.set1)
gene.set1 <- read_delim("./genesets/genes.deg.Leu.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
gene.set1 <- gene.set1[1]
gene.set1 <- as.list(gene.set1)
names(gene.set1) <- "Leu"
gene.set <- c(gene.set,gene.set1)

# genes.leu <- read_excel("./genesets/40425_2017_215_MOESM1_ESM.xlsx",sheet="S3. Candidate markers")
# leu <- as.list(unique(genes.leu[,2]))$Cell
# leu <- leu[-c(9,17,20,24,26,28)]
# #leu.l <- leu[-c(1,3:4,7:8,13:18,21,20:30)]
# #leu.lin <- c("Myeloid","Lymphoid","Lymphoid","Lymphoid","Myeloid","Myeloid","Lymphoid","Myeloid","Myeloid","Myeloid","Myeloid","Myeloid","Lymphoid","Lymphoid","Lymphoid","Myeloid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid","Lymphoid")
# #leu.lin <- c("Lymphoid","Myeloid","Myeloid","Myeloid","Myeloid","Myeloid","Myeloid","Lymphoid")  
# genes.leu <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu),]
# #genes.leu.l <- genes.leu[unlist(genes.leu[,2]) %in% unlist(leu.l),]
# #genes.leu.l[nrow(genes.leu.l)+1,] <- list("MS4A7","Macrophages",0)
# #genes.leu.l[nrow(genes.leu.l)+1,] <- list("CD14","Macrophages",0)
# # genes.leu.l[nrow(genes.leu.l)+1,] <- list("CD86","Macrophages",0)
# # genes.leu.l[nrow(genes.leu.l)+1,] <- list("S100A9","Neutrophils",0)
# # genes.leu.l[nrow(genes.leu.l)+1,] <- list("EREG","Neutrophils",0)
# # genes.leu.l[nrow(genes.leu.l)+1,] <- list("S100A8","Neutrophils",0)
# # genes.leu.l[nrow(genes.leu.l)+1,] <- list("FCN1","Neutrophils",0)
# gene.set1 <- split(genes.leu[,1], genes.leu[,2])
# #gene.set1 <- split(genes.leu.l[,1], genes.leu.l[,2])
# gene.set1 <- lapply(gene.set1,unname)
# gene.set1 <- lapply(gene.set1,unlist)
# #names(gene.set1) <- leu.lin
# names(gene.set1) <- leu
# gene.set <- c(gene.set,gene.set1)
rm(gene.set1)

min <- min(table(sc10x$integrated_snn_res.0.5))
results <- scQuSAGE(sc10x,gs=gene.set,save=TRUE,type="lg",id="integrated_snn_res.0.5",ds=0,nm="all.pops",print="2")
sc10x <- results[[1]]
results.cor.all.pops <- results[[2]]
results.clust.all.pops.id <- results[[3]]
rm(results)
rm(gene.set)

save(sc0x,file=paste0("./analysis/",project.name,".RData"))
save.image(file="./analysis/Data.RData")