Harmony integration in Seurat V3 pipeline

最新推荐文章于 2026-06-26 17:11:23 发布

原创最新推荐文章于 2026-06-26 17:11:23 发布 · 772 阅读

2 ·

本内容遵循CC 4.0 BY-SA版权协议

收录于

单细胞pipeline

在Seurat V3的工作流程中，通常将每个样本的矩阵分别读入为独立的Seurat对象进行整合。然而，使用Harmony时，需要先合并所有样本的矩阵再创建一个Seurat对象。确保在元数据中包含分组信息，如'stim'变量。整合前数据的异质性显著，整合过程可以采用Harmony而不是PCA作为降维方法。整合后的效果可以通过下游分析进行评估。

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


load('data/pbmc_stim.RData')

创建seurat对象

在原有的seurat CCA整合中，每个样本的matrix都被单独读入成为一个seurat object，整合是对若干个seurat object进行的。但在Harmony中应将所有样本matrix合并后直接读入成一个seurat object.

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)

应确保分组信息在seurat metadata中，此处用‘stim’变量表示

pbmc@meta.data$stim <- c(rep("STIM", ncol(stim.sparse)), rep("CTRL", ncol(ctrl.sparse)))

可以看出，未整合的数据异质性非常大

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)

整合数据

#注意，此例中为stim control两个样本，故参数选择“stim”,样本较多时注意更改

options(repr.plot.height = 2.5, repr.plot.width = 6)
pbmc <- pbmc %>% 
    RunHarmony("stim", plot_convergence = TRUE)

观察一下整合效果，注意reduction要换成harmony而非PCA

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)

下游分析

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

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

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

标签

#r语言