Harmony integrates spatially resolved transcriptomic with dissociated scRNAseq datasets

#R包导入

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

#数据

• 数据下载pbmc_stim

• 数据导入

load('data/pbmc_stim.RData')

#初始化Seurat 对象

使用Harmony之前，需要构建一个Seurat对象， 进行一个标准的Seurat分析（包括PCA）。

Seurat中使用Harmony的流程与常规流程的区别是：可以所有细胞创建一个Seurat 对象，不需要为每个数据集创建一个Seurat 对象（Seurat list）。

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)

在object的metadata中定义细胞ID信息，变量名为stim.

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

现在还没有使用Harmony矫正主成分分析结果的数据，数据集之间还是有很大的差异。

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)

Before Harmony

#运行Harmony进行数据整合(矫正批次效应)

• 输入：使用Harmony，需要一个Seurat 对象和指定metadata信息中需要整合的变量名。
• 输出：返回一个Seurat对象，以及矫正之后的Harmony 信息
• plot_convergenc参数设置为TRUE，保证Harmony 在运行中每一次迭代都在使矫正越累越好。

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

RunHarmony

##获取Harmony 矫正之后的信息，使用Embeddings()函数

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

harmony embeddings

##查看数据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)

After harmony

#下游分析

下游分析都是基于Harmony矫正之后的PCA降维数据，不是基因表达数据和直接的PCA降维数据。设置reduction = 'harmony'，后续分析就会调用Harmony矫正之后的PCA降维数据。

使用Harmony 矫正之后的数据，UMAP 和 Nearest Neighbor分析。

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

##UMAP 结果

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

Umap

##聚类分析

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

Cluster analysis

#参考：

Seurat V3 with harmony