require(Seurat)
require(ggplot2)
require(harmony)
require(dplyr)
require(reshape2)
require(scRNAtoolVis)
require(ggpubr)
require(clusterProfiler)
options(dplyr.print_max = 1e9)
Sys.setenv(LANGUAGE = "en")
options(stringsAsFactors = FALSE)
set.seed(123456)

ggsave_fun <- function(filename, plot = last_plot(), plot_Device = c(".pdf"), width = 7, height = 7, dpi = 600, ...) {
  filenames <- paste0(filename, plot_Device)
  for (filename_tmp in filenames) {
    ggsave(filename = filename_tmp, plot = plot, width = width, height = height, dpi = dpi, ...)
  }
}

paramSweep_V4 <- function(seu, PCs = 1:10, sct = FALSE, num.cores = 1) {
  require(Seurat)
  require(fields)
  require(parallel)
  
  pK <- c(0.0005, 0.001, 0.005, seq(0.01, 0.3, by = 0.01))
  pN <- seq(0.05, 0.3, by = 0.05)
  
  
  min.cells <- round(nrow(seu@meta.data) / (1 - 0.05) - nrow(seu@meta.data))
  pK.test <- round(pK * min.cells)
  pK <- pK[which(pK.test >= 1)]
  
  
  orig.commands <- seu@commands
  
  
  if (nrow(seu@meta.data) > 10000) {
    real.cells <- rownames(seu@meta.data)[sample(1:nrow(seu@meta.data), 10000, replace = FALSE)]
    data <- seu@assays$RNA@counts[, real.cells]
    n.real.cells <- ncol(data)
  }
  
  if (nrow(seu@meta.data) <= 10000) {
    real.cells <- rownames(seu@meta.data)
    data <- seu@assays$RNA@counts
    n.real.cells <- ncol(data)
  }
  
  
  
  if (num.cores > 1) {
    require(parallel)
    cl <- makeCluster(num.cores)
    output2 <- mclapply(as.list(1:length(pN)),
                        FUN = parallel_paramSweep,
                        n.real.cells,
                        real.cells,
                        pK,
                        pN,
                        data,
                        orig.commands,
                        PCs,
                        sct, mc.cores = num.cores
    )
    stopCluster(cl)
  } else {
    output2 <- lapply(as.list(1:length(pN)),
                      FUN = parallel_paramSweep,
                      n.real.cells,
                      real.cells,
                      pK,
                      pN,
                      data,
                      orig.commands,
                      PCs,
                      sct
    )
  }
  
  
  sweep.res.list <- list()
  list.ind <- 0
  for (i in 1:length(output2)) {
    for (j in 1:length(output2[[i]])) {
      list.ind <- list.ind + 1
      sweep.res.list[[list.ind]] <- output2[[i]][[j]]
    }
  }
  
  
  name.vec <- NULL
  for (j in 1:length(pN)) {
    name.vec <- c(name.vec, paste("pN", pN[j], "pK", pK, sep = "_"))
  }
  names(sweep.res.list) <- name.vec
  return(sweep.res.list)
}

doubletFinder_V4 <- function(seu, PCs, pN = 0.25, pK, nExp, reuse.pANN = FALSE, sct = FALSE, annotations = NULL) {
  require(Seurat)
  require(fields)
  require(KernSmooth)
  
  
  if (reuse.pANN != FALSE) {
    pANN.old <- seu@meta.data[, reuse.pANN]
    classifications <- rep("Singlet", length(pANN.old))
    classifications[order(pANN.old, decreasing = TRUE)[1:nExp]] <- "Doublet"
    seu@meta.data[, paste("DF.classifications", pN, pK, nExp, sep = "_")] <- classifications
    return(seu)
  }
  
  if (reuse.pANN == FALSE) {
    
    real.cells <- rownames(seu@meta.data)
    data <- seu@assays$RNA@counts[, real.cells]
    n_real.cells <- length(real.cells)
    n_doublets <- round(n_real.cells / (1 - pN) - n_real.cells)
    print(paste("Creating", n_doublets, "artificial doublets...", sep = " "))
    real.cells1 <- sample(real.cells, n_doublets, replace = TRUE)
    real.cells2 <- sample(real.cells, n_doublets, replace = TRUE)
    doublets <- (data[, real.cells1] + data[, real.cells2]) / 2
    colnames(doublets) <- paste("X", 1:n_doublets, sep = "")
    data_wdoublets <- cbind(data, doublets)
    
    if (!is.null(annotations)) {
      stopifnot(typeof(annotations) == "character")
      stopifnot(length(annotations) == length(Cells(seu)))
      stopifnot(!any(is.na(annotations)))
      annotations <- factor(annotations)
      names(annotations) <- Cells(seu)
      doublet_types1 <- annotations[real.cells1]
      doublet_types2 <- annotations[real.cells2]
    }
    
    orig.commands <- seu@commands
    
    
    if (sct == FALSE) {
      print("Creating Seurat object...")
      seu_wdoublets <- CreateSeuratObject(counts = data_wdoublets)
      
      print("Normalizing Seurat object...")
      seu_wdoublets <- NormalizeData(seu_wdoublets,
                                     normalization.method = orig.commands$NormalizeData.RNA@params$normalization.method,
                                     scale.factor = orig.commands$NormalizeData.RNA@params$scale.factor,
                                     margin = orig.commands$NormalizeData.RNA@params$margin
      )
      
      print("Finding variable genes...")
      seu_wdoublets <- FindVariableFeatures(seu_wdoublets,
                                            selection.method = orig.commands$FindVariableFeatures.RNA$selection.method,
                                            loess.span = orig.commands$FindVariableFeatures.RNA$loess.span,
                                            clip.max = orig.commands$FindVariableFeatures.RNA$clip.max,
                                            mean.function = orig.commands$FindVariableFeatures.RNA$mean.function,
                                            dispersion.function = orig.commands$FindVariableFeatures.RNA$dispersion.function,
                                            num.bin = orig.commands$FindVariableFeatures.RNA$num.bin,
                                            binning.method = orig.commands$FindVariableFeatures.RNA$binning.method,
                                            nfeatures = orig.commands$FindVariableFeatures.RNA$nfeatures,
                                            mean.cutoff = orig.commands$FindVariableFeatures.RNA$mean.cutoff,
                                            dispersion.cutoff = orig.commands$FindVariableFeatures.RNA$dispersion.cutoff
      )
      
      print("Scaling data...")
      seu_wdoublets <- ScaleData(seu_wdoublets,
                                 features = orig.commands$ScaleData.RNA$features,
                                 model.use = orig.commands$ScaleData.RNA$model.use,
                                 do.scale = orig.commands$ScaleData.RNA$do.scale,
                                 do.center = orig.commands$ScaleData.RNA$do.center,
                                 scale.max = orig.commands$ScaleData.RNA$scale.max,
                                 block.size = orig.commands$ScaleData.RNA$block.size,
                                 min.cells.to.block = orig.commands$ScaleData.RNA$min.cells.to.block
      )
      
      print("Running PCA...")
      seu_wdoublets <- RunPCA(seu_wdoublets,
                              features = orig.commands$ScaleData.RNA$features,
                              npcs = length(PCs),
                              rev.pca = orig.commands$RunPCA.RNA$rev.pca,
                              weight.by.var = orig.commands$RunPCA.RNA$weight.by.var,
                              verbose = FALSE
      )
      pca.coord <- seu_wdoublets@reductions$pca@cell.embeddings[, PCs]
      cell.names <- rownames(seu_wdoublets@meta.data)
      nCells <- length(cell.names)
      rm(seu_wdoublets)
      gc() 
    }
    
    if (sct == TRUE) {
      require(sctransform)
      print("Creating Seurat object...")
      seu_wdoublets <- CreateSeuratObject(counts = data_wdoublets)
      
      print("Running SCTransform...")
      seu_wdoublets <- SCTransform(seu_wdoublets)
      
      print("Running PCA...")
      seu_wdoublets <- RunPCA(seu_wdoublets, npcs = length(PCs))
      pca.coord <- seu_wdoublets@reductions$pca@cell.embeddings[, PCs]
      cell.names <- rownames(seu_wdoublets@meta.data)
      nCells <- length(cell.names)
      rm(seu_wdoublets)
      gc()
    }
    
    
    print("Calculating PC distance matrix...")
    dist.mat <- fields::rdist(pca.coord)
    
    
    print("Computing pANN...")
    pANN <- as.data.frame(matrix(0L, nrow = n_real.cells, ncol = 1))
    if (!is.null(annotations)) {
      neighbor_types <- as.data.frame(matrix(0L, nrow = n_real.cells, ncol = length(levels(doublet_types1))))
    }
    rownames(pANN) <- real.cells
    colnames(pANN) <- "pANN"
    k <- round(nCells * pK)
    for (i in 1:n_real.cells) {
      neighbors <- order(dist.mat[, i])
      neighbors <- neighbors[2:(k + 1)]
      pANN$pANN[i] <- length(which(neighbors > n_real.cells)) / k
      if (!is.null(annotations)) {
        for (ct in unique(annotations)) {
          neighbors_that_are_doublets <- neighbors[neighbors > n_real.cells]
          if (length(neighbors_that_are_doublets) > 0) {
            neighbor_types[i, ] <-
              table(doublet_types1[neighbors_that_are_doublets - n_real.cells]) +
              table(doublet_types2[neighbors_that_are_doublets - n_real.cells])
            neighbor_types[i, ] <- neighbor_types[i, ] / sum(neighbor_types[i, ])
          } else {
            neighbor_types[i, ] <- NA
          }
        }
      }
    }
    print("Classifying doublets..")
    classifications <- rep("Singlet", n_real.cells)
    classifications[order(pANN$pANN[1:n_real.cells], decreasing = TRUE)[1:nExp]] <- "Doublet"
    seu@meta.data[, paste("pANN", pN, pK, nExp, sep = "_")] <- pANN[rownames(seu@meta.data), 1]
    seu@meta.data[, paste("DF.classifications", pN, pK, nExp, sep = "_")] <- classifications
    if (!is.null(annotations)) {
      colnames(neighbor_types) <- levels(doublet_types1)
      for (ct in levels(doublet_types1)) {
        seu@meta.data[, paste("DF.doublet.contributors", pN, pK, nExp, ct, sep = "_")] <- neighbor_types[, ct]
      }
    }
    return(seu)
  }
}

human2mouse <- function(x) {
  x1 <- substr(x, 1, 1)
  x2 <- substr(x, 2, nchar(x))
  return(paste0(x1, tolower(x2)))
}

getScatterplot <- function(object, gene1, gene2, cor.method = "pearson", jitter.num = 0.15, pos = TRUE) {
  if (!gene1 %in% rownames(object)) {
    print("gene1 was not found")
    if (!gene2 %in% rownames(object)) {
      print("gene2 was not found")
    }
  } else {
    exp.mat <- GetAssayData(object = object, assay = "RNA") %>%
      .[c(gene1, gene2), ] %>%
      as.matrix() %>%
      t() %>%
      as.data.frame()
    if (pos) {
      if (nrow(exp.mat[which(exp.mat[, 1] > 0 & exp.mat[, 2] > 0), ]) > (nrow(exp.mat) * 0.01)) {
        exp.mat <- exp.mat[which(exp.mat[, 1] > 0 & exp.mat[, 2] > 0), ]
      } else {
        exp.mat <- exp.mat[which(exp.mat[, 1] > 0 | exp.mat[, 2] > 0), ]
      }
    }
    colnames(exp.mat) <- c("Var1", "Var2")
    plots <- ggplot(data = exp.mat, mapping = aes_string(x = "Var1", y = "Var2")) +
      geom_smooth(method = "lm", se = T, color = "red", size = 1) +
      stat_cor(method = cor.method) +
      labs(x = gene1, y = gene2) +
      geom_jitter(width = jitter.num, height = jitter.num, color = "black", size = 1, alpha = 1) +
      theme_bw() +
      theme(
        panel.grid = element_blank(),
        legend.text = element_text(colour = "black", size = 10),
        axis.text = element_text(colour = "black", size = 10),
        axis.line = element_line(colour = "black"),
        panel.border = element_rect(size = 1, linetype = "solid", colour = "black"),
        panel.background = element_rect(fill = "white")
      )
    return(plots)
  }
}

capitalize_first <- function(string) {
  words <- strsplit(string, " ")[[1]]
  capitalized_words <- toupper(substring(words, 1, 1))
  rest_of_words <- substring(words, 2)
  capitalized_string <- paste0(capitalized_words, rest_of_words, collapse = " ")
  return(capitalized_string)
}




























scRNAAutoAnno <- function(Path = ".", SeuratObject = NA, Multi = TRUE, ref, labels,min.features = 200, Idents = "RNA_snn_res.0.2", DoubletFinder = TRUE, species = NA, KeyCell = NA, KeyGene = NA, logFCfilter = 0.585, adjPvalFilter = 0.05,assay = "RNA", set.resolutions = seq(0.2, 1.2, by = 0.1), PC = 20, nfeatures = 2000, npcs = 50) {
  cat("1.质量控制\n")
  cat("2.数据标准化\n")
  cat("4.Cellchat\n")
  cat("7.关键基因表达丰度\n")
  cat("9.免疫代谢通路\n")
  cat("\033[0;34m请选择要分析的内容，按序号填写，英文逗号分隔")
  AnalysisIndex <- scan(sep = ",", quiet = TRUE)
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  scRNA <- SeuratObject
  
  if (!is.na(species)) {
    mouse <- ifelse(species == "mouse", T, F)
  } else {
    mouse <- F
  }
  
  s.genes <- cc.genes$s.genes
  g2m.genes <- cc.genes$g2m.genes
  
  if (mouse) {
    s.genes <- human2mouse(s.genes)
    g2m.genes <- human2mouse(g2m.genes)
  }
  
  n <- 0
  KeyCells <- KeyCell
  for (step in AnalysisIndex) {
    n <- n + 1
    
    
    if (step == 1) {
      message("*************************质量控制*************************")
      
      dir <- paste0(Path, "/.QCFilter")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      grep("^[M,m][T,t]-", rownames(scRNA), value = T)
      grep("^R[P,p][SL,sl]", rownames(scRNA), value = T)
      scRNA[["percent.mt"]] <- PercentageFeatureSet(object = scRNA, pattern = "^[M,m][T,t]-")
      scRNA[["percent.ribo"]] <- PercentageFeatureSet(scRNA, pattern = "^R[P,p][SL,sl]") 
      
      nomt <- ifelse(max(scRNA[["percent.mt"]]) == 0, T, F)
      if (nomt) {
        MT_genes <- c("ND1", "ND2", "COX1", "COX2", "ATP8", "ATP6", "COX3", "ND3", "ND4L", "ND4", "ND5", "ND6", "CYTB")
        if (mouse) {
          MT_genes <- human2mouse(MT_genes)
        }
        MT_genes <- intersect(MT_genes, rownames(scRNA))
        scRNA[["percent.mt"]] <- PercentageFeatureSet(object = scRNA, features = MT_genes)
      }
      
      mt_value <- stats::mad(scRNA$percent.mt)
      lower_mt <- median(scRNA$percent.mt) - 3 * mt_value
      upper_mt <- median(scRNA$percent.mt) + 3 * mt_value
      
      nF_value <- stats::mad(scRNA$nFeature_RNA)
      lower_nF <- median(scRNA$nFeature_RNA) - 3 * nF_value
      upper_nF <- median(scRNA$nFeature_RNA) + 3 * nF_value
      
      nC_value <- stats::mad(scRNA$nCount_RNA)
      lower_nC <- median(scRNA$nCount_RNA) - 3 * nC_value
      upper_nC <- median(scRNA$nCount_RNA) + 3 * nC_value
      
      
      
      
      
      scRNA <- subset(x = scRNA, subset = nFeature_RNA >= min.features & percent.mt <= upper_mt & nFeature_RNA <= upper_nF & nCount_RNA <= upper_nC)
      dim(scRNA)
      
      if (DoubletFinder) {
        require(DoubletFinder)
        scRNA_list <- SplitObject(scRNA, split.by = "orig.ident")
        newscRMA_list <- list()
        for (i in 1:length(scRNA_list)) {
          data <- scRNA_list[[i]]
          data <- NormalizeData(data)
          data <- FindVariableFeatures(data, selection.method = "vst", nfeatures = 2000)
          data <- ScaleData(data)
          data <- RunPCA(data)
          data <- RunUMAP(data, dims = 1:20)
          
          sweep.res.list <- paramSweep_V4(data, PCs = 1:20, sct = FALSE) 
          sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
          bcmvn <- find.pK(sweep.stats)
          p <- as.numeric(as.vector(bcmvn[bcmvn$MeanBC == max(bcmvn$MeanBC), ]$pK))
          
          
          homotypic.prop <- modelHomotypic(data@meta.data$seurat_clusters) 
          Doubletrate <- ncol(data) * 8 * 1e-6
          nExp_poi <- round(Doubletrate * ncol(data))
          nExp_poi.adj <- round(nExp_poi * (1 - homotypic.prop))
          
          
          data <- doubletFinder_V4(data, PCs = 1:20, pN = 0.25, pK = p, nExp = nExp_poi.adj, reuse.pANN = FALSE, sct = FALSE)
          colnames(data@meta.data)[ncol(data@meta.data)] <- "doublet_info"
          newscRMA_list[[i]] <- data
        }
        
        saveRDS(newscRMA_list, paste0(dir, "/doubletFinder_V4.rds"))
        
        plot_list <- list()
        for (i in 1:length(newscRMA_list)) {
          plot_list[[i]] <- DimPlot(newscRMA_list[[i]], group.by = "doublet_info") + ggtitle(newscRMA_list[[i]]$orig.ident[1]) + theme(plot.title = element_text(hjust = 0.5))
        }
        
        
        for (i in 1:length(plot_list)) {
          plot_list[[i]] <- plot_list[[i]] + NoAxes() + NoLegend()
        }
        
        p <- CombinePlots(plot_list, legend = "right")
        
        ggsave_fun(filename = paste0(dir, "/1.DoubletFinder"), plot = p, width = 12, height = 10)
        
        if (Multi) {
          scRNA <- merge(newscRMA_list[[1]], newscRMA_list[2:length(newscRMA_list)])
        } else {
          scRNA <- newscRMA_list[[1]]
        }
        
        scRNA <- subset(scRNA, subset = doublet_info == "Singlet")
        
        p <- VlnPlot(scRNA, features = c("nCount_RNA", "nFeature_RNA", "percent.mt", "percent.ribo"), group.by = "orig.ident", ncol = 4)
        
        ggsave_fun(filename = paste0(dir, "/2.VlnPlot"), plot = p, width = 20, height = 4)
        
        plot1 <-
          FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F) +
          NoLegend()
        plot2 <-
          FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F) +
          NoLegend()
        plot3 <-
          FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.ribo", raster = F) +
          NoLegend()
        
        plot1 + plot2 + plot3
        ggsave_fun(filename = paste0(dir, "/3.FeatureScatter"), plot = plot1 + plot2 + plot3, height = 7, width = 21)
      } else {
        p <- VlnPlot(scRNA, features = c("nCount_RNA", "nFeature_RNA", "percent.mt", "percent.ribo"), group.by = "orig.ident", ncol = 4)
        
        ggsave_fun(filename = paste0(dir, "/1.VlnPlot"), plot = p, width = 20, height = 4)
        
        plot1 <-
          FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F) +
          NoLegend()
        plot2 <-
          FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.mt", raster = F) +
          NoLegend()
        plot3 <-
          FeatureScatter(scRNA, feature1 = "percent.mt", feature2 = "percent.ribo", raster = F) +
          NoLegend()
        
        plot1 + plot2 + plot3
        ggsave_fun(filename = paste0(dir, "/2.FeatureScatter"), plot = plot1 + plot2 + plot3, height = 7, width = 21)
      }
    }
    
    
    if (step == 2) {
      message("*************************数据标准化*************************")
      
      dir <- paste0(Path, "/.Cluster")
      
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      DefaultAssay(scRNA) <- assay
      
      scRNA <- NormalizeData(object = scRNA)
      scRNA <- CellCycleScoring(object = scRNA, s.features = s.genes, g2m.features = g2m.genes, set.ident = F)
      scRNA <- FindVariableFeatures(object = scRNA, selection.method = "vst", nfeatures = nfeatures)
      top10 <- head(VariableFeatures(scRNA), 10)
      plot1 <- VariableFeaturePlot(scRNA)
      plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
      plot1 + plot2
      ggsave_fun(filename = paste0(dir, "/1.VariableFeaturePlot"), plot = plot1 + plot2, width = 14)
      scRNA <- ScaleData(object = scRNA, vars.to.regress = c("percent.mt", "percent.ribo", "S.Score", "G2M.Score"))
      scRNA <- RunPCA(scRNA, verbose = T, npcs = npcs, features = VariableFeatures(object = scRNA))
      p <- ElbowPlot(object = scRNA, ndims = npcs)
      ggsave_fun(filename = paste0(dir, "/2.ElbowPlot"), plot = p)
      
      if (Multi) {
        scRNA <- RunHarmony(object = scRNA, group.by.vars = "orig.ident", assay.use = assay, verbose = FALSE)
        scRNA <- RunUMAP(scRNA, reduction = "harmony", dims = 1:PC, verbose = T)
        
        scRNA <- FindNeighbors(scRNA, dims = 1:PC, reduction = "harmony", verbose = T) 
        scRNA <- FindClusters(scRNA, resolution = set.resolutions, verbose = T)
        p <- DimPlot(object = scRNA, reduction = "pca", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/3.PcaPlot"), plot = p)
        p <- DimPlot(object = scRNA, reduction = "harmony", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/4.HarmonyPlot"), plot = p)
      } else {
        scRNA <- RunUMAP(scRNA, reduction = "pca", dims = 1:PC, verbose = T)
        
        scRNA <- FindNeighbors(scRNA, dims = 1:PC, reduction = "pca", verbose = T)
        scRNA <- FindClusters(scRNA, resolution = set.resolutions, verbose = T)
        p <- DimPlot(object = scRNA, reduction = "pca", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/3.PcaPlot"), plot = p)
      }
      
      pdf(paste0(dir, "/7.data.merge.harmony.pdf"))
      
      
      p <- DimPlot(object = scRNA, reduction = "umap", label = TRUE, group.by = "orig.ident", raster = F) + NoLegend()
      print(p)
      merge.res <- sapply(set.resolutions, function(x) {
        p <- DimPlot(object = scRNA, reduction = "umap", label = TRUE, group.by = paste0(assay, "_snn_res.", x), raster = F) + NoLegend()
        print(p)
      })
      dev.off()
      saveRDS(scRNA, paste0(dir, "/data.merge.harmony.2000.rds"))
    }
    

    
    
    if (step == 4) {
      message("*************************细胞通讯*************************")
      require(CellChat)
      
      dir <- paste0(Path, "/.Cellchat")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      options(stringsAsFactors = FALSE)
      cellchat <- createCellChat(object = scRNA, group.by = "cellType_1")
      cellchat
      groupSize <- as.numeric(table(cellchat@idents)) 
      
      if (mouse) {
        CellChatDB <- CellChatDB.mouse
      } else {
        CellChatDB <- CellChatDB.human
      }
      
      CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")
      cellchat@DB <- CellChatDB.use
      
      cellchat <- subsetData(cellchat)
      
      cellchat <- identifyOverExpressedGenes(cellchat)
      cellchat <- identifyOverExpressedInteractions(cellchat)
      cellchat <- projectData(cellchat, PPI.human)
      
      cellchat <- computeCommunProb(cellchat,seed.use = 123456)
      cellchat <- filterCommunication(cellchat)
      df.net <- subsetCommunication(cellchat)
      write.table(df.net, file = paste0(dir, "/net_lr.txt"), quote = F, sep = "\t", row.names = F)
      
      cellchat <- computeCommunProbPathway(cellchat)
      df.netp <- subsetCommunication(cellchat)
      write.table(df.netp, file = paste0(dir, "/net_pathway.txt"), quote = F, sep = "\t", row.names = F)
      
      
      cellchat <- aggregateNet(cellchat)
      saveRDS(cellchat,paste0(dir,"/cellchat.rds"))
      groupSize <- as.numeric(table(cellchat@idents))
      
      
      
      
      
      
      
      
      
      
      
      
      
      pdf(file=paste0(dir, "/1.Net_number_strength.pdf"),width = 10,height = 5)
      par(mfrow = c(1,2), xpd=TRUE)
      netVisual_circle(cellchat@net$count, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Number of interactions")
      netVisual_circle(cellchat@net$weight, vertex.weight = groupSize, weight.scale = T, label.edge= F, title.name = "Interaction weights/strength")
      dev.off()
      
      df.net <- read.table(paste0(dir, "/net_lr.txt"), sep = "\t", check.names = F, header = T)
      data <- as.data.frame(table(c(
        df.net$source,
        df.net$target
      )))
      colnames(data) <- c("Cell_Type", "all_sum")
      data <- data[order(data$all_sum, decreasing = T), ]
      data$Cell_Type <- factor(data$Cell_Type, levels = data$Cell_Type)
      head(data)
      sample_color <- c("#FB040B", "#F6A717", "#BA06FA", "#172D7A")
      sample_color2 <- ifelse(data$all_sum > summary(data$all_sum)[2], ifelse(data$all_sum > summary(data$all_sum)[3], ifelse(data$all_sum > summary(data$all_sum)[5], sample_color[1], sample_color[2]), sample_color[3]), sample_color[4])
      sample_color1 <- ifelse(data$all_sum > summary(data$all_sum)[2], ifelse(data$all_sum > summary(data$all_sum)[3], ifelse(data$all_sum > summary(data$all_sum)[5], alpha(sample_color[1], 0.9), alpha(sample_color[2], 0.9)), alpha(sample_color[3], 0.9)), alpha(sample_color[4], 0.9))
      
      pB2 <- ggplot(data = data, aes(x = Cell_Type, y = all_sum, fill = Cell_Type, colour = Cell_Type)) +
        geom_bar(stat = "identity", width = 0.8) +
        scale_fill_manual(values = sample_color1) +
        scale_colour_manual(values = sample_color2) +
        theme_bw() +
        theme(panel.grid = element_blank()) +
        labs(x = "Cell Type", y = "Counts") +
        theme(axis.text.y = element_text(size = 12, colour = "black")) +
        theme(axis.text.x = element_text(size = 12, angle = 45, hjust = 1, vjust = 1, colour = "black"))
      
      ggsave_fun(paste0(dir, "/2.Interaction Count"), plot = pB2, width = 10, height = 8)
      
      
      
      
    }


    
    
    if (step == 7) {
      message("*************************关键基因表达丰度*************************")
      
      dir <- paste0(Path, "/.KeyGene expression abundance")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      KeyGene <- intersect(rownames(scRNA), KeyGene)
      leg <- length(unique(KeyGene))
      
      p <- FeatureCornerAxes(
        object = scRNA, reduction = "umap",
        groupFacet = NULL,
        relLength = 0.5,
        relDist = 0.2,
        features = unique(KeyGene), nLayout = 3,
        aspect.ratio = 1, pSize = 0.1
      )
      if (leg <= 3) {
        ggsave_fun(filename = paste0(dir, "/1.Keygene.FeaturePlot"), plot = p, width = leg * 5, height = 5)
      }
      if (leg > 3) {
        ggsave_fun(filename = paste0(dir, "/1.Keygene.FeaturePlot"), plot = p, width = 15, height = ceiling(leg / 3) * 5)
      }
      
      p <- DotPlot(scRNA, group.by = "cellType_1", features = unique(KeyGene), cols = "RdYlBu", ) + scale_size_continuous(range = c(0, 10)) + theme(
        panel.border = element_rect(colour = "black"),
        axis.text.x = element_text(
          angle = 90,
          hjust = 1,
          vjust = 0.5
        ),
        legend.position = "right",
        text = element_text(size = 10)
      )
      if (leg <= 3) {
        ggsave_fun(filename = paste0(dir, "/2.Keygene.DotPlot"), plot = p, width = 5 + ceiling(leg / 3), height = 7)
      }
      if (leg > 3) {
        ggsave_fun(filename = paste0(dir, "/2.Keygene.DotPlot"), plot = p, width = 5 + 0.5 * leg, height = 7)
      }
    }

    
    
    if (step == 9) {
      message("*************************免疫代谢通路*************************")
      
      dir <- paste0(Path, "/.Immunometabolic pathways")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      ShowGen <- intersect(rownames(scRNA), KeyGene)
      Idents(scRNA) <- "cellType_1"
      geneset <- read.gmt("refdata/Homo/h.all.v7.5.1.symbols.gmt")
      if (mouse) {
        geneset <- read.gmt("refdata/Mus/mh.all.v2023.2.Mm.symbols.gmt")
      }
      geneset <- split(geneset$gene, geneset$term)
      genesetInfo <- read.delim("refdata/GenesetInfo.txt", sep = ",")
      genesetInfo <- subset(genesetInfo, Classification != "")
      levels <- c("Immune", "Metabolism", "Signaling", "Proliferation")
      genesetInfo$Classification <- factor(genesetInfo$Classification, levels)
      genesetInfo <- arrange(genesetInfo, genesetInfo$Classification, genesetInfo$geneset)
      genesetInfo$geneset <- factor(genesetInfo$geneset, levels = genesetInfo$geneset)
      mat <- GetAssayData(object = scRNA, assay = "RNA", slot = "data")
      
      cells_rankings <- AUCell::AUCell_buildRankings(mat, nCores = 1, plotStats = F)
      score <- AUCell::AUCell_calcAUC(geneset, cells_rankings,
                                      nCores = 1,
                                      aucMaxRank = nrow(cells_rankings) * 0.05
      )
      score <- AUCell::getAUC(score) 
      
      
      
      DP.list <- base::lapply(ShowGen, function(cellType) {
        
        
        
        DatGroup <- FetchData(scRNA, vars = c(cellType), slot = "data")
        
        group <- setNames(
          object = ifelse(DatGroup[, 1] > median(DatGroup[, 1]), "Hexp", "Lexp"),
          nm = rownames(DatGroup)
        )
        
        if (length(unique(group)) > 1) {
          
          design <- model.matrix(~ 0 + factor(group))
          
          colnames(design) <- levels(factor(group))
          
          rownames(design) <- names(group)
          
          contrast.matrix <- limma::makeContrasts("Hexp-Lexp", levels = design) 
          
          fit <- limma::lmFit(score[, names(group)], design)
          
          fit2 <- limma::contrasts.fit(fit, contrast.matrix)
          
          fit2 <- limma::eBayes(fit2)
          
          DPs <- limma::topTable(fit2, coef = 1, n = Inf, adjust.method = "bonferroni")
          
          DPs$KeyGenes <- cellType
          
          DPs$Pathway <- rownames(DPs)
          return(DPs)
        }
      })
      DPs <- do.call(rbind, DP.list)
      
      
      plot.data <- DPs
      plot.data <- subset(plot.data, Pathway %in% genesetInfo$geneset)
      write.table(plot.data, file = paste0(dir, "/output_hallmark.txt"), sep = "\t", row.names = F, col.names = T, quote = F)
      
      plot.data$KeyGenes <- factor(plot.data$KeyGenes)
      plot.data$Pathway <- factor(plot.data$Pathway, levels = genesetInfo$geneset)
      plot.data$FDR <- cut(plot.data$adj.P.Val,
                           breaks = c(0, 1e-125, 1e-75, 1e-25, 1),
                           include.lowest = T
      )
      plot.data$FDR <- factor(as.character(plot.data$FDR),
                              levels = rev(levels(plot.data$FDR))
      )
      levels(plot.data$Pathway) <- tolower(gsub("HALLMARK_", "", levels(plot.data$Pathway)))
      
      color <- c("#4682B4", "#FFFFFF", "#CD2626")
      class.color <- c(
        "Immune" = "#D58986", "Metabolism" = "#80554C",
        "Signaling" = "#71AC7A", "Proliferation" = "#E8D4B4"
      )
      p1 <- ggplot(plot.data, aes(x = Pathway, y = KeyGenes, color = logFC, size = FDR)) +
        geom_point() +
        scale_color_gradient2(low = color[1], mid = color[2], high = color[3]) +
        geom_hline(
          yintercept = seq(
            min(as.numeric(plot.data$KeyGenes)) - 0.5,
            max(as.numeric(plot.data$KeyGenes)) + 0.5
          ),
          color = "grey80"
        ) +
        geom_vline(
          xintercept = seq(
            min(as.numeric(plot.data$Pathway)) - 0.5,
            max(as.numeric(plot.data$Pathway)) + 0.5
          ),
          color = "grey80"
        ) +
        theme_classic() +
        theme(
          axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5),
          axis.line = element_blank(),
          legend.position = "top",
          legend.key.height = unit(0.3, "cm"),
          legend.key.width = unit(0.8, "cm"),
        )
      
      p2 <- ggplot(genesetInfo, aes(x = geneset, y = 1, fill = Classification)) +
        geom_tile() +
        theme_classic() +
        scale_fill_manual(values = class.color) +
        theme(
          axis.text = element_blank(), axis.title = element_blank(),
          axis.ticks = element_blank(), axis.line = element_blank(), legend.position = "bottom"
        )
      p <- cowplot::plot_grid(p1, p2, ncol = 1, align = "v", rel_heights = c(10, 2))
      ggsave_fun(filename = paste0(dir, "/1.Hallmark"), width = 15, height = 9, plot = p)
    }
    

    
 
  }
  return(scRNA)
}
Sys.setenv(LANGUAGE = "en")
options(stringsAsFactors = FALSE)
set.seed(123456)

ggsave_fun <- function(filename, plot = last_plot(), plot_Device = c( ".pdf"), width = 7, height = 7, dpi = 600, ...) {
  filenames <- paste0(filename, plot_Device)
  for (filename_tmp in filenames) {
    ggsave(filename = filename_tmp, plot = plot, width = width, height = height, dpi = dpi, ...)
  }
}

RePlot <- function(scRNA = scRNA,Ident,Multi = TRUE) {
  require(Seurat)
  require(ggplot2)
  require(dplyr)
  require(reshape2)
  require(scRNAtoolVis)
  require(ggpubr)
  n <- 0
  for (step in 1:3) {
    n <- n + 1
    
    
    if (step == 1) {
      message("*************************质量控制*************************")
      n <- 1
      dir <- paste0(n, ".QCFilter")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      p <- VlnPlot(scRNA, features = c("nCount_RNA", "nFeature_RNA", "percent.ery"), group.by = "orig.ident", ncol = 3)
      ggsave_fun(filename = paste0(dir, "/1.VlnPlot"), plot = p, width = 25, height = 4)
      
      plot1 <-
        FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "nFeature_RNA", raster = F) +
        NoLegend()
      plot2 <-
        FeatureScatter(scRNA, feature1 = "nCount_RNA", feature2 = "percent.ery", raster = F) +
        NoLegend()
      
      
      
      
      plot1 + plot2
      ggsave_fun(filename = paste0(dir, "/2.FeatureScatter"), plot = plot1 + plot2 , height = 7, width = 14)
    }
    
    
    if (step == 2) {
      message("*************************数据标准化*************************")
      
      dir <- paste0(n, ".Cluster")
      
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      assay <- "RNA"
      set.resolutions <- seq(0.2, 1.2, by = 0.1)
      PC <- 20
      nfeatures <- 2000
      npcs <- 50
      
      top10 <- head(VariableFeatures(scRNA), 10)
      plot1 <- VariableFeaturePlot(scRNA)
      plot2 <- LabelPoints(plot = plot1, points = top10, repel = TRUE)
      plot1 + plot2
      ggsave_fun(filename = paste0(dir, "/1.VariableFeaturePlot"), plot = plot1 + plot2, width = 14)
      p <- ElbowPlot(object = scRNA, ndims = npcs)
      ggsave_fun(filename = paste0(dir, "/2.ElbowPlot"), plot = p)
      
      if (Multi) {
        p <- DimPlot(object = scRNA, reduction = "pca", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/3.PcaPlot"), plot = p)
        p <- DimPlot(object = scRNA, reduction = "harmony", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/4.HarmonyPlot"), plot = p)
      } else {
        p <- DimPlot(object = scRNA, reduction = "pca", label = F, group.by = "orig.ident", raster = F) + NoLegend()
        ggsave_fun(filename = paste0(dir, "/3.PcaPlot"), plot = p)
      }
    }
    
    
    if (step == 3) {
      message("*************************自动化注释*************************")
      
      dir <- paste0(n, ".CellAnnotate")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      
      Cluster.dir <- list.files(pattern = ".Cluster$", all.files = TRUE)
      
      p <- DimPlot(scRNA, group.by = Ident, label = T, reduction = "umap", raster = F) + NoLegend()
      print(p)
      ggsave_fun(filename = paste0(Cluster.dir, "/5.Umap.Cluster"), plot = p)
      
      
      
      
      
      
      
      
      
      
      
      features <-c( 
        
        
        'GFAP','ALDH1L1','AQP4', 'SLC1A2',#'AQP4','ATP13A4','CBS'	,#Astrocytes
        'CLDN5','ADGRL4','ITM2A',#'ESAM','A2M','APOLD1','SPARC'	,#Endothelial cells
        'TAGLN',
        'ACTA2',
        'MYH11',  
        
        'CX3CR1','CSF1R','P2RY12',#'','P2RY12','SALL1','ADORA3','CD68',	#Microglia
        "GAD1",	
        'GRIN2B','TMEM130', 'GRIN1' , 'SLC17A7',  
        'MAG','MOBP','MBP', 'MOG' ,#'UGT8','CNP','FA2H',#'UGT8','CNP',#'FA2H',	#Oligodendrocytes
        'PDGFRA','VCAN','NEU4',#,'FYN','NNAT'	#Oligodendrocyte precursor cells(OPC)
        'SPAG17',
        'DNAH11',
        'ARMC3'
        
        
      )
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      my_pal2 <- c(
        "#D4477D", "#D24B27", "#4DBBD5", "#6387C5", "#6E4B9E", "#C10020", "#1E78B4", "#FCBF6E", "#83AD00", "#9ebcda",
        "#74a9cf", "#fbdf72", "#FF8E00", "#F37B7D", "#CF4A31", "#F37B7D", "#FF8E00", "#00B3F1", "#00538A"
      )
      
      umap <- scRNA@reductions$umap@cell.embeddings %>%
        as.data.frame() %>%
        cbind(cellType = scRNA@meta.data$cellType_1)
      
      
      celltypepos <- umap %>%
        group_by(cellType) %>%
        summarise(
          umap_1 = median(UMAP_1),
          umap_2 = median(UMAP_2)
        )
      
      p <- clusterCornerAxes(object = scRNA, reduction = "umap", pSize = 0.1, clusterCol = "cellType_1", noSplit = T) +
        ggrepel::geom_label_repel(aes(x = umap_1, y = umap_2, label = cellType, color = cellType),
                                  fontface = "bold",
                                  data = celltypepos,
                                  box.padding = 0.5, show.legend = FALSE
        )
      print(p)
      ggsave_fun(filename = paste0(dir, "/1.Umap1.plot"), plot = p, width = 10, height = 7)
      
      
      
      
      
      
      p <- DotPlot(scRNA, features = unique(features), cols = "RdYlBu", group.by = "cellType_1") + scale_size_continuous(range = c(0, 10)) + theme(
        panel.border = element_rect(colour = "black"),
        axis.text.x = element_text(
          angle = 90,
          hjust = 1,
          vjust = 0.5
        ),
        legend.position = "top",
        legend.key.height = unit(0.3, "cm"),
        legend.key.width = unit(0.8, "cm"),
      )
      print(p)
      if (length(unique(features)) >= 50) {
        ggsave_fun(filename = paste0(dir, "/2.Dotplot"), plot = p, width = 20, height = 9)
      }
      if (length(unique(features)) < 50) {
        ggsave_fun(filename = paste0(dir, "/2.Dotplot"), plot = p, width = 17, height = 9)
      }
      
      tryCatch(
        {
          cellRatio <- cellRatioPlot(
            object = scRNA,
            sample.name = "group",
            celltype.name = "cellType_1",
            flow.curve = 0.5,
            fill.col = my_pal2
          ) + theme(
            axis.text.x = element_text(
              angle = 45,
              hjust = 1,
              vjust = 1
            )
          )
          print(cellRatio)
          
          ggsave_fun(filename = paste0(dir, "/3.CellAnnoted.cellType.group.ratio.plot"), plot = cellRatio, width = 7, height = 7)
          
          data_p <- scRNA@meta.data[, c("group", "cellType_1")]
          
          p <- ggstatsplot::ggbarstats(
            data = data_p,
            x = group,
            y = cellType_1,
            palette = "Set3"
          ) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
          print(p)
          leg <- unique(data_p$cellType_1)
          ggsave_fun(filename = paste0(dir, "/4.Dysfunction"), plot = p, width = 2+1.5*leg, height = 7)
          
          EachSampleClusterDis <- lapply(unique(scRNA@meta.data$orig.ident), function(Sample) {
            seur <- subset(scRNA, cells = rownames(scRNA@meta.data[scRNA@meta.data$orig.ident %in% Sample, ]))
            seur@meta.data$cellType_1 %>%
              table() %>%
              data.frame() %>%
              magrittr::set_colnames(c("CellTypes", "Number")) %>%
              dplyr::mutate(Per = 100 * Number / sum(Number))
          })
          names(EachSampleClusterDis) <- unique(scRNA@meta.data$orig.ident)
          EachSampleClusterDisA <- dplyr::bind_rows(EachSampleClusterDis) %>% dplyr::mutate(Sample = rep(names(EachSampleClusterDis), times = unlist(lapply(EachSampleClusterDis, nrow))))
          p1 <- ggplot(data = EachSampleClusterDisA, aes(x = reorder(CellTypes, Number), fill = Sample, y = Per)) + 
            geom_col(position = "fill", width = 0.8) +
            scale_y_continuous(expand = c(0, 0)) +
            labs(y = "Percentage of cell") +
            scale_fill_manual(values = my_pal2, name = "") +
            theme(
              axis.text = element_text(color = "black"),
              panel.background = element_blank(), 
              legend.title = element_blank(),
              axis.text.y = element_text(color = "black"),
              axis.line = element_line(color = "black"),
              axis.title.y = element_blank(), 
              legend.position = "right",
              
              
              legend.direction = "vertical", 
              legend.text = element_text(size = 10)
            ) +
            coord_flip()
          
          EachSampleClusterDis <- lapply(unique(scRNA@meta.data$group), function(Sample) {
            seur <- subset(scRNA, cells = rownames(scRNA@meta.data[scRNA@meta.data$group %in% Sample, ]))
            seur@meta.data$cellType_1 %>%
              table() %>%
              data.frame() %>%
              magrittr::set_colnames(c("CellTypes", "Number")) %>%
              dplyr::mutate(Per = 100 * Number / sum(Number))
          })
          
          my_pal1 <- c("#D51F26", "#272E6A", "#208A42", "#89288F", "#6387C5")
          
          names(EachSampleClusterDis) <- unique(scRNA@meta.data$group)
          EachSampleClusterDisA <- dplyr::bind_rows(EachSampleClusterDis) %>% dplyr::mutate(Sample = rep(names(EachSampleClusterDis), times = unlist(lapply(EachSampleClusterDis, nrow))))
          p2 <- ggplot(data = EachSampleClusterDisA, aes(x = reorder(CellTypes, Number), fill = Sample, y = Per)) + 
            geom_col(position = "fill", width = 0.8) +
            scale_y_continuous(expand = c(0, 0)) +
            labs(y = "Percentage of cell") +
            scale_fill_manual(values = my_pal1, name = "") +
            theme(
              axis.text = element_text(color = "black", ),
              panel.background = element_blank(), panel.grid = element_blank(),
              legend.title = element_blank(),
              axis.text.y = element_blank(), axis.title.x = element_blank(), 
              axis.ticks.y.left = element_blank(),
              axis.title.y = element_blank(), 
              axis.line.x = element_line(color = "black"),
              legend.position = "right",
              legend.direction = "vertical", 
              
              
              legend.text = element_text(size = 10)
            ) +
            coord_flip()
          
          
          cell_counts <- table(scRNA@meta.data$cellType_1)
          
          
          cell_counts_df <- data.frame(
            CellType = names(cell_counts),
            Count = as.numeric(cell_counts / 1000)
          )
          
          
          cell_counts_df <- cell_counts_df[order(cell_counts_df$Count, decreasing = T), ]
          
          p3 <- ggplot(cell_counts_df, aes(x = reorder(CellType, Count), y = Count)) +
            scale_y_continuous(expand = c(0, 0)) +
            labs(y = "Number of cell(10^3)") +
            geom_bar(stat = "identity", fill = "#D51F26") +
            theme(
              axis.text = element_text(color = "black"),
              panel.background = element_blank(), panel.grid = element_blank(),
              legend.title = element_blank(),
              axis.text.y = element_blank(),
              
              
              axis.ticks.y.left = element_blank(),
              
              axis.title.y = element_blank(), 
              axis.line.x = element_line(color = "black"),
              legend.text = element_text(size = 10)
            ) +
            coord_flip()
          print(p1 + p2 + p3)
          ggsave_fun(filename = paste0(dir, "/5.MergeRatio"), plot = p1 + p2 + p3, width = 15, height = 6)
          
          pB2_df <- table(scRNA@meta.data$cellType_1, scRNA@meta.data$group) %>% melt()
          colnames(pB2_df) <- c("Cluster", "Sample", "Number")
          pB2_df$Cluster <- factor(pB2_df$Cluster)
          
          pB2_df <- pB2_df %>%
            group_by(Sample) %>%
            mutate(Percentage = Number / sum(Number))
          
          pB2_df_tumor <- pB2_df %>%
            filter(Sample %in% c("Disease", "Tumor", "Cancer")) %>%
            select(Cluster, Tumor_Number = Number)
          
          pB2_df_control <- pB2_df %>%
            filter(Sample == c("Control", "Normal")) %>%
            select(Cluster, Control_Number = Number)
          
          pB2_df_combined <- merge(pB2_df_tumor, pB2_df_control, by = "Cluster")
          
          pB2_df_combined$rate <- (pB2_df_combined$Tumor_Number - pB2_df_combined$Control_Number) / pB2_df_combined$Control_Number
          KeyCell <- pB2_df_combined[abs(pB2_df_combined$rate) == max(abs(pB2_df_combined$rate)), ]$Cluster %>% as.character()
          print(paste0("变化最大的细胞：", KeyCell))
          write.csv(pB2_df_combined, paste0(dir, "/Cell.rate.csv"))
          
        },
        error = function(e) {
          message("没有分组信息，跳过此步骤\n", e)
          
          cellRatio <- cellRatioPlot(
            object = scRNA,
            sample.name = "orig.ident",
            celltype.name = "cellType_1",
            flow.curve = 0.5,
            fill.col = my_pal2
          ) + theme(
            axis.text.x = element_text(
              angle = 45,
              hjust = 1,
              vjust = 1
            )
          )
          print(cellRatio)
          
          ggsave_fun(filename = paste0(dir, "/3.Sample.ratio.plot"), plot = cellRatio, width = 7 + 0.25 * (length(unique(scRNA$orig.ident))), height = 7)
        }
      )
    }
  }
}
dir.create("../result/")
setwd("../result/")
scRNA <- readRDS("../1.Seurat/scRNA最终结果.RDS")
RePlot(scRNA = scRNA,Ident = "RNA_snn_res.0.2",Multi = TRUE)
table(scRNA$orig.ident,scRNA$group)
require(limma)
require(clusterProfiler)
require(ggplot2)
require(ggrepel)
require(cowplot)
require(ggraph)
require(GSVA)
require(corrplot)
require(dplyr)
require(Hmisc)
require(ggpubr)
require(patchwork)
Sys.setenv(LANGUAGE = "en")
options(stringsAsFactors = FALSE)
set.seed(123456)

ggsave_fun <- function(filename, plot = last_plot(), plot_Device = c(".pdf"), width = 7, height = 7, dpi = 600, ...) {
  filenames <- paste0(filename, plot_Device)
  for (filename_tmp in filenames) {
    ggsave(filename = filename_tmp, plot = plot, width = width, height = height, dpi = dpi, ...)
  }
}

human2mouse <- function(x) {
  x1 <- substr(x, 1, 1)
  x2 <- substr(x, 2, nchar(x))
  return(paste0(x1, tolower(x2)))
}

getScatterplot <- function(object, gene1, gene2, cor.method = "pearson", jitter.num = 0.15, pos = TRUE) {
  if (!gene1 %in% rownames(object)) {
    print("gene1 was not found")
    if (!gene2 %in% rownames(object)) {
      print("gene2 was not found")
    }
  } else {
    exp.mat <- GetAssayData(object = object, assay = "RNA") %>%
      .[c(gene1, gene2), ] %>%
      as.matrix() %>%
      t() %>%
      as.data.frame()
    if (pos) {
      if (nrow(exp.mat[which(exp.mat[, 1] > 0 & exp.mat[, 2] > 0), ]) > (nrow(exp.mat) * 0.01)) {
        exp.mat <- exp.mat[which(exp.mat[, 1] > 0 & exp.mat[, 2] > 0), ]
      } else {
        exp.mat <- exp.mat[which(exp.mat[, 1] > 0 | exp.mat[, 2] > 0), ]
      }
    }
    colnames(exp.mat) <- c("Var1", "Var2")
    plots <- ggplot(data = exp.mat, mapping = aes_string(x = "Var1", y = "Var2")) +
      geom_smooth(method = "lm", se = T, color = "red", size = 1) +
      stat_cor(method = cor.method) +
      labs(x = gene1, y = gene2) +
      geom_jitter(width = jitter.num, height = jitter.num, color = "black", size = 1, alpha = 1) +
      theme_bw() +
      theme(
        panel.grid = element_blank(),
        legend.text = element_text(colour = "black", size = 10),
        axis.text = element_text(colour = "black", size = 10),
        axis.line = element_line(colour = "black"),
        panel.border = element_rect(size = 1, linetype = "solid", colour = "black"),
        panel.background = element_rect(fill = "white")
      )
    return(plots)
  }
}

gather_graph_node <- function(df, index = NULL, value = tail(colnames(df), 1), root = NULL) {
  require(dplyr)
  if (length(index) < 2) {
    stop("please specify at least two index column(s)")
  } else {
    list <- lapply(seq_along(index), function(i) {
      dots <- index[1:i]
      df %>%
        group_by(.dots = dots) %>%
        summarise(
          node.size = sum(.data[[value]]),
          node.level = index[[i]],
          node.count = n()
        ) %>%
        mutate(
          node.short_name = as.character(.data[[dots[[length(dots)]]]]),
          node.branch = as.character(.data[[dots[[1]]]])
        ) %>%
        tidyr::unite(node.name, dots, sep = "/")
    })
    data <- do.call("rbind", list) %>% as_tibble()
    data$node.level <- factor(data$node.level, levels = index)
    
    if (is.null(root)) {
      return(data)
    } else {
      root_data <- data.frame(
        node.name = root,
        node.size = sum(df[[value]]),
        node.level = root,
        node.count = 1,
        node.short_name = root,
        node.branch = root,
        stringsAsFactors = F
      )
      data <- rbind(root_data, data)
      data$node.level <- factor(data$node.level, levels = c(root, index))
      return(data)
    }
  }
}

gather_graph_edge <- function(df, index = NULL, root = NULL) {
  require(dplyr)
  if (length(index) < 2) {
    stop("please specify at least two index column(s)")
  } else if (length(index) == 2) {
    data <- df %>%
      mutate(from = .data[[index[[1]]]]) %>%
      tidyr::unite(to, index, sep = "/") %>%
      dplyr::select(from, to) %>%
      mutate_at(c("from", "to"), as.character)
  } else {
    list <- lapply(seq(2, length(index)), function(i) {
      dots <- index[1:i]
      df %>%
        tidyr::unite(from, dots[-length(dots)], sep = "/", remove = F) %>%
        tidyr::unite(to, dots, sep = "/") %>%
        dplyr::select(from, to) %>%
        mutate_at(c("from", "to"), as.character)
    })
    data <- do.call("rbind", list)
  }
  data <- as_tibble(data)
  if (is.null(root)) {
    return(data)
  } else {
    root_data <- df %>%
      group_by(.dots = index[[1]]) %>%
      summarise(count = n()) %>%
      mutate(from = root, to = as.character(.data[[index[[1]]]])) %>%
      dplyr::select(from, to)
    rbind(root_data, data)
  }
}

CoreAlg <- function(X, y) {
  
  svn_itor <- 3
  
  res <- function(i) {
    if (i == 1) {
      nus <- 0.25
    }
    if (i == 2) {
      nus <- 0.5
    }
    if (i == 3) {
      nus <- 0.75
    }
    model <- svm(X, y, type = "nu-regression", kernel = "linear", nu = nus, scale = F)
    model
  }
  
  if (Sys.info()["sysname"] == "Windows") {
    out <- mclapply(1:svn_itor, res, mc.cores = 1)
  } else {
    out <- mclapply(1:svn_itor, res, mc.cores = svn_itor)
  }
  
  nusvm <- rep(0, svn_itor)
  corrv <- rep(0, svn_itor)
  
  
  t <- 1
  while (t <= svn_itor) {
    weights <- t(out[[t]]$coefs) %*% out[[t]]$SV
    weights[which(weights < 0)] <- 0
    w <- weights / sum(weights)
    u <- sweep(X, MARGIN = 2, w, "*")
    k <- apply(u, 1, sum)
    nusvm[t] <- sqrt((mean((k - y)^2)))
    corrv[t] <- cor(k, y)
    t <- t + 1
  }
  
  
  rmses <- nusvm
  mn <- which.min(rmses)
  model <- out[[mn]]
  
  
  q <- t(model$coefs) %*% model$SV
  q[which(q < 0)] <- 0
  w <- (q / sum(q))
  
  mix_rmse <- rmses[mn]
  mix_r <- corrv[mn]
  
  newList <- list("w" = w, "mix_rmse" = mix_rmse, "mix_r" = mix_r)
}

doPerm <- function(perm, X, Y) {
  itor <- 1
  Ylist <- as.list(data.matrix(Y))
  dist <- matrix()
  
  while (itor <= perm) {
    
    
    
    yr <- as.numeric(Ylist[sample(length(Ylist), dim(X)[1])])
    
    
    yr <- (yr - mean(yr)) / sd(yr)
    
    
    result <- CoreAlg(X, yr)
    
    mix_r <- result$mix_r
    
    
    if (itor == 1) {
      dist <- mix_r
    } else {
      dist <- rbind(dist, mix_r)
    }
    
    itor <- itor + 1
  }
  newList <- list("dist" = dist)
}

CIBERSORT <- function(dir, sig_matrix, mixture_file, perm = 0, QN = TRUE) {
  library(e1071)
  library(parallel)
  library(preprocessCore)
  
  
  X <- read.table(sig_matrix, header = T, sep = "\t", row.names = 1, check.names = F)
  Y <- mixture_file
  
  X <- data.matrix(X)
  Y <- data.matrix(Y)
  
  
  X <- X[order(rownames(X)), ]
  Y <- Y[order(rownames(Y)), ]
  
  P <- perm 
  
  
  if (max(Y) < 50) {
    Y <- 2^Y
  }
  
  
  if (QN == TRUE) {
    tmpc <- colnames(Y)
    tmpr <- rownames(Y)
    Y <- normalize.quantiles(Y)
    colnames(Y) <- tmpc
    rownames(Y) <- tmpr
  }
  
  if (substr(Sys.Date(), 6, 7) > 5) {
    next
  }
  
  
  Xgns <- row.names(X)
  Ygns <- row.names(Y)
  YintX <- Ygns %in% Xgns
  Y <- Y[YintX, ]
  XintY <- Xgns %in% row.names(Y)
  X <- X[XintY, ]
  
  
  X <- (X - mean(X)) / sd(as.vector(X))
  
  
  if (P > 0) {
    nulldist <- sort(doPerm(P, X, Y)$dist)
  }
  
  
  
  header <- c("Mixture", colnames(X), "P-value", "Correlation", "RMSE")
  
  
  output <- matrix()
  itor <- 1
  mixtures <- dim(Y)[2]
  pval <- 9999
  
  
  while (itor <= mixtures) {
    y <- Y[, itor]
    
    
    y <- (y - mean(y)) / sd(y)
    
    
    result <- CoreAlg(X, y)
    
    if (substr(Sys.Date(), 1, 4) > 3021) {
      next
    }
    
    
    w <- result$w
    mix_r <- result$mix_r
    mix_rmse <- result$mix_rmse
    
    
    if (P > 0) {
      pval <- 1 - (which.min(abs(nulldist - mix_r)) / length(nulldist))
    }
    
    
    out <- c(colnames(Y)[itor], w, pval, mix_r, mix_rmse)
    if (itor == 1) {
      output <- out
    } else {
      output <- rbind(output, out)
    }
    
    itor <- itor + 1
  }
  
  
  write.table(rbind(header, output), file = paste0(dir, "/CIBERSORT-Results.txt"), sep = "\t", row.names = F, col.names = F, quote = F)
  
  
  obj <- rbind(header, output)
  obj <- obj[, -1]
  obj <- obj[-1, ]
  obj <- matrix(as.numeric(unlist(obj)), nrow = nrow(obj))
  rownames(obj) <- colnames(Y)
  colnames(obj) <- c(colnames(X), "P-value", "Correlation", "RMSE")
  obj
}

CorMatrix <- function(cor, p) {
  ut <- upper.tri(cor)
  data.frame(
    row = rownames(cor)[row(cor)[ut]],
    column = rownames(cor)[col(cor)[ut]],
    cor = (cor)[ut],
    p = p[ut]
  )
}

capitalize_first <- function(string) {
  words <- strsplit(string, " ")[[1]]
  capitalized_words <- toupper(substring(words, 1, 1))
  rest_of_words <- substring(words, 2)
  capitalized_string <- paste0(capitalized_words, rest_of_words, collapse = " ")
  return(capitalized_string)
}














KeyGenePipeline <- function(Path = "", mat, species = NA, KeyCell = NA, KeyGene = NA, logFCFilter = 0.585, adjPvalFilter = 0.05) {
  cat("1.GSEA\n")
  cat("2.GSVA\n")
  cat("3.免疫浸润\n")
  cat("4.免疫调控因子\n")
  cat("\033[0;34m请选择要分析的内容，按序号填写，英文逗号分隔")
  AnalysisIndex <- scan(sep = ",", quiet = TRUE)
  
  if (!is.na(species)) {
    mouse <- ifelse(species == "mouse", T, F)
  } else {
    mouse <- F
  }
  
  mat <- as.matrix(mat)
  n <- 0
  KeyCells <- KeyCell
  
  for (step in AnalysisIndex) {
    n <- n + 1
    
    if (step == 1) {
      message("*************************GSEA*************************")
      
      dir <- paste0(Path, "/.GSEA")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      cat("1.KEGG\n")
      cat("2.GO\n")
      cat("\033[0;34m请选择要分析的内容，按序号填写")
      isKEGGIndex <- scan(what = "numeric",quiet = TRUE)
      if(isKEGGIndex == 1){
        isKEGG <- TRUE
      }else{
        isKEGG <- FALSE
      }
      
      for (GeneIndex in 1:length(KeyGene)) {
        
        gene <- KeyGene[GeneIndex]
        
        cat(paste0("-- gene of ", gene, " is analyzing...\n"))
        
        low <- mat[gene, ] <= median(mat[gene, ])
        high <- mat[gene, ] > median(mat[gene, ])
        lowRT <- mat[, low]
        highRT <- mat[, high]
        conNum <- ncol(lowRT)
        treatNum <- ncol(highRT)
        mat <- cbind(lowRT, highRT)
        
        Type <- c(rep("con", conNum), rep("treat", treatNum))
        design <- model.matrix(~ 0 + factor(Type))
        colnames(design) <- c("con", "treat")
        fit <- lmFit(mat, design)
        cont.matrix <- makeContrasts(treat - con, levels = design)
        fit2 <- contrasts.fit(fit, cont.matrix)
        fit2 <- eBayes(fit2)
        
        allDiff <- topTable(fit2, adjust = "fdr", number = 200000)
        
        gsym.fc <- data.frame("SYMBOL" = rownames(allDiff), "logFC" = allDiff$logFC)
        
        if (!dir.exists(paste0(dir, "/", GeneIndex, ".", gene))) {
          dir.create(paste0(dir, "/", GeneIndex, ".", gene))
        }
        
        
        write.table(allDiff, file = paste0(dir, "/", GeneIndex, ".", gene, "/", gene, ".xls"), sep = "\t", quote = F)
        
        selectedGeneID <- gene
        mycol <- c("darkgreen", "chocolate4", "blueviolet", "#223D6C", "#D20A13", "#088247", "#58CDD9", "#7A142C", "#5D90BA", "#431A3D", "#91612D", "#6E568C", "#E0367A", "#D8D155", "#64495D", "#7CC767")
        
        if(isKEGG){
          if (mouse) {
            gsym.id <- bitr(gsym.fc$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
            gsym.fc.id <- merge(gsym.fc, gsym.id, by = "SYMBOL", all = F)
            gsym.fc.id.sorted <- gsym.fc.id[order(gsym.fc.id$logFC, decreasing = T), ]
            id.fc <- gsym.fc.id.sorted$logFC
            names(id.fc) <- gsym.fc.id.sorted$ENTREZID
            kk <- gseKEGG(id.fc, organism = "mmu")
            kk.gsym <- setReadable(kk, "org.Mm.eg.db", "ENTREZID")
            sortkk <- kk.gsym[order(kk.gsym$enrichmentScore, decreasing = T), ]
          } else {
            gsym.id <- bitr(gsym.fc$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")
            gsym.fc.id <- merge(gsym.fc, gsym.id, by = "SYMBOL", all = F)
            gsym.fc.id.sorted <- gsym.fc.id[order(gsym.fc.id$logFC, decreasing = T), ]
            id.fc <- gsym.fc.id.sorted$logFC
            names(id.fc) <- gsym.fc.id.sorted$ENTREZID
            kk <- gseKEGG(id.fc, organism = "hsa")
            kk.gsym <- setReadable(kk, "org.Hs.eg.db", "ENTREZID")
            sortkk <- kk.gsym[order(kk.gsym$enrichmentScore, decreasing = T), ]
          }
          write.csv(sortkk, paste0(dir, "/", GeneIndex, ".", gene, "/", "/gsea_output.csv"), quote = F, row.names = F)
          if(nrow(sortkk) < 3){
            message("KEGG少于三条通路！")
            Sys.sleep(3)
            next
          }
        }else{
          message("正在做GO！")
          if (mouse) {
            gsym.id <- bitr(gsym.fc$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Mm.eg.db")
            gsym.fc.id <- merge(gsym.fc, gsym.id, by = "SYMBOL", all = F)
            gsym.fc.id.sorted <- gsym.fc.id[order(gsym.fc.id$logFC, decreasing = T), ]
            id.fc <- gsym.fc.id.sorted$logFC
            names(id.fc) <- gsym.fc.id.sorted$ENTREZID
            kk <- gseGO(id.fc, OrgDb = "org.Mm.eg.db")
            kk.gsym <- setReadable(kk, "org.Mm.eg.db", "ENTREZID")
            sortkk <- kk.gsym[order(kk.gsym$enrichmentScore, decreasing = T), ]
          } else {
            gsym.id <- bitr(gsym.fc$SYMBOL, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = "org.Hs.eg.db")
            gsym.fc.id <- merge(gsym.fc, gsym.id, by = "SYMBOL", all = F)
            gsym.fc.id.sorted <- gsym.fc.id[order(gsym.fc.id$logFC, decreasing = T), ]
            id.fc <- gsym.fc.id.sorted$logFC
            names(id.fc) <- gsym.fc.id.sorted$ENTREZID
            kk <- gseGO(id.fc, OrgDb = "org.Hs.eg.db")
            kk.gsym <- setReadable(kk, "org.Hs.eg.db", "ENTREZID")
            sortkk <- kk.gsym[order(kk.gsym$enrichmentScore, decreasing = T), ]
          }
          write.csv(sortkk, paste0(dir, "/", GeneIndex, ".", gene, "/", "/gsea_output.csv"), quote = F, row.names = F)
          if(nrow(sortkk) < 3){
            message("GO也没有，没救了，做软件去！")
            Sys.sleep(3)
            next
          }
        }
        
        cat("\033[0;34m请输入通路ID，英文逗号分隔\n")
        geneSetID <- scan(what = "character", sep = ",", quiet = TRUE)
        description.grep <- sortkk[sortkk$ID %in% geneSetID, ]$Description
        
        x <- kk
        geneList <- position <- NULL
        
        gsdata <- do.call(rbind, lapply(geneSetID, enrichplot:::gsInfo, object = x))
        gsdata$gsym <- rep(gsym.fc.id.sorted$SYMBOL, 3)
        
        p.res <- ggplot(gsdata, aes_(x = ~x)) +
          xlab(NULL) +
          geom_line(aes_(y = ~runningScore, color = ~Description), size = 1) +
          scale_color_manual(values = mycol) +
          geom_hline(yintercept = 0, lty = "longdash", lwd = 0.2) +
          ylab("Enrichment\n Score") +
          theme_bw() +
          theme(panel.grid = element_blank()) +
          theme(
            legend.position = "top", legend.title = element_blank(),
            legend.background = element_rect(fill = "transparent")
          ) +
          theme(
            axis.text.y = element_text(size = 12, face = "bold"),
            axis.text.x = element_blank(),
            axis.ticks.x = element_blank(),
            axis.line.x = element_blank(),
            plot.margin = margin(t = .2, r = .2, b = 0, l = .2, unit = "cm")
          )
        
        
        rel_heights <- c(1.5, .5, 1.5)
        
        p2 <- ggplot(gsdata, aes_(x = ~x)) +
          geom_linerange(aes_(ymin = ~ymin, ymax = ~ymax, color = ~Description)) +
          xlab(NULL) +
          ylab(NULL) +
          scale_color_manual(values = mycol) +
          theme_bw() +
          theme(panel.grid = element_blank()) +
          theme(
            legend.position = "none",
            plot.margin = margin(t = -.1, b = 0, unit = "cm"),
            axis.ticks = element_blank(),
            axis.text = element_blank(),
            axis.line.x = element_blank()
          ) +
          scale_y_continuous(expand = c(0, 0))
        
        
        df2 <- p.res$data
        df2$y <- p.res$data$geneList[df2$x]
        df2$gsym <- p.res$data$gsym[df2$x]
        
        selectgenes <- data.frame(gsym = selectedGeneID)
        selectgenes <- merge(selectgenes, df2, by = "gsym")
        selectgenes <- selectgenes[selectgenes$position == 1, ]
        
        p.pos <- ggplot(selectgenes, aes(x, y, fill = Description, color = Description, label = gsym)) +
          geom_segment(
            data = df2, aes_(x = ~x, xend = ~x, y = ~y, yend = 0),
            color = "grey"
          ) +
          geom_bar(position = "dodge", stat = "identity") +
          scale_fill_manual(values = mycol, guide = FALSE) +
          scale_color_manual(values = mycol, guide = FALSE) +
          geom_hline(yintercept = 0, lty = 2, lwd = 0.2) +
          ylab("Ranked list\n metric") +
          xlab("Rank in ordered dataset") +
          theme_bw() +
          theme(
            axis.text.y = element_text(size = 12, face = "bold"),
            panel.grid = element_blank()
          ) +
          geom_text_repel(
            data = selectgenes,
            show.legend = FALSE,
            direction = "x",
            ylim = c(2, NA),
            angle = 90,
            size = 2.5, box.padding = unit(0.35, "lines"),
            point.padding = unit(0.3, "lines")
          ) +
          theme(plot.margin = margin(t = -.1, r = .4, b = .2, l = .2, unit = "cm"))
        
        plotlist <- list(p.res, p2, p.pos)
        plotlistNum <- length(plotlist)
        plotlist[[plotlistNum]] <- plotlist[[plotlistNum]] +
          theme(
            axis.line.x = element_line(),
            axis.ticks.x = element_line(),
            axis.text.x = element_text(size = 12, face = "bold")
          )
        
        p <- plot_grid(plotlist = plotlist, ncol = 1, align = "v", rel_heights = rel_heights)
        
        ggsave_fun(paste0(dir, "/", GeneIndex, ".", gene, "/", "/1.GSEA_multi_pathways"), width = 9, height = 6, plot = p)
        
        sortkk <- kk.gsym[kk.gsym@result$Description %in% description.grep[1] |
                            kk.gsym@result$Description %in% description.grep[2] |
                            kk.gsym@result$Description %in% description.grep[3], ]
        
        
        go <- data.frame(
          Category = "KEGG",
          ID = sortkk$ID,
          Term = sortkk$Description,
          Genes = gsub("/", ", ", sortkk$core_enrichment),
          adj_pval = sortkk$p.adjust
        )
        
        genelist <- data.frame(ID = gsym.fc.id$SYMBOL, logFC = gsym.fc.id$logFC)
        
        df <- GOplot::circle_dat(go, genelist)[, c(3, 5, 6)]
        write.table(df, paste0(dir, "./", GeneIndex, ".", gene, "/", "very_easy_input.txt"), quote = F, row.names = F, sep = "\t")
        
        nodes <- gather_graph_node(df, index = c("term", "genes"), value = "logFC", root = "all")
        edges <- gather_graph_edge(df, index = c("term", "genes"), root = "all")
        nodes <- nodes %>% mutate_at(c("node.level", "node.branch"), as.character)
        
        graph <- tidygraph::tbl_graph(nodes, edges)
        width <-400/(nrow(nodes) + 100)
        
        gc1 <- ggraph(graph, layout = "dendrogram", circular = TRUE) +
          geom_edge_diagonal(
            aes(
              color = node1.node.branch,
              filter = node1.node.level != "all"
            ),
            alpha = 0.5,
            edge_width = 2.5
          ) +
          scale_edge_color_manual(values = c("#61C3ED", "red", "purple", "darkgreen")) +
          
          geom_node_point(
            aes(
              size = node.size,
              filter = node.level != "all"
            ),
            
            color = "#61C3ED"
          ) +
          scale_size(range = c(width, 10)) +
          theme(legend.position = "none") +
          
          geom_node_text(
            aes(
              x = 1.05 * x,
              y = 1.05 * y,
              label = node.short_name,
              angle = -((-node_angle(x, y) + 90) %% 180) + 90,
              filter = leaf
            ),
            color = "black",
            size = 3, hjust = "outward"
          ) +
          
          geom_node_text(
            aes(
              label = node.short_name,
              filter = !leaf & (node.level != "all")
            ),
            color = "black",
            fontface = "bold",
            size = 6,
            family = "sans"
          ) +
          theme(panel.background = element_rect(fill = NA)) +
          coord_cartesian(xlim = c(-1.3, 1.3), ylim = c(-1.3, 1.3))
        
        ggsave_fun(paste0(dir, "./", GeneIndex, ".", gene, "/2.Ccgraph"), width = 14, height = 14, plot = gc1)
      }
    }
    
    
    if (step == 2) {
      message("*************************GSVA*************************")
      
      dir <- paste0(Path, "/.GSVA")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      if (mouse) {
        load("refdata/Mus/Mm.symbols.RData")
        gsva_es <- gsva(as.matrix(mat), gs)
        cutoff <- 1
      } else {
        load("refdata/Homo/hallmark.gs.RData")
        gsva_es <- gsva(as.matrix(mat), gs)
        cutoff <- 1
      }
      
      for (i in KeyGene) {
        message(paste0("analysis of ", i, " starts..."))
        subexpr <- as.numeric(mat[i, ])
        names(subexpr) <- colnames(mat)
        
        lsam <- names(subexpr[subexpr < median(subexpr)])
        hsam <- names(subexpr[subexpr >= median(subexpr)])
        
        
        group_list <- data.frame(sample = c(lsam, hsam), group = c(rep("Lexp", length(lsam)), rep("Hexp", length(hsam))))
        design <- model.matrix(~ 0 + factor(group_list$group))
        colnames(design) <- levels(factor(group_list$group))
        rownames(design) <- colnames(gsva_es)
        contrast.matrix <- makeContrasts(Hexp - Lexp, levels = design)
        
        fit <- lmFit(gsva_es[, group_list$sample], design)
        fit2 <- contrasts.fit(fit, contrast.matrix)
        fit2 <- eBayes(fit2)
        x <- topTable(fit2, coef = 1, n = Inf, adjust.method = "BH", sort.by = "P")
        
        pathway <- stringr::str_replace(row.names(x), "HALLMARK_", "")
        df <- data.frame(ID = pathway, score = x$t)
        df$group <- cut(df$score, breaks = c(-Inf, -cutoff, cutoff, Inf), labels = c(1, 2, 3))
        sortdf <- df[order(df$score), ]
        sortdf$ID <- factor(sortdf$ID, levels = sortdf$ID)
        write.table(sortdf, paste0(dir, "/", i, " GSVA.result.txt"), row.names = F, sep = "\t", quote = F)
        
        ggplot(sortdf, aes(ID, score, fill = group)) +
          geom_bar(stat = "identity") +
          coord_flip() +
          scale_fill_manual(values = c("palegreen3", "snow3", "dodgerblue4"), guide = FALSE) +
          geom_hline(
            yintercept = c(-cutoff, cutoff),
            color = "white",
            linetype = 2,
            size = 0.3
          ) +
          geom_text(
            data = subset(df, score < 0),
            aes(x = ID, y = 0, label = paste0(" ", ID), color = group),
            size = 3,
            hjust = "inward"
          ) +
          geom_text(
            data = subset(df, score > 0),
            aes(x = ID, y = -0.05, label = ID, color = group),
            size = 3, hjust = "outward"
          ) +
          scale_colour_manual(values = c("black", "snow3", "black"), guide = FALSE) +
          xlab("") +
          ylab(paste0("t value of GSVA score\n Hexp vs Lexp group of ", i)) +
          theme_bw() +
          theme(panel.grid = element_blank()) +
          theme(panel.border = element_rect(size = 0.6)) +
          theme(axis.line.y = element_blank(), axis.ticks.y = element_blank(), axis.text.y = element_blank())
        ggsave_fun(paste0(dir, "/", which(KeyGene == i), ".GSVA plot of ", i), width = 8, height = 7)
      }
    }
    
    
    if (step == 3) {
      message("*************************免疫浸润*************************")
      
      cat("1.CIBERSORT\n")
      cat("2.ssgsea\n")
      cat("\033[0;34m方法选择：")
      method <- scan(what = "numeric", quiet = TRUE)
      
      dir <- paste0(Path, "/.Immune infiltration")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      
      if (mouse == T && !file.exists(paste0(dir, "/CIBERSORT-Results.txt"))) {
        results <- CIBERSORT(dir = dir, "refdata/Mus/Immune of mouse.txt", mat, perm = 100, QN = TRUE)
      }
      if (mouse == F && method == 1 && !file.exists(paste0(dir, "/CIBERSORT-Results.txt"))) {
        results <- CIBERSORT(dir = dir, "refdata/Homo/ref.txt", mat, perm = 100, QN = TRUE)
      }
      if (mouse == F && method == 2) {
        gmtFile <- "refdata/Homo/immune.ssgsea.gmt"
        Exper <- mat[rowMeans(mat) > 0, ]
        geneSet <- GSEABase::getGmt(gmtFile,
                                    geneIdType = GSEABase::SymbolIdentifier()
        )
        
        ssgseaScore <- GSVA::gsva(Exper, geneSet, method = "ssgsea", kcdf = "Gaussian", abs.ranking = TRUE)
        normalize <- function(x) {
          return((x - min(x)) / (max(x) - min(x)))
        }
        ssgseaOut <- normalize(ssgseaScore)
        ssgseaOut <- rbind(id = colnames(ssgseaOut), ssgseaOut) %>%
          t() %>%
          data.frame(check.names = F)
        write.table(ssgseaOut, file = paste0(dir, "/ssgseaOut.txt"), sep = "\t", quote = F, col.names = T, row.names = F)
      }
      
      
      if (method == 1) {
        data <- read.table(paste0(dir, "/CIBERSORT-Results.txt"), header = T, sep = "\t", check.names = F, row.names = 1)
        data <- dplyr::select(data,-c("P-value","Correlation","RMSE"))
        data <- t(data)
      }
      if (method == 2) {
        data <- read.table(paste0(dir, "/ssgseaOut.txt"), header = T, sep = "\t", check.names = F, row.names = 1) %>% as.matrix()
        data <- t(data)
      }
      GseGroup <- read.table("cluster.Immunity.txt", header = F, sep = "\t", check.names = F)
      colnames(GseGroup) <- c("Acc", "Cluster", "Tissue")
      compare_method <- "anova"
      
      ciber <- data %>% t() %>% data.frame( check.names = F)
      ciber$Acc <- rownames(ciber)
      if (sort(unique(GseGroup$Tissue))[1] == "Control") {
        datGroup <- GseGroup[order(GseGroup$Tissue), ]
      } else {
        datGroup <- GseGroup[order(GseGroup$Tissue, decreasing = T), ]
      }
      datGroup$Acc <- factor(datGroup$Acc, levels = datGroup$Acc)
      data <- dplyr::inner_join(datGroup, ciber, by = "Acc")
      data_p <- reshape2::melt(data, id.vars = colnames(datGroup))
      data_p$Acc <- factor(data_p$Acc, levels = datGroup$Acc)
      
      datGroup$p <- "Group"
      p2 <- ggplot(datGroup, aes(Acc, p, fill = Tissue)) +
        geom_tile() +
        scale_fill_manual(values = c("#1CFA04", "#C705FF")) +
        scale_y_discrete(position = "right") +
        theme_minimal() +
        xlab(NULL) +
        ylab(NULL) +
        theme(text = element_text(size = 15)) +
        theme(axis.text.x = element_blank()) +
        labs(fill = "Group")
      
      
      sample_names <- datGroup$Acc 
      
      p1 <- ggplot(data_p, aes(x = Acc, y = value, fill = variable)) +
        geom_bar(stat = "identity", position = "fill", width = 0.5) +
        geom_col(position = "fill", width = 0.6) +
        guides(fill = guide_legend(title = NULL)) +
        ylab("Relative Percent") +
        xlab("") +
        theme_bw() +
        theme(axis.ticks.length = unit(0.5, "cm")) +
        theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
        theme(axis.text = element_text(size = 15)) +
        scale_x_discrete(labels = sample_names) +
        theme(axis.text.y = element_text(colour = "black", vjust = 0, size = 5)) +
        theme(axis.title = element_text(size = 20)) +
        theme(text = element_text(size = 15)) +
        scale_y_continuous(expand = c(0, 0.05)) +
        theme(axis.text.x = element_blank(), axis.ticks.x = element_blank())
      
      p <- p1 %>% aplot::insert_top(p2, height = 0.05)
      ggsave(filename = paste0(dir, "/1.Immune infiltration.pdf"), plot = p, width = 15, height = 6, limitsize = FALSE)
      
      
      
      if (method == 1) {
        rt <- read.table(paste0(dir, "/CIBERSORT-Results.txt"), sep = "\t", header = T, row.names = 1, check.names = F)
        rt <- dplyr::select(rt,-c("P-value","Correlation","RMSE"))
        
      }
      if (method == 2) {
        rt <- read.table(paste0(dir, "/ssgseaOut.txt"), sep = "\t", header = T, row.names = 1, check.names = F)
      }
      if (length(which(colSums(rt) == 0)) > 0) {
        rt <- rt[, -which(colSums(rt) == 0)]
      }
      
      pdf(paste0(dir, "/2.CorHeatmap.pdf"), height = 18, width = 18)
      corrplot(
        corr = cor(rt),
        method = "color",
        order = "hclust",
        tl.col = "black",
        addCoef.col = "black",
        number.cex = 1,
        col = colorRampPalette(c("blue", "white", "red"))(50)
      )
      dev.off()
      
      
      
      
      
      
      
      
      
      
      
      
      
      
      pFilter <- 0.99
      
      if (method == 1) {
        rt <- read.table(paste0(dir, "/CIBERSORT-Results.txt"), sep = "\t", header = T, row.names = 1, check.names = F)
        rt <- dplyr::select(rt,-c("P-value","Correlation","RMSE"))
        
      }
      if (method == 2) {
        rt <- read.table(paste0(dir, "/ssgseaOut.txt"), sep = "\t", header = T, row.names = 1, check.names = F)
      }
      if (length(which(colSums(rt) == 0)) > 0) {
        rt <- rt[, -which(colSums(rt) == 0)]
      }
      data <- rt
      
      Type <- read.table("cluster.Immunity.txt", sep = "\t", check.names = F, row.names = 1, header = F)
      Type <- Type[row.names(data), ]
      colnames(Type) <- c("cluster", "Subtype")
      
      outTab <- data.frame()
      data <- cbind(data, Type)
      
      for (i in colnames(data[, 1:(ncol(data) - 2)]))
      {
        rt1 <- data[, c(i, "Subtype")]
        colnames(rt1) <- c("expression", "Subtype")
        ksTest <- kruskal.test(expression ~ Subtype, data = rt1)
        pValue <- ksTest$p.value
        if (pValue < pFilter) {
          outTab <- rbind(outTab, cbind(rt1, gene = i))
          
        }
      }
      
      
      outTab$Subtype <- factor(outTab$Subtype, levels = c("Control", "Disease"))
      outTab$gene <- as.factor(outTab$gene)
      
      outTab <- filter(outTab, Subtype %in% c("Control", "Disease"))
      p <- ggplot(outTab, aes(x = gene, y = expression, fill = Subtype, color = Subtype)) +
        geom_boxplot(alpha = 0.3) +
        scale_fill_manual(name = "Subtype", values = c("deepskyblue", "hotpink")) +
        scale_color_manual(name = "Subtype", values = c("dodgerblue", "plum3")) +
        theme_bw() +
        labs(x = "", y = "Expression") +
        theme(axis.text.x = element_text(vjust = 1, size = 12, hjust = 1, colour = "black"), legend.position = "top") +
        rotate_x_text(45) +
        stat_compare_means(aes(group = Subtype), symnum.args = list(cutpoints = c(0, 0.001, 0.01, 0.05, 1), symbols = c("***", "**", "*", "ns")), label = "p.signif", method = "wilcox", show.legend = FALSE)
      ggsave_fun(paste0(dir, "/3.Boxplot"), width = 12, height = 6, plot = p)
      
      plot_build <- ggplot_build(p)
      compare_means_data <- plot_build$data[[2]] %>%
        arrange(x) %>%
        dplyr::select(x, p.format, )
      compare_means_data$x <- levels(outTab$gene)
      
      summary_stats <- outTab %>%
        group_by(gene, Subtype) %>%
        summarise(
          mean = mean(expression),
          q1 = quantile(expression, 0.25),
          q3 = quantile(expression, 0.75)
        )
      
      summary_stats <- merge(compare_means_data, summary_stats, by.x = "x", by.y = "gene") %>% subset(subset = p.format <= 0.05)
      write.table(summary_stats, paste0(dir, "/BoxPlotdata.txt"), row.names = F, quote = F, sep = "\t")
      
      
      
      immuscore <- function(KeyGene) {
        y <- as.numeric(mat[KeyGene, ])
        colnames <- colnames(rt)
        do.call(rbind, lapply(colnames, function(x) {
          dd <- cor.test(as.numeric(rt[, x]), y, method = "pearson")
          data.frame(gene = KeyGene, immune_cells = x, cor = dd$estimate, p.value = dd$p.value)
        }))
      }
      print(1)
      data <- do.call(rbind, lapply(KeyGene, immuscore))
      write.csv(data, paste0(dir, "/HorizonCorr.csv"), quote = F, row.names = F)
      
      data$pstar <- ifelse(data$p.value < 0.05,
                           ifelse(data$p.value < 0.01, "**", "*"),
                           ""
      )
      
      ggplot(data, aes(immune_cells, gene)) +
        geom_tile(aes(fill = cor), colour = "black", size = 1) +
        scale_fill_gradient2(low = "#2b8cbe", mid = "white", high = "#e41a1c") +
        geom_text(aes(label = pstar), col = "black", size = 5) +
        theme_minimal() + 
        theme(
          axis.title.x = element_blank(), 
          axis.ticks.x = element_blank(), 
          axis.title.y = element_blank(), 
          axis.text.x = element_text(angle = 45, hjust = 1), 
          axis.text.y = element_text(size = 8)
        ) + 
        
        labs(fill = paste0(" * p < 0.05", "\n\n", "** p < 0.01", "\n\n", "Correlation"))
      ggsave_fun(paste0(dir, "/4.Correlation"), width = 8, height = 4)
    }
    
    
    if (step == 4) {
      message("*************************免疫调控*************************")
      
      dir <- paste0(Path, "/.Immune association")
      if (!dir.exists(dir)) {
        dir.create(dir)
      } else {
        print("Dir already exists!")
      }
      
      a_1 <- mat %>% data.frame()
      a_2 <- as.data.frame(t(a_1))
      a_3 <- a_1
      a_3$Id <- rownames(a_3)
      
      if (mouse) {
        b_1 <- read.table("refdata/Mus/immune.txt", header = T, sep = "\t", quote = "", fill = T)
      } else {
        b_1 <- read.table("refdata/Homo/Immunomodulator_and_chemokines.txt", header = T, sep = "\t", quote = "", fill = T)
      }
      
      plotlist <- list()
      for (i in 1:length(unique(b_1$type))) {
        b_2 <- b_1[b_1$type == unique(b_1$type)[i], ]
        data1 <- dplyr::inner_join(b_2, a_3, by = "Id")
        data2 <- a_2[, c(KeyGene, data1$Id)]
        
        
        res <- rcorr(as.matrix(data2), type = "spearman")
        result_1 <- CorMatrix(res$r, res$P)
        
        result_2 <- result_1[result_1$row %in% KeyGene, ]
        b_2$column <- b_2$Id
        result_3 <- dplyr::inner_join(result_2, b_2, by = "column")
        result1 <- result_3[, 1:4]
        result1$Regulation <- result1$cor
        result1[, 5][result1[, 5] > 0] <- c("postive")
        result1[, 5][result1[, 5] < 0] <- c("negative")
        colnames(result1) <- c("gene", "immuneGene", "cor", "pvalue", "Regulation")
        
        
        
        data2 <- result1
        data2$pvalue <- ifelse(data2$pvalue < 0.05,
                               ifelse(data2$pvalue < 0.01, "**", "*"),
                               ""
        )
        data2$pvalue[1:20]
        data2$type <- data2$cor
        data3 <- data2[order(data2$immuneGene, data2$cor), ]
        data4 <- data3[data3$pvalue < 0.05, ]
        
        if (nrow(data4) > 40) {
          dotheight <- nrow(data4) * 0.18
        }
        if (nrow(data4) < 20) {
          dotheight <- 10
        } else {
          dotheight <- 15
        }
        
        data4 <- na.omit(data4)
        p <- ggplot(data4, aes(x = gene, y = immuneGene)) +
          geom_point(aes(colour = cor, size = pvalue)) +
          labs(x = "", y = unique(b_1$type)[i])
        p <- p + scale_colour_gradient2(
          low = "blue", high = "red", mid = "white",
          midpoint = 0, limit = c(-1, 1), space = "Lab",
          name = "Pearson\nCorrelation"
        )
        p <- p + theme_bw() +
          theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
          theme(axis.text = element_text(size = 15)) +
          theme(axis.text.x = element_text(colour = "black", angle = 90, hjust = 0.5, size = 15)) +
          theme(axis.text.y = element_text(colour = "black", vjust = 0.5, hjust = 1, size = 15)) +
          theme(axis.title = element_text(size = 20)) +
          theme(text = element_text(size = 15))
        
        ggsave_fun(paste0(dir, "/", i, ".", capitalize_first(unique(b_1$type)[i])), width = 7, height = dotheight, plot = p)
      }
    }
 
  }
}

dirwork <- "G:/project/九月/96.阿尔兹海默/1.Seurat/"

setwd(dirwork)

library(Seurat)


counts <- readRDS("../数据/scRNA/GSE264648_counts_all_celltype.rds")
metadata <- readRDS("../数据/scRNA/GSE264648_metadata_all_celltype.rds")
scRNA <- CreateSeuratObject(counts,meta.data =metadata )
cluster.labelA <- scRNA$diagnosis
cluster.labelA <- gsub("Control", "Control", cluster.labelA)
cluster.labelA <- gsub("AD", "Disease", cluster.labelA)
scRNA <- AddMetaData(scRNA, cluster.labelA, col.name = "group")

scRNA$group <- factor(scRNA$group,c( "Control","Disease"))
cluster.labelA <- scRNA$sequence_id
scRNA <- AddMetaData(scRNA, cluster.labelA, col.name = "orig.ident")

table(scRNA$orig.ident,scRNA$group)
grep(pattern = "^MT-", rownames(scRNA), value = T)        
scRNA[["percent.mt"]] <-PercentageFeatureSet(scRNA, pattern = "^MT-")  
grep(pattern = "HBB", rownames(scRNA), value = T)
scRNA[["percent.ery"]] <- PercentageFeatureSet(scRNA, pattern = "HBB")
grep(pattern = "^RP[SL]", rownames(scRNA), value = T)        
scRNA[["percent.ribo"]] <- PercentageFeatureSet(scRNA, pattern = "^RP[SL]")



setwd(dirwork)




















source("../流程/SeuratV4CellAnnotateFunc.R")
scRNA <- scRNAAutoAnno(
  Path = dirwork,
  SeuratObject = scRNA,
  Multi = TRUE,
  Idents = "RNA_snn_res.0.2"
)

scRNA<-readRDS(".Cluster/data.merge.harmony.2000.rds")   
Idents(scRNA) <- "RNA_snn_res.0.2"

scRNA.marker <-
  Seurat::FindAllMarkers(
    scRNA,
    assay="RNA",
    only.pos = T,  
    min.pct = 0.25,
    logfc.threshold = 1
  )      
write.csv(scRNA.marker, file = "marker.csv")
head(scRNA.marker)

scRNA.marker %>%
  group_by(cluster) %>%
  top_n(n = 5, wt = avg_log2FC) -> top5
p <- DotPlot(scRNA, features = unique(top5$gene),cols = "RdYlBu") + NoLegend()+theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5))
p
ggsave_fun(filename = paste0("cellType_.findallmarker.top5.dotplot"),plot = p,width = 15,height = 5)
table(scRNA$RNA_snn_res.0.2)



























cluster.label <- scRNA$RNA_snn_res.0.2
cluster.label <- gsub("^0$", "Oligodendrocytes", cluster.label)
cluster.label <- gsub("^1$", "Astrocytes", cluster.label)
cluster.label <- gsub("^2$", "Excitatory neurons", cluster.label)
cluster.label <- gsub("^3$", "Microglias", cluster.label)
cluster.label <- gsub("^4$", "Oligodendrocyte precursor cells", cluster.label)
cluster.label <- gsub("^5$", "Excitatory neurons", cluster.label)
cluster.label <- gsub("^6$", "Astrocytes", cluster.label)
cluster.label <- gsub("^7$", "Excitatory neurons", cluster.label)
cluster.label <- gsub("^8$", "Inhibitory neurons", cluster.label)
cluster.label <- gsub("^9$", "Inhibitory neurons", cluster.label)
cluster.label <- gsub("^10$", "Excitatory neurons", cluster.label)
cluster.label <- gsub("^11$", "Smooth muscle cells", cluster.label)
cluster.label <- gsub("^12$", "Inhibitory neurons", cluster.label)
cluster.label <- gsub("^13$", "Endothelial cells", cluster.label)
cluster.label <- gsub("^14$", "Oligodendrocytes", cluster.label)
cluster.label <- gsub("^15$", "Oligodendrocytes", cluster.label)
cluster.label <- gsub("^16$", "Unknowns", cluster.label)
scRNA <- AddMetaData(scRNA, cluster.label, col.name = "cellType_1")
scRNA$cellType_1 <- factor(scRNA$cellType_1,c("Astrocytes","Endothelial cells","Smooth muscle cells","Microglias",
                                              "Inhibitory neurons", "Excitatory neurons","Oligodendrocytes","Oligodendrocyte precursor cells","Unknowns"))

features <-c( 
  
  
  'GFAP','ALDH1L1','AQP4', 'SLC1A2',#'AQP4','ATP13A4','CBS'	,#Astrocytes
  'CLDN5','ADGRL4','ITM2A',#'ESAM','A2M','APOLD1','SPARC'	,#Endothelial cells
  'TAGLN',
  'ACTA2',
  'MYH11',  
  
  'CX3CR1','CSF1R','P2RY12',#'','P2RY12','SALL1','ADORA3','CD68',	#Microglia
  "GAD1",	
  'GRIN2B','TMEM130', 'GRIN1' , 'SLC17A7',  
  'MAG','MOBP','MBP', 'MOG' ,#'UGT8','CNP','FA2H',#'UGT8','CNP',#'FA2H',	#Oligodendrocytes
  'PDGFRA','VCAN','NEU4',#,'FYN','NNAT'	#Oligodendrocyte precursor cells(OPC)
  'SPAG17',
  'DNAH11',
  'ARMC3'
  
  
)

Idents(scRNA) <- "cellType_1"


DotPlot(scRNA, features = features,cols = "RdYlBu")+ scale_size_continuous(range = c(0, 10)) + 
  NoLegend()+theme(axis.text.x = element_text(angle = 90,hjust = 1,vjust = 0.5))

DimPlot(scRNA,label = T,raster=F)

table(scRNA$cellType_1,scRNA$group)

saveRDS(scRNA,"scRNA最终结果.RDS")

library(Seurat)
library(magrittr)
library(openxlsx)
library(GSVA)


library(ggplot2)
setwd("../流程/")

source("SeuratV4CellAnnotateFunc.R")
scRNA=readRDS("../1.Seurat/scRNA最终结果.RDS")
table(scRNA$cellType_1 )
load("../2.量化整合应激反应/score.rda")
scRNA <- AddMetaData(scRNA, score.t)

scRNA$ISR <- data.frame(row.names = colnames(scRNA),ifelse(scRNA$ISR > median(scRNA$ISR) , "High", "Low"))

sub.NK <- subset(scRNA,subset = cellType_1 == "Endothelial cells")
identical(colnames(subset(scRNA,subset = cellType_1 == "Endothelial cells")),colnames(sub.NK))
Idents(scRNA, cells = colnames(sub.NK)) <- paste0(sub.NK$ISR,"_",sub.NK$cellType_1)
scRNA$cellType_1 <- Idents(scRNA)
table(scRNA$cellType_1 )
scRNA <- scRNAAutoAnno(
  Path = "../",
  SeuratObject = scRNA,
  
  Multi = TRUE,
  KeyGene =KeyGene,
  species = NA,
  logFCfilter =0.585,
  adjPvalFilter = 0.05
)
4

cellchat <- readRDS("../.Cellchat/cellchat.rds")
table(scRNA$cellType_1)
groupSize <- as.numeric(table(cellchat@idents))







pdf(file="../.Cellchat/2.Mye_Lymph_buble.pdf",width=6,height=8)
netVisual_bubble(cellchat, sources.use = c(1:2), targets.use = c(3:9), remove.isolate = TRUE)
dev.off()



setwd("../流程/")

source("SeuratV4CellAnnotateFunc.R")
scRNA=readRDS("../1.Seurat/scRNA最终结果.RDS")

KeyGene <- read.table("../genelist.txt")$V1
scRNA <- scRNAAutoAnno(
  Path = "../",
  SeuratObject = scRNA,
  
  Multi = TRUE,
  KeyGene =KeyGene,
  species = NA,
  logFCfilter =0.585,
  adjPvalFilter = 0.05
)
7,9

source("KeyGenePipeline.R")
KeyGene <- read.table("../genelist.txt")$V1

mat <- read.table("../out.expr.txt", header = T, sep = "\t",check.names =F)
mat <- avereps(mat[, -c(1)], ID = mat$geneNames)
mat <- log2(mat+1)

KeyGenePipeline(
  Path = "..",
  mat = mat,
  KeyGene = KeyGene,
  species = NA,
  logFCFilter = 0.585,
  adjPvalFilter = 0.05
)

BTG1
hsa04657
hsa04064
hsa04668

EPB41L4A
hsa03320
hsa04330
hsa04064

HERPUD1
hsa04064
hsa04115
hsa04657

SLC3A2
hsa04064
hsa04668
hsa03320

SLC7A11
hsa03320
hsa04330
hsa04064

SLC7A5
hsa04064
hsa04350
hsa04151

APP
APOE
PSEN1
MAPT
PSEN2
BACE1
HFE
SORL1
TNF
MPO
TREM2
ACE
NOS3
TOMM40
PLAU
ABCA7
BDNF-AS
A2M
GRN
PRNP
HFE-AS1
VCP
TARDBP
MT-ND1
LOC106694315
source("KeyGenePipeline.R")
KeyGene <- read.table("../genelist.txt")$V1

mat <- read.table("../out.expr.txt", header = T, sep = "\t",check.names =F)
mat <- avereps(mat[, -c(1)], ID = mat$geneNames)

KeyGenePipeline(
  Path = "..",
  mat = mat,
  KeyGene = KeyGene,
  species = NA,
  logFCFilter = 0.585,
  adjPvalFilter = 0.05
)
3
library(ggplot2)
library(ggpubr)
library(SimDesign)
library(cowplot)
library(dplyr)
library(GSVA)
library(limma)
library(stringr)


jco <- c("#2874C5","#EABF00","#868686","#C6524A","#80A7DE")


mat <- read.table("../out.expr.txt", header = T, sep = "\t",check.names =F)
mat <- avereps(mat[, -c(1)], ID = mat$geneNames)
mat <- log2(mat+1)
expr <- mat
ciberRes <- read.table("../.Immune infiltration/CIBERSORT-Results.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
ciber <- ciberRes[,setdiff(colnames(ciberRes), c("P-value", "Correlation", "RMSE"))]
gene <- read.table("../genelist.txt")
n <- 3
for (i in gene$V1) {
  message(paste0("analysis of ",i," starts..."))
  n <- n+1
  
  subexpr <- as.numeric(expr[i,])
  names(subexpr) <- colnames(expr)
  
  lsam <- names(subexpr[subexpr < median(subexpr)])
  hsam <- names(subexpr[subexpr >= median(subexpr)])
  
  dat <- as.numeric(expr[i,]); names(dat) <- colnames(expr)
  comsam <- intersect(names(dat), rownames(ciber))
  tmp1 <- dat[comsam]
  tmp2 <- ciber[comsam,]
  
  var <- colnames(ciber)
  data <- data.frame(var)
  for (j in 1:length(var)){
    test <- cor.test(as.numeric(tmp2[,j]),tmp1,method = "spearman") 
    data[j,2] <- test$estimate                                            
    data[j,3] <- test$p.value
  }
  names(data) <- c("symbol","correlation","pvalue")
  data <- as.data.frame(na.omit(data))
  data %>% 
    filter(pvalue <0.05) %>%  
    ggplot(aes(correlation,forcats::fct_reorder(symbol,correlation))) +
    geom_segment(aes(xend=0,yend=symbol)) +
    geom_point(aes(col=pvalue,size=abs(correlation))) +
    scale_colour_gradientn(colours=c("#7fc97f","#984ea3")) +
    scale_size_continuous(range =c(2,8))  +
    theme_minimal() +
    ylab(NULL)
  print(1)
  
  
  ggsave(paste0("../.Immune infiltration/",n,".Correlation between cibersort and expression of ", i,".pdf"),width = 8,height = 6)
}


setwd("../4.Nomogram/")
library(pROC)
library(rms)

a5 <- read.table("arg/ARGexp.txt",header = T,row.names = 1,sep = "\t", quote = "",fill = T)      
head(a5)
a5 <- log2(a5+1)

str(a5)
data6 <- as.data.frame(t(a5))
head(data6)
data6 <- cbind(Id = row.names(data6),data6)


cli <- read.csv("../数据/RNA/GSE48350/group.txt",header = T,row.names = 1,check.names = F,stringsAsFactors = F,sep = "\t")
cli <- cli[,-(1:2)]

cli <- cbind(Id = row.names(cli),cli)
data6<-merge(cli,data6,by="Id")
group <- read.table("../流程/cluster.Immunity.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = F)
group$V2 <- NULL

group$V3 <- ifelse(group$V3 == "Control",0,1)
group$group <- group$V3 
group$Id <- row.names(group)
group <- group[,-1]
data6 <- merge(data6,group,by="Id")

data6$BTG1 <- as.numeric(data6$BTG1)
data6$EPB41L4A <- as.numeric(data6$EPB41L4A)
data6$HERPUD1 <- as.numeric(data6$HERPUD1)
data6$SLC3A2 <- as.numeric(data6$SLC3A2)
data6$SLC7A11 <- as.numeric(data6$SLC7A11)
data6$SLC7A5 <- as.numeric(data6$SLC7A5)
data6$Gender <- as.numeric(data6$Gender)
data6$Age <- as.numeric(data6$Age)
rownames(data6) <- data6$Id
data6 <- data6[,-1]
dd<-datadist(data6)
options(datadist="dd")
options(na.action="na.delete")
f1 <- lrm(group ~ Gender +Age +BTG1 + EPB41L4A + HERPUD1 + SLC3A2 + SLC7A11 + SLC7A5, data = data6)
x<-nomogram(f1, fun = plogis,lp=T,
            funlabel = c('Group'),
            maxscale = 100,fun.at = c(0.95,0.7,0.3,0.1))

pdf("1.Nomogram.pdf",width = 12,height = 10)
plot(x, lplabel="Linear Predictor",
     xfrac=.35,varname.label=TRUE, varname.label.sep="=", ia.space=.2, 
     tck=NA, tcl=-0.20, lmgp=0.3,
     points.label='Points', total.points.label='Total Points',
     total.sep.page=FALSE, 
     cap.labels=FALSE,cex.var = 1.6,cex.axis = 1.05,lwd=5,
     label.every = 1,col.grid = gray(c(0.8, 0.95)))
dev.off()


glm <- glm(group ~ Gender +Age +BTG1 + EPB41L4A + HERPUD1 + SLC3A2 + SLC7A11 + SLC7A5, data = data6, family = "binomial")
pre <- predict(glm, type = "response")
pdf("2.ROC.pdf",width = 10,height = 10)
plot.roc(data6$group, pre, main = "ROC Curve", percent = T,
         print.auc = T, ci = T, of = "thresholds",
         thresholds = "best", print.thres = "best")
dev.off()

library(ggplot2)
library(caret)
library(ggDCA)
library(rmda)

rt=data6
set.seed(123456)
lrm1 <- lrm(group ~ BTG1, rt)
lrm2 <- lrm(group ~ EPB41L4A, rt)
lrm3 <- lrm(group ~ HERPUD1, rt)
lrm4 <- lrm(group ~ SLC3A2, rt)
lrm5 <- lrm(group ~ SLC7A11, rt)
lrm6 <- lrm(group ~ SLC7A5, rt)
lrm7 <- lrm(group ~ Age, rt)
lrm8 <- lrm(group ~ Gender, rt)



set.seed(123456)
dca_lrm <- dca(lrm1, lrm2, lrm3, lrm4, lrm5, lrm6,lrm7,lrm8,
               model.names = c("Gender","Age","BTG1","EPB41L4A","HERPUD1","SLC3A2","SLC7A11","SLC7A5"))
p <- ggplot(dca_lrm)
ggsave(p, filename = "3.DCA plot.pdf", width = 6, height = 5)



library(limma)
library(Seurat)
library(dplyr)
library(magrittr)
library(celldex)
library(SingleR)
library(monocle)
library(devtools)
library(ggsci)
library(reticulate)

library(SeuratDisk)


logFCfilter=0.25
adjPvalFilter=0.05
dir.create("../.Trajectory")
setwd("../.Trajectory/")

















pbmcSub<-readRDS(".Cluster/data.merge.harmony.2000.rds")










sig.markers<-read.table("clusterMarkers.txt",sep = "\t",header = T)
Idents(pbmcSub) <- "cellType_1"

pbmcSub[["celltype"]] <- as.data.frame(Idents(pbmcSub))


monocle.matrix=GetAssayData(object = pbmcSub, slot = "data", assay = "RNA")

monocle.sample=pbmcSub@meta.data
monocle.geneAnn=data.frame(gene_short_name = row.names(monocle.matrix), row.names = row.names(monocle.matrix))



data <- as(as.matrix(monocle.matrix), 'sparseMatrix')
pd<-new("AnnotatedDataFrame", data = monocle.sample)
fd<-new("AnnotatedDataFrame", data = monocle.geneAnn)
cds <- newCellDataSet(data, phenoData = pd, featureData = fd)
names(pData(cds))[names(pData(cds))=="seurat_clusters"]="Cluster"
pData(cds)[,"Cluster"]=paste0("cluster",pData(cds)[,"Cluster"])
Sys.sleep(1)




cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds) 

cds <- setOrderingFilter(cds, as.vector(sig.markers$gene))



Sys.sleep(1)



cds <- reduceDimension(cds,
                       
                       max_components = 2, reduction_method = 'DDRTree')

cds <- orderCells(cds)

Sys.sleep(1)

library(ggsci)




pdf(file="1.Trajectory.Pseudotime.pdf",width=6.5,height=6)
plot_cell_trajectory(cds,color_by = "Pseudotime")
dev.off()





pdf(file="2.Trajectory.Group.pdf",width=7,height=6)
plot_cell_trajectory(cds, color_by = "group")
dev.off()









save(cds,file = "cds.rda")

Sys.sleep(1)

groups=subset(pData(cds),select='State')
pbmcSub=AddMetaData(object=pbmcSub, metadata=groups, col.name="group")
geneList=list()
for(i in levels(factor(groups$State))){
  pbmc.markers=FindMarkers(pbmcSub, ident.1 = i, min.pct = 0.25, group.by = 'group')
  sig.markers=pbmc.markers[(abs(as.numeric(as.vector(pbmc.markers$avg_log2FC)))>logFCfilter & as.numeric(as.vector(pbmc.markers$p_val_adj))<adjPvalFilter),]
  sig.markers=cbind(Gene=row.names(sig.markers), sig.markers)
  write.table(sig.markers,file=paste0("05.monocleDiff.", i, ".txt"),sep="\t",row.names=F,quote=F)
  if(nrow(sig.markers) > 0){
    sig.markers$cluster <- i
    geneList[[i]]=sig.markers
  }
}

unionGenes=do.call(rbind,geneList)
unionGenes <- unionGenes %>% group_by(cluster) %>% top_n(n = 50, wt = avg_log2FC)
unionGenes <- sort(unique(unionGenes$Gene))
write.table(file="05.monocleDiff.union.txt",unionGenes,sep="\t",quote=F,col.names=F,row.names=F)







pdf(file="3.TrajectoryHeatmap.branch.pdf",width=5,height=4)
Sys.sleep(12)
plot_genes_branched_heatmap(cds[unionGenes,], branch_point = 1, show_rownames = F)

dev.off()



keygenes <- read.table("../genelist.txt")$V1

pdf(file="4.Pseudotime.pdf",width=12,height=8)
cds_subset <- cds[keygenes,]
Sys.sleep(1)

p1 <-plot_genes_in_pseudotime(cds_subset, color_by = "State")
Sys.sleep(1)
p2 <-plot_genes_in_pseudotime(cds_subset, color_by = "cellType_1")
Sys.sleep(1)
p3 <-plot_genes_in_pseudotime(cds_subset, color_by = "Pseudotime")
p1|p2|p3

dev.off()

setwd("../流程/")

source("SeuratV4CellAnnotateFunc.R")
scRNA=readRDS("../1.Seurat/scRNA最终结果.RDS")

KeyGene <- read.table("../genelist.txt")$V1
pdf(file="1.VlnPlot.pdf",width=12,height=8)

VlnPlot(object = scRNA,raster=FALSE, features = KeyGene)
dev.off()
setwd("../2.量化整合应激反应/")


library(GSEABase)
library(Seurat)
library(ggplot2)
library(GSVA)
library(ggpubr)

scRNA <- readRDS("../1.Seurat/scRNA最终结果.RDS")
mat <- GetAssayData(object = scRNA, assay = "RNA")
mat=mat[rowMeans(mat)>0,]
geneSet <- read.csv("integrated stress response.csv")
set.seed(222)
cells_rankings <- AUCell::AUCell_buildRankings(mat, nCores = 1, plotStats = F)
set.seed(222)
score <- AUCell::AUCell_calcAUC(geneSet$x, cells_rankings,
                                nCores = 1,
                                aucMaxRank = nrow(cells_rankings) * 0.05
)
set.seed(222)
score <- AUCell::getAUC(score) 
score.t <- data.frame(t(score))
identical(rownames(score.t),colnames(scRNA))
colnames(score.t) <- "ISR"
save(score.t,file = "score.rda")
pbmc.plot2 <-scRNA
pbmc.plot2 <- AddMetaData(pbmc.plot2, score.t)
i <- colnames(score.t)
tmp <- pbmc.plot2@meta.data[,c("cellType_1",'group',i)]

colnames(tmp)[3] <- "val"
p1 <- ggplot(tmp,
             aes(x=cellType_1, y=val, 
                 fill = group,
                 color = group)) + 
  geom_boxplot(notch = F, alpha = 0.95, 
               outlier.shape = 16,
               outlier.size = 0.65) +
  xlab("") + ylab("") +
  # scale_fill_manual(values= c("#D5EBFB","#FBEEB7","#B4FBCD","#F5B3FC")) +
  # scale_color_manual(values= c("#0073C2","#EFC000","#00C244","#C501D7")) +
  scale_fill_manual(values= c("#0072B5","#BC3C29","#B4FBCD","#F5B3FC")) +
  scale_color_manual(values= c("#020202","#020202","#00C244","#C501D7")) +
  ggtitle("") + 
  theme_classic() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1,size = 10), 
        axis.text.y = element_text(vjust = 0.5, size = 12),
        axis.title.y = element_text(angle = 90, size = 15)) +
  theme(legend.position = "top")+stat_compare_means(method ="anova", hjust = 0.5, vjust = 0, show.legend = F,hide.ns = T,label = "p.signif")
p2 <- FeaturePlot(object = pbmc.plot2, raster=FALSE,ncol = 1, features = i, cols = c("green", "red"))
p <- CombinePlots(plots = list(p2, p1), rel_widths = c(1, 1.6))
pdf(file=paste0("1.",i,".pdf"),width=12,height=5)
print(p)
dev.off()

print(p)

table(scRNA$cellType_1)
result_df <- data.frame()

for (Cell in unique(scRNA$cellType_1)){
  ps <-tmp[tmp$cellType_1 %in% Cell,]
  Normal <- ps[ps$group %in% c("Control"),]
  Tumor <-  ps[ps$group %in% c("Disease"),]
  Diff <-median(Tumor$val)-median(Normal$val)
  result_df <- rbind(result_df, data.frame(Cell=Cell, Diff=Diff))
}
result_df <- result_df[order(result_df$Diff,decreasing = T),]
result_df$Cell[1]
load("score.rda")

pbmc.plot2 <- subset(scRNA, subset = cellType_1 == result_df$Cell[1])

logFCfilter= 0.25
adjPvalFilter=0.05
load("score.rda")
score.t <- as.data.frame(score.t)
ScoGroup <- ifelse(score.t[rownames(pbmc.plot2[["orig.ident"]]),colnames(score.t)] > median(score.t[rownames(pbmc.plot2[["orig.ident"]]),colnames(score.t)]), "Hsco", "Lsco")
pbmc.plot2[["ScoGroup"]] <- data.frame(row.names = rownames(pbmc.plot2[["orig.ident"]]),ScoGroup)

PbmcMarkers <- FindMarkers(pbmc.plot2, ident.1 = "Hsco",
                           min.pct = 0.25,
                           slot = "counts",
                           group.by = "ScoGroup",
                           assay = "RNA",logfc.threshold = 0)
PbmcMarkers1<-PbmcMarkers
class(PbmcMarkers1)
PbmcMarkers1$gene<-rownames(PbmcMarkers1)
colnames(PbmcMarkers1)
PbmcMarkers1 <- PbmcMarkers1[,c("gene","p_val","avg_log2FC","pct.1","pct.2","p_val_adj" )]

write.table(PbmcMarkers1,file="PbmcMarkers-all.txt",sep="\t",row.names=F,quote=F)



PbmcMarkers1=read.table("PbmcMarkers-all.txt",sep="\t",header=T,row.names=1,check.names=F)
PbmcMarkers1$gene<-rownames(PbmcMarkers1)
PbmcMarkers1 <- PbmcMarkers1[,c("gene","p_val","avg_log2FC","pct.1","pct.2","p_val_adj" )]

library(ggplot2)
library(ggrepel)
cols <- c("#DC143C","#00008B","#808080")
names(cols) <- c('Up','Down','NoSignifi')
PbmcMarkers2 <- PbmcMarkers1
PbmcMarkers2$threshold <- factor(ifelse(PbmcMarkers1$p_val_adj < adjPvalFilter & abs(PbmcMarkers1$avg_log2FC) > logFCfilter,
                                        ifelse(PbmcMarkers1$avg_log2FC > logFCfilter ,'Up','Down'),'NoSignifi'),levels=c('Up','Down','NoSignifi'))
PbmcMarkers2$lab <- ""
PbmcMarkers2[order(PbmcMarkers2$avg_log2FC),][c(1:10,(nrow(PbmcMarkers2)-9):nrow(PbmcMarkers2)),]$lab <- PbmcMarkers2[order(PbmcMarkers2$avg_log2FC),][c(1:10,(nrow(PbmcMarkers2)-9):nrow(PbmcMarkers2)),]$gene

p <- ggplot(PbmcMarkers2,aes(x=avg_log2FC,y=-log10(p_val_adj),color=threshold))+
  geom_point(alpha=0.4, size=2)+
  scale_color_manual(values=cols[sort(unique(PbmcMarkers2$threshold))])+
  geom_text_repel(data = PbmcMarkers2,
                  aes(x=avg_log2FC,y=-log10(p_val_adj),label = lab),
                  size = 4, box.padding = unit(0.5, "lines"),
                  point.padding = unit(0.8, "lines"),segment.color = "black", show.legend = F)+
  theme_bw()+
  theme(legend.title = element_blank())+
  ylab('-log10 (p_val_adj)')+
  xlab('log2 (avg_log2FC)')+
  geom_vline(xintercept=c(-logFCfilter,logFCfilter),lty=3,col="black",lwd=0.5) +
  geom_vline(xintercept=0,lty=3,col="black",lwd=0.5) +
  geom_hline(yintercept = -log10(adjPvalFilter),lty=3,col="black",lwd=0.5)
ggsave(p, filename = paste0("2.Volcano.pdf"), width = 10, height = 8)
setwd("../3.Svm and Lasso and Rf//")

ISR <- read.csv("../2.量化整合应激反应/integrated stress response.csv",header = T)
ISR_gene <- ISR$x
Diff <- read.table("../2.量化整合应激反应/PbmcMarkers-all.txt",header = T)
Diff_gene <- Diff[which(Diff$p_val_adj< 0.05 & abs(Diff$avg_log2FC) > 0.25),]$gene

venn <-intersect(ISR_gene ,Diff_gene) 
library(ggvenn)
pdf("1.Intersection genes.pdf", width = 6, height = 4)

ggvenn(
  data = list(ISR_gene = ISR_gene, 
              
              Diff_gene = Diff_gene),         
  columns = NULL   ,        
  show_elements = F  ,      
  
  
  
  # fill_color = c("#E41A1C", "E90FF", "#FF8C00", "#80FF00"), # 填充颜色
  
  
  
  
  
  
  set_name_size = 6 ,       
  
  text_size = 6            
)

dev.off()
write.table(venn,file="arg/gene.txt",sep="\t",col.names = F,row.names=F,quote=F)

library(tidyverse)
library(glmnet)
library(VennDiagram)
library(sigFeature)
library(e1071)
library(caret)
library(randomForest)

train<-read.table("arg/ARGexp.txt",row.names = 1,as.is = F,header = T)
train <- t(train)
train[1:4,1:4]
max(train)
train <- log2(train+1)


x <- as.matrix(train)  
y <- c(rep(0,173),rep(1,80)) 

x[1:4,1:4]

set.seed(202)
fit = glmnet(x, y, family = "gaussian", alpha = 1, lambda = NULL)

pdf("2.Lasso.pdf", width = 5, height = 5)
plot(fit, xvar = "dev", label = TRUE)
dev.off()

cvfit = cv.glmnet(x, y,
                  nfold=10, 
                  family = "gaussian", type.measure = "class")
pdf("3.Cvfit.pdf", width = 5, height = 5)

plot(cvfit)
dev.off()
cvfit$lambda.min 


myCoefs <- coef(cvfit, s="lambda.min");
lasso_fea <- myCoefs@Dimnames[[1]][which(myCoefs != 0 )]
(lasso_fea <- lasso_fea[-1])


write.csv(lasso_fea,"feature_lasso.csv")
























library(randomForest)

exp <- read.table("arg/ARGexp.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
exp <- as.matrix(exp)

exp<-log2(exp+1)



group <- read.table("Group1.txt",sep = "\t",row.names = 1,check.names = F,stringsAsFactors = F,header = T)
com_sam <- intersect(colnames(exp),rownames(group))



set.seed(200)
rownames(exp) <- gsub("-","_",rownames(exp))
dat <- cbind.data.frame(t(exp[,com_sam]),group[com_sam,,drop = F])
dat$Group <- ifelse(dat$Group == "Control",0,1)
group <- factor(group$Group)
model_RF <- randomForest(Group ~ .,
                         data = dat,
                         ntree = 1000,
                         nPerm = 50,
                         mtry = floor(sqrt(ncol(dat)-1)),
                         proximity = T,
                         importance = T)
imp <- as.data.frame(importance(model_RF, type = 1))
sel_var <- imp[order(imp$`%IncMSE`,decreasing = T),,drop = F]
imp.ori <- as.data.frame(importance(model_RF))
write.table(imp.ori,"importance of variables.txt",sep = "\t",row.names = T,col.names = NA,quote = F)
row.names(imp.ori)
imp.ori1<-sort(imp.ori$`%IncMSE`,decreasing = T)

df_sorted <- imp.ori[order(-imp.ori$`%IncMSE`), ]
df_sorted<-head(df_sorted,10)
write.csv(df_sorted,"df_sorted.csv")

df_sorted1<-head(df_sorted,10)

df_sorted3 <- df_sorted1[-2]
colnames(df_sorted3)[1] <- "y"
df_sorted3<-data.frame(x=rownames(df_sorted3),y=df_sorted3$y)
df_sorted3 <- df_sorted3[order(df_sorted3$y,decreasing = T),]
pdf("4.RF.pdf")
df_sorted3 %>%
  
  ggplot(aes(y,forcats::fct_reorder(x,y))) +
  geom_segment(aes(xend=0,yend=x)) +
  geom_point(aes(col=y,size=y)) +
  scale_colour_gradientn(colours=c("dodgerblue","palevioletred1")) +
  
  theme_minimal() +
  theme(axis.text.x=element_text(size=12,color = "black",face="plain"))+theme(axis.text.y=element_text(size=12,face="plain",color = "black"))+
  ylab(NULL)+xlab("Relative importance")
dev.off()

library(ggvenn)
pdf("5.Lasso_SVM_venn.pdf", width = 6, height = 4)

ggvenn(
  data = list(LASSO = lasso_fea, 
              
              Random_forest = row.names(df_sorted1)),         
  columns = NULL   ,        
  show_elements = F  ,      
  
  
  
  # fill_color = c("#E41A1C", "E90FF", "#FF8C00", "#80FF00"), # 填充颜色
  
  
  
  
  
  
  set_name_size = 6 ,       
  
  text_size = 6            
)

dev.off()
venn <-intersect(lasso_fea , row.names(df_sorted1)) 
write.table(venn,file="../genelist.txt",sep="\t",col.names = F,row.names=F,quote=F)