重生之我在剑桥大学学习单细胞RNA-seq分析——8. scRNA-seq数据整合（3）

8.6 Harmony, 3’ vs 5’ 10k PBMC
使用harmony比任何其他方法都要快得多，并且在最近的标准测试中发现其表现相当好，还可以方便地与Seurat交互。让我们首先合并对象（不进行整合）。

> pbmc_harmony    <- merge(srat_3p,srat_5p)
Warning message:
In CheckDuplicateCellNames(object.list = objects) :
  Some cell names are duplicated across objects provided. Renaming to enforce unique cell names.

进行和以前一样的分析：

> pbmc_harmony <- NormalizeData(pbmc_harmony, verbose = F)
> pbmc_harmony <- FindVariableFeatures(pbmc_harmony, selection.method = "vst", nfeatures = 2000, verbose = F)
> pbmc_harmony <- ScaleData(pbmc_harmony, verbose = F)
> pbmc_harmony <- RunPCA(pbmc_harmony, npcs = 30, verbose = F)
> pbmc_harmony <- RunUMAP(pbmc_harmony, reduction = "pca", dims = 1:30, verbose = F)

再次绘制整合前的UMAP图：

> DimPlot(pbmc_harmony,reduction = "umap") + 
    plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, before integration")

在合并的Seurat对象上运行RunHarmony，使用orig.ident作为批处理：

> pbmc_harmony <- pbmc_harmony %>% RunHarmony("orig.ident", plot_convergence = T)
Harmony 1/10
Harmony 2/10
Harmony 3/10
Harmony 4/10
Harmony 5/10
Harmony converged after 5 iterations
Warning: Invalid name supplied, making object name syntactically valid. New object name is Seurat..ProjectDim.RNA.harmony; see ?make.names for more details on syntax validity

检查生成的降维结果：

> harmony_embeddings <- Embeddings(pbmc_harmony, 'harmony')
> harmony_embeddings[1:5, 1:5]
                      harmony_1 harmony_2 harmony_3  harmony_4 harmony_5
AAACCCACATAACTCG-1_1  -9.206607 -2.351619 -2.374652  -1.897186 -1.011885
AAACCCACATGTAACC-1_1   7.124223 21.600131 -0.292039   1.530283 -5.792142
AAACCCAGTGAGTCAG-1_1 -18.134725  3.405369  5.256459   4.220001  3.961466
AAACGAACAGTCAGTT-1_1 -18.103262 15.279955 12.301681 -18.115094 31.785955
AAACGAACATTCGGGC-1_1  11.097966 -2.330278 -2.723953   1.546468  1.552332

查看PCA图：

> p1 <- DimPlot(object = pbmc_harmony, reduction = "harmony", pt.size = .1, group.by = "orig.ident") + NoLegend()
> p2 <- VlnPlot(object = pbmc_harmony, features = "harmony_1", group.by = "orig.ident", pt.size = .1) + NoLegend()
> plot_grid(p1,p2)

进行UMAP和聚类：

> pbmc_harmony <- pbmc_harmony %>% 
    RunUMAP(reduction = "harmony", dims = 1:30, verbose = F) %>% 
    FindNeighbors(reduction = "harmony", k.param = 10, dims = 1:30) %>% 
    FindClusters() %>% 
    identity()

绘制与上章相似的UMAP图：

> pbmc_harmony <- SetIdent(pbmc_harmony,value = "orig.ident")
> DimPlot(pbmc_harmony,reduction = "umap") + 
    plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, after integration (Harmony)")

> DimPlot(pbmc_harmony, reduction = "umap", group.by = "orig.ident", pt.size = .1, split.by = 'orig.ident') + NoLegend()

harmony的结果看起来比Seurat的要差一点：

> pbmc_harmony <- SetIdent(pbmc_harmony,value = "seurat_clusters")
> DimPlot(pbmc_harmony,label = T) + NoLegend()

最后来看一下cluster内容：

> plot_integrated_clusters(pbmc_harmony)

> rm(pbmc_harmony)

Cluster及其内容与我们在Seurat整合后获得的结果非常相似。为了进行更详细的分析，需要进行细胞类型注释。
8.7 LIGER, 3’ vs 5’ 10k PBMC
与其他方法类似，我们创建一个统一的对象并对其进行标准化/HVG/scale。LIGER在scale时不会将数据中心化，因此ScaleData中指定了do.center选项。最后两个函数是使用orig.ident作为批处理变量运行rliger的包装器。

> pbmc_liger    <- merge(srat_3p,srat_5p)
> pbmc_liger    <- NormalizeData(pbmc_liger)
> pbmc_liger    <- FindVariableFeatures(pbmc_liger)
> pbmc_liger    <- ScaleData(pbmc_liger, split.by = "orig.ident", do.center = F)
> pbmc_liger    <- RunOptimizeALS(pbmc_liger, k = 30, lambda = 5, split.by = "orig.ident")
> pbmc_liger    <- RunQuantileNorm(pbmc_liger, split.by = "orig.ident")

可以选择在RunQuantileNorm之后执行Louvain聚类（FindNeighbors和FindClusters），以将结果与之前的整合方法进行比较。这里使用相同的参数（k=10，Louvain聚类的默认分辨率）。

> pbmc_liger    <- FindNeighbors(pbmc_liger,reduction = "iNMF",k.param = 10,dims = 1:30)
> pbmc_liger    <- FindClusters(pbmc_liger)

和前面一样，我们进行降维并绘图：

> pbmc_liger    <- RunUMAP(pbmc_liger, dims = 1:ncol(pbmc_liger[["iNMF"]]), reduction = "iNMF", verbose = F)
> pbmc_liger    <- SetIdent(pbmc_liger,value = "orig.ident")
> DimPlot(pbmc_liger,reduction = "umap") + plot_annotation(title = "10k 3' PBMC and 10k 5' PBMC cells, after integration (LIGER)")

> DimPlot(pbmc_liger, reduction = "umap", group.by = "orig.ident", pt.size = .1, split.by = 'orig.ident') + NoLegend()

使用LIGER整合数据进行聚类似乎更加精细：

这些cluster看起来非常不同，并且每个cluster的分布似乎证实了这一点（两个cluster分别被认为是3'和5'数据集所特有的）：

> plot_integrated_clusters(pbmc_liger)

> rm(pbmc_liger)
> rm(srat_3p)
> rm(srat_5p)

往期内容：
重生之我在剑桥大学学习单细胞RNA-seq分析——8. scRNA-seq数据整合（1）
重生之我在剑桥大学学习单细胞RNA-seq分析——8. scRNA-seq数据整合（2）