单细胞的批次效应

批次：细胞以不同的分组处理时，可能发生批次效应。不同单细胞样本的数据合并通常有多种方式，包含锚定点法、Harmony、RPCA和SCTransform等。可以根据实际情况选择需要的合并方式。下面分别介绍几种消除批次效应的合并方式。

一、锚定点法

1. 总体代码

###加载所需要的包
library(Seurat)
library(tidyverse)
library(dplyr)
library(patchwork)

x =list.files()

dir = c('BC2/', "BC21/")
names(dir) = c('BC2',  'BC21')      
dir
counts <- Read10X(data.dir =dir)
scRNA_m = CreateSeuratObject(counts,min.cells = 3, min.features = 200)
table(scRNA_m@meta.data$orig.ident)
dir[1]
scRNA_list <- SplitObject(scRNA_m, split.by = "orig.ident")
?SplitObject

scRNAlist <- list()
for(i in 1:length(dir)){
  counts <- Read10X(data.dir = dir[i])
  scRNAlist[[i]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)
}

counts <- Read10X(data.dir = "BC2/")
scRNAlist[[1]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)

counts <- Read10X(data.dir = "BC3/")
scRNAlist[[2]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)

counts <- Read10X(data.dir = "BC5/")
scRNAlist[[3]] <- CreateSeuratObject(counts, min.cells = 3, min.features =300)

##归一化数据
for (i in 1:length(scRNA_list)) {
  scRNA_list[[i]] <- NormalizeData(scRNA_list[[i]])
  scRNA_list[[i]] <- FindVariableFeatures(scRNA_list[[i]], selection.method = "vst",nfeatures = 3000)
}

library(future)
plan("multicore", workers =128)
options(future.globals.maxSize = 2000 * 1024^2)
?SelectIntegrationFeatures

features <- SelectIntegrationFeatures(object.list = scRNA_list)

?FindIntegrationAnchors

scRNA_anchors <- FindIntegrationAnchors(object.list = scRNA_list,anchor.features = features)

scRNA_inte <- IntegrateData(anchorset = scRNA_anchors)

save(scRNA1,file = "IntegrateData.Rdata")
scRNA_inte@

x=scRNA1@meta.data
z1=x$orig.ident
z2=x[,1]
z3=x[1,]
class(x)
character()
data.frame()
matrix()
list()

?DefaultAssay

DefaultAssay(scRNA_inte) <- "integrated"

scRNA_inte_scale <- ScaleData(scRNA_inte)

scRNA_inte_PCA <- RunPCA(scRNA_inte_scale, npcs = 30, verbose = T)

# t-SNE and Clustering

scRNA_inte_FN <- FindNeighbors(scRNA_inte_PCA, reduction = "pca", dims = 1:20)
scRNA_inte_FC <- FindClusters(scRNA_inte_FN, resolution = 0.8)
scRNA_UMAP <- RunUMAP(scRNA_inte_FC, reduction = "pca", dims = 1:20)
colnames(scRNA1@meta.data)
scRNA_TSNE <- RunTSNE(scRNA_inte_FC, dims = 1:20)
colnames(scRNA1@meta.data)
DimPlot(scRNA_UMAP, reduction = "umap", group.by = "orig.ident")
DimPlot(scRNA_UMAP, reduction = "umap", label = TRUE)

DimPlot(scRNA_TSNE, reduction = "tsne", group.by = "orig.ident")
DimPlot(scRNA_TSNE, reduction = "tsne", label = TRUE)

2. SplitObject()

确定各个数据集之间的 "锚点"。首先，将组合对象拆分成一个列表，每个数据集都是其中的一个元素（这是必要的，因为数据是捆绑在一起的，便于分发）。

在寻找锚点之前，会进行标准的预处理（对数正态化），并为每个锚点单独识别变量特征。请注意，Seurat 基于方差稳定变换（"vst"）实现了一种改进的变量特征选择方法

3.SelectIntegrationFeatures()

选择在不同数据集中反复变化的特征进行整合.

在整合多个数据集时选择要使用的特征。此函数根据被认为是变量的数据集数量对特征进行排序，并根据各数据集变量特征排序的中位数打破平局。它会根据这一排名返回得分最高的特征。

4. FindIntegrationAnchors()

使用 FindIntegrationAnchors() 函数（该函数将 Seurat 对象列表作为输入）识别锚点，并使用 IntegrateData() 将两个数据集整合在一起。

5. DefaultAssay()

指定我们将对修正后的数据进行下游分析。原始未经修改的数据仍保留在 "RNA "检测中。

二、Harmony

1.Harmony总体优势

官方网站：https://github.com/immunogenomics/harmony

（1）整合数据的同时对稀有细胞存在一定的敏感性；

（2）省内存；

（3）适合于复杂的单细胞分析设计，可以比较不同平台、组织和供体的细胞。

2. 总体代码

DefaultAssay(scRNA_inte) <- "RNA"
scRNA_harmony <- ScaleData(scRNA_inte)
##安装Harmony包

install.packages("devtools")
library(devtools)
install_github("immunogenomics/harmony")
library(harmony)
BiocManager::install("SingleCellExperiment")

scRNA_1 <- Read10X("BC2/")
scRNA_2 <- Read10X("BC21/")
scRNA_1 <- CreateSeuratObject(scRNA_1 ,project="sample_1",min.cells = 3, min.features = 200)
scRNA_2 <- CreateSeuratObject(scRNA_2 ,project="sample_2",min.cells = 3, min.features = 200)
?CreateSeuratObject
?merge
scRNA_harmony <- merge(scRNA_1, y=c(scRNA_2 ))
scRNA_harmony <- NormalizeData(scRNA_harmony) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA(verbose=FALSE)
?system.time
?RunHarmony
system.time({scRNA_harmony <- RunHarmony(scRNA_harmony, group.by.vars = "orig.ident")})

#指定harmony
scRNA_harmony <- FindNeighbors(scRNA_harmony, reduction = "harmony", dims = 1:15) %>% FindClusters(resolution = 0.5)

scRNA_harmony <- RunUMAP(scRNA_harmony, reduction = "harmony", dims = 1:16)

plot1 =DimPlot(scRNA_harmony, reduction = "umap",label = T) 
plot2 = DimPlot(scRNA_harmony, reduction = "umap", group.by='orig.ident') 
#combinate
plot_c <- plot1+plot2
plot_c

3. merge()

通过共同列名或行名合并两个数据帧，或执行其他版本的数据库 join 操作]。

4. 管道符

管道操作符（Pipe Operator）是一个特定的符号，它可以将前一行代码的输出传递给后一行代码作为输入，从而将原本相互独立的两行代码连接在一起。而通过不断地使用管道操作符，最终可以将多行代码写成“流”的形式。使用管道操作符既可以简化代码，又可以使代码间的逻辑关系更加清晰，还可以省去中间变量的输出。

%>%

如果一行代码需要输入的参数值刚好是它前一行的输出结果，可以使用%>%省略中间的输入过程。如果参数是位于第一的位置可以直接省略（大多数都是这种情况），其他位置的参数照常书写。参数不是位于第一的位置，则需要额外使用占位符。

当有多个参数的输入值依赖于前面代码的输出结果时，需要结合大括号{}进行使用

# 不使用管道操作符
data <- mtcars[,1]
n <- length(data)
min <- min(data)
max <- max(data)
rdn <- runif(n = n, min = min, max = max)

# 使用管道操作符
mtcars %>%
  pluck(1) %>%
  {
    n = min(.)
    min = min(.)
    max = max(.)
    runif(n = n, min = min, max = max)
  } -> rdn

原文链接：https://blog.csdn.net/weixin_54000907/article/details/114481943

5. system.time()

system.time() 函数是我们可以用来估计代码执行所需时间的工具之一。它的输出给出了三个值：用户、系统和经过的时间（以秒为单位）。用户时间是处理用户应用程序代码（R 代码）所花费的 CPU 时间。如果用户应用程序访问系统资源，则该处理时间将报告为系统时间。

6. RunHarmony()

使用 Seurat 和 SingleCellAnalysis 流程运行Harmony算法。

三、RPCA

一个稍作修改的工作流程，用于整合 scRNA-seq 数据集。不使用典型相关分析（CCA）来确定锚点，而是使用相互PCA（RPCA）。在使用 RPCA 确定任意两个数据集之间的锚点时，我们会将每个数据集投影到另一个 PCA 空间，并根据相同的互邻要求对锚点进行约束。这两种工作流程的命令基本相同，但这两种方法可以在不同的环境中应用。

1. 总体代码

# split the dataset into a list of two seurat objects (stim and CTRL)
ifnb.list <- SplitObject(ifnb, split.by = "stim")

# normalize and identify variable features for each dataset independently
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    x <- NormalizeData(x)
    x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

# select features that are repeatedly variable across datasets for integration run PCA on each
# dataset using these features
features <- SelectIntegrationFeatures(object.list = ifnb.list)
ifnb.list <- lapply(X = ifnb.list, FUN = function(x) {
    x <- ScaleData(x, features = features, verbose = FALSE)
    x <- RunPCA(x, features = features, verbose = FALSE)
})
##先缩放数据、降维，后寻找锚定点
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features, reduction = "rpca")
#整合数据，后面的代码均与锚定点法相同
immune.combined <- IntegrateData(anchorset = immune.anchors)
# specify that we will perform downstream analysis on the corrected data note that the
# original unmodified data still resides in the 'RNA' assay
DefaultAssay(immune.combined) <- "integrated"

# Run the standard workflow for visualization and clustering
immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)
# Visualization
p1 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
    repel = TRUE)
p1 + p2

2. lapply()

lapply是主力的函数。他的主要用途是，对列表（list）对象而言，你想在其内部做一个循环，并对列表中的每一个元素运用函数。

sapply是lapply的一个变体，简化了lapply的结果

apply是一个对数组进行行或列运算的函数。如果你想对矩阵或其他高维数组求和，这个函数会非常好用。

tapply是`table apply()的缩写，将函数应用于向量的子集。

mapply是`apply的多变量版本。

3. 更改整合强度

基于 rpca 的整合更为保守，在这种情况下，细胞子集（即幼稚和记忆 T 细胞）在不同实验中的配准并不完美。您可以通过增加默认设置为 5 的 k.anchor 参数来提高配准强度。将该参数提高到 20 将有助于这些细胞群的配准。

immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, anchor.features = features, reduction = "rpca",
    k.anchor = 20)
immune.combined <- IntegrateData(anchorset = immune.anchors)

immune.combined <- ScaleData(immune.combined, verbose = FALSE)
immune.combined <- RunPCA(immune.combined, npcs = 30, verbose = FALSE)
immune.combined <- RunUMAP(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindNeighbors(immune.combined, reduction = "pca", dims = 1:30)
immune.combined <- FindClusters(immune.combined, resolution = 0.5)

四、SCTtransform

SCTransform可以用于数据标准化，代替三个函数（NormalizeData, ScaleDate, FindVariableFeatures）的运行。其对测序深度的较正效果好于log标准化。（10万以内的细胞均建议使用SCT）。不能用于批次较正。

SCTranform可以选择更多的主成分（PC）。使用SCTranform在算法上更加“自信”，可以提取更多生物差异信息。

高变基因选择：默认选择3000个基因，可以涵盖之前未选择的基因。

1. 单个样本的归一化

scRNA.counts=Read10X("BC21/")

###创建Seurat对象
scRNA = CreateSeuratObject(scRNA.counts ,min.cells = 3,project="os", min.features = 300)
scRNA <- PercentageFeatureSet(scRNA, pattern = "^MT-", col.name = "percent.mt")
scRNA <- SCTransform(scRNA, vars.to.regress = "percent.mt", verbose = FALSE)

######别看SCTransform只有一个单独的函数，其实它做了：NormalizeData 、ScaleData、FindVariableFeatures 的事情
########并且也支持ScaleData的vars.to.regress

##运行的结果存储在：
scRNA@assays$SCT) 
#或者
scRNA[["SCT"]]

##降维，然后聚类

scRNA <- RunPCA(scRNA, verbose = FALSE)
scRNA <- RunUMAP(scRNA, dims = 1:30, verbose = FALSE)

scRNA <- FindNeighbors(scRNA, dims = 1:30, verbose = FALSE)
scRNA <- FindClusters(scRNA, verbose = FALSE)
DimPlot(scRNA, label = TRUE) + NoLegend()

2. 多个样本利用SCT整合（非RPCA法）

在 Hafemeister 和 Satija 2019 年的文章中，介绍了一种基于正则化负二项式回归的 scRNA-seq 归一化改进方法。该方法被命名为 "sctransform"，避免了标准归一化工作流程的一些缺陷，包括添加伪计数和对数变换。

在整合之前，通过 SCTransform() 而不是 NormalizeData() 对数据集进行单独的规范化处理
正如我们在 SCTransform vignette 中进一步讨论的，我们通常使用 3,000 个或更多特征进行 sctransform 下游分析。
在识别锚点前运行 PrepSCTIntegration() 函数
运行 FindIntegrationAnchors() 和 IntegrateData() 时，将 normalization.method 参数设置为 SCT 值。
运行基于 sctransform 的工作流（包括集成）时，请勿运行 ScaleData() 函数。

LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform)
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
    anchor.features = features)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT")
immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)
p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
    repel = TRUE)
p1 + p2

3. 多个样本利用SCT整合（RPCA法）

SCTransform 对数据集进行归一化处理。我们可以选择将方法参数设置为 glmGamPoi，以便在 SCTransform() 中更快地估计回归参数。

LoadData("ifnb")
ifnb.list <- SplitObject(ifnb, split.by = "stim")
ifnb.list <- lapply(X = ifnb.list, FUN = SCTransform, method = "glmGamPoi")
features <- SelectIntegrationFeatures(object.list = ifnb.list, nfeatures = 3000)
ifnb.list <- PrepSCTIntegration(object.list = ifnb.list, anchor.features = features)
ifnb.list <- lapply(X = ifnb.list, FUN = RunPCA, features = features)
immune.anchors <- FindIntegrationAnchors(object.list = ifnb.list, normalization.method = "SCT",
    anchor.features = features, dims = 1:30, reduction = "rpca", k.anchor = 20)
immune.combined.sct <- IntegrateData(anchorset = immune.anchors, normalization.method = "SCT", dims = 1:30)
immune.combined.sct <- RunPCA(immune.combined.sct, verbose = FALSE)
immune.combined.sct <- RunUMAP(immune.combined.sct, reduction = "pca", dims = 1:30)
# Visualization
p1 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "stim")
p2 <- DimPlot(immune.combined.sct, reduction = "umap", group.by = "seurat_annotations", label = TRUE,
    repel = TRUE)
p1 + p2

本文遵循知识共享许可 CC-BY-SA 4.0 协议，转载请标明转自“吃了吃了”博客，地址：https://blog.ineuro.net。

THE END