【原】兩種不同的方法實現(xiàn)harmony的多個單細胞整合

健明 2024-01-08 發(fā)布于廣東

展開全文

本筆記會被收錄于《生信技能樹》公眾號的《單細胞2024》專輯，而且我們從2024開始的教程都是基于Seurat的V5版本啦，之前已經(jīng)演示了如何讀取不同格式的單細胞轉(zhuǎn)錄組數(shù)據(jù)文件，如下所示：

初試Seurat的V5版本

使用Seurat的v5來讀取多個10x的單細胞轉(zhuǎn)錄組矩陣

使用Seurat的v5來讀取多個不是10x標準文件的單細胞項目

但是其它代碼基本上就跟Seurat早期的v4沒有區(qū)別，比如harmony整合多個單細胞樣品。

但實際上Seurat有了自己的創(chuàng)新，官網(wǎng)文檔：

https:///seurat/articles/seurat5_integration

https:///seurat/articles/integration_introduction

里面提到了它內(nèi)置了多種整合多個單細胞樣品的算法，可以 Perform streamlined (one-line) integrative analysis

Anchor-based CCA integration (method=CCAIntegration)

Anchor-based RPCA integration (method=RPCAIntegration)

Harmony (method=HarmonyIntegration)

FastMNN (method= FastMNNIntegration)

scVI (method=scVIIntegration)

如果是一次性讀取了多個10x樣品

這個時候，因為函數(shù)Read10X可以一次性讀取多個合理的路徑，所以我們會把多個樣品就被統(tǒng)一讀取成為了一個稀疏矩陣而不是每個樣品獨立的稀疏矩陣，如下所示；

統(tǒng)一讀取成為了一個稀疏矩陣

詳見：使用Seurat的v5來讀取多個10x的單細胞轉(zhuǎn)錄組矩陣，它就不適合走Seurat的v5的內(nèi)置的多個單細胞樣品的整合算法，所以我們會先split它，代碼如下所示：table(sce.all$orig.ident)

obj = sce.all

obj[["RNA"]] <- split(obj[["RNA"]], f = obj$orig.ident)

效果如下所示，可以看到每個樣品的矩陣這個時候被上面的split函數(shù)拆開了：

split函數(shù)拆開

接下來，如下所示走內(nèi)置的harmony流程，就是IntegrateLayers函數(shù)里面的：obj <- NormalizeData(obj)

obj <- FindVariableFeatures(obj)

obj <- ScaleData(obj)

obj <- RunPCA(obj)

# 運行的harmony需要基于上面的PCA結(jié)果哦

sce.all <- IntegrateLayers(

object = obj, method = HarmonyIntegration,

orig.reduction = "pca", new.reduction = "harmony",

verbose = FALSE

)

sce.all@reductions$harmony

sce.all <- FindNeighbors(sce.all, reduction = "harmony", dims = 1:30)

res = c(0.01, 0.05, 0.1, 0.2, 0.3, 0.5,0.8 )

sce.all <- FindClusters(sce.all, resolution = res )

sce.all <- RunUMAP(sce.all, reduction = "harmony", dims = 1:30)

p1 <- DimPlot( sce.all,label = T)

是很容易看到harmony被成功運行啦。

如果是先自己跑RunHarmony函數(shù)

這個時候又不能用split函數(shù)拆開了的Seurat對象哦，需要最開始那個多個樣品就被統(tǒng)一讀取成為了一個稀疏矩陣的Seurat對象，這個時候可以命名為 input_sce 變量，然后走下面的代碼：print(dim(input_sce))

input_sce <- NormalizeData(input_sce,

normalization.method = "LogNormalize",

scale.factor = 1e4)

input_sce <- FindVariableFeatures(input_sce)

input_sce <- ScaleData(input_sce)

input_sce <- RunPCA(input_sce, features = VariableFeatures(object = input_sce))

seuratObj <- RunHarmony(input_sce, "orig.ident")

names(seuratObj@reductions)

seuratObj <- RunUMAP(seuratObj, dims = 1:15,

reduction = "harmony")

# p = DimPlot(seuratObj,reduction = "umap",label=T )

# ggsave(filename='umap-by-orig.ident-after-harmony',plot = p)

input_sce=seuratObj

input_sce <- FindNeighbors(input_sce, reduction = "harmony",

dims = 1:15)

這個效果跟前面的IntegrateLayers函數(shù)是一回事，如果你這個時候的Seurat對象里面的矩陣是按照樣品拆開的，就需要先合并才能走RunHarmony函數(shù)，它來自于harmony這個r包啦。

通常是不建議走內(nèi)置的harmony流程

其實我就不應該介紹這個IntegrateLayers函數(shù)的，因為它需要split那個矩陣，這樣的話后面的很多分析都會有問題，比如我們跑 cosg 函數(shù)針對那個矩陣去找marker就會遇到報錯：# remotes::install_github('genecell/COSGR')

# genexcell <- Seurat::GetAssayData(object = object[[assay]],slot = slot)

marker_cosg <- cosg(

sce.all.int,

groups='all',

assay='RNA',

slot='data',

mu=1,

n_genes_user=100)

如下所示：Error in `Seurat::GetAssayData()`:

! GetAssayData doesn't work for multiple layers in v5 assay.

Run `rlang::last_trace()` to see where the error occurred.

> rlang::last_trace()