#sc-TissueMapper
#Author: Gervaise H. Henry
#Email: gervaise.henry@utsouthwestern.edu
#Lab: Strand Lab, Deparment of Urology, University of Texas Southwestern Medical Center


scFolders <- function(){
  if (!dir.exists("./analysis/qc/")){
    dir.create("./analysis/qc/")
  }
  if (!dir.exists("./analysis/qc/cellcycle")){
    dir.create("./analysis/qc/cellcycle")
  }
  if (!dir.exists("./analysis/vis")){
    dir.create("./analysis/vis")
  }
  if (!dir.exists("./analysis/score_id")){
    dir.create("./analysis/score_id")
  }
  if (!dir.exists("./analysis/cor")){
    dir.create("./analysis/cor")
  }
}


scLoad <- function(p,mc=3,mg=200,sub=FALSE){
  #Load and prefilter filtered_gene_bc_matrices_mex output from cellranger
  
  #Inputs:
  #p = project name
  #mc = minimum cells (filter out genes which occur in fewer cells)
  #mg = minimum genes (filter out cells which conntain fewer genes)
  #sub = multi-dataset experiment  
  
  #Outputs:
  #Seurat object
  
  #Load filtered_gene_bc_matrices_mex output from cellranger
  if (sub==FALSE){
    sc10x.data <- Read10X(data.dir="./analysis/DATA/10x/filtered_gene_bc_matrices_mex/GRCh38/")
  } else {
    sc10x.data <- Read10X(data.dir=paste0("./analysis/DATA/",p,"/10x/filtered_gene_bc_matrices_mex/GRCh38/"))
  }
  sc10x <- new("seurat",raw.data=sc10x.data)
  
  #Label sample names from aggregation_csv.csv
  if (sub==FALSE){
    sc10x.aggr <- read_csv("./analysis/DATA/10x/aggregation_csv.csv")
  } else {
    sc10x.aggr <- read_csv(paste0("./analysis/DATA/",p,"/10x/aggregation_csv.csv"))
  }
  cell.codes <- as.data.frame(sc10x@raw.data@Dimnames[[2]])
  colnames(cell.codes) <- "barcodes"
  rownames(cell.codes) <- cell.codes$barcodes
  cell.codes$lib.codes <- as.factor(gsub(pattern=".+-",replacement="",cell.codes$barcodes))
  cell.codes$samples <- sc10x.aggr$library_id[match(cell.codes$lib.codes,as.numeric(rownames(sc10x.aggr)))]
  sc10x <- CreateSeuratObject(counts=sc10x.data,project=p,assay="RNA",min.cells=mc,min.features=mg,meta.data=cell.codes["samples"])
  
  #Create groups found in demultiplex.csv
  if (sub==FALSE){
    sc10x.demultiplex <- read_csv(paste0("./analysis/DATA/",p,"-demultiplex.csv"))
  } else {
    sc10x.demultiplex <- read_csv(paste0("./analysis/DATA/",p,"/",p,"-demultiplex.csv"))
  }
  for (i in 2:ncol(sc10x.demultiplex)){
    Idents(sc10x) <- "samples"
    merge.cluster <- apply(sc10x.demultiplex[,i],1,as.character)
    merge.cluster[merge.cluster==1] <- colnames(sc10x.demultiplex[,i])
    
    Idents(sc10x) <- plyr::mapvalues(x=Idents(sc10x),from=sc10x.demultiplex$Samples,to=merge.cluster)
    sc10x@meta.data[,colnames(sc10x.demultiplex[,i])] <- Idents(sc10x)
  }
  
  return(sc10x)
}


scSubset <- function(sc10x,i="ALL",g="ALL"){
  #Subset cells based on an identity
  
  #Inputs:
  #sc10x = seruat object
  #i = identity to use
  #g = group to subset by
  
  #Outputs:
  #Seurat object
  
  Idents(sc10x) <- i
  
  sc10x.sub <- subset(x=sc10x,idents=g)
  
  return(sc10x.sub)
}


scQC <- function(sc10x,lg=500,hg=2500,hm=0.1,sub=FALSE){
  #QC and filter Seurat object
  
  #Inputs:
  #sc10x = Seruat object
  #lg = threshold for cells with minimum genes
  #hg = threshold for cells with maximum genes
  #hm = threshold for cells with maximum %mito genes
  #sub = Subfolder to save output files
  
  #Outputs:
  #result[1] = Seurat object
  #result[2] = raw cell count
  #result[3] = raw gene count
  #result[4] = filtered cell count
  #result[5] = filtered gene count
  
  #Make sub-folders if necessary
  if (sub==FALSE){
    folder <- "./analysis/qc/"
  } else {
    folder <- paste0("./analysis/qc/",sub,"/")
    if (!dir.exists(folder)){
      dir.create(folder)
    }}
  
  #Calculate stats
  mito.genes <- grep(pattern="^MT-",x=rownames(x=GetAssayData(object=sc10x)),value=TRUE)
  percent.mito <- Matrix::colSums(GetAssayData(object=sc10x,slot="counts")[mito.genes,])/Matrix::colSums(GetAssayData(object=sc10x,slot="counts"))
  sc10x$percent.mito <- percent.mito
  
  #Plot raw stats
  for(x in c("ALL","samples")){
    Idents(object=sc10x) <- x
    postscript(paste0(folder,"Violin_qc.raw.",x,".eps"))
    plot1 <- VlnPlot(object=sc10x,features="nCount_RNA",pt.size=0.1)
    plot2 <- VlnPlot(object=sc10x,features="nFeature_RNA",pt.size=0.1)+geom_hline(yintercept=lg,size=1,color="red")+geom_hline(yintercept=hg,size=1,color="red")
    plot3 <- VlnPlot(object=sc10x,features="percent.mito",pt.size=0.1)+geom_hline(yintercept=hm,size=1,color="red")
    grid.arrange(plot1,plot2,plot3,nrow=1)
    dev.off()
  }
  postscript(paste0(folder,"Plot_qc.raw.eps"))
  density1 <- density(sc10x$nCount_RNA,sc10x$nFeature_RNA,n=1000)
  plot1 <- ggplot(data.frame(cbind(sc10x$nCount_RNA,sc10x$nFeature_RNA)))+geom_point(aes(x=X1,y=X2,color=density1),size=0.1)+scale_color_viridis(option="inferno")+geom_hline(yintercept=lg,size=1,color="red")+geom_hline(yintercept=hg,size=1,color="red")+labs(x="nUMI",y="nGenes",color="Count")+cowplot::theme_cowplot()
  density2 <- density(sc10x$nCount_RNA,sc10x$percent.mito,n=1000)
  plot2 <- ggplot(data.frame(cbind(sc10x$nCount_RNA,sc10x$percent.mito)))+geom_point(aes(x=X1,y=X2,color=density2),size=0.1)+scale_color_viridis(option="inferno")+geom_hline(yintercept=hm,size=1,color="red")+labs(x="nUMI",y="Percent Mitochondrial",color="Count")+cowplot::theme_cowplot()
  grid.arrange(plot1,plot2,nrow=1)
  dev.off()
  counts.cell.raw <- ncol(GetAssayData(object=sc10x,slot="counts"))
  counts.gene.raw <- nrow(GetAssayData(object=sc10x,slot="counts"))
  
  #Filter/normalize data
  sc10x.sub <- subset(x=sc10x,subset=nFeature_RNA > lg)
  sc10x.sub <- subset(x=sc10x.sub,subset=nFeature_RNA < hg)
  sc10x.sub <- subset(x=sc10x.sub,subset=percent.mito < hm)
  sc10x.sub <- NormalizeData(object=sc10x.sub,normalization.method="LogNormalize",scale.factor=10000,verbose=FALSE)
  
  #Plot filtered stats
  for(x in c("ALL","samples")){
    Idents(object=sc10x.sub) <- x
    postscript(paste0(folder,"Violin_qc.filtered.",x,".eps"))
    plot1 <- VlnPlot(object=sc10x.sub,features="nCount_RNA",pt.size=0.1)
    plot2 <- VlnPlot(object=sc10x.sub,features="nFeature_RNA",pt.size=0.1)
    plot3 <- VlnPlot(object=sc10x.sub,features="percent.mito",pt.size=0.1)
    grid.arrange(plot1,plot2,plot3,nrow=1)
    dev.off()
  }
  postscript(paste0(folder,"Plot_qc.filtered.eps"))
  density1 <- density(sc10x.sub$nCount_RNA,sc10x.sub$nFeature_RNA,n=1000)
  plot1 <- ggplot(data.frame(cbind(sc10x.sub$nCount_RNA,sc10x.sub$nFeature_RNA)))+geom_point(aes(x=X1,y=X2,color=density1),size=0.1)+scale_color_viridis(option="inferno")+labs(x="nUMI",y="nGenes",color="Count")+cowplot::theme_cowplot()
  density2 <- density(sc10x.sub$nCount_RNA,sc10x.sub$percent.mito,n=1000)
  plot2 <- ggplot(data.frame(cbind(sc10x.sub$nCount_RNA,sc10x.sub$percent.mito)))+geom_point(aes(x=X1,y=X2,color=density2),size=0.1)+scale_color_viridis(option="inferno")+labs(x="nUMI",y="Percent Mitochondrial",color="Count")+cowplot::theme_cowplot()
  grid.arrange(plot1,plot2,nrow=1)
  dev.off()
  counts.cell.filtered <- ncol(GetAssayData(object=sc10x.sub,slot="counts"))
  counts.gene.filtered <- nrow(GetAssayData(object=sc10x.sub,slot="counts"))
  
  results <- list(
    sc10x=sc10x.sub,
    counts.cell.raw=counts.cell.raw,
    counts.gene.raw=counts.gene.raw,
    counts.cell.filtered=counts.cell.filtered,
    counts.gene.filtered=counts.gene.filtered
  )
  return(results)
}


scCellCycle <- function(sc10x,sub=FALSE){
  #Runs Seurat based PCA analysis for cell cycle ID
  
  #Inputs:
  #sc10x = Seruat object
  #sub = Subfolder to save output files
  
  #Outputs:
  #results[1] = Seurat object
  #results[2] = s genes
  #results[3] = g2m genes
  
  #Make sub-folders if necessary
  if (sub==FALSE){
    folder <- "./analysis/qc/cellcycle/"
  } else {
    folder <- paste0("./analysis/qc/cellcycle/",sub,"/")
    if (!dir.exists(folder)){
      dir.create(folder)
    }}
  
  #score cell cycle
  genes.cc <- readLines(con="./genesets/regev_lab_cell_cycle_genes.txt")
  genes.s <- genes.cc[1:43]
  genes.g2m <- genes.cc[44:97]
  sc10x <- NormalizeData(object=sc10x,verbose=FALSE)
  sc10x <- ScaleData(object=sc10x,do.par=TRUE,num.cores=45,verbose=FALSE)
  sc10x <- CellCycleScoring(object=sc10x,s.features=genes.s,g2m.features=genes.g2m,set.ident=TRUE)
  
  #plot cell cycle specific genes
  genes=c("PCNA","TOP2A","MCM6","MKI67")
  postscript(paste0(folder,"Violin_cc.Raw.eps"))
  plot <- VlnPlot(object=sc10x,features=genes,ncol=2,pt.size=1)
  plot(plot)
  dev.off()
  
  # sc10x <- RunPCA(object=sc10x,features=c(genes.s,genes.g2m),npcs=2,verbose=FALSE)
  # postscript(paste0(folder,"PCA_cc.Raw.eps"))
  # plot <- DimPlot(object=sc10x,reduction="pca")
  # plot(plot)
  # dev.off()
  # gc()
  # sc10x <- ScaleData(object=sc10x,vars.to.regress=c("S.Score","G2M.Score"),do.par=TRUE,num.cores=45,verbose=TRUE)
  # gc()
  # sc10x <- RunPCA(object=sc10x,features=c(genes.s,genes.g2m),npcs=2,verbose=FALSE)
  # postscript(paste0(folder,"PCA_cc.Norm.eps"))
  # plot <- DimPlot(object=sc10x,reduction="pca")
  # plot(plot)
  # dev.off()
  
  results <- list(
    sc10x=sc10x,
    genes.s=genes.s,
    genes.g2m=genes.g2m
  )
  return(results)
}


scPC <- function(sc10x,pc=50,hpc=0.9,file="pre.stress",print="tsne",cca=FALSE){
  #Scale Seurat object & calculate # of PCs to use
  
  #Inputs:
  #sc10x = Seruat object
  #pc = number of PCs to cacluate
  #hpc = max variance cutoff for PCs to use"
  #file = file for output
  
  #Outputs:
  #result[1] = Seurat object
  #result[2] = # of PCs to use
  
  #Find variable genes
  if (cca==FALSE){sc10x <- FindVariableFeatures(object=sc10x,verbose=FALSE)}
  
  #Run PCA
  Idents(object=sc10x) <- "ALL"
  sc10x <- RunPCA(object=sc10x,npcs=pc,verbose=FALSE)
  
  #Calculate PCs to use
  pc.use <- sc10x[["pca"]]@stdev^2
  pc.use <- pc.use/sum(pc.use)
  pc.use <- cumsum(pc.use)
  pc.use <- min(which(pc.use>=hpc))
  
  postscript(paste0("./analysis/qc/Plot_PCElbow_",file,".eps"))
  plot <- ElbowPlot(object=sc10x,ndims=pc)
  plot <- plot+geom_vline(xintercept=pc.use,size=1,color="red")
  plot(plot)
  dev.off()
  
  results <- list(
    sc10x=sc10x,
    pc.use=pc.use
  )
  return(results)
}


scCCA <-  function(sc10x.l){
  for (i in 1:length(sc10x.l)){
    sc10x.l[[i]] <- NormalizeData(object=sc10x.l[[i]],verbose=FALSE)
    sc10x.l[[i]] <- FindVariableFeatures(object=sc10x.l[[i]],verbose=FALSE)
  }
  
  sc10x <- FindIntegrationAnchors(object.list=sc10x.l,dims=1:30)
  sc10x <- IntegrateData(anchorset=sc10x,dims=1:30)
  
  DefaultAssay(object=sc10x) <- "integrated"
  
  return(sc10x)
}


scCluster <- function(sc10x,res=0.1,red="pca",dim,print="tsne",folder=FALSE){
  #Cluster Seurat object and produce visualizations
  
  #Inputs:
  #sc10x = Seruat object
  #res = resolution to calculate clustering
  #red = rediction type to use for clustering
  #dim = number of dimentions to use for display
  #print = dimentionality reduction to use for display
  #folder = folder for output
  
  #Outputs:
  #result = Seurat object
  
  #Create subfolder if necessary
  if (folder==FALSE){
    sub <- ""
  } else {
    if (!dir.exists(paste0("./analysis/vis/",folder))){
      dir.create(paste0("./analysis/vis/",folder))
    }
    sub <- paste0(folder,"/")
    
  }
  
  #Calculste Vis
  sc10x <- RunTSNE(sc10x,dims=1:dim,reduction="pca")
  sc10x <- RunUMAP(sc10x,dims=1:dim,reduction="pca")
  
  sc10x <- FindNeighbors(sc10x,reduction=red,verbose=FALSE)
  
  for (i in res){
    sc10x <- FindClusters(sc10x,resolution=i,verbose=FALSE)
    
    plot1 <- DimPlot(sc10x,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
    plot2 <- DimPlot(sc10x,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
    plot3 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
    
    if (print=="tsne"){
      postscript(paste0("./analysis/vis/",sub,"tSNE_",i,".eps"))
      print(print2)
      dev.off()
    } else if (print=="umap"){
      postscript(paste0("./analysis/vis/",sub,"UMAP_",i,".eps"))
      print(print3)
      dev.off()
    } else if (print=="2"){
      plot2 <- plot2+theme(legend.position="none")
      plot3 <- plot3+theme(legend.position="none")
      postscript(paste0("./analysis/vis/",sub,"2Vis_",i,".eps"))
      grid.arrange(plot2,plot3,ncol=1)
      dev.off()
    } else if (print=="3"){
      plot1 <- plot1+theme(legend.position="none")
      plot2 <- plot2+theme(legend.position="none")
      plot3 <- plot3+theme(legend.position="none")
      postscript(paste0("./analysis/vis/",sub,"3Vis_",i,".eps"))
      grid.arrange(plot1,plot2,plot3,ncol=1)
      dev.off()
    }}
  
  plot1 <- DimPlot(sc10x,reduction="pca",group.by="samples")
  plot2 <- DimPlot(sc10x,reduction="tsne",group.by="samples")
  plot3 <- DimPlot(sc10x,reduction="umap",group.by="samples")
  legend <- cowplot::get_legend(plot1)

  if (print=="tsne"){
    postscript(paste0("./analysis/vis/",sub,"tSNE_samples.eps"))
    grid.arrange(plot2,legend,ncol=1)
    dev.off()
  } else if (print=="umap"){
    postscript(paste0("./analysis/vis/",sub,"UMAP_samples.eps"))
    grid.arrange(plot3,legend,ncol=1)
    dev.off()
  } else if (print=="2"){
    plot2 <- plot2+theme(legend.position="none")
    plot3 <- plot3+theme(legend.position="none")
    postscript(paste0("./analysis/vis/",sub,"2Vis_samples.eps"))
    grid.arrange(plot2,plot3,legend,ncol=1)
    dev.off()
  } else if (print=="3"){
    plot1 <- plot1+theme(legend.position="none")
    plot2 <- plot2+theme(legend.position="none")
    plot3 <- plot3+theme(legend.position="none")
    postscript(paste0("./analysis/vis/",sub,"3Vis_samples.eps"))
    grid.arrange(plot1,plot2,plot3,legend,ncol=1)
    dev.off()
  }
  
  return(sc10x)
}


scScore <- function(sc10x,score,geneset,cut.pt=0.9,anchor=FALSE){
  #Runs custom PCA analysis for stress ID
  
  #Inputs:
  #sc10x = Seruat object
  #score = name of geneset for scoring
  #geneset = geneset to use for ID
  #cut.pt = % of cells to keep
  
  #Outputs:
  #results[1] = Seurat object (original + score)
  #results[2] = Seurat object (negatively filtered)
  #results[3] = Seurat object (positively filtered)
  
  #Make subdirectory
  if (!dir.exists(paste0("./analysis/score_id/",score))){
    dir.create(paste0("./analysis/score_id/",score))
  }
  if (!dir.exists(paste0("./analysis/vis/",score))){
    dir.create(paste0("./analysis/vis/",score))
  }
  
  #Score geneset
  sc10x <- AddModuleScore(object=sc10x,features=genes.stress,name=score,assay="RNA")
  Idents(object=sc10x) <- paste0(score,"1")
  
  #CDF
  cdf <- ecdf(as.numeric(levels(sc10x)))
  cut.x <- quantile(cdf,probs=cut.pt)
  cut.x <- unname(cut.x)
  postscript(paste0("./analysis/score_id/",score,"/CDF_",score,".eps"))
  plot(cdf,main=paste0("Cumulative Distribution of ",score," Score"),xlab=paste0(score," Score"),ylab="CDF")
  abline(v=cut.x,col="red")
  dev.off()  
  
  #KDE
  postscript(paste0("./analysis/score_id/",score,"/Histo_",score,".eps"))
  plot(ggplot(data.frame(Score=as.numeric(levels(sc10x))),aes(x=Score))+geom_histogram(bins=100,aes(y=..density..))+geom_density()+geom_vline(xintercept=cut.x,size=1,color="red")+ggtitle(paste0(score," Score"))+cowplot::theme_cowplot())
  dev.off()
  
  Idents(object=sc10x) <- "ALL"
  predicate <- paste0(score,"1 >= ",cut.x)
  Idents(object=sc10x, cells = WhichCells(object=sc10x,expression= predicate)) <- score
  sc10x[[score]] <- Idents(object=sc10x)
  sc10x.negative <- subset(x=sc10x,idents="ALL")
  sc10x.positive <- subset(x=sc10x,idents=score)
  
  #Generate plots
  postscript(paste0("./analysis/vis/",score,"/3Vis_",score,".eps"))
  plot1 <- DimPlot(sc10x,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("ALL")
  plot2 <- DimPlot(sc10x.negative,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot3 <- DimPlot(sc10x.positive,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot4 <- DimPlot(sc10x,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Negative")
  plot5 <- DimPlot(sc10x.negative,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot6 <- DimPlot(sc10x.positive,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot7 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")+ggtitle("Positive")
  plot8 <- DimPlot(sc10x.negative,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot9 <- DimPlot(sc10x.positive,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
  grid.arrange(plot1,plot2,plot3,plot4,plot5,plot6,plot7,plot8,plot9,ncol=3)
  dev.off()
  
  #Generate violin plot of gene exvpression
  if (anchor!=FALSE){
    postscript(paste0("./analysis/score_id/",score,"/Violin_",score,".eps"))
    plot <- VlnPlot(object=sc10x,features=anchor,pt.size=0.1)
    plot(plot)
    dev.off()
  }
  
  results <- list(
    sc10x <- sc10x,
    sc10x.negative <- sc10x.negative,
    sc10x.positive <- sc10x.positive
  )
  return(results)
}


scQuSAGE <- function(sc10x,gs,save=FALSE,type,id,ds=0,nm="pops",print="tsne"){
  #Runs QuSAGE
  
  #Inputs:
  #sc10x = Seruat object
  #gs = geneset to use for correlation
  #save = save ID
  #type = type of qusage to run (id: create ID's based on run, sm: cor only using small genesets, lg: cor only with large genesets)
  #id = ident to use
  #nm = name of test
  #print = dimentionality reduction to use for display
  
  #Outputs:
  #results[1] = Seurat object
  #results[2] = correlation table
  #results[3] = correlation results
  
  if (!dir.exists(paste0("./analysis/cor/",nm))){
    dir.create(paste0("./analysis/cor/",nm))
  }
  if (!dir.exists(paste0("./analysis/cor/",nm,"/geneset"))){
    dir.create(paste0("./analysis/cor/",nm,"/geneset"))
  }
  if (!dir.exists(paste0("./analysis/cor/",nm,"/cluster"))){
    dir.create(paste0("./analysis/cor/",nm,"/cluster"))
  }
  if (!dir.exists(paste0("./analysis/vis/",nm))){
    dir.create(paste0("./analysis/vis/",nm))
  }
  
  Idents(object=sc10x) <- id
  number.clusters <- length(unique(levels(x=sc10x)))
  
  labels <- paste0("Cluster_",as.vector(Idents(object=sc10x)))
  
  cell.sample <- NULL
  for (i in unique(labels)){
    cell <- WhichCells(sc10x,ident=sub("Cluster_","",i))
    if (length(cell)>ds & ds!=0){
      rnd <- sample(1:length(cell),ds)
      cell <- cell[rnd]
    }
    cell.sample <- c(cell.sample,cell)
  }
  data <- as.data.frame(as.matrix(GetAssayData(sc10x[,colnames(sc10x) %in% cell.sample])))
  labels <- labels[colnames(sc10x) %in% cell.sample]
  groups <- sort(unique(labels))
  
  col <- hcl(h=(seq(15,375-375/length(groups),length=length(groups))),c=100,l=65)
  
  #Make labels for QuSAGE
  clust <- list()
  clust.comp <- list()
  for (i in groups){
    t <- labels
    t[t != i] <- "REST"
    clust[i] <- list(i=t)
    rm(t)
    clust.comp[i] <- paste0(i,"-REST")
  }
  
  #Run QuSAGE
  for (i in groups){
    assign(paste0("results.",i),qusage(data,unlist(clust[i]),unlist(clust.comp[i]),gs))
  }
  
  #Generate ID table
  results.cor <- NULL
  results.cor <- qsTable(get(paste0("results.",groups[1])),number=length(gs))
  results.cor$Cluster <- groups[1]
  for (i in groups[-1]){
    qs <- qsTable(get(paste0("results.",i)),number=length(gs))
    qs$Cluster <- i
    results.cor <- rbind(results.cor,qs)
  }
  results.cor <- results.cor[,-3]
  rownames(results.cor) <- NULL
  
  results.clust.id <- NULL
  if (max(results.cor[results.cor[,4]==groups[1] & results.cor[,3]<=0.05,][,2],na.rm=TRUE)>=0){
    results.clust.id <- results.cor[results.cor[,4]==groups[1] & results.cor[,3]<=0.05,][which.max(results.cor[results.cor[,4]==groups[1] & results.cor[,3]<=0.05,][,2]),]
  } else {
    results.clust.id$pathway.name <- "Unknown"
    results.clust.id$log.fold.change <- 0
    results.clust.id$FDR <- 0
    results.clust.id$Cluster <- groups[1]
    results.clust.id <- as.data.frame(results.clust.id)
  }
  for (i in groups[-1]){
    if (max(results.cor[results.cor[,4]==i & results.cor[,3]<=0.05,][,2],na.rm=TRUE)>=0){
      results.clust.id <- rbind(results.clust.id,results.cor[results.cor[,4]==i & results.cor[,3]<=0.05,][which.max(results.cor[results.cor[,4]==i & results.cor[,3]<=0.05,][,2]),])
    } else {
      results.clust.id <- rbind(results.clust.id,data.frame(pathway.name="Unknown",log.fold.change=0,FDR=0,Cluster=i))
    }}
  rownames(results.clust.id) <- NULL
  
  #Determine axes for correlation plots
  max.x.rg <- 0
  min.x.rg <- 0
  max.y.rg <- 0
  for (i in groups){
    qs <- get(paste0("results.",i))
    if (max(qs$path.mean)>max.x.rg){
      max.x.rg <- max(qs$path.mean)
    }
    if (min(qs$path.mean)<min.x.rg){
      min.x.rg <- min(qs$path.mean)
    }
    if (max(qs$path.PDF)>max.y.rg){
      max.y.rg <- max(qs$path.PDF)
    }}
  
  #Plot correlation plots by geneset
  for (i in 1:length(gs)){
    postscript(paste0("./analysis/cor/",nm,"/geneset/QuSAGE_",nm,".",names(gs)[i],".eps"))
    for (j in groups){
      qs <- get(paste0("results.",j))
      if (j==groups[1]){
        plotDensityCurves(qs,path.index=i,col=col[match(j,groups)],main=names(gs)[i],xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
      } else {
        plotDensityCurves(qs,path.index=i,add=TRUE,col=col[match(j,groups)],lwd=5)
      }}
    legend("topright",col=col,legend=groups,lty=1,lwd=5,cex=2,ncol=1,bty="n",pt.cex=2)
    box(lwd=5)
    dev.off()
  }
  
  if (type=="sm"){
    #Plot correlation plots by cluster
    for (i in groups){
      qs <- get(paste0("results.",i))
      postscript(paste0("./analysis/cor/",nm,"/cluster/QuSAGE_",nm,"_",i,".eps"))
      for (j in 1:length(gs)){
        if (j==1){
          plotDensityCurves(qs,path.index=j,col=viridis(length(gs))[j],main=i,xlim=c(min.x.rg-0.05,max.x.rg+0.05),ylim=c(0,50*ceiling(max.y.rg/50)),xlab="Gene Set Activation",lwd=5,cex.main=2.5,cex.axis=1.5,cex.lab=2)
        } else {
          plotDensityCurves(qs,path.index=j,add=TRUE,col=viridis(length(gs))[j],lwd=5)
        }}
      legend("topright",col=viridis(length(gs)),legend=names(gs),lty=1,lwd=5,cex=1,ncol=2,bty="n",pt.cex=2)
      box(lwd=5)
      dev.off()
    }} else {
      for (i in groups){
        qs <- get(paste0("results.",i))
        postscript(paste0("./analysis/cor/",nm,"/cluster/QuSAGE_",nm,"_",i,".eps"))
        plotCIs(qs,path.index=1:numPathways(qs),cex.lab=1.5)
        dev.off()
      }}
  
  if (save==TRUE){
    merge.cluster <- NULL
    for (i in groups){
      if (max(qsTable(get(paste0("results.",i)),number=length(gs))[qsTable(get(paste0("results.",i)),number=length(gs))[,4]<=0.05,][,2],na.rm=TRUE)>=0){
        sc10x<-eval(parse(text=paste0("RenameIdents(object=sc10x,'",sub("Cluster_","",i),"' = '",qsTable(get(paste0("results.",i)),number=length(gs))[qsTable(get(paste0("results.",i)),number=length(gs))[2]==max(qsTable(get(paste0("results.",i)),number=length(gs))[qsTable(get(paste0("results.",i)),number=length(gs))[,4]<=0.05,][,2])][1],"')")))
      } else {
        sc10x<-eval(parse(text=paste0("RenameIdents(object=sc10x,'",sub("Cluster_","",i),"' = 'Unknown')")))
      }}
    sc10x[[nm]] <- Idents(object=sc10x)
  }
  
  plot1 <- DimPlot(sc10x,reduction="pca",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot2 <- DimPlot(sc10x,reduction="tsne",label=TRUE,repel=TRUE)+theme(legend.position="none")
  plot3 <- DimPlot(sc10x,reduction="umap",label=TRUE,repel=TRUE)+theme(legend.position="none")
  if (print=="tsne"){
    postscript(paste0("./analysis/vis/",nm,"/tSNE_",id,"_",nm,".eps"))
    print(print2)
    dev.off()
  } else if (print=="umap"){
    postscript(paste0("./analysis/vis/",nm,"/UMAP_",id,"_",nm,".eps"))
    print(print3)
    dev.off()
  } else if (print=="2"){
    postscript(paste0("./analysis/vis/",nm,"/2Vis_",id,"_",nm,".eps"))
    grid.arrange(plot2,plot3,ncol=1)
    dev.off()
  } else if (print=="3"){
    postscript(paste0("./analysis/vis/",nm,"/3Vis_",id,"_",nm,".eps"))
    grid.arrange(plot1,plot2,plot3,ncol=1)
    dev.off()
  }
  
  results <- list(
    sc10x=sc10x,
    results.cor=results.cor,
    results.clust.id=results.clust.id
  )
  names(results)=c("sc10x",paste0("results.cor.",nm),paste0("results.clust.",nm,".id"))
  return(results)
}