使用背景
我们一般拿到单细胞矩阵需要做的一个预处理步骤是去双胞，而常用的一个R包是DoubletFinder，详细教程可查阅其GitHub主页，本文主要展示如何使用这个包，以及批量化处理的脚本示例。因为实际分析情况中会有多个矩阵需要单独去除双胞，样本多起来就比较麻烦，因此封装成函数是一个办法。这篇文章会详细展示使用DoubletFinder去除双胞的流程。
双胞形成的原因以及为什么要去双胞？
双胞是指由于实验原因，在单细胞微流控等流程中某个液滴中含有2个及以上的细胞，而在同一个液滴中的细胞在后续分析中带有相同的Cell Barcode，从而会被认为是一个细胞的伪细胞，因此，实际上双胞不仅是只两个细胞的混合情况，还可能是多个细胞混淆为一个细胞。
双胞的存在，对研究结论会产生极大的影响。这类伪细胞的主要特点是检测的UMI数和基因数往往比正常细胞要多一倍及以上，另外可能会带有不同细胞类型的经典marker基因，这一方面会给细胞类型鉴定待了阻碍，另一方面会使得某些结果不显著或者呈现与正常趋势相反的结果，因此去除双胞的步骤相当重要。去双胞的软件很多，目前公认效果比较好的是DoubletFinder包。
Part 1
加载所需的R包
library(DoubletFinder)
library(Seurat)
library(dplyr)
Part 2
封装函数
预设参数简介：name 样本名称，input 矩阵路径，dim.usage 设置主成分PCA个数，auto 与矩阵路径结构有关。（因为我们得到的矩阵路径一般有两种结构，所以我在函数中加入了auto参数以兼容多种结构）。
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/")#这是另一种矩阵路径结构
}
##去双胞需要先聚类，均使用默认参数即可，这一块可以用管道符%>%写得简洁些，因为我懒就不额外写了
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去双胞的标准流程，封装成一个函数
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)
}
##调用上面写好的函数，返回的是一个Seurat对象，meta.data信息里会有双胞信息，需要自己手动删除
EC<-Find_doublet(EC)
EC<-subset(EC,subset=doublet_info=="Singlet")
EC@meta.data$library = name #顺便打上这个样本的label
##生成的Seurat对象有个问题，会在meta.data里多了很多pANN_开头的列，需要手动删除
c <- grep("pANN_",colnames(EC@meta.data))
EC@meta.data <- EC@meta.data[,-c]
##输出此样本的细胞数
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
对单个样本取双胞。
返回的是一个Seurat对象。
obj <- rm_doublet(name,input)
可保存为RDS格式。
saveRDS(obj,paste0(name,".rds"))
建议写一个R脚本传参批量投递脚本得到RDS（这部分就不写了，如果有需求可留言），之后也可批量读取RDS，具体可参考下面的函数：
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是给样本打上标签。
以上方式适合多样本（样本>=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")
小结与补充
以上参数都是实践后确定的，当然也可以根据官方文档和具体研究课题进行调试。DoubletFinder去双胞根据不同矩阵质量耗时长短差异很大，所以很多样本的建议多个任务并行跑，不然要等很久很久很久（亲测）。有时候用DoubletFinder也不能完全去掉所有双胞，但是双胞一般会自动聚成一个亚群，可以手动删除这种亚群。
