单细胞数据整合-去除批次效应harmony和CCA (学习)

目录

单细胞批次效应学习

定义

理解

常用的去批次方法-基于Seurat

1） Seurat-integration（CCA）

2） Seurat-harmony

去批次代码

①Seurat-integration（CCA）

②Seurat-harmony

---

单细胞批次效应学习

定义

在单细胞RNA测序(scRNA-seq)数据分析中，去除批次效应是至关重要的步骤，因为技术噪声和批次效应可能掩盖真正的生物学差异。批次效应来源于实验过程中的多种因素，如不同的实验批次、实验人员、仪器或实验条件的变化等。这些效应可能导致数据不一致性、增加噪声、降低实验的灵敏度和准确性，以及降低可重复性。

scRNA-seq实验成本的降低鼓励建立大型项目，如人类细胞图谱，这些项目对数千至数百万个细胞的转录组进行分析。对于如此大规模的研究，物流限制不可避免地要求数据单独生成，即在不同的时间和不同的操作员生成。数据也可以在多个实验室中使用不同的细胞解离和处理方案、文库制备技术和/或测序平台生成。所有这些因素都会导致批量效应，其中一个批次中的基因表达与另一个批次中的基因表达存在系统性差异。这种差异可能会掩盖潜在的生物学特征，或在数据中引入虚假结构;因此，为了避免产生误导性的结论，在进一步分析之前必须纠正它们。

---

批次校正-单细胞转录组分析批次效应的校正scRNA-02_scrna 去批次效应-CSDN博客

理解

描述不同的特征：

同学A [身高：180，体重：50，肤色：白色...]
同学B [身高：170，体重：55，肤色：黄色...]
同学C [身高：175，体重:60，肤色：白色...]
同学D [身高：165，体重：50，肤色：黑色...]

单细胞转录组分析来说，同一类细胞其基因表达具有相似性。

#'单细胞转录组的counts矩阵细胞A 细胞B 细胞C 细胞D
基因1   0    1    1    0
基因2   1    0    0    1
基因3   1   1     0    0
...

聚类：根据每个细胞的所有基因表达情况对细胞进行基因表达模式的分析，将相似的基因表达模式的细胞聚类在一起。

批次效应：批次效应，就是因为每个批次的实验方法等不同造成同一个类型的细胞的某些（非细胞类型的关键基因）基因表达量存在较大的差异，使得在聚类可视化的时候，明明是同一个类型的细胞，却因为批次没有聚在一起。

图红框中的细胞，根据细胞类型注释它们都是CD14+细胞，然而在最左边的可视化图中，却因为batch不同形成了分群。

样本的批次无法消除，要数据量多就需要尽可能使用不同实验室、不同机构测出来的数据。去除批次效应，只能通过各种算法，重新计算细胞之间的相似度（根据细胞类型的关键基因），将同一个类型的细胞的基因表达量校正回来。例如在肿瘤分析中，尽量减少不同病人之间同一细胞类型的误差。

---

常用的去批次方法-基于Seurat

关于去除多组数据中的批次效应，有多种算法，如Seurat包的CCA(canonical correlation analysis)、LIGER的NMF(non-negative matrix factorization)、Scran包的mnnCorrect、Seurat包的sctransform

1） Seurat-integration（CCA）

Seurat-integration是Seurat包带有的单细胞批次效应处理方法；函数：FindIntegrationAnchors

Seurat-integration方法为CCA（Canonical Correlation Analysis）分析，其原理为MNN算法。

单细胞笔记7-scRNA-seq去除批次效应的方法总结 - 简书 (jianshu.com)

源于江湾青年

• 使用CCA分析将两个数据集降维到同一个低维空间，因为CCA降维之后的空间距离不是相似性而是相关性，所以相同类型与状态的细胞可以克服技术偏倚重叠在一起。

• CCA降维之后细胞在低维空间有了可以度量的“距离”，MNN(mutual nearest neighbor)算法以此找到两个数据集之间互相“距离”最近的细胞，Seurat将这些相互最近邻细胞称为“锚点细胞”。我们用两个数据集A和B来说明锚点，假设：

  • A样本中的细胞A3与B样本中距离最近的细胞有3个（B1,B2,B3）
  • B样本中的细胞B1与A样本中距离最近的细胞有4个（A1,A2,A3,A4）
  • B样本中的细胞B2与A样本中距离最近的细胞有2个（A5,A6）
  • B样本中的细胞B3与A样本中距离最近的细胞有3个（A1,A2,A7）

从图中可以看出，A3与B1是“相互“最近邻的细胞（双箭头），而A3与B2、B3不是相互最近邻细胞，A3+B1就是A、B两个数据集中的锚点之一。实际数据中，两个数据集之间的锚点可能有几百上千个。

---

• 理想情况下，相同类型和状态的细胞才能构成配对锚点细胞，但是异常的情况也会出现，如C图中query数据集中黑色的细胞团。它在reference数据集没有相同类型的细胞，但是它也找到了锚点配对细胞（红色连线）。Seurat会通过两步过滤这些不正确的锚点：

  • 在CCA低维空间找到的锚点，返回到基因表达数据构建的高维空间中验证，如果它们的转录特征相似性高则保留，否则过滤此锚点。
  • 检查锚点细胞所在数据集最邻近的30个细胞，查看它们重叠的锚点配对细胞的数量，重叠越多分值越高，代表锚点可靠性更高。

• 经过层层过滤剩下的锚点细胞对，可以认为它们是相同类型和状态的细胞，它们之间的基因表达差异是技术偏倚引起的。Seurat计算它们的差异向量，然后用此向量校正这个锚点锚定的细胞子集的基因表达值。校正后的基因表达值即消除了技术偏倚，实现了两个单细胞数据集的整合。

Tutorial：https://satijalab.org/seurat/articles/integration_mapping.html
Paper：https://www.cell.com/cell/fulltext/S0092-8674(19)30559-8

---

2） Seurat-harmony

如何使用Harmony整合单细胞数据_哔哩哔哩_bilibili

Harmony原理

• Harmony需要输入低维空间的坐标值（embedding），一般使用PCA的降维结果。Harmony导入PCA的降维数据后，会采用soft k-means clustering算法将细胞聚类。常用的聚类算法仅考虑细胞在低维空间的距离，但是soft clustering算法会考虑我们提供的校正因素。
• 例如，细胞c2距离cluster1有点远，本来不能算作cluster1的一份子；但是c2和cluster1的细胞来自不同的数据集，因为我们期望不同的数据集融合，所以破例让它加入cluster1了。
• 聚类之后，先计算每个cluster内各个数据集的细胞的中心点，然后根据这些中心点计算各个cluster的中心点。最后通过算法让cluster内的细胞向中心聚集，实在收敛不了的离群细胞就过滤掉。调整之后的数据重复：聚类—计算cluster中心点—收敛细胞—聚类的过程，不断迭代直至聚类效果趋于稳定。

Tutorial：http://htmlpreview.github.io/?https://github.com/immunogenomics/harmony/blob/master/docs/SeuratV3.html
Paper：https://www.nature.com/articles/s41592-019-0619-0

---

去批次代码

①Seurat-integration（CCA）

library(Seurat)
sc<-"" #'Seurat对象
split.by<-"" #'需要去除批次的批次标签名字
method<-"" #'标准化的方法 log或者sct或者harmony
nfeatures<-3500 #'高变基因保留的数目，默认3500
dims<-30 #'降维的维度
resolution<-0.5 #'聚类分辨率
sc.list<-SplitObject(sc, split.by = split.by)
if (method == "log") {sc.list = lapply(sc.list, function(x) {x = NormalizeData(x) %>% FindVariableFeatures(nfeatures = nfeatures) %>% ScaleData()return(x)})sc.anchors <- FindIntegrationAnchors(object.list = sc.list, dims = 1:30)sc.combined <- IntegrateData(anchorset = sc.anchors, dims = 1:30)DefaultAssay(sc.combined) <- "integrated"sc.combined <- sc.combined %>% ScaleData() %>% RunPCA() %>% FindNeighbors(dims = 1:dims) %>% RunUMAP(dims = 1:dims) %>% FindClusters(resolution = resolution)data <- sc.combined
}

pdf("umap_seurat_clusters.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="seurat_clusters")
dev.off()
pdf("umap_groups.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="groups")
dev.off()
pdf("umap_samples.pdf",height=7,width=10)
DimPlot(data,reduction = "umap",group.by="orig.ident")
dev.off()
#'data就是去除批次后的数据

 Tutorial：https://satijalab.org/seurat/articles/integration_mapping.html

---

②Seurat-harmony

library(Seurat)
library(harmony)
sc<-"" #'Seurat对象
split.by<-"" #'需要去除批次的批次标签名字
method<-"" #'标准化的方法 log或者sct或者harmony
nfeatures<-3500 #'高变基因保留的数目，默认3500
dims<-30 #'降维的维度
resolution<-0.5 #'聚类分辨率
sc.list<-SplitObject(sc, split.by = split.by)
if(method=="harmony"){sc <- NormalizeData(sc, normalization.method = "LogNormalize", scale.factor = 10000)sc <- FindVariableFeatures(sc, selection.method = "vst", nfeatures = nfeatures)all.genes <- rownames(sc) #'ScaleData使用的基因，否则后面可能报错sc <- ScaleData(sc , features = all.genes, vars.to.regress = "nCount_RNA")sc <- RunPCA(sc, features = VariableFeatures(object = sc))sc <- RunHarmony(sc, split.by , plot_convergence = F,dims.use = 1:dims)sc.combine <- scsc.combine = RunTSNE(sc.combine, reduction = "harmony", dims = 1:dims)sc.combine = RunUMAP(sc.combine, reduction = "harmony", dims = 1:dims)sc.combine = FindNeighbors(sc.combine, reduction = "harmony",dims = 1:dims)sc.combine = FindClusters(sc.combine, resolution = resolution)data<-sc.combine
}

Tutorial：http://htmlpreview.github.io/?https://github.com/immunogenomics/harmony/blob/master/docs/SeuratV3.html

参考：

1，批次校正-单细胞转录组分析批次效应的校正scRNA-02_scrna 去批次效应-CSDN博客

2，单细胞笔记7-scRNA-seq去除批次效应的方法总结 - 简书 (jianshu.com)

3，如何使用Harmony整合单细胞数据_哔哩哔哩_bilibili

4，SCS【23】单细胞转录组之数据整合 (Harmony)-CSDN博客