这篇文章从GSE132465下载，https://www.jianshu.com/p/44d30d3194e4这篇文章详细的解释了SCTransform是干什么的，所以我用这个方法试试如何去除批次效应

常规方法跑一遍

本文数据经过过滤过的，所以不需要过滤

library(Seurat)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)
library(ggsci)

library(tidyverse)
library(SingleR)
library(harmony)
library(patchwork)

exp1 <- data.table::fread('GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix.txt',data.table = F)
test <- exp1[1:10,1:10]
rownames(exp1) <- exp1$Index
exp1 <-exp1[,-1]
experiment.aggregate <- CreateSeuratObject(
  exp1,
  project = "multi", 
  min.cells = 10,
  min.features = 200,
  names.field = 1,
  names.delim = "_")
View(experiment.aggregate@meta.data)

experiment.aggregate[["percent.mt"]] <- PercentageFeatureSet(experiment.aggregate, 
                                                             pattern = "^MT-")
#save(experiment.aggregate, file = "experiment.aggregate.Rdata")
experiment.aggregate <- NormalizeData(experiment.aggregate, 
                                      normalization.method = "LogNormalize",
                                      scale.factor = 10000)

#计算表达量变化显著的基因FindVariableFeatures
experiment.aggregate <- FindVariableFeatures(experiment.aggregate, 
                                             selection.method = "vst",
                                             nfeatures = 1000)
all.genes <- rownames(experiment.aggregate)
experiment.aggregate <- ScaleData(experiment.aggregate,
                                  features = all.genes, verbose = FALSE
)
experiment.aggregate <- RunPCA(object = experiment.aggregate, 
                               features = VariableFeatures(experiment.aggregate),
                               verbose = F,npcs = 50)
ElbowPlot(experiment.aggregate, ndims = 50)
dim.use <- 1:20
experiment.aggregate <- FindNeighbors(experiment.aggregate, dims = dim.use)
#基于细胞之间的相似性计算具体的cluster
experiment.aggregate <- FindClusters(experiment.aggregate, resolution = 0.6)
#TSNE算法聚类（降维）
experiment.aggregate <- RunTSNE(experiment.aggregate, dims = dim.use, do.fast = TRUE)
experiment.aggregate <- RunUMAP(experiment.aggregate, dims = dim.use, do.fast = TRUE)
dir.create("Overview")
source("./function.R")
refdata <- get(load("Resource/SingleR_ref/ref_BE_259s.RData"))
library(SingleR)
experiment.aggregate <- cell_identify(experiment.aggregate, refdata, output = "Overview")

##tSNE图展示结果
p4 <- DimPlot(experiment.aggregate, reduction = "tsne", label = T) + NoLegend()

p5 <- DimPlot(experiment.aggregate, reduction = "tsne", label = T,
              group.by = "SingleR") + NoLegend()
pc1 = p4|p5
pc1

p6 <- DimPlot(experiment.aggregate, reduction = "tsne", group.by = "orig.ident", 
              split.by = "orig.ident", ncol = 4) + NoLegend()
ggsave("tsne.png", p6, width = 9, height = 12)
p1 <- DimPlot(experiment.aggregate, reduction = "umap", label = T) + NoLegend()

p2 <- DimPlot(experiment.aggregate, reduction = "umap", label = T,
              group.by = "SingleR") + NoLegend()
pc = p1|p2
pc

p3 <- DimPlot(experiment.aggregate, reduction = "umap", group.by = "orig.ident", 
              split.by = "orig.ident", ncol = 4) + NoLegend()
ggsave("umap.png", p3, width = 9, height = 12)

图片.png
图片.png
图片.png
图片.png

感觉没啥批次效应

导入数据入seraut对象，常规操作一下

library(Seurat)
library(ggplot2)
library(cowplot)
library(Matrix)
library(dplyr)
library(ggsci)
library(tidyverse)
library(SingleR)
library(harmony)
library(patchwork)

exp1 <- data.table::fread('GSE132465_GEO_processed_CRC_10X_raw_UMI_count_matrix.txt',data.table = F)
test <- exp1[1:10,1:10]
rownames(exp1) <- exp1$Index
exp1 <-exp1[,-1]
experiment.aggregate <- CreateSeuratObject(
  exp1,
  project = "multi", 
  min.cells = 10,
  min.features = 200,
  names.field = 1,
  names.delim = "_")
View(experiment.aggregate@meta.data)

experiment.aggregate[["percent.mt"]] <- PercentageFeatureSet(experiment.aggregate, 
                                                             pattern = "^MT-")
save(experiment.aggregate, file = "experiment.aggregate.Rdata")
experiment.aggregate <- NormalizeData(experiment.aggregate, 
                                      normalization.method = "LogNormalize",
                                      scale.factor = 10000)

#计算表达量变化显著的基因FindVariableFeatures
experiment.aggregate <- FindVariableFeatures(experiment.aggregate, 
                                             selection.method = "vst",
                                             nfeatures = 1000)
all.genes <- rownames(experiment.aggregate)
experiment.aggregate <- ScaleData(experiment.aggregate,
                                  features = all.genes, verbose = FALSE
)

细胞周期基因

if(T){
  g2m_genes <- cc.genes$g2m.genes
  g2m_genes <- CaseMatch(search=g2m_genes, match=rownames(experiment.aggregate))
  s_genes <- cc.genes$s.genes    
  s_genes <- CaseMatch(search=s_genes, match=rownames(experiment.aggregate))
  experiment.aggregate <- CellCycleScoring(experiment.aggregate, g2m.features=g2m_genes, s.features=s_genes)
  tmp <- RunPCA(experiment.aggregate, features = c(g2m_genes, s_genes), verbose = F)
  p <- DimPlot(tmp, reduction = "pca", group.by = "orig.ident")
  ggsave("QC/CellCycle_pca.png", p, width = 8, height = 6)
  rm(tmp)
}

去除批次

experiment.aggregate <- SCTransform(experiment.aggregate, return.only.var.genes = F, 
                     vars.to.regress = c( vars.to.regress = c("S.Score", "G2M.Score")))
memory.limit()
##使用harmony整合数据
experiment.aggregate <- RunPCA(experiment.aggregate, npcs=50, verbose=FALSE)
experiment.aggregate <- RunHarmony(experiment.aggregate, group.by.vars="orig.ident", assay.use="SCT")

ElbowPlot(experiment.aggregate, ndims = 40)
pc.num=1:20
scRNA <- RunTSNE(experiment.aggregate, reduction="harmony", dims=pc.num) %>% 
  RunUMAP(reduction="harmony", dims=pc.num) %>%
  FindNeighbors(reduction="harmony", dims=pc.num) %>% 
  FindClusters(resolution=0.6)

细胞注释

dir.create("Overview")
source("./function.R")
refdata <- get(load("Resource/SingleR_ref/ref_BE_259s.RData"))
library(SingleR)
scRNA <- cell_identify(scRNA, refdata, output = "Overview")

作Umap图

p1 <- DimPlot(scRNA, reduction = "umap", label = T) + NoLegend()

p2 <- DimPlot(scRNA, reduction = "umap", label = T,
              group.by = "SingleR") + NoLegend()
pc = p1|p2
pc

图片.png

查看批次效应

p3 <- DimPlot(scRNA, reduction = "umap", group.by = "orig.ident", 
              split.by = "orig.ident", ncol = 4) + NoLegend()
ggsave("umap.png", p3, width = 9, height = 12)

图片.png

作tsne图

p4 <- DimPlot(scRNA, reduction = "tsne", label = T) + NoLegend()

p5 <- DimPlot(scRNA, reduction = "tsne", label = T,
              group.by = "SingleR") + NoLegend()
pc1 = p4|p5
pc1

图片.png

批次效应

p6 <- DimPlot(scRNA, reduction = "tsne", group.by = "orig.ident", 
              split.by = "orig.ident", ncol = 4) + NoLegend()
ggsave("tsne.png", p6, width = 9, height = 12)

图片.png