（未测试）“harmony”整合不同平台的单细胞数据之旅

其实在Seurat v3官方网站的Vignettes中就曾见过该算法，但并没有太多关注，直到看了北大张泽民团队在2019年10月31日发表于Cell的《Landscap and Dynamics of Single Immune Cells in Hepatocellular Carcinoma》，为了同时整合两类数据（包括SMART-seq2和10X）（Hemberg-lab单细胞转录组数据分析（七）- 导入10X和SmartSeq2数据Tabula Muris），使不同平台的数据可以整合一起进行非监督聚类（基因共表达聚类分析和可视化），作者使用了harmony算法。

其实该算法于2018年就已经发表于bioRxiv(https://www.biorxiv.org/content/early/2018/11/04/461954) ，其算法逻辑如下图所示：

harmony算法与其他整合算法相比的优势：

（1）整合数据的同时对稀有细胞的敏感性依然很好；
（2）省内存；
（3）适合于更复杂的单细胞分析实验设计，可以比较来自不同供体，组织和技术平台的细胞。

基本原理：我们用不同颜色表示不同数据集，用形状表示不同的细胞类型。首先，Harmony应用主成分分析（一文看懂PCA主成分分析）将转录组表达谱嵌入到低维空间中，然后应用迭代过程去除数据集特有的影响。

（A）Harmony概率性地将细胞分配给cluster，从而使每个cluster内数据集的多样性最大化。
（B）Harmony计算每个cluster的所有数据集的全局中心，以及特定数据集的中心。
（C）在每个cluster中，Harmony基于中心为每个数据集计算校正因子。
（D）最后，Harmony使用基于C的特定于细胞的因子校正每个细胞。由于Harmony使用软聚类，因此可以通过多个因子的线性组合对其A中进行的软聚类分配进行线性校正，来修正每个单细胞。
重复步骤A到D，直到收敛为止。聚类分配和数据集之间的依赖性随着每一轮的减少而减小。

安装

library(devtools)
install_github("immunogenomics/harmony")

流程

我们以Seurat v3为例，使用harmony进行数据整合：

library(Seurat)
library(cowplot)
library(harmony)

首先，下载稀疏矩阵示例(https://www.dropbox.com/s/t06tptwbyn7arb6/pbmc_stim.RData?dl=1)并将其移动到文件夹下（例如data/）。

load('data/pbmc_stim.RData') #加载矩阵数据

Initialize Seurat Object

在运行Harmony之前，创建一个Seurat对象并按照标准PCA（用了这么多年的PCA可视化竟然是错的！！！）进行分析。

pbmc <- CreateSeuratObject(counts = cbind(stim.sparse, ctrl.sparse), project = "PBMC", min.cells = 5) %>%
    Seurat::NormalizeData(verbose = FALSE) %>%
    FindVariableFeatures(selection.method = "vst", nfeatures = 2000) %>%
    ScaleData(verbose = FALSE) %>%
    RunPCA(pc.genes = pbmc@var.genes, npcs = 20, verbose = FALSE)

R语言中%>%的含义是什么呢，管道函数啦，就是把左件的值发送给右件的表达式，并作为右件表达式函数的第一个参数。

pbmc@meta.data$stim <- c(rep("STIM", ncol(stim.sparse)), rep("CTRL", ncol(ctrl.sparse)))#赋值条件变量

未经校正的PC中的数据集之间存在明显差异：

options(repr.plot.height = 5, repr.plot.width = 12)
p1 <- DimPlot(object = pbmc, reduction = "pca", pt.size = .1, group.by = "stim", do.return = TRUE)
p2 <- VlnPlot(object = pbmc, features = "PC_1", group.by = "stim", do.return = TRUE, pt.size = .1)
plot_grid(p1,p2)

Run Harmony

运行Harmony的最简单方法是传递Seurat对象并指定要集成的变量。RunHarmony返回Seurat对象，并使用更正后的Harmony坐标。让我们将plot_convergence设置为TRUE，这样我们就可以确保Harmony目标函数在每一轮中都变得更好。

options(repr.plot.height = 2.5, repr.plot.width = 6)
pbmc <- pbmc %>%
RunHarmony("stim", plot_convergence = TRUE)
Harmony 1/10
Harmony 2/10
Harmony 3/10
Harmony 4/10
Harmony 5/10
Harmony 6/10
Harmony 7/10
Harmony 8/10
Harmony converged after 8 iterations

要直接访问新的Harmony embeddings，请使用Embeddings命令。

harmony_embeddings <- Embeddings(pbmc, 'harmony')
harmony_embeddings[1:5, 1:5]

让我们查看确认数据集在Harmony运行之后的前两个维度中得到很好的整合。

options(repr.plot.height = 5, repr.plot.width = 12)
p1 <- DimPlot(object = pbmc, reduction = "harmony", pt.size = .1, group.by = "stim", do.return = TRUE)
p2 <- VlnPlot(object = pbmc, features = "harmony_1", group.by = "stim", do.return = TRUE, pt.size = .1)
plot_grid(p1,p2)

Downstream analysis

许多下游分析是在低维嵌入而不是基因表达上进行的。要使用校正后的Harmony embeddings而不是PC（还在用PCA降维？快学学大牛最爱的t-SNE算法吧, 附Python/R代码），请设置reduction ='harmony'。例如，让我们使用Harmony降维后的数据执行UMAP和Nearest Neighbor分析。

pbmc <- pbmc %>%
    RunUMAP(reduction = "harmony", dims = 1:20) %>%
    FindNeighbors(reduction = "harmony", dims = 1:20) %>%
    FindClusters(resolution = 0.5) %>%
    identity()

在UMAP embedding中，我们可以看到更复杂的结构。由于我们使用harmony embeddings，因此UMAP embeddings混合得很好。

options(repr.plot.height = 4, repr.plot.width = 10)
DimPlot(pbmc, reduction = "umap", group.by = "stim", pt.size = .1, split.by = 'stim')

在这种充分混合的嵌入中，我们可以开始使用聚类分析来识别细胞类型（Celaref | 单细胞测序细胞类型注释工具）。

options(repr.plot.height = 4, repr.plot.width = 6)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = .1)

快来试一试：https://github.com/immunogenomics/harmony