1.寫個(gè)函數(shù)

scRNA_workflow <- function(scRNA_object, scRNA_direction, ncolum){

  library(Seurat)
  library(patchwork)
  library(SingleCellExperiment)
  library(msigdbr)
  library(fgsea)
  library(dplyr)
  library(ggplot2)
  library(stringr)
  library(SingleR)
  library(psych)
  library(clustree)
  library(reshape2)
  library(paletteer)
  library(harmony)
  library(RColorBrewer) 
  library(viridis)
  library(wesanderson)
  
  setwd(scRNA_direction)
  
  scRNA_object <- Macrophage
  if(("final_cluster" %in% colnames(scRNA_object@meta.data)) & ("Group" %in% colnames(scRNA_object@meta.data))){
    cellinfo <- subset(scRNA_object@meta.data, select = c("orig.ident", "percent.mt", 
                                                          "percent.rb", "percent.HB",
                                                          "Group","final_cluster"))
    
    scRNA_new <- CreateSeuratObject(scRNA_object@assays$RNA@counts, meta.data = cellinfo)
  }else{
    
    cellinfo <- subset(scRNA_object@meta.data, select = c("orig.ident", "percent.mt"
                                                          # ,"percent.rb", "percent.HB"
                                                          #,"Sample","Group"
    ))
    
    scRNA_new <- CreateSeuratObject(scRNA_object@assays$RNA@counts, meta.data = cellinfo)
  }
  
  all.genes <- rownames(scRNA_new)
  scRNA_new <- NormalizeData(scRNA_new) %>% FindVariableFeatures(nfeatures = 3000) %>% ScaleData(features = all.genes)
  scRNA_new <- RunPCA(scRNA_new, features = VariableFeatures(object = scRNA_new),
                      verbose = F,npcs = 50)
  ElbowPlot(scRNA_new, ndims = 50)
  
  # 然后準(zhǔn)備harmony去批次
  #首先SCT整合
  scRNA_new <- SCTransform(scRNA_new)
  scRNA_new <- RunPCA(scRNA_new, npcs=50, verbose=FALSE)
  
  DefaultAssay(scRNA_new)
  scRNA_new <- RunHarmony(scRNA_new, group.by.vars="orig.ident", assay.use="SCT", max.iter.harmony = 20)
  
  
  pc.num=1:40
  
  scRNA_new <- RunTSNE(scRNA_new, reduction="harmony", dims=pc.num) %>% RunUMAP(reduction="harmony", dims=pc.num)
  
  
  p <- DimPlot(scRNA_new, group.by = "orig.ident")
  
  ggsave(paste0(scRNA_direction,"/","UMAP_Samples_harmony.pdf"), p, width = 8, height = 6)
  
  p <- DimPlot(scRNA_new, group.by = "orig.ident", split.by = "orig.ident", ncol = ncolum)
  
  ggsave(paste0(scRNA_direction,"/","UMAP_Samples_Split_harmony.pdf"), p, width = 18, height = 12)
  
  
  # 找鄰居
  scRNA_new <- FindNeighbors(scRNA_new, dims = 1:40)
  
  # 在這里對分辨率進(jìn)行循環(huán)
  scRNA_new <- FindClusters(scRNA_new, resolution = seq(0.4, 4, 0.2))
  
  p_tree = clustree(scRNA_new@meta.data, prefix = "SCT_snn_res.",node_label = "SCT_snn_res.")
  
  ggsave(plot = p_tree, filename="Tree_diff_resolution.pdf",width = 20,height = 14)
  
  saveRDS(scRNA_new, "scRNA_SCT_harmony.rds")
}

2.具體用法

這里只針對細(xì)胞亞群分析好后裤唠，需要對其中的某一亞群進(jìn)行細(xì)致分析花椭。例如成纖維細(xì)胞、巨噬細(xì)胞疲陕、淋巴細(xì)胞方淤、髓系細(xì)胞等等。

2.1 首先將上面的函數(shù)保存成一個(gè)R文件

2.2 運(yùn)行以下代碼

source("~/Function/Harmony.R")

# 創(chuàng)建目錄
dir.create("~/Endometrium/Macrophage_harmony")

#給變量賦值
scRNA_direction = "/home/yourID/Endometrium/Macrophage_harmony"

scRNA_workflow(scRNA_object = Macrophage, scRNA_direction = scRNA_direction, ncolum = 40)

2.3 等待結(jié)果