其實在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) ，其算法邏輯如下圖所示：

圖1. Harmony算法概述

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

總結

以上是生活随笔為你收集整理的“harmony”整合不同平台的单细胞数据之旅的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。