单细胞转录组，因为其技术原因，在很多时候，会造成双细胞的情况，一般我们进行基础的质控之后，好点的数据是没有问题的。目前也有很多算法进行双细胞的预测，但是一定要注意，无论是哪种方法，都是有风险的，可能能够去除双细胞，但也有可能造成你正常细胞的删除，所以在使用的时候需要慎重。这也就是我们之前一直没有进行这个分析的原因，因为这是一个可选的内容。目前基于R的去除双细胞的方法有DoubletFinder、DoubletDecon，python的有scrublet、DoubletDetection。考虑到很多人使用R以及文献中出现最多的，我们这里讲解DoubletFinder。用一句通俗的语言解释它的原理就是：人为模拟双细胞，然后用你的细胞与模拟双细胞对比，那么更靠近的就有可能是双细胞。
DoubletFinder原文和github可以查看学习更深入的原理：https://github.com/chris-mcginnis-ucsf/DoubletFinder。至于DoubletFinder的使用，在其官网上其实有详细的教程，也很简单。但是需要注意的是，如果你目前使用，最好更新成最新的版本，DoubletFinder也支持seurat V5，此外还需要注意，更新后的DoubletFinder有些函数的名称改变了，所以如果你参考之前的一些网上的教程可能会有问题。DoubletFinder的使用需要注意，只能单个样本使用，不能用整合的data。
那么这里我们的改变只是，我们让这个分析过程更加简单，包装了一个函数，快捷的进行分析，不用一步步的写。
首先我们单个的构建seurat object，要跑完整个seurat流程才可以进行分析：
setwd("D:\\KS项目\\公众号文章\\单细胞ATAC-signac分析教程")
library(Seurat)
dirs = list.files('./',pattern='[scRNA]$')
seurat_object = lapply(dirs,function(dirs){ CreateSeuratObject(counts = Read10X(dirs),
project = dirs,min.cells = 3, min.features = 200)})
#qc
for (i in seq_along(seurat_object)) {
seurat_object[[i]][["percent.mt"]] <- PercentageFeatureSet(seurat_object[[i]], pattern = "^MT-")
seurat_object[[i]][["percent.rp"]] <- PercentageFeatureSet(seurat_object[[i]],pattern = "^RP[SL]")
seurat_object[[i]] <- subset(seurat_object[[i]],
subset = nCount_RNA > 1000 &
nFeature_RNA < 6000 &
percent.mt < 20 & percent.rp < 40)
}
names(seurat_object) <- c("AA","HC","SD")
for (i in seq_along(seurat_object)) {
seurat_object[[i]] <- NormalizeData(seurat_object[[i]]) %>% FindVariableFeatures(., nfeatures = 4000)%>%
ScaleData(., vars.to.regress = c("percent.mt", "percent.rp"), verbose = FALSE)
}
#run PCA
for (i in seq_along(seurat_object)) {
seurat_object[[i]] <- RunPCA(seurat_object[[i]], npcs = 50, verbose = FALSE)
}
# Seurat::ElbowPlot(seurat_object[[1]], ndims = 50)
#cluster
for (i in seq_along(seurat_object)) {
seurat_object[[i]] <- RunUMAP(seurat_object[[i]], reduction = 'pca', dims = 1:20)%>%
FindNeighbors(., reduction = 'pca', dims = 1:20)%>%
FindClusters(., resolution = 0.4)
}
上面就是简单的seurat流程，这是按照seurat V4跑的。安装下最新的DoubletFinder包。看看函数的参数：我们将整个流程包装在了函数里面：
#之前的卸载了，咱们用最新的测试
# remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')
# library(DoubletFinder)
image.png
去除双细胞：
seurat_object[[1]] <- ks_detectDoublet(seurat_object[[1]],dims = 1:30,estDubRate=0.065,
ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")
seurat_object[[2]] <- ks_detectDoublet(seurat_object[[2]],dims = 1:30,estDubRate=0.025,
ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")
seurat_object[[3]] <- ks_detectDoublet(seurat_object[[3]],dims = 1:30,estDubRate=0.060,
ncores = 2, SCTransform=F, Homotypic=F, annotation="seurat_clusters")
#plot
duPlot <- list()
for (i in seq_along(seurat_object)) {
p = DimPlot(seurat_object[[i]], group.by = "DF.classify")+ggtitle(unique(seurat_object[[i]]$orig.ident))
duPlot[[i]] <- p
}
library(cowplot)
plot_grid(duPlot[[1]],duPlot[[2]],duPlot[[3]], ncol = 3)
###remove Doublet
for (i in seq_along(seurat_object)) {
seurat_object[[i]] <- subset(seurat_object[[i]], subset = (DF.classify == "Singlet"))
}
最后我们重新合并样本，跑流程即可：
###recluster
scRNA_ha <- merge(seurat_object[[1]], y=c(seurat_object[[2]],seurat_object[[3]]),
add.cell.ids =c("AA","HC","SD"))
scRNA_ha <- NormalizeData(scRNA_ha, normalization.method = "LogNormalize") %>%
FindVariableFeatures(., selection.method = "vst", nfeatures = 3000)
scRNA_ha <- ScaleData(scRNA_ha, vars.to.regress = c("percent.mt", "percent.rp"), verbose = T)
scRNA_ha <- RunPCA(scRNA_ha, npcs = 50, verbose = FALSE)
Seurat::ElbowPlot(scRNA_ha, ndims = 50)
#去批次
library(harmony)
scRNA_ha <- RunHarmony(scRNA_ha, "orig.ident")
scRNA_ha <- RunUMAP(scRNA_ha, dims = 1:20, reduction = "harmony")
scRNA_ha <- FindNeighbors(scRNA_ha, reduction = "harmony",dims = 1:20)
scRNA_ha <- FindClusters(object = scRNA_ha , resolution = seq(from = 0.1, to = 1.0, by = 0.1))
Idents(scRNA_ha) <- "RNA_snn_res.0.4"
scRNA_ha$seurat_clusters <- scRNA_ha@active.ident
DimPlot(scRNA_ha,label = T)
这就是今天所有内容了，感兴趣可以学习一下，更多参考：
https://mp.weixin.qq.com/s/xtO08NXWaMN5BB3CuBr-Xg
网友评论