單細胞數據分析:批次效應與Seurat Integration分析詳解

前言

在經過了單細胞定序簡介以及Seurat套件的介紹後,相信大家都已經對單細胞有了初步的認識。但是在分析單細胞數據時,樣本往往並非只有單一樣本,而是數個,甚至十數個。如果實驗室再有錢點,或是使用資料庫的單細胞數據來進行分析的話,那樣本數會突破數十個也並不奇怪。因此,在分析多個單細胞數據時,我們首要面對的問題就是數據整合(integration)問題。

甚麼是批次效應(Batch Effect)

在單細胞數據分析中,批次效應 (batch effect)是指不同批次,或不同實驗條件下所產生的技術差異,這些差異可能與樣本處理、不同分析流程、不同定序過程或平台、儀器差異等相關。批次效應會對數據的分析和解釋產生顯著的影響,但這些差異卻不是由所研究的生物學變異所引起的。換句話說,批次效應可能掩蓋或乾擾了研究人員所關註的生物學信號,從而導致對樣本間的差異產生錯誤的解讀。

上圖的結果是根據其來源樣本著色,在進行批次效應校正之前,批次之間的分離明顯可見,但在校正之後則不那麼明顯,詳細內容請參考文獻[1]。

批次效應的矯正方法

單細胞數據分析已經日漸成熟,現如今已經有許多方法可以矯正單細胞數據的批次效應,例如: MNN Correct,fastMNN,MultiCCA Seurat 2,Seurat 3,MMD-ResNet,Harmony,Scanorama,BBKNN,scGen,ComBat,LIGER,limma等等,這裡就不一一介紹所有矯正方法,有興趣的朋友可以參考文獻[2]。

Seurat Integration

Seurat Integration [3] 是Seurat團隊在2019年6月發表在Cell期刊上的方法,主要旨在解決”樣本整合”和”數據轉移”的問題。他們引入了錨 (Anchor) 的概念,通過在共用空間中尋找共有的細胞,實現跨不同模態之間的信息整合和標簽轉移。這個方法不僅結合了傳統的CCA (canonical correlation analysis) 演算法,還引入了MNN (Mutual Nearest Neighbor) 的概念。

簡略來說,他們的方法使用了兩個或多個具有不同模態的單細胞數據集,通過計算每個數據集中細胞之間的相似性,找到兩個數據集之間的共有細胞,稱為”Anchor”細胞。這些”Anchor”細胞在不同數據集之間具有高度相似的表達模式,接下來進一步使用CCA演算法來構建共用的低維空間,將不同數據集中的細胞投影到該共用空間中。在這個共用空間中,細胞的表達模式將得到整合,同時保留了不同數據集之間的差異。

為了更好地保持數據的一致性和減少批次效應,Seurat團隊還引入了MNN的概念。MNN是指在不同數據集中尋找相互最近鄰的細胞對,並將它們作為額外的對齊錨點。通過使用MNN,可以更準確地對齊細胞,減少數據整合過程中的技術差異和偏差。今天本篇主要是以Seurat Integration來為大家示範如何進行單細胞數據的整合

範例數據

在這次的範例中,我們將分析來自10X Genomics的周邊血(PBMC)的兩組單細胞數據集,分別包含了2700和5025顆單細胞的數據集,原始數據可以在我的雲端下載,如果有疑慮可以到10X Genomics數據庫自行下載。

數據集下載位置:

3k PBMCs from a Healthy Donor

5k Peripheral Blood Mononuclear Cells (PBMCs) from a Healthy Donor (v3 chemistry)

1. 加載套件和載入PBMC數據集

將PBMC數據解壓縮,加載所需套件進入R環境,然後使用CreateSeuratObject()函數導入數據,其中counts參數是輸入的數據,project參數是用於標識Seurat對象的名稱,min.cells和min.features參數是為了過濾掉一些低質量細胞和基因。如果對Seurat的分析步驟不甚了解的朋友,可以參考我之前所發布的文章。

# 加載套件
library(Seurat)
library(dplyr)
library(patchwork)

# 載入PBMC數據集,需用電腦中的正確路徑,須將"\"更改為"/"
pbmc3K.data <- Read10X(data.dir = "C:/Users/Administrator/Desktop/Seurat integration/pbmc3k")
pbmc5K.data <- Read10X(data.dir = "C:/Users/Administrator/Desktop/Seurat integration/pbmc5k")
# 創建Seurat對象
pbmc3K <- CreateSeuratObject(counts = pbmc3K.data, project = "pbmc3k", min.cells = 3, min.features = 10)
pbmc5K <- CreateSeuratObject(counts = pbmc5K.data, project = "pbmc5K", min.cells = 3, min.features = 10)
# 檢查數據
pbmc3K
pbmc5K

2. 過濾掉低品質的細胞

低質量或者瀕死的細胞通常具有廣泛的粒線體污染,因此我們可以使用PercentageFeatureSet()函數來計算粒線體的指標,並進一步將異常的細胞移除,以免影響後續資料分析。至於移除參數的cut off值,則取決於分析人員

# 過濾掉低品質的細胞 (pbmc3K)
pbmc3K[["percent.mt"]] <- PercentageFeatureSet(pbmc3K, pattern = "^MT-")
# 使用小提琴圖將QC指標視覺化
VlnPlot(pbmc3K, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# 過濾細胞
pbmc3K <- subset(pbmc3K, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 20)

# 過濾掉低品質的細胞 (pbmc5K)
pbmc5K[["percent.mt"]] <- PercentageFeatureSet(pbmc5K, pattern = "^MT-")
# 使用小提琴圖將QC指標視覺化
VlnPlot(pbmc5K, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
# 過濾細胞
pbmc5K <- subset(pbmc5K, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 20)

# 檢查數據
pbmc3K
pbmc5K

從結果中我們可以看出,當我們將percent.mt()的cut off值統一設定在20時,過濾後分別可以得到pbmc3K有2693顆細胞,而pbmc5K則剩下3038顆細胞。從上述兩個單細胞數據例子來看,pbmc5K的數據集可能有較多的低質量細胞,為了配合pbmc5K數據集,我也將pbmc3K數據集的percent.mt()參數調整至與pbmc5K數據一致。這種調整參數的做法,主要是為了保持不同數據集之間的一致性,讓所有參數標準化,這樣就能減少數據整合過程中的偏差。再重申一次,參數cut off值的設定,完全取決於分析人員對自身數據的考量以及判斷

3. 標準化以及選擇高度變異的特徵

# 設定合併檔案
pbmc.merge <- merge (pbmc3K, y = c(pbmc5K), add.cell.ids = c("pbmc3K", "pbmc5K"), project = "groups")
table(pbmc.merge$orig.ident)

# 設定數據分割列表,進行標準化以及選擇高度變異的特徵的分析
pbmc.list<- SplitObject(pbmc.merge, split.by = "orig.ident")
pbmc.list
pbmc.list <- lapply(X = pbmc.list, FUN = function(x) {
  x <- NormalizeData(x,normalization.method = "LogNormalize", scale.factor = 10000)
  x <- FindVariableFeatures(x, selection.method = "vst", nfeatures = 2000)
})

在合併數據部分,將pbmc3K和pbmc5K兩個數據集合並為一個新的數據集pbmc.merge。add.cell.ids()參數用於為每個數據集添加唯一的細胞標識,project()參數則用於指定項目名稱。table(pbmc.merge$orig.ident)用於顯示合並後數據集中每個細胞的原始標識符號,並進一步確認合併檔案的正確性。

在分割數據部分,根據原始標識符號(orig.ident)將合並後的數據集pbmc.merge進行分割,生成一個數據分割列表pbmc.list。每個分割後的數據子集將按照原始標識符號進行命名。

在標準化以及選擇高度變異特徵的部分,使用lapply()函數對pbmc.list中的每個數據子集進行操作。首先,使用NormalizeData()函數對數據進行標準化處理,採用對數標準化方法(normalization.method = “LogNormalize”),並設置比例因數為10000。接著,我們使用FindVariableFeatures()函數選擇高度變異的特徵,這裡採用Variance Stabilizing Transformation (VST)方法,選擇2000個高變異特徵進行後續分析。

4. 進行數據整合(Seurat integration)

# 選擇跨數據集重複變化的特徵進行整合
features <- SelectIntegrationFeatures(object.list = pbmc.list)
# 整合分析
pbmc.anchors <- FindIntegrationAnchors(object.list = pbmc.list, dims = 1:20, anchor.features = features)
# 創建一個“integrated”數據分析
pbmc.combined <- IntegrateData(anchorset = pbmc.anchors, dims = 1:20)

使用FindIntegrationAnchors()函數來找到數據集之間的錨點。這個函數接受數據集列表pbmc.list作為輸入,以及使用維度(這里以1到20個dimension為例子)和選定的特徵作為參數。

使用IntegrateData()函數將數據集進行整合,以錨點信息pbmc.anchors和維度範圍(這里以1到20個dimension為例子)作為參數,並將數據集整合到一個共同的整合空間中,生成一個名為pbmc.combined的整合數據分析對象。

5. 進行整合分析

接下來的分析與Seurat的基本流程一樣,進行數據縮放(scaling data)、PCA降維分析、確定數據維度以及進行細胞聚集分析(cell clustering)。

# 進行整合數據分析
# 設定使用先前整合的integrated assay進行分析。
# 原始未修改的數據仍然存在於RNA assay中,可以用DefaultAssay(pbmc.combined) <- "RNA"進行轉換
DefaultAssay(pbmc.combined) <- "integrated"
pbmc.combined <- ScaleData(pbmc.combined, verbose = TRUE)
pbmc.combined <- RunPCA(pbmc.combined, npcs = 50, verbose = TRUE)
#Dertermine the dimensionality of the dataset
ElbowPlot(pbmc.combined, ndims= 50)
# Cluster the cells (1:20)
pbmc.combined <- FindNeighbors(pbmc.combined, reduction = "pca", dims = 1:20)
pbmc.combined <- FindClusters(pbmc.combined, resolution = 0.5)
table(pbmc.combined@active.ident)
table(Idents(pbmc.combined), pbmc.combined$orig.ident)

6. 非線性降維分析法(UMAP/tSNE)

最後一步便是進行可視化非線性降維分析了,圖示以UMAP圖為例子。

# UMAP (1:20)
pbmc.combined <- RunUMAP(pbmc.combined, reduction = "pca", dims = 1:20)
DimPlot(pbmc.combined, reduction = "umap", group.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "umap")
DimPlot(pbmc.combined, reduction = "umap", split.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "umap",label = TRUE, split.by = "orig.ident")
# tSNE (1:20)
pbmc.combined <- RunTSNE(pbmc.combined, reduction = "pca", dims = 1:20)
DimPlot(pbmc.combined, reduction = "tsne", group.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "tsne")
DimPlot(pbmc.combined, reduction = "tsne", split.by = "orig.ident")
DimPlot(pbmc.combined, reduction = "tsne",label = TRUE, split.by = "orig.ident")

得到UMAP和tSNE圖後,接下來就要花時間將所有的細胞群進行差異表達基因分析(DEGs, Differentially expression genes analysis)和鑑定每個群集的細胞類型,鑑定細胞類型的方法可以參考SingleR (2.2.0)以及細胞特徵數據庫的文章,兩篇文章都有很詳細的說明,歡迎大家前往閱讀。

結語

在本篇文章中,我們重點討論了多組單細胞數據集的整合方法和步驟。分析時的各項函數數值的設定,都取決於分析人員對自身數據的考量以及判斷,分析人員應該根據數據的特點、實驗設計和研究目的來選擇合適的函數和參數設置,避免盲目跟從教學參數。

外部連結

以下是我在youtube搜尋到的教學視頻,我認為她教得很好,有興趣的朋友可以點進去看看。

參考文獻

1. Malte D Luecken and Fabian J Theis, Current best practices in single‐cell RNA‐seq analysis: a tutorial, Mol Syst Biol. 2019 Jun; 15(6): e8746.

2.  Hoa Thi Nhu Tran, Kok Siong Ang, Marion Chevrier, Xiaomeng Zhang et al., A benchmark of batch-effect correction methods for single-cell RNA sequencing data, Genome Biology volume 21, Article number: 12 (2020).

3. Tim Stuart, Andrew Butler, Paul Hoffman, Christoph Hafemeister et al., Comprehensive integration of single-cell data, Cell. 2019 Jun 13;177(7):1888-1902.e21.

4. Seurat 官網: https://satijalab.org/seurat/articles/integration_introduction.html

您的分享,我非常感激!!!
MillionQuesn
MillionQuesn

一個旅居台灣的異鄉人,分享突然靈光一現的亮點。

文章: 46

發佈留言

發佈留言必須填寫的電子郵件地址不會公開。 必填欄位標示為 *

error: Alert: Content selection is disabled!!