重生之我在劍橋大學(xué)學(xué)習(xí)單細胞RNA-seq分析——8. scRNA-seq數(shù)據(jù)整合(3)

8.6 Harmony, 3’ vs 5’ 10k PBMC
使用harmony比任何其他方法都要快得多,并且在最近的標準測試中發(fā)現(xiàn)其表現(xiàn)相當好,還可以方便地與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

檢查生成的降維結(jié)果:

> 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的結(jié)果看起來比Seurat的要差一點:

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

最后來看一下cluster內(nèi)容:

> plot_integrated_clusters(pbmc_harmony)
> rm(pbmc_harmony)

Cluster及其內(nèi)容與我們在Seurat整合后獲得的結(jié)果非常相似。為了進行更詳細的分析,需要進行細胞類型注釋。
8.7 LIGER, 3’ vs 5’ 10k PBMC
與其他方法類似,我們創(chuàng)建一個統(tǒng)一的對象并對其進行標準化/HVG/scale。LIGER在scale時不會將數(shù)據(jù)中心化,因此ScaleData中指定了do.center選項。最后兩個函數(shù)是使用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之后執(zhí)行Louvain聚類(FindNeighborsFindClusters),以將結(jié)果與之前的整合方法進行比較。這里使用相同的參數(shù)(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整合數(shù)據(jù)進行聚類似乎更加精細:

這些cluster看起來非常不同,并且每個cluster的分布似乎證實了這一點(兩個cluster分別被認為是3'和5'數(shù)據(jù)集所特有的):

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

往期內(nèi)容:
重生之我在劍橋大學(xué)學(xué)習(xí)單細胞RNA-seq分析——8. scRNA-seq數(shù)據(jù)整合(1)
重生之我在劍橋大學(xué)學(xué)習(xí)單細胞RNA-seq分析——8. scRNA-seq數(shù)據(jù)整合(2)

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

友情鏈接更多精彩內(nèi)容