使用背景
我們一般拿到單細(xì)胞矩陣需要做的一個(gè)預(yù)處理步驟是去雙胞鲁猩，而常用的一個(gè)R包是DoubletFinder，詳細(xì)教程可查閱其GitHub主頁(yè)谤职，本文主要展示如何使用這個(gè)包导街，以及批量化處理的腳本示例。因?yàn)閷?shí)際分析情況中會(huì)有多個(gè)矩陣需要單獨(dú)去除雙胞，樣本多起來就比較麻煩便脊，因此封裝成函數(shù)是一個(gè)辦法。這篇文章會(huì)詳細(xì)展示使用DoubletFinder去除雙胞的流程光戈。
雙胞形成的原因以及為什么要去雙胞哪痰？
雙胞是指由于實(shí)驗(yàn)原因，在單細(xì)胞微流控等流程中某個(gè)液滴中含有2個(gè)及以上的細(xì)胞久妆，而在同一個(gè)液滴中的細(xì)胞在后續(xù)分析中帶有相同的Cell Barcode卿城，從而會(huì)被認(rèn)為是一個(gè)細(xì)胞的偽細(xì)胞拣凹，因此，實(shí)際上雙胞不僅是只兩個(gè)細(xì)胞的混合情況，還可能是多個(gè)細(xì)胞混淆為一個(gè)細(xì)胞家厌。
雙胞的存在谍椅，對(duì)研究結(jié)論會(huì)產(chǎn)生極大的影響弧械。這類偽細(xì)胞的主要特點(diǎn)是檢測(cè)的UMI數(shù)和基因數(shù)往往比正常細(xì)胞要多一倍及以上粗仓，另外可能會(huì)帶有不同細(xì)胞類型的經(jīng)典marker基因，這一方面會(huì)給細(xì)胞類型鑒定待了阻礙奸绷，另一方面會(huì)使得某些結(jié)果不顯著或者呈現(xiàn)與正常趨勢(shì)相反的結(jié)果边灭，因此去除雙胞的步驟相當(dāng)重要。去雙胞的軟件很多健盒，目前公認(rèn)效果比較好的是DoubletFinder包绒瘦。
Part 1
加載所需的R包
library(DoubletFinder)
library(Seurat)
library(dplyr)
Part 2
封裝函數(shù)
預(yù)設(shè)參數(shù)簡(jiǎn)介：name 樣本名稱称簿，input 矩陣路徑，dim.usage 設(shè)置主成分PCA個(gè)數(shù)惰帽，auto 與矩陣路徑結(jié)構(gòu)有關(guān)憨降。（因?yàn)槲覀兊玫降木仃嚶窂揭话阌袃煞N結(jié)構(gòu)，所以我在函數(shù)中加入了auto參數(shù)以兼容多種結(jié)構(gòu)）该酗。
rm_doublet <- function(name=NULL,input=NULL,dim.usage=30,auto="true") {
##獲取矩陣路徑授药，這里傳入的矩陣是10X的格式
inpath <- list.files(path = input,pattern = name,full.names = T)
if (auto=="true") {
inpath <- paste0(inpath,"/","04.Matrix/")#這是另一種矩陣路徑結(jié)構(gòu)
}
##去雙胞需要先聚類，均使用默認(rèn)參數(shù)即可呜魄，這一塊可以用管道符%>%寫得簡(jiǎn)潔些悔叽，因?yàn)槲覒芯筒活~外寫了
EC <- Read10X(inpath,gene.column=1)
EC <- CreateSeuratObject(EC)
EC <- NormalizeData(EC)
EC <- FindVariableFeatures(EC, selection.method = "vst", nfeatures = 3000)
EC <- ScaleData(EC)
EC <- RunPCA(EC)
EC <- RunUMAP(EC, dims = 1:dim.usage)
##DoubletFinder去雙胞的標(biāo)準(zhǔn)流程，封裝成一個(gè)函數(shù)
Find_doublet <- function(data){
sweep.res.list <- paramSweep_v3(data, PCs = 1:dim.usage, sct = FALSE)
sweep.stats <- summarizeSweep(sweep.res.list, GT = FALSE)
bcmvn <- find.pK(sweep.stats)
nExp_poi <- round(0.05*ncol(data))
p<-as.numeric(as.vector(bcmvn[bcmvn$MeanBC==max(bcmvn$MeanBC),]$pK))
data <- doubletFinder_v3(data, PCs = 1:dim.usage, pN = 0.25, pK = p, nExp = nExp_poi, reuse.pANN = FALSE, sct = FALSE)
colnames(data@meta.data)[ncol(data@meta.data)] = "doublet_info"
return(data)
}
##調(diào)用上面寫好的函數(shù)爵嗅，返回的是一個(gè)Seurat對(duì)象娇澎，meta.data信息里會(huì)有雙胞信息，需要自己手動(dòng)刪除
EC<-Find_doublet(EC)
EC<-subset(EC,subset=doublet_info=="Singlet")
EC@meta.data$library = name #順便打上這個(gè)樣本的label
##生成的Seurat對(duì)象有個(gè)問題睹晒，會(huì)在meta.data里多了很多pANN_開頭的列趟庄，需要手動(dòng)刪除
c <- grep("pANN_",colnames(EC@meta.data))
EC@meta.data <- EC@meta.data[,-c]
##輸出此樣本的細(xì)胞數(shù)
print(paste0(name," cells: ", length(EC@meta.data$orig.ident)," is read!"," Time:",format(Sys.time(), "%Y%m%d %X"),sep = " "))
return(EC)
}
Part 3
對(duì)單個(gè)樣本取雙胞。
返回的是一個(gè)Seurat對(duì)象伪很。
obj <- rm_doublet(name,input)
可保存為RDS格式戚啥。
saveRDS(obj,paste0(name,".rds"))
建議寫一個(gè)R腳本傳參批量投遞腳本得到RDS（這部分就不寫了，如果有需求可留言）锉试，之后也可批量讀取RDS猫十，具體可參考下面的函數(shù)：
merge_seob_rds <- function(infile,input,batch=NULL){
filelist <- readLines(infile)
inpath <- input
filelist <- unlist(lapply(filelist,list.files, path=inpath,full.names = T))
seob_list <- lapply(filelist, readRDS)
n <- length(seob_list)
all <- merge(seob_list[[1]], seob_list[c(2:n)])
if (!is.null(batch)) {
all$batch <- batch
}
return(all)
}
obj <- merge_seob_rds(infile="datasets_list",input="rds_inpath")
#infile 是樣本名稱列表txt文件，input是去雙胞后保留的RDS文件所在路徑呆盖，batch是給樣本打上標(biāo)簽拖云。
以上方式適合多樣本（樣本>=20）的情況。
Part 4
如果樣本不多絮短，可以參考以下腳本，建議樣本量<=20使用昨忆。
merge_seob <- function(infile,input,batch=NULL){
filelist <- readLines(infile)
seob_list <- lapply(filelist, rm_doublet, input=input)
n <- length(seob_list)
all <- merge(seob_list[[1]], seob_list[c(2:n)])
if (!is.null(batch)) {
all$batch <- batch
}
return(all)
}
ob1 <- merge_seob(infile="datasets_list",input="rds_inpath")
小結(jié)與補(bǔ)充
以上參數(shù)都是實(shí)踐后確定的丁频，當(dāng)然也可以根據(jù)官方文檔和具體研究課題進(jìn)行調(diào)試。DoubletFinder去雙胞根據(jù)不同矩陣質(zhì)量耗時(shí)長(zhǎng)短差異很大邑贴，所以很多樣本的建議多個(gè)任務(wù)并行跑席里，不然要等很久很久很久（親測(cè)）。有時(shí)候用DoubletFinder也不能完全去掉所有雙胞拢驾，但是雙胞一般會(huì)自動(dòng)聚成一個(gè)亞群奖磁，可以手動(dòng)刪除這種亞群。
