其實(shí)在Seurat v3官方網(wǎng)站的Vignettes中就曾見(jiàn)過(guò)該算法韧衣，但并沒(méi)有太多關(guān)注喇勋，直到看了北大張澤民團(tuán)隊(duì)在2019年10月31日發(fā)表于Cell的《Landscap and Dynamics of Single Immune Cells in Hepatocellular Carcinoma》纵竖，為了同時(shí)整合兩類數(shù)據(jù)（包括SMART-seq2和10X）（Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（七）- 導(dǎo)入10X和SmartSeq2數(shù)據(jù)Tabula Muris）葱她，使不同平臺(tái)的數(shù)據(jù)可以整合一起進(jìn)行非監(jiān)督聚類（基因共表達(dá)聚類分析和可視化）涵亏，作者使用了harmony算法晰韵。
其實(shí)該算法于2018年就已經(jīng)發(fā)表于bioRxiv(https://www.biorxiv.org/content/early/2018/11/04/461954) 链瓦，其算法邏輯如下圖所示：
圖1. Harmony算法概述
harmony算法與其他整合算法相比的優(yōu)勢(shì)：
（1）整合數(shù)據(jù)的同時(shí)對(duì)稀有細(xì)胞的敏感性依然很好拆魏；
（2）省內(nèi)存；
（3）適合于更復(fù)雜的單細(xì)胞分析實(shí)驗(yàn)設(shè)計(jì)慈俯，可以比較來(lái)自不同供體渤刃，組織和技術(shù)平臺(tái)的細(xì)胞。
基本原理：我們用不同顏色表示不同數(shù)據(jù)集贴膘，用形狀表示不同的細(xì)胞類型卖子。首先，Harmony應(yīng)用主成分分析（一文看懂PCA主成分分析）將轉(zhuǎn)錄組表達(dá)譜嵌入到低維空間中刑峡，然后應(yīng)用迭代過(guò)程去除數(shù)據(jù)集特有的影響洋闽。
（A）Harmony概率性地將細(xì)胞分配給cluster玄柠，從而使每個(gè)cluster內(nèi)數(shù)據(jù)集的多樣性最大化。
（B）Harmony計(jì)算每個(gè)cluster的所有數(shù)據(jù)集的全局中心诫舅，以及特定數(shù)據(jù)集的中心羽利。
（C）在每個(gè)cluster中，Harmony基于中心為每個(gè)數(shù)據(jù)集計(jì)算校正因子刊懈。
（D）最后这弧，Harmony使用基于C的特定于細(xì)胞的因子校正每個(gè)細(xì)胞。由于Harmony使用軟聚類虚汛，因此可以通過(guò)多個(gè)因子的線性組合對(duì)其A中進(jìn)行的軟聚類分配進(jìn)行線性校正匾浪，來(lái)修正每個(gè)單細(xì)胞。
重復(fù)步驟A到D卷哩，直到收斂為止蛋辈。聚類分配和數(shù)據(jù)集之間的依賴性隨著每一輪的減少而減小。
安裝
library(devtools)
install_github("immunogenomics/harmony")
流程
我們以Seurat v3為例将谊，使用harmony進(jìn)行數(shù)據(jù)整合：
library(Seurat)
library(cowplot)
library(harmony)
首先梯浪，下載稀疏矩陣示例(https://www.dropbox.com/s/t06tptwbyn7arb6/pbmc_stim.RData?dl=1)并將其移動(dòng)到文件夾下（例如data/）。
load('data/pbmc_stim.RData') #加載矩陣數(shù)據(jù)
Initialize Seurat Object
在運(yùn)行Harmony之前瓢娜，創(chuàng)建一個(gè)Seurat對(duì)象并按照標(biāo)準(zhǔn)PCA（用了這么多年的PCA可視化竟然是錯(cuò)的９衣濉！眠砾！）進(jìn)行分析虏劲。
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語(yǔ)言中%>%的含義是什么呢，管道函數(shù)啦褒颈，就是把左件的值發(fā)送給右件的表達(dá)式柒巫，并作為右件表達(dá)式函數(shù)的第一個(gè)參數(shù)。
pbmc@meta.data$stim <- c(rep("STIM", ncol(stim.sparse)), rep("CTRL", ncol(ctrl.sparse)))#賦值條件變量
未經(jīng)校正的PC中的數(shù)據(jù)集之間存在明顯差異：
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
運(yùn)行Harmony的最簡(jiǎn)單方法是傳遞Seurat對(duì)象并指定要集成的變量谷丸。RunHarmony返回Seurat對(duì)象堡掏，并使用更正后的Harmony坐標(biāo)。讓我們將plot_convergence設(shè)置為TRUE刨疼，這樣我們就可以確保Harmony目標(biāo)函數(shù)在每一輪中都變得更好泉唁。
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
要直接訪問(wèn)新的Harmony embeddings，請(qǐng)使用Embeddings命令揩慕。
harmony_embeddings <- Embeddings(pbmc, 'harmony')
harmony_embeddings[1:5, 1:5]
讓我們查看確認(rèn)數(shù)據(jù)集在Harmony運(yùn)行之后的前兩個(gè)維度中得到很好的整合亭畜。
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
許多下游分析是在低維嵌入而不是基因表達(dá)上進(jìn)行的。要使用校正后的Harmony embeddings而不是PC（還在用PCA降維迎卤？快學(xué)學(xué)大牛最愛(ài)的t-SNE算法吧, 附Python/R代碼）拴鸵，請(qǐng)?jiān)O(shè)置reduction ='harmony'。例如，讓我們使用Harmony降維后的數(shù)據(jù)執(zhí)行UMAP和Nearest Neighbor分析劲藐。
pbmc <- pbmc %>%
RunUMAP(reduction = "harmony", dims = 1:20) %>%
FindNeighbors(reduction = "harmony", dims = 1:20) %>%
FindClusters(resolution = 0.5) %>%
identity()
在UMAP embedding中八堡，我們可以看到更復(fù)雜的結(jié)構(gòu)。由于我們使用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')
在這種充分混合的嵌入中浴韭，我們可以開(kāi)始使用聚類分析來(lái)識(shí)別細(xì)胞類型（Celaref | 單細(xì)胞測(cè)序細(xì)胞類型注釋工具）。
options(repr.plot.height = 4, repr.plot.width = 6)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = .1)
快來(lái)試一試：https://github.com/immunogenomics/harmony
