gc()
library(optparse)
library(Seurat)
library(readr)
library(dplyr)
library(autothresholdr)
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("--m",action="store",default="Patient",type='character',help="Group to merge on")
)
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")

project.name=opt$p

scFolders()

results <- scLoad(p=project.name,cellranger=4,aggr=FALSE,ncell=0,nfeat=0,seurat=TRUE)
sc10x <- results[[1]]
sc10x.groups <- results[[2]]
rm(results)

sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$ALL==1,1]))]

results <- scQC(sc10x,sp=opt$s,feature=c("nFeature_RNA"))
sc10x <- results[[1]]
counts.cell.raw <- results[[2]]
counts.gene.raw <- results[[3]]
counts.cell.filtered.1umi <- results[[4]]
counts.gene.filtered.1umi <- results[[5]]
rm(results)

results <- scQC(sc10x,sp=opt$s,feature="percent.mito")
sc10x <- results[[1]]
counts.cell.filtered.2mito <- results[[4]]
counts.gene.filtered.2mito <- results[[5]]
rm(results)

#results <- scQC(sc10x,sp=opt$s,feature="percent.ribo")
#sc10x <- results[[1]]
#counts.cell.filtered.3ribo <- results[[4]]
#counts.gene.filtered.3ribo <- results[[5]]
#rm(results)

#results <- scQC(sc10x,sp=opt$s,feature="nCount_RNA")
#sc10x <- results[[1]]
#counts.cell.filtered.3gene <- results[[4]]
#counts.gene.filtered.3gene <- results[[5]]
#rm(results)

counts.cell.filtered <- counts.cell.filtered.2mito
counts.gene.filtered <- counts.gene.filtered.2mito

sc10x <- sc10x[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))]

# pc.use.prestress <- list()
# for (i in names(sc10x)){
#   sc10x.temp <- sc10x[[i]]
#   sc10x.temp <- SCTransform(sc10x.temp,vars.to.regress=c("nFeature_RNA","percent.mito"),verbose=FALSE,assay="RNA")
#   if (ncol(sc10x.temp) > 100) {
#     pc.calc <- 100
#   } else if (ncol(sc10x.temp) <= 100) {
#     pc.calc <- ncol(sc10x.temp)-1
#   }
#   results <- scPC(sc10x.temp,pc=pc.calc,hpc=0.9,file=paste0(i,".pre.stress"),print="umap",assay="SCT")
#   sc10x.temp <- results[[1]]
#   pc.use.prestress.temp <- results[[2]]
#   rm(results)
#   sc10x.temp <- scCluster(sc10x.temp,res=0.5,red="pca",dim=pc.use.prestress.temp,print="umap",folder=paste0("preStress/",i),assay="SCT")
#   sc10x[i] <- sc10x.temp
#   pc.use.prestress[i] <- pc.use.prestress.temp
#   rm(sc10x.temp)
#   rm(pc.calc)
#   rm(pc.use.prestress.temp)
# }
# rm(i)

merges <- NULL
# for (i in names(counts.cell.destress[as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1]))])){
#   if (counts.cell.destress[[i]]<750){
#     merges <- c(merges,i)
#   }
# }
# rm(i)
# if (length(merges)>1){
#   sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
#   sc10x[merges] <- NULL
#   sc10x$merge <- sc10x.merge
#   rm(sc10x.merge)
# } else if (length(merges)==1){
#   merges <- c(merges,names(counts.cell.destress[counts.cell.destress == min(sapply(counts.cell.destress[setdiff(as.character(unlist(sc10x.groups[sc10x.groups$Keep==1,1])),merges)],min))]))
#   sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
#   sc10x[merges] <- NULL
#   sc10x$merge <- sc10x.merge
# }
# rm(merges)
merges <- NULL
for (i in unique(sc10x.groups$Patient[sc10x.groups$Keep==1])){
  for (j in sc10x.groups$Samples[sc10x.groups$Keep==1]){
    if (sc10x.groups[sc10x.groups$Samples==j,][opt$merge]==i){
      merges <- c(merges,j)
    }
  }
  if (length(merges)>1){
  sc10x.merge <- merge(x=sc10x[[merges[1]]],y=sc10x[merges[-1]])
  } else {
    sc10x.merge <- sc10x[[merges[1]]]
  }
  sc10x[merges] <- NULL
  sc10x[i] <- sc10x.merge
  merges <- NULL
  rm(sc10x.merge)
}
rm(merges)
rm(i)
rm(j)

gc()
if (length(sc10x)>1){
  sc10x <- scAlign(sc10x)
  sc10x@project.name <- project.name
  sc10x$ALL <- "ALL"
}
gc()

if (ncol(sc10x) > 1000) {
  pc.calc <- 1000
} else if (ncol(sc10x) > 500) {
  pc.calc <- 500
} else if (ncol(sc10x) > 100) {
  pc.calc <- ncol(sc10x)-1
}
results <- scPC(sc10x,pc=pc.calc,hpc=0.9,file="ALL",print="umap")
sc10x <- results[[1]]
pc.use.poststress <- results[[2]]
rm(results)
rm(pc.calc)

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="umap",folder="ALL")

if (opt$s == "hu"){
  genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/genes.deg.Stress.csv",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
  genes.stress <- genes.stress[1]
  colnames(genes.stress) <- "scDWS.Stress"
  anchor <- c("EGR1","FOS","JUN")
} else if (opt$s == "mu"){
  genes.stress <- read_delim("/work/urology/ghenry/RNA-Seq/SingleCell/ANALYSIS/REF/stress/van.den.Brink1.txt",",",escape_double=FALSE,trim_ws=TRUE,col_names=TRUE)
  genes.stress <- genes.stress[1]
  colnames(genes.stress) <- "van.den.Brink1.Stress"
  anchor <- c("Egr1","Fos","Jun")
}

results <- scScore(list(ALL=sc10x),score="Stress",geneset=as.list(genes.stress),cut.pt="renyi",anchor=anchor)
sc10x.preStress <- results[[1]]
sc10x <- results[[2]]
#sc10x <- results[[1]]$ALL
rm(results)
#counts.cell.destress <- lapply(sc10x,ncol)
counts.cell.destress <- ncol(sc10x)
rm(anchor)

DefaultAssay(object=sc10x) <- "SCT"

sc10x@meta.data <- sc10x@meta.data[,c("samples","integrated_snn_res.0.1","integrated_snn_res.0.2","integrated_snn_res.0.3","integrated_snn_res.0.4","integrated_snn_res.0.5","integrated_snn_res.0.75","integrated_snn_res.1","integrated_snn_res.2","integrated_snn_res.3","integrated_snn_res.4","integrated_snn_res.5","nCount_RNA","nFeature_RNA","percent.mito","percent.ribo","Stress1")]

saveRDS(sc10x,paste0("./analysis/",project.name,"_raw.rds"))

library(sceasy)
library(reticulate)
use_condaenv('sceasy')
convertFormat(sc10x,from="seurat",to="anndata",outFile=paste0("/project/urology/Strand_lab/shared/cellxgene/anndata/",project.name,"_raw.h5ad"),assay="SCT",main_layer="scale.data")
saveRDS(sc10x,paste0("/project/urology/Strand_lab/shared/cellxgene/seurat/",project.name,"_raw.rds"))

rm(list=ls(pattern="sc10x"))
save.image(file="./analysis/sc10x.raw.RData")