精准分析单细胞数据：使用DoubletFinder去除双胞的详细步骤

天意生信云

发布于 2025-02-25 19:53:08

94010

代码可运行

运行总次数：0

代码可运行

在单细胞RNA测序中，我们希望对每个单独的细胞进行测序分析。然而，在样本制备过程中，可能会出现多个细胞聚集在一起的现象，形成所谓的“双细胞”或“多细胞”。双胞的存在可能会对下游分析（如聚类、差异表达等）产生误导，因此在分析过程中识别和过滤双胞是一个重要的步骤。

准备工作

为了满足本篇笔记分析的需求，我从10X官网上下载了一组小鼠心脏的单细胞测序下机数据，并按照之前笔记的教程单细胞数据分析 | 单细胞计数矩阵（Seurat)使用cellranger v8.0.0和小鼠的GRCm39版本的参考基因集对其进行了比对定量工作。


# 下载数据并解压
# 下载数据
$ cat  download.sh
# mouse_heart_5k_1
wget https://cf.10xgenomics.com/samples/cell-exp/7.0.0/5k_mouse_heart_CNIK_3pv3/5k_mouse_heart_CNIK_3pv3_fastqs.tar
# mouse_heart_10k_2
wget  https://s3-us-west-2.amazonaws.com/10x.files/samples/cell-exp/3.0.0/heart_10k_v3/heart_10k_v3_fastqs.tar
# mouse_heart_1k_3
wget  https://cf.10xgenomics.com/samples/cell-exp/3.0.0/heart_1k_v3/heart_1k_v3_fastqs.tar

# 解压
for  line in  $(ls ./\*.tar)
do
tar  -xvf  ${line}
done

批量比对定量：


# 先创建一个与样品编号对应的分组编号表格
$ cat sample_list.txt
# sample_list.txt 内容
5k_mouse_heart_CNIK_3pv3 mouse1
heart_10k_v3 mouse2
heart_1k_v3 mouse3
# 编写批量比对定量脚本
$ cat run_qc.sh
reference=/file/to/reference/refdata-gex-GRCm39-2024-A
fq_path=/file/to/data/0_rawdata
cat sample_list.txt | while  read  col1  col2
do
        cellranger  count --id= $col2  # 输出文件的名字，即分析时的分组编号
        --transcriptome = ${reference} 
        --fastqs = ${fq_path}/${col1}_fastqs/  # 需要质控的下机fastq文件路径
        --sample = $col1  # sample是输入文件的名字，即样品编号，通常是两种，一种是SRR开头的编号，一种是10X的文件名去掉最后的_fastqs
        --create-bam = true 
        --localcores = 20 
        --localmem = 200
done

定量完成之后就会得到标准分析的三文件:

Doubletfinder去除双胞

本文按照seurat v5的处理方式，选择Doubletfinder去除双胞。


# DoubletFinder v2.0.4 及配套seurat V5安装
# R命令行安装：DoubletFinder
remotes::install_github('chris-mcginnis-ucsf/DoubletFinder')
# shell命令行安装下述软件
conda  install  conda-forge::r-seurat
conda  install  bioconda::bioconductor-dropletutil

根据官网工具（https://github.com/chris-mcginnis-ucsf/DoubletFinder）的使用要求：doubleFinde是对单个样本处理的，并且输入数据不能包含低质量的细胞簇。那么在使用doublefinder之前，先进行过滤。


# 加载需要的R包
library(Seurat)
library(cowplot)
library(DropletUtils)
library(DoubletFinder)
# 读入单细胞矩阵数据，并按照官网要求进行过滤预处理
object <- Read10X('your/path/to/filtered_feature_bc_matrix', gene.column = 2)
# 创建Seurat对象
object <- CreateSeuratObject(counts = object.data, project = "brain_1", min.cells = 3, min.features = 200)
# 计算线粒体基因百分比，0值补成非0的极小数
object[["percent.mt"]] <- PercentageFeatureSet(object, pattern = "ATP6|ATP8|COX1|COX2|COX3|CYTB|ND1|ND2|ND3|ND4|ND4L|ND5|ND6")
if (sum(object@meta.data$percent.mt) == 0) {
    object@meta.data$percent.mt[1] <- 0.000001
}
# 按照相对宽泛的指标过滤
object <- subset(object, subset = nFeature_RNA > 200 & nFeature_RNA < 6000 & nCount_RNA < 20000 & percent.mt < 5)

对于seurat的预处理流程，有NormalizeData标准化和SCT标准化两种方法，对于本文中的例子，选择哪个都可以。


# 三步标准化流程
object <- NormalizeData(object)
object <- FindVariableFeatures(object, selection.method = "vst", nfeatures = 2000)
object <- ScaleData(object)
object <- RunPCA(object)
object <- RunUMAP(object, dims = 1:20)
# SCT标准化流程
object <- SCTransform(object)
object <- RunPCA(object)
object <- RunUMAP(object, dims = 1:20)
# 用分群结果代替注释信息
object <- FindNeighbors(object, reduction = "pca", dims = 1:20)
object <- FindClusters(object, resolution = 0.5)

DoubletFinder去除双胞步骤：

1）从现有的 scRNA-seq 数据生成人工双联体 2）预处理合并的真实-人工数据 3）进行PCA，利用PC距离矩阵求出每个细胞人工k近邻的比例（pANN） 4）根据预期双联体数量排序和阈值 pANN 值。


# 对"object"这个单细胞对象进行pN-pK参数扫描，以生成人工双细胞并计算每个细胞的pANN值
sweep.res.list <- paramSweep_v3(object, PCs = 1:20, sct = T)
# 对参数扫描的结果进行汇总，计算每种pN和pK组合下的双细胞检测指标
sweep.stats_object <- summarizeSweep(sweep.res.list, GT = FALSE)
# 找到最优的pK参数
bcmvn <- find.pK(sweep.stats_object)
# 提取出全局最优的pK值，储存于"pK_bcmvn"
pK_bcmvn <- as.numeric(bcmvn$pK[which.max(bcmvn$BCmetric)])
# 把聚类分群信息当作注释信息
annotations <- object@meta.data$SCT_snn_res.0.5
# 估计的同源双细胞比例
homotypic.prop <- modelHomotypic(annotation)
# 计算总的双细胞数量（假设双细胞形成率为 7.5%）
nExp_poi <- round(0.075 * nrow(object@meta.data))
nExp_poi.adj <- round(nExp_poi * (1 - homotypic.prop))
# 使用确定好的参数鉴定doublets
object <- doubletFinder_v3(object, PCs = 1:30, pN = 0.25, pK = pK_bcmvn,
                                  nExp = nExp_poi.adj, reuse.pANN = FALSE, sct = T)

以上就是doubletfinder去除双胞的关键步骤，可视化去除双胞结果。


# 可视化
DimPlot(object, reduction = "UMAP", group.by = "DF.classifications_0.25_33_9717", raster = FALSE)
# 将双细胞剔除后生成新的对象scRNA_harmony.singlet，便于后续分析
object.singelt <- subset(object, subset = DF.classifications_0.25_33_9717 == "Singlet")

# 输出结果保存为标准格式
DropletUtils::::write10xCounts("singlet_matrix", object_Singlet@assays$RNA@counts, version = "3")

批量化脚本：


# 设置R环境配置
options(future.globals.maxSize = 2000 * 1024^2)
# 加载需要的R包
library(Matrix)
library(Seurat)
library(DropletUtils)
library(DoubletFinder)
library(optparse)
# 解析命令行参数
option_list <- list(
  make_option(c("-i", "--input"), type="character", help="输入矩阵目录路径"),
  make_option(c("-s", "--sample"), type="character", help="样本名称"),
  make_option(c("-f", "--max_features"), type="integer", default=6000, help="nFeature_RNA最大值阈值"),
  make_option(c("-c", "--max_counts"), type="integer", default=20000, help="nCount_RNA最大值阈值"),
  make_option(c("-m", "--max_mt"), type="numeric", default=5, help="线粒体百分比最大值阈值"),
  make_option(c("-o", "--output"), type="character", help="输出目录路径")
)
opt <- parse_args(OptionParser(option_list=option_list))
# 创建输出目录
if (!dir.exists(opt$output)) dir.create(opt$output, recursive = TRUE)
# 数据读取和预处理
data <- Read10X(opt$input, gene.column = 2)
seurat_obj <- CreateSeuratObject(counts = data, project = opt$sample, min.cells = 3, min.features = 200)
# 线粒体基因计算
seurat_obj[["percent.mt"]] <- PercentageFeatureSet(seurat_obj, pattern = "^mt-")
if (sum(seurat_obj$percent.mt) == 0) seurat_obj$percent.mt[1] <- 1e-6  # 处理零值情况
# 细胞过滤
seurat_obj <- subset(seurat_obj, subset = nFeature_RNA > 200 & nFeature_RNA < opt$max_features &
                    nCount_RNA < opt$max_counts & percent.mt < opt$max_mt)
# Seurat标准流程
seurat_obj <- SCTransform(seurat_obj)
seurat_obj <- RunPCA(seurat_obj)
seurat_obj <- RunUMAP(seurat_obj, dims = 1:20)
seurat_obj <- FindNeighbors(seurat_obj, dims = 1:20)
seurat_obj <- FindClusters(seurat_obj, resolution = 0.5)

# 双胞检测
sweep.res <- paramSweep(seurat_obj, PCs = 1:20, sct = TRUE)
sweep.stats <- summarizeSweep(sweep.res, GT = FALSE)
bcmvn <- find.pK(sweep.stats)
pK <- as.numeric(as.character(bcmvn$pK[bcmvn$BCmetric == max(bcmvn$BCmetric)]))
# 双胞比例计算
homotypic.prop <- modelHomotypic(seurat_obj$SCT_snn_res.0.5)
nExp_poi <- round(0.075 * ncol(seurat_obj))
nExp_poi.adj <- round(nExp_poi * (1 - homotypic.prop))
# 双胞检测和分类
seurat_obj <- doubletFinder(seurat_obj, PCs = 1:20, pN = 0.25, pK = pK,
                                  nExp = nExp_poi, reuse.pANN = FALSE, sct = TRUE)
seurat_obj <- doubletFinder(seurat_obj, PCs = 1:20, pN = 0.25, pK = pK,
                           nExp = nExp_poi.adj, reuse.pANN = FALSE, sct = TRUE)
# 结果整合
df_cols <- grep("DF.classifications", names(seurat_obj@meta.data), value = TRUE)
seurat_obj$Doublet <- ifelse(seurat_obj@meta.data[[df_cols[1]]]=="Doublet" &
  seurat_obj@meta.data[[df_cols[2]]]=="Singlet", "Doublet_low",
  ifelse(seurat_obj@meta.data[[df_cols[1]]]=="Doublet", "Doublet_high", "Singlet"))
# 保存单细胞矩阵
singlets <- subset(seurat_obj, subset = Doublet == "Singlet")
counts_matrix <- GetAssayData(singlets, assay = "RNA", layer = "counts")
write10xCounts(opt$output, counts_matrix, version = "3", overwrite = TRUE)

此时得到的文件就是去除双胞的标准三文件！

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2025-02-25，如有侵权请联系 cloudcommunity@tencent.com 删除

数据