单细胞测序并不一定需要去除批次效应

# load Seurat package 
library(Seurat)

dir.create("~/gzh/harmony_sct",recursive = TRUE)
setwd("~/gzh/harmony_sct/")
getwd()
list.files()


# 1 不去除批次效应，教程的步骤----
{
  pfc2.data <- Read10X(data.dir = "raw-data/pfc-sample2")
  pfc3.data <- Read10X(data.dir = "raw-data/pfc-sample3")
  pfc5.data <- Read10X(data.dir = "raw-data/pfc-sample5")
  pfc7.data <- Read10X(data.dir = "raw-data/pfc-sample7")
  
  # create a new Seurat object for each sample
  # min.cells = 3, only genes detected in at least 3 cells will be included
  # min.features = 200, only cells with at least 200 genes detected will be included
  
  pfc2 <- CreateSeuratObject(counts = pfc2.data, project = "pfc-demo", min.cells = 3, min.features = 200)
  pfc3 <- CreateSeuratObject(counts = pfc3.data, project = "pfc-demo", min.cells = 3, min.features = 200)
  pfc5 <- CreateSeuratObject(counts = pfc5.data, project = "pfc-demo", min.cells = 3, min.features = 200)
  pfc7 <- CreateSeuratObject(counts = pfc7.data, project = "pfc-demo", min.cells = 3, min.features = 200)
  
  pfc2
  
  # remove the raw data to save processing space
  rm(pfc2.data, pfc3.data, pfc5.data, pfc7.data)
  
  # inspect the metadata for one of our objects using the 'head' function
  head(pfc2@meta.data)
  
  
  pfc2@meta.data$sample_number <- "2" 
  pfc3@meta.data$sample_number <- "3" 
  pfc5@meta.data$sample_number <- "5" 
  pfc7@meta.data$sample_number <- "7" 
  
  # merge multiple Seurat objects
  pfc <- merge(x = pfc2, y = list(pfc3, pfc5, pfc7))
  
  
  # remove individual objects to save processing space
  rm(pfc2, pfc3, pfc5, pfc7)
  
  
  # inpect our new combined object
  pfc
  # an important metadata slot to add in every experiment is the ratio of mitochondrial genes
  # detected in each cell - this can be used as a proxy for cell quality in most preparations
  pfc[["percent.mt"]] <- PercentageFeatureSet(object = pfc, pattern = "^mt-")
  
  # percent.mt cutoffs typically range from 5-10% depending on the sample
  
  VlnPlot(pfc, features = c("percent.mt"), pt.size=0, y.max=15) 
  
  pfc <- subset(pfc, subset = percent.mt < 5)
  
  
  VlnPlot(pfc, features = c("nFeature_RNA"), pt.size=0, y.max=2000)
  
  pfc <- subset(pfc, subset = nFeature_RNA > 600)
  # inspect our QC metrics again
  VlnPlot(pfc, features = c("nFeature_RNA","nCount_RNA","percent.mt"), pt.size=0)
  
  
  # this may take several minutes to execute, and progress will display in the console
  pfc <- SCTransform(pfc)
  
  
  pfc <- RunPCA(pfc, npcs = 60)
  
  
  pfc <- FindNeighbors(pfc, dims = 1:60) 
  pfc <- FindClusters(pfc, resolution = 0.7)
  pfc <- RunUMAP(pfc, dims = 1:60)
  
  DimPlot(pfc, label=T)
  DimPlot(pfc, group.by="sample_number")
}

DefaultAssay(pfc)
Assays(pfc)

2# 2 harmony去除批次效应 流程---------

subset_data=pfc

DefaultAssay(subset_data)='RNA'

{
library(dplyr)


DimPlot(subset_data)


subset_data[["percent.mt"]] <- PercentageFeatureSet(subset_data, pattern = "^mt-")
VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)


subset_data = subset_data %>%
  Seurat::NormalizeData(verbose = FALSE) %>%  
  FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
  ScaleData(verbose = FALSE) %>%
  RunPCA(npcs = 30, verbose = FALSE)
ElbowPlot(subset_data, ndims = 30)
VlnPlot(subset_data, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

head(subset_data@meta.data)
library(stringr)
table(str_split(colnames(subset_data),pattern = "_",simplify = TRUE)[,2])
subset_data@meta.data$stim <-paste0('mice',str_split(colnames(subset_data),pattern = "_",simplify = TRUE)[,2])
#table(subset_data$stim)

library('harmony')

subset_data <- subset_data %>% RunHarmony("stim", plot_convergence = TRUE)
harmony_embeddings <- Embeddings(subset_data, 'harmony') 
#######################cluster
dims = 1:30
subset_data <- subset_data %>% 
  RunUMAP(reduction = "harmony", dims = dims) %>% 
  RunTSNE(reduction = "harmony", dims = dims) %>% 
  FindNeighbors(reduction = "harmony", dims = dims)

subset_data=FindClusters(subset_data,resolution =0.7)
DimPlot(subset_data,group)
head(subset_data@meta.data)
  head(pfc@meta.data)
  
  
}

3#3 标准流程-----
head(subset_data@meta.data)
All=CreateSeuratObject(counts = subset_data@assays$RNA@counts,meta.data = subset_data@meta.data)

{
  
  VlnPlot(All, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
  All = All%>%Seurat::NormalizeData(verbose = FALSE) %>%  
    FindVariableFeatures(selection.method = "vst"    ) %>%
    ScaleData(verbose = FALSE)
  All = RunPCA(All, npcs = 30, verbose = FALSE)
  ElbowPlot(All, ndims = 30)
  
  
  
#######################cluster
  All <- All %>% 
    RunUMAP(reduction = "pca", dims = 1:30) %>% 
    RunTSNE(reduction = "pca", dims = 1:30) %>% 
    FindNeighbors(reduction = "pca", dims = 1:30)
  All<-All%>% FindClusters(resolution = 0.7) %>% identity()
  
  DimPlot(All,group.by ='stim' )+ggtitle("standard_pipeline")
  
}
head(All@meta.data)