相信各位做單細(xì)胞的童鞋們囊嘉，都遇到過去除多細(xì)胞的情況，而Cellranger默認(rèn)的是前1%（3000cell），對于低質(zhì)量的細(xì)胞有了算法來確認(rèn)是否保留柱宦，但是對于多細(xì)胞而言，Cellranger和Seurat并沒有相應(yīng)的算法進(jìn)行去除播瞳，今天給大家介紹一個去除多細(xì)胞的方法-DoubletFinder掸刊，一個R包，現(xiàn)在我們開回顧一下這個包的用法

DoubletFinder概述

DoubletFinder工作流程是從處理過的Seurat對象開始赢乓。

seu_kidney <- CreateSeuratObject(kidney.data)
seu_kidney <- NormalizeData(seu_kidney)
seu_kidney <- ScaleData(seu_kidney)
seu_kidney <- FindVariableFeatures(seu_kidney, selection.method = "vst", nfeatures = 2000)
seu_kidney <- RunPCA(seu_kidney)
seu_kidney <- RunUMAP(seu_kidney, dims = 1:10)

## Pre-process Seurat object (sctransform) -----------------------------------------------------------------------------------
seu_kidney <- CreateSeuratObject(kidney.data)
seu_kidney <- SCTransform(seu_kidney)
seu_kidney <- RunPCA(seu_kidney)
seu_kidney <- RunUMAP(seu_kidney, dims = 1:10)

至于更多的Seurat的細(xì)節(jié)處理忧侧，這里不再詳述，大家可以根據(jù)自己的需求編寫Seurat代碼牌芋，我們這里從構(gòu)建好Seurat對象開始分析DoubletFinder蚓炬。
然后，通過對通過隨機(jī)抽樣與替換選擇的細(xì)胞對的基因表達(dá)譜求平均躺屁，從原始UMI計數(shù)矩陣中生成Artificial doublets肯夏。然后生成足夠的Artificial doublets，占合成數(shù)據(jù)的25％（pN = 0.25）犀暑。接下來驯击，使用原始數(shù)據(jù)分析工作流程中使用的相同scale，縮放和可變基因定義參數(shù)耐亏，對真實數(shù)據(jù)和人工數(shù)據(jù)進(jìn)行合并和預(yù)處理余耽。值得注意的是，為了保留singlets and doublets之間的差異苹熏，在合并的人工數(shù)據(jù)集預(yù)處理期間不執(zhí)行nUMI回歸碟贾。使用在原始數(shù)據(jù)預(yù)處理期間選擇的具有統(tǒng)計意義的相同數(shù)量的PC，然后使用“ fields” R包中的“ rdist”函數(shù)將PC單元嵌入轉(zhuǎn)換為歐幾里得距離矩陣（Nychka等人轨域，2017）袱耽。
然后，從此距離矩陣中定義每個單元格的最近鄰居干发，然后通過將每個實際單元格的人工鄰居數(shù)除以鄰域大兄炀蕖（pK），為每個實際單元格計算人工最近鄰居的比例枉长。然后冀续，將最終的doublet分類分配給具有最高pANN的單元格，其中將n設(shè)置為具有或不具有同型doublet調(diào)整的期望doublet的總數(shù)必峰。
在這里我們只用單細(xì)胞數(shù)據(jù)來進(jìn)行模擬
DoubletFinder可以分為4個步驟：
一洪唐、 從現(xiàn)有的scRNA-seq數(shù)據(jù)生成 artificial doublets
二、預(yù)處理合并的real-artificial data
三吼蚁、執(zhí)行PCA并使用PC距離矩陣來找到每個單元的人工k最近鄰（pANN）的比例
四凭需、根據(jù) expected number of doublets排列順序和閾值pANN值
DoubletFinder采用以下參數(shù)：
（1）seu?這是一個經(jīng)過充分處理的Seurat對象（ after NormalizeData, FindVariableGenes, ScaleData, RunPCA, and RunTSNE have all been run）
（2）PC?具有統(tǒng)計意義的主成分?jǐn)?shù)，指定為范圍（e.g., PCs = 1:10）
（3）pN?定義生成的人工雙峰的數(shù)量，表示為合并的真實人工數(shù)據(jù)的一部分粒蜈。 基于DoubletFinder性能在很大程度上是pN是不變的顺献，默認(rèn)設(shè)置為25％
（4）pK?定義用于計算pANN的PC鄰域大小，表示為合并的真實real-artificial數(shù)據(jù)的一部分枯怖。 未設(shè)置默認(rèn)值注整，因為應(yīng)該為每個scRNA-seq數(shù)據(jù)集調(diào)整pK。 最佳pK值應(yīng)使用以下描述的策略進(jìn)行估算度硝。
（5）nExp?這定義了用于進(jìn)行最終doublet/singlet預(yù)測的pANN閾值肿轨。 可以從10X / Drop-Seq裝置中的細(xì)胞裝載密度來最好地估計該值，并根據(jù)homotypic doublets的估計比例進(jìn)行調(diào)整塘淑。

'Best-Practices' for scRNA-seq data generated without sample multiplexing

Input scRNA-seq Data

不要將DoubletFinder應(yīng)用于代表多個不同樣本的匯總scRNA-seq數(shù)據(jù)。 例如蚂斤，如果您對表示跨10Xlane測序的WT和突變細(xì)胞系的匯總數(shù)據(jù)運(yùn)行DoubletFinder存捺，則將從WT和突變細(xì)胞生成artificial doublet，它們在您的數(shù)據(jù)中不存在曙蒸。 這些artificial doublets 會歪曲結(jié)果捌治。 值得注意的是，可以對通過在多個10Xlane上分割單個樣本而生成的數(shù)據(jù)運(yùn)行DoubletFinder纽窟。
確保清除輸入數(shù)據(jù)中的低質(zhì)量cluster

pK Selection

通過對模擬的scRNAseq數(shù)據(jù)進(jìn)行pN-pK參數(shù)掃描的ROC分析肖油，其中（I）細(xì)胞狀態(tài)數(shù)目可變，并且（II）轉(zhuǎn)錄異質(zhì)性的變量大小表明（I）最佳pK值選擇取決于細(xì)胞狀態(tài)總數(shù)臂港，并且（II ）DoubletFinder性能在應(yīng)用于轉(zhuǎn)錄同源數(shù)據(jù)時會受到影響森枪。
使用BCMVN最大化優(yōu)化pK
雙峰系數(shù)（BC）測量數(shù)據(jù)分布中與單峰的偏差.對于DoubletFinder參數(shù)擬合，我們認(rèn)為產(chǎn)生非單峰pANN分布的參數(shù)集將最佳地將單重態(tài)與雙態(tài)分離审孽，結(jié)果將表現(xiàn)最佳县袱。 Thus scRNA-seq datasets,我們測試了在pN-pK參數(shù)掃描期間生成的每個pANN分布，以發(fā)現(xiàn)具有較高BC值的pANN分布佑力。 具體來說式散，我們按照modes R包中的“ bimodality_coefficient”函數(shù)中實現(xiàn)的方式計算了BC，其形式為：

$9W7KS9~PDTGZ8ZUG}J3E9D.png

2.png

sweep.res.list_kidney <- paramSweep_v3(pbmc, PCs = 1:20, sct = FALSE)
#head(sweep.res.list_kidney)
sweep.stats_kidney <- summarizeSweep(sweep.res.list_kidney, GT = FALSE)
#head(sweep.res.list_kidney)
bcmvn_kidney <- find.pK(sweep.stats_kidney)

根據(jù)上面的描述什往，最佳的pk值應(yīng)該是

mpK<-as.numeric(as.vector(bcmvn_kidney$pK[which.max(bcmvn_kidney$BCmetric)]))

annotations <- pbmc@meta.data$seurat_clusters
## Homotypic Doublet Proportion Estimate -------------------------------------------------------------------------------------
homotypic.prop <- modelHomotypic(annotations)           ## ex: annotations <- seu_kidney@meta.data$ClusteringResults
nExp_poi <- round(0.075*length(seu_kidney@cell.names))  ## Assuming 7.5% doublet formation rate - tailor for your dataset
nExp_poi.adj <- round(nExp_poi*(1-homotypic.prop))

最后，我們運(yùn)用

 pbmc <- doubletFinder_v3(pbmc, PCs = 1:20, pN = 0.25, pK = mpK, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)

就大功告成了~~~~~