单细胞测序数据整合
一、单细胞测序数据整合的原理
单细胞测序数据的整合过程具有很强的相似性，在基因序列的拼接与比较的过程中表现出高度的一致性或相似性。这种整合反映的是在两个不同实验条件下进行采集的数据集中存在一部分具有相同生物学特性的细胞群体（尽管这些簇中的基因表达强度可能不同但它们表现出高度的一致性或相似性）。
A 两组分别来自不同实验研究的单细胞测序数据中包含了具有相似生物学状态的细胞群体，在该查询数据集中存在独特的特定细胞群体
B 进行常规的相关性分析，并进行Log2标准化处理。
在同一个共享空间下，识别两个 dataset 之间彼此最近的 MNN 即可作为这两个 dataset 之间的 anchor 点，从而促进 dataset 的整合。
D 对于每一个anchor对，会给予一致性给出打分值
E 以这些 anchor 细胞和评分值为基础，并通过计算得到矫正向量进而完成数据集的整合
来源：[单细胞测序数据整合 - 知乎 (zhihu.com)
https://zhuanlan.zhihu.com/p/158974557?ivk_sa=1024320u](https://zhuanlan.zhihu.com/p/158974557?ivk_sa=1024320u "单个细胞水平的测序数据分析整合方法在知乎平台的具体应用")
二、一般步骤
来源：[单细胞数据整合 Comprehensive Integration - 知乎 (zhihu.com)
https://zhuanlan.zhihu.com/p/465227964](https://zhuanlan.zhihu.com/p/465227964 "Single-Cell Omics Data Integration - 知乎 (zhihu.com)")
https://zhuanlan.zhihu.com/p/465227964](https://zhuanlan.zhihu.com/p/465227964 "Single-Cell Omics Data Integration - 知乎 (zhihu.com)")
（1）Normalization
所有分析均需对单细胞RNA-Seq数据集执行标准化的预处理步骤。除非另有特别说明，在后续分析中我们默认对所有数据集进行对数正态化（lognormalization）。该方法的主要目的是防止数据分析结果受到高度表达基因的影响，在经过log-normalization处理后，该方法确保了各基因表达水平与其所引起的差异相互独立；同时显著降低了那些低表达基因的噪声干扰。
接着，在完成标准化转换后（它是执行主成分分析等降维操作之前的必要步骤）。标准化转换的目的是可见于下图所示的技术基础下的结果差异处理方法。基于技术特性，在各细胞检测到的转录水平存在差异的情况下（即数据存在噪声干扰），经过标准化转换后（使得各细胞间的转录水平具有可比性）。
（2）Feature selection for individual datasets
对于每一个数据集而言，在其中我们需要从各个维度中筛选出一系列关键基因作为候选特征集合，并通过统计学方法验证其生物学意义和临床价值
对于每一个基因，在所有细胞中计算其标准化值的方差，并利用这一结果对特征进行排序以确定前2000个具有最大标准差异性的基因，并将其选为"高度变异性的基因"。通过Seurat v3中的FindVariableFeatures函数来完成这一过程。
(3) Feature selection for integrated analysis of multiple datasets
当我们整合所有数据集时, 我们要特别关注高变基因的重要性. 为此, 我们将分别对每个独立的数据集进行特征选择. 接着统计每个基因在各独立分析中被识别为高度变异的次数, 并据此进行筛选排序, 最终选取前2000个高变基因用于 下游分析. 这些方法步骤均可通过Seurat v3软件中的SelectIntegrationFeatures函数轻松实现.
(4)Identification of anchor correspondences between two datasets
确定整合分析的关键步骤在于识别两个数据集之间的 锚点（anchors）。这些 锚点 代表每个数据集中的一个独特细胞，在这两个数据集中彼此最为高度相关。我们推测这些细胞可能源自同一个生物状态。一旦确定了这些 锚点 ，我们便能够推断出不同数据集之间的细胞关联关系。对于参考组装或迁移学习场景中 锚点 的计算方法而言，在 Seurat v3 包中可以分别调用 FindIntegrationAnchors 和 FindTransferAnchors 函数进行操作。
三、代码实现
方法一：Seurat v3的标准整合流程
（1）数据获取
#清空环境，养成好习惯
rm(list = ls())
setwd("")
library(Seurat)
library(multtest)
library(dplyr)
library(ggplot2)
library(patchwork)
library(SeuratData)#涵盖函数InstallData和LoadData
library(tidyverse)该测试数据集选自Seurat库中的ifnb数据集，在其中包含了两个来自PBMC的样本中进行分析：其中一个样本是经过干扰素刺激处理的细胞群，在另一个则作为对照组进行比较研究
#install.packages('devtools')
library(devtools)
devtools::install_github('satijalab/seurat-data',force = TRUE)
InstallData("ifnb")
LoadData("ifnb")将ifnb依据stim状态分成两个列表（stim和CTRL），其中涉及两个PBMC样本（一个为干扰素刺激处理的样本、另一个为对照样本）。
ifnb.list <- SplitObject(ifnb, split.by = "stim")
C57 <- ifnb.list$CTRL
AS1 <- ifnb.list$STIM（2）数据处理
分别进行数据标准化＋选择高变基因
myfunction1 <- function(testA.seu){
testA.seu <- NormalizeData(testA.seu, normalization.method = "LogNormalize", scale.factor = 10000)
testA.seu <- FindVariableFeatures(testA.seu, selection.method = "vst", nfeatures = 2000)
return(testA.seu)
}
C57 <- myfunction1(C57)
AS1 <- myfunction1(AS1)The integration process involves IdentifyIntegrationAnchors to locate the data point set (anchors) within the dataset.
testAB.anchors <- FindIntegrationAnchors(object.list = list(C57,AS1), dims = 1:20)
testAB.integrated <- IntegrateData(anchorset = testAB.anchors, dims = 1:20)需要注意的是：该整合步骤相较于harmony整合方法，在处理规模较大的数据集（约几万个细胞）时会占用大量内存并耗费时间，在测试案例中发现超过9GB的数据会导致32GB内存不足以支持处理过程；此外，在某些情况下（如某个seurat对象中的细胞数量极少）会出现错误提示信息，则建议采用harmony整合方法进行处理
（3）数据分析及可视化
切换至整合后的试验设计，并将其配置为integrated类型；此数据集代表了整合后的结果，请您理解为经过批处理优化的参数配置；后续分析将基于此整合数据集进行
DfaultAssay(testAB.integrated) <- "integrated"
# # Run the standard workflow for visualization and clustering
testAB.integrated <- ScaleData(testAB.integrated, features = rownames(testAB.integrated))
testAB.integrated <- RunPCA(testAB.integrated, npcs = 50, verbose = FALSE)
testAB.integrated <- FindNeighbors(testAB.integrated, dims = 1:30)
testAB.integrated <- FindClusters(testAB.integrated, resolution = 0.5)
testAB.integrated <- RunUMAP(testAB.integrated, dims = 1:30)
testAB.integrated <- RunTSNE(testAB.integrated, dims = 1:30)
p1<- DimPlot(testAB.integrated,label = T,split.by = 'stim')#integrated
# Visualization
p2 <- DimPlot(immune.combined, reduction = "umap", group.by = "stim")
p3 <- DimPlot(immune.combined, reduction = "umap", label = TRUE, repel = TRUE)
将DefaultAssay设为RNA，则预示着后续分析将基于原始数据展开，并用于识别保守细胞类型标记。同时适合利用整合数据进行聚类分析；然而针对差异基因采用整合方法不合适，并且大部分工具仅接受原始的DE计数结果。
DefaultAssay(testAB.integrated) <- "RNA"
testAB.integrated <- ScaleData(testAB.integrated, features = rownames(testAB.integrated))
testAB.integrated <- RunPCA(testAB.integrated, npcs = 50, verbose = FALSE)
testAB.integrated <- FindNeighbors(testAB.integrated, dims = 1:30)
testAB.integrated <- FindClusters(testAB.integrated, resolution = 0.5)
testAB.integrated <- RunUMAP(testAB.integrated, dims = 1:30)
testAB.integrated <- RunTSNE(testAB.integrated, dims = 1:30)
p4 <- DimPlot(testAB.integrated,label = T,split.by = 'stim')
p1|p4
观察到基于原始值的数据分析时会引入所谓的批次效应；其对应的细胞群在不同实验批次中几乎没有任何重叠；通过Integration方法能够有效消除这种干扰；这些方法首先发现能够在不同数据集之间建立匹配生物学状态（生物锚点）的细胞对；进而利用这些生物锚点对各实验批次间的数据差异进行校正；经过上述批效消除后；同一类型的细胞群则会聚集在一起；从而使得对于不同类型的细胞群进行鉴别更加容易。
方法二、harmony的方法
主要优点是：速度快，内存小
基本原理：我们用不同颜色代表不同的数据集，并用形状标记不同的细胞类型。首先采用Harmony方法通过主成分分析将转录组表达谱嵌入到低维空间中，并随后采用迭代过程来去除各数据集特有的影响因素。
（A）该算法以概率方式将细胞分配至各聚类中，并通过此过程最大化各聚类内部的数据多样性。
（B）该算法确定所有数据集在各自聚类中的全局中心位置，并确定其各自对应的局部中心位置。
（C）通过上述步骤C所得出的结果，在每个聚类中通过确定各核心位置来计算出针对各个数据集的独特修正系数。
（D）最终该算法利用这些修正系数对单个细胞进行精确调整。由于该算法采用软聚类方法，在每一轮迭代过程中会生成新的分类结果，并利用这些结果不断优化修正系数直至收敛完成。
聚类分配和数据间关系随着迭代次数减少而逐步减弱。
来源：融合多平台数据以探索单细胞数据分析深度与广度 - 云端社群活动 - 腾讯云（tencent.com）
https://cloud.tencent.com/developer/article/1650648](https://cloud.tencent.com/developer/article/1650648 "通过‘harmony’实现了不同平台的单细胞数据融合之旅 - 云+社区 - 腾讯云 (tencent.com)")
1.数据处理
在启动HARBU软件包前，在开始处理数据时，在开始处理数据时（注：此处可能存在笔误或表述不够清晰的情况），创建一个新的Seurat对象，并应用基于主成分分析的标准PCA模型来提取特征信息。
if(!require(harmony))devtools::install_github("immunogenomics/harmony")
test.seu <- pbmc
test.seu <- test.seu%>%
Seurat::NormalizeData() %>%
FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
ScaleData()
test.seu <- RunPCA(test.seu, npcs = 50, verbose = FALSE)2.Run Harmony
通过调用RunHarmony是最简单的方法来传递Seurat对象并指定需要整合的变量。该函数后会返回一个完整的Seurat对象，并利用修正后的harmony坐标进行后续操作。建议我们将plot_convergence参数配置为True值以便追踪每个迭代中的harmony目标函数优化情况。
#####先运⾏Seurat标准流程到PCA这⼀步，然后就是Harmony整合，可以简单把这⼀步理解为⼀种新的降维########
test.seu=test.seu %>% RunHarmony("stim", plot_convergence = TRUE)
#设置reduction = 'harmony'，后续分析就会调用Harmony矫正之后的PCA降维数据。
test.seu <- test.seu %>%
RunUMAP(reduction = "harmony", dims = 1:30) %>%
FindNeighbors(reduction = "harmony", dims = 1:30) %>%
FindClusters(resolution = 0.5) %>%
identity()
test.seu <- test.seu %>%
RunTSNE(reduction = "harmony", dims = 1:30)3.可视化
#按照合并的分组进行tsne作图，group.by = "stim"
p5 <- DimPlot(test.seu, reduction = "tsne", group.by = "stim", pt.size=0.5)+theme(
axis.line = element_blank(),
axis.ticks = element_blank(),axis.text = element_blank()
)
#根据细胞类型进行tsne作图，group.by = "ident"
p6 <- DimPlot(test.seu, reduction = "tsne", group.by = "ident", pt.size=0.5, label = TRUE,repel = TRUE)+theme(
axis.line = element_blank(),
axis.ticks = element_blank(),axis.text = element_blank()
)
