單細胞分析實錄(4): doublet檢測

最近Cell Systems雜誌發表了一篇針對現有幾種檢測單細胞測序doublet的工具的評估文章，系統比較了常見的例如Scrublet、DoubletFinder等工具在檢測準確性、計算效率等方面的優劣，以及比較了使用不一樣方法去除doublet後對下游DE分析、軌跡分析的影響。

現有的檢測方法，基本都會先構造出虛擬doublet，而後將候選droplet與這些虛擬doublet比較，很類似的那些就定義爲doublet。這裏的虛擬doublet是經過隨機組合兩個（類）細胞的表達值獲得的虛擬的doublet，能夠做爲檢測時的參照。

在現有的9種方法中(Scrublet、doubletCells、cxds、bcds、Hybrid、DoubletDetection、DoubletFinder、Solo、DoubletDecon)，文章的結論是DoubletFinder的準確率最高。

裏面的方法我用過三種：Scrublet、DoubletFinder和DoubletDecon。前面兩個用法相似，須要提供一個參數表示doublet的佔比。DoubletDecon的原文我看過，算法比較簡單，不須要提供doublet的佔比，減小了用戶額外的輸入，不過也形成了一些麻煩，有時候會報告出不少doublet，多得驚人。實際分析中，我採用「兩步走」的策略：選取了Scrublet和DoubletFinder的共同結果做爲doublet去除掉，此外在後續聚類分亞羣的過程當中，根據一些已知的經典的細胞marker來判斷doublet，好比CD45和EPCAM都高表達的亞羣極有多是doublet。

接下來，我會簡單介紹DoubletDecon、DoubletFinder、Scrublet三個工具的使用。

1. DoubletDecon

該方法中間有一步用到了相似bulk RNA-seq裏面deconvolution的思路來評估每個細胞的表達模式，所以叫DoubletDecon。
這裏用含有500個細胞的真實數據做爲例子。在使用DoubletDecon以前，須要先用seurat對數據進行初聚類，seurat的使用我後面會詳細講，這裏先把標準流程放上來。

library(tidyverse)
library(Seurat)
library(DoubletDecon)

test=read.table("test.count.txt",header = T,row.names = 1)
test.seu <- CreateSeuratObject(counts = test)
#Normalize
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
#FindVariableFeatures
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
#Scale
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
#PCA
test.seu <- RunPCA(test.seu, features = VariableFeatures(test.seu),npcs = 50)
#cluster
test.seu <- FindNeighbors(test.seu, dims = 1:20)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:20)
test.seu <- RunTSNE(test.seu, dims = 1:20)

而後纔是DoubletDecon的代碼

#Improved_Seurat_Pre_Process()
newFiles=Improved_Seurat_Pre_Process(test.seu, num_genes=50, write_files=FALSE)

這一步主要是找差別基因，會返回三個表格，分別表示marker基因的表達矩陣、全部基因的表達矩陣、細胞的seurat_cluster註釋，前兩個文件的第一行第一列有相應的註釋。

而後就是找doublet的主要步驟了

#Main_Doublet_Decon
results=Main_Doublet_Decon(rawDataFile=newFiles$newExpressionFile, 
                           groupsFile=newFiles$newGroupsFile, 
                           filename="tmp", 
                           location="./",
                           species="hsa",
                           rhop=1,
                           num_doubs=80,
                           write=FALSE,
                           heatmap=TRUE,
                           centroids=TRUE,
                           nCores=2)

這裏面有幾個很重要的參數，rhop默認值爲1，用它來調節皮爾森相關係數的閾值（以下圖）。在seurat聚類以後，這個軟件會進一步將很類似的cluster合併，利用的就是最初cluster之間表達的相關性。這個值也很麻煩，前面提到的DoubletDecon會檢測出不少doublet的問題可能就是這個參數設置不對致使的。那這個參數應該如何設置，可能軟件做者也意識到了這個問題，後面又發了一篇Protocol，專門講參數如何選擇，核心思想就是多試幾回。（事兒真多，準備放棄這個軟件了）

接下來把DoubletDecon的檢測結果保存成單獨的文件，方便後面使用

doublet_df <- as.data.frame(results$Final_doublets_groups)
doublet_df$DoubletDecon="Doublet"
singlet_df <- as.data.frame(results$Final_nondoublets_groups)
singlet_df$DoubletDecon="Singlet"
DD_df <- rbind(doublet_df,singlet_df)
DD_df <- DD_df[colnames(test),]
DD_df$CB = colnames(test)
DD_df <- DD_df[,c("CB","DoubletDecon")]
write.table(DD_df, file = "DoubletDecon_result.txt", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)

也能夠將doublet的結果投到tsne圖上看看效果

test.seu@meta.data$CB=rownames(test.seu@meta.data)
test.seu@meta.data=inner_join(test.seu@meta.data,DD_df,by="CB")
rownames(test.seu@meta.data)=test.seu@meta.data$CB
DimPlot(test.seu,reduction = "tsne",pt.size = 2,group.by = "DoubletDecon")

看上去還不戳，符合doublet單獨成羣的預期

2. DoubletFinder

DoubletFinder找doublet的原理也比較簡單，看細胞羣裏面虛擬doublet的佔比，超過某個閾值就認定這一羣的真實細胞是doublet。在運行DoubletFinder以前，須要對細胞進行初聚類，這和上一種方法是同樣的。

library(Seurat)
library(DoubletFinder)
test.seu <- Createtest.seuratObject(counts = test)
test.seu <- NormalizeData(test.seu, normalization.method = "LogNormalize", scale.factor = 10000)
test.seu <- FindVariableFeatures(test.seu, selection.method = "vst", nfeatures = 2000)
test.seu <- ScaleData(test.seu, features = rownames(test.seu))
test.seu <- RunPCA(test.seu, features = VariableFeatures(test.seu),npcs = 50)
test.seu <- FindNeighbors(test.seu, dims = 1:20)
test.seu <- FindClusters(test.seu, resolution = 0.5)
test.seu <- RunUMAP(test.seu, dims = 1:20)
test.seu <- RunTSNE(test.seu, dims = 1:20)

接下來選擇一個重要參數pK，這個參數定義了PC的neighborhood size

sweep.res.list <- paramSweep_v3(test.seu, PCs = 1:10, sct = FALSE)
for(i in 1:length(sweep.res.list)){
  if(length(sweep.res.list[[i]]$pANN[is.nan(sweep.res.list[[i]]$pANN)]) != 0){
    if(i != 1){
      sweep.res.list[[i]] <- sweep.res.list[[i - 1]]
    }else{
      sweep.res.list[[i]] <- sweep.res.list[[i + 1]]
    }
  }
}
sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <- find.pK(sweep.stats)
pk_v <- as.numeric(as.character(bcmvn$pK))
pk_good <- pk_v[bcmvn$BCmetric==max(bcmvn$BCmetric)]

nExp_poi <- round(0.1*length(colnames(test.seu)))

指按期望的doublet數

test.seu <- doubletFinder_v3(test.seu, PCs = 1:10, pN = 0.25, pK = pk_good, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)

這一行是主要代碼，會在test.seu@meta.data數據框的基礎上加上兩列

colnames(test.seu@meta.data)[ncol(test.seu@meta.data)]="DoubletFinder"

第二列換一個列名

DF_df <- test.seu@meta.data[,c("CB","DoubletFinder")]
write.table(DF_df, file = "DoubletFinder_result.txt", quote = FALSE, sep = '\t', row.names = FALSE, col.names = TRUE)
DimPlot(test.seu,reduction = "tsne",pt.size = 2,group.by = "DoubletFinder")

最終的效果以下圖所示：

3. Scrublet

是一個Python包，安裝能夠參考：https://github.com/AllonKleinLab/scrublet

>>> import scrublet as scr
>>> import numpy as np
>>> infile = "test.count.txt"
>>> outfile = "Scrublet_result.txt"

下面的代碼對輸入文件作預處理，包括提取出CB，提取count矩陣並轉置

>>> finallist = []
>>> with open(infile, 'r') as f:
...     header = next(f)
...     cell_barcodes = header.rstrip().split('\t')
...     for line in f:
...             tmpline = line.rstrip().split('\t')[1: ]
...             tmplist = [int(s) for s in tmpline]
...             finallist.append(tmplist)

>>> finalarray = np.array(finallist)
>>> count_matrix = np.transpose(finalarray)

Scrublet檢測doublet主要代碼以下：

>>> scrub = scr.Scrublet(count_matrix, expected_doublet_rate = 0.1)
>>> doublet_scores, predicted_doublets = scrub.scrub_doublets()
>>> predicted_doublets_final = scrub.call_doublets(threshold = 0.2)

第三行換閾值，更新註釋結果，接下來保存結果

>>> with open(outfile, 'w') as f:
...     f.write('\t'.join(['CB', 'Scrublet', 'Scrublet_Score']) + '\n')
...     for i in range(len(doublet_scores)):
...             if predicted_doublets_final[i] == 0:
...                     result = 'Singlet'
...             else:
...                     result = 'Doublet'
...             f.write('\t'.join([cell_barcodes[i], result, str(doublet_scores[i])]) + '\n')

切換到R環境，在tsne上看看效果

SR_df=read.table("Scrublet_result.txt",header = T,sep = "\t",stringsAsFactors = F)
test.seu@meta.data=inner_join(test.seu@meta.data,SR_df,by="CB")
rownames(test.seu@meta.data)=test.seu@meta.data$CB
DimPlot(test.seu,reduction = "tsne",pt.size = 2,group.by = "Scrublet")

---

到這裏就把三個軟件的基本使用講完了，我只使用了一個實際數據演示，結果並不足以反映這幾個軟件誰好誰壞，小夥伴們須要結合本身的數據選擇合適的軟件。開篇提到的文獻可信度仍是挺高的，你們有須要的話能夠認真學一下DoubletFinder這個軟件。
另外，能夠在github上看到這幾個軟件是相互推薦的，在生信圈子還挺少見~

因水平有限，有錯誤的地方，歡迎批評指正！