利用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)
4. 辨識pK值
這段參數是針對處理過的 scRNA-seq 數據 pbmc，進行 pK（PC neighborhood size） 的識別，並使用參數搜索來估算最佳的 pK 值。解釋如下：
paramSweep_v3(: 此函數用於進行參數搜索，尋找在指定的主成分數量範圍內（1到20個主成分），使用經過預處理的 scRNA-seq 數據pbmc, PCs = 1:20, sct = FALSE)pbmc進行 PCA 的不同參數組合。參數sct = FALSE表示不使用 SCTransform來進行分析。summarizeSweep(: 這個函數用於總結參數搜索的結果，返回每個主成分數量下的不同參數組合的統計信息。參數sweep.res.list_pbmc, GT = FALSE)GT = FALSE表示沒有基準真實（ground-truth）的數據用於比較，即沒有已知的雙細胞信息。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套件的注意事項
- 不要將 DoubletFinder用於多個不同樣本的整合式 scRNA-seq 數據。在單個樣本上運行DoubletFinder套件，然後再將所有數據整合在一起才是正確的分析流程。
- 確保輸入數據中不包含低質量的細胞群，可以根據以下方法進行數據的過濾和篩選：
- 使用標準分析流程對數據進行預處理。
- 過濾具有以下特徵的細胞群：(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教學官網
