利用DoubletFinder移除單細胞數據中的雙細胞

前言

單細胞RNA定序(Single Cell RNA Sequencing,scRNA-seq)能夠揭示單細胞基因表達的異質性。然而,這種技術面臨著一個常見的問題,那就是存在雙細胞(Doublets)的干擾。雙細胞是指兩個或更多不同細胞的RNA混合體,其存在會導致基因表達的不準確解讀。因此,我們可以使用 DoubletFinder套件來去除這些雙細胞,從而獲得更精確以及更乾淨的scRNA-seq數據。

DoubletFinder的基本原理

DoubletFinder 算法的四個主要步驟:

(1) 生成人工雙細胞:DoubletFinder通過結合隨機選擇的真實單細胞的基因表達概況來創建人工雙細胞。這一步旨在模擬潛在雙細胞的基因表達模式。

(2) 合併真實和人工數據:將人工雙細胞與真實單細胞數據結合,創建一個合併的數據集,然後對這個合併數據集進行預處理,包括正規化和標準化,以準備進行後續分析。

(3) 執行PCA並計算 pANN:對合併數據集應用主成分分析(PCA),降低數據的維度。使用主成分距離矩陣來找到每個細胞的artificial k-nearest neighbors的比例(pANN)。這個指標會根據它與人工雙細胞的相似性,來估計細胞為雙細胞的可能性。

(4) 排名和閾值化 pANN 值:將 pANN 值進行排名和閾值化,pANN 值超過閾值的細胞被認為是潛在的雙細胞,而低於閾值的細胞則被認為是單細胞。

如何下載DoubletFinder套件?

# 下載DoubletFinder套件
remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')

# 檢查版本
packageVersion("DoubletFinder")

本次教學使用的是DoubletFinder v2.0.3版本,沒有remotes套件的朋友可以使用install.packages("remotes")先行下載套件。

範例數據

在這次範例中,我們分析來自10X Genomics的周邊血(PBMC)單細胞數據集,這是一個包含11,439顆單細胞的數據集,採用Illumina NovaSeq 6000進行测序,原始數據可以在這裡下載Feature / cell matrix (filtered))。如果想要了解詳細的Seurat單細胞分析步驟,可以參考Seurat V4.9.9 – 一個強大的單細胞分析R套件教學。

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

# 加載套件
library(Seurat)
library(ggplot2)
library(tidyverse)
library(DoubletFinder)

# 載入PBMC數據集,需用電腦中的正確路徑,須將"\"更改為"/"
pbmc.data <- Read10X(data.dir = "C:/Users/USER/Desktop/filtered_feature_bc_matrix")

# 創建Seurat對象
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc", min.cells = 3, min.features = 100)
pbmc

在這裡使用CreateSeuratObject()函數導入數據,其中counts參數是輸入的數據,project參數是用於標識Seurat對象的名稱,min.cells和min.features參數是為了過濾掉一些低質量細胞和基因。接著我們便可以繼續進一步分析,這裡使用的分析方式是標準分析法。

2. 過濾掉低品質的細胞

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

# 過濾低品質細胞
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# 使用小提琴圖將QC指標視覺化
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)

# 過濾細胞
pbmc <- subset(pbmc, subset = nCount_RNA > 800 & nFeature_RNA > 500 & percent.mt < 10)
pbmc

執行完篩選條件後,用於後續分析的細胞還剩下9,983顆細胞。

3. 進行標準單細胞分析流程

在進行DoubletFinder 分析之前,先要進行標準單細胞分析流程:

# 進行標準分析流程
pbmc <- NormalizeData(object = pbmc)
pbmc <- FindVariableFeatures(object = pbmc)
pbmc <- ScaleData(object = pbmc)
pbmc <- RunPCA(object = pbmc)
ElbowPlot(pbmc)
pbmc <- FindNeighbors(object = pbmc, dims = 1:20)
pbmc <- FindClusters(object = pbmc)
pbmc <- RunUMAP(object = pbmc, dims = 1:20)
DimPlot(pbmc, reduction = "umap", label = T)
圖例為UMAP結果圖

4. 辨識pK值

這段參數是針對處理過的 scRNA-seq 數據 pbmc,進行 pK(PC neighborhood size) 的識別,並使用參數搜索來估算最佳的 pK 值。解釋如下:

  1. paramSweep_v3(pbmc, PCs = 1:20, sct = FALSE): 此函數用於進行參數搜索,尋找在指定的主成分數量範圍內(1到20個主成分),使用經過預處理的 scRNA-seq 數據pbmc進行 PCA 的不同參數組合。參數 sct = FALSE 表示不使用 SCTransform來進行分析
  2. summarizeSweep(sweep.res.list_pbmc, GT = FALSE): 這個函數用於總結參數搜索的結果,返回每個主成分數量下的不同參數組合的統計信息。參數 GT = FALSE 表示沒有基準真實(ground-truth)的數據用於比較,即沒有已知的雙細胞信息
  3. find.pK(sweep.stats_pbmc): 這個函數用於根據參數搜索的統計結果來估算最佳的 pK 值。pK 是 DoubletFinder 中用於計算每個細胞的 pANN的一個重要參數。find.pK()函數通過估算不同 pK 值對應的 BCmvn(雙細胞檢測性能指標)值,來找到在不同主成分數量下的最優 pK 值。
# pK值辨識 (no ground-truth) 
sweep.res.list_pbmc <- paramSweep_v3(pbmc, PCs = 1:20, sct = FALSE)
sweep.stats_pbmc <- summarizeSweep(sweep.res.list_pbmc, GT = FALSE)
bcmvn_pbmc <- find.pK(sweep.stats_pbmc)

# 選擇具有最大BCmvn(雙細胞檢測性能指標)值的pK
ggplot(bcmvn_pbmc, aes(pK, BCmetric, group = 1)) + geom_point() + geom_line()
pK <- bcmvn_pbmc %>% filter(BCmetric == max(BCmetric)) %>% select(pK) 
pK <- as.numeric(as.character(pK[[1]]))
head(pK)

從上圖我們可以看到,具有最大BCmvn(雙細胞檢測性能指標)值的pK落在了0.14的數值,因此我們接下來會使用pK=0.14的參數進行下一步的分析。

5. 估算同型雙細胞(Homotypic Doublet)的比例

# 估算同型雙細胞的比例
annotations <- pbmc@meta.data$seurat_clusters
homotypic.prop <- modelHomotypic(annotations)  
   
# 假設雙細胞形成率為7.5%,這個數值可以從10X genomics的protocol中找到
nExp_poi <- round(0.075*nrow(pbmc@meta.data))  
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))
head(nExp_poi.adj)

套件估算出來的同型雙細胞(Homotypic Doublet)為680顆,接下來就能進行雙細胞鑑定分析了。

6. 鑑定雙細胞

# 運行DoubletFinder 
pbmc <- doubletFinder_v3(pbmc, PCs = 1:20, pN = 0.25, pK = pK, nExp = nExp_poi.adj, reuse.pANN = FALSE, sct = FALSE)
head(pbmc)

# 可視化雙細胞
DimPlot(pbmc, reduction = 'umap', group.by = "DF.classifications_0.25_0.14_680")

# 檢視單細胞和雙細胞的結果
table(pbmc@meta.data$DF.classifications_0.25_0.14_680)

head()函數查看pbmc資料中前幾行的結果,選取DF.classifications_0.25_0.14_680結果來進行可視化。

7. 移除雙細胞

# 移除雙細胞
pbmc_NoDoublets <- subset(pbmc, subset = DF.classifications_0.25_0.14_680 == "Singlet")
DimPlot(pbmc_NoDoublets, reduction = "umap", group.by="DF.classifications_0.25_0.14_680")
pbmc_NoDoublets

移除雙細胞後,就能繼續往下做後續的單細胞數據分析了!

8. 使用DoubletFinder套件的注意事項

  1. 不要將 DoubletFinder用於多個不同樣本的整合式 scRNA-seq 數據。在單個樣本上運行DoubletFinder套件,然後再將所有數據整合在一起才是正確的分析流程。
  2. 確保輸入數據中不包含低質量的細胞群,可以根據以下方法進行數據的過濾和篩選:
    • 使用標準分析流程對數據進行預處理。
    • 過濾具有以下特徵的細胞群:(A) RNA UMIs 很低的細胞,(B) 高比例的線粒體讀數,和/或 (C) 不具有信息價值的標記基因。

外部連結

參考文獻

1. Christopher S McGinnis, Lyndsay M Murrow, Zev J Gartner. DoubletFinder: Doublet Detection in Single-Cell RNA Sequencing Data Using Artificial Nearest Neighbors. Cell Syst. 2019 Apr 24;8(4):329-337.e4.

2. DoubletFinder – GitHub教學官網

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

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

文章: 46

發佈留言

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

error: Alert: Content selection is disabled!!