从 Barcode Rank Plot 到细胞条形码过滤的原理与 R 实操

图 1  Barcode Rank 曲线（来自10X官网）。Cliff（悬崖，陡降的一段），Knee（膝盖，变平的一段）

核心结论： raw 矩阵（图上蓝色线+灰色线）保存的是所有检测到的有效 barcode，其中绝大多数并不对应完整细胞。只有 raw 时，应先做 cell calling，生成候选细胞矩阵（图上蓝色线），再进入 Seurat 或 Scanpy。

先看你拿到的是哪种矩阵

根据文件名判断：有 filtered 就从 filtered 开始；只有 raw，必须先做细胞识别（cell calling）。raw 不是更完整的“细胞矩阵”，而是尚未区分含细胞液滴与背景液滴的计数空间。

图 2  GEO中常见的矩阵分流方式

先把四个概念说清楚

10x 微流控把细胞悬液、试剂和带 barcode 的凝胶珠封装成 GEM 液滴。同一液滴中的 RNA 共享 cell barcode；UMI 用于区分原始 RNA 分子并去除 PCR 重复。游离在悬液中的环境RNA（ambient RNA） 也会进入空液滴，形成低水平计数。

术语	可以怎样理解	分析时代表什么
cell barcode	同一液滴内分子的来源标签	矩阵中的一列，不一定就是一个细胞
UMI	单个RNA 分子的去重标签	同一 barcode 的总 UMI 反映捕获到的分子量
ambient RNA	漂在悬液中的游离 RNA	空液滴也可能因此出现少量计数
cell calling	判断液滴是否含有细胞相关信号	把 raw 中的候选细胞挑出来

 

液滴不等于细胞

微流控不能保证每个液滴恰好包裹一个细胞。实际会同时出现空液滴、单细胞液滴、双细胞或多细胞液滴，以及只含碎片或裂解细胞 RNA 的液滴。因此，液滴 ≠ 细胞，barcode 也不能直接等同于一个细胞。

图 3  微流控液滴的四种常见状态  液滴与细胞不是一一对应关系

拿到的文件	可以做什么	不能据此做什么
filtered	创建Seurat 或 Scanpy 对象，并继续细胞级 QC	无法恢复已经被排除的 barcode
raw 和 filtered	默认从 filtered 开始；必要时用 raw 重新 cell calling	不要为了细胞更多就无条件改用 raw
只有 raw	先做 cell calling，再使用新写出的 filtered 矩阵	不能把约 100 万个barcodes全部当成细胞

raw 为什么不能直接用于下游分析

10x 实验产生的 barcode 数远多于实际细胞数。以 GEO 样品 GSM6395363 为例，下载下来的barcodes.tsv文件有约 600 万行；多数来自空液滴或低水平背景。细胞裂解释放的 环境RNA会被空液滴捕获，因此 raw 中常见大量彼此相似的低 UMI barcode。

直接把 raw 导入 Seurat，即使设置 min.features，也只是用固定门槛粗筛，不能替代 cell calling。门槛太低会保留背景，太高又会误删低 RNA 细胞，例如中性粒细胞 RNA 含量和可检测转录本通常较低，设置较高阈值，则中性粒细胞会在 10x 数据中缺失或低估；判断时应结合细胞标记、样本组成和后续 QC，而不是只看 UMI 总数。

图 4  raw 矩阵直接进入下游分析时的常见后果

Barcode Rank Plot 到底在看什么

Barcode rank plot 先汇总每个 barcode 的 UMI，再按总 UMI 从高到低排序。横轴是 rank，纵轴是 UMI，通常使用双对数坐标。左上方多为含细胞液滴，右下方多为背景；两者之间常出现 cliff（悬崖，陡降的一段）或 knee（膝盖，陡降到平缓的转折点）。

一条理想的 Barcode Rank Plot 曲线，顺序是：平缓（细胞）→ cliff 悬崖（陡降）→ knee 拐点（变平）→ 长尾（背景）。曲线用于观察细胞信号与背景是否分开，不能单独给每个 barcode 定性。混合高、低 RNA 细胞时可出现多段平台；碎片多或细胞裂解严重时，背景尾部也会抬高。

图 5 Barcode Rank Plot (来自10X官网)   蓝色为 Cells  灰色为 Background

几种需要警惕的曲线形态

图 6  四种需要进一步排查的 barcode rank 形态示意图

这些形态可能与细胞裂解、环境 RNA 偏高、捕获效率低或实际回收量不足有关。曲线只能提示问题方向，还要结合 web summary、活率、reads in cells、每细胞 UMI 和样本组成判断。

Cell Ranger 如何从 raw 得到 filtered

Cell Ranger 的 Gene Expression cell calling 可概括为两步：OrdMag 先识别高 RNA 细胞，EmptyDrops 再利用 ambient RNA 模型寻找低 RNA 但表达谱不像背景的液滴。filtered_feature_bc_matrix 是两步结果的合并，不是固定 UMI 阈值的产物。

图 7  cell calling 与后续细胞级 QC 的边界

第一步 OrdMag 先抓住高 RNA 细胞

给定预期细胞数 N 时，OrdMag 根据前 N 个 barcode 的 UMI 分布建立初始阈值；Cell Ranger 7.0 起也可自动估计 expect-cells。阈值随数据变化，因此不应跨项目固定为 500 或 1,000 UMI。

第二步 EmptyDrops 救回低 RNA 细胞

EmptyDrops 先用低计数 barcode 建立 ambient RNA 背景，再检验候选 barcode 的表达组成是否显著偏离背景，从而找回部分低 RNA 细胞。

按当前 10x 说明，EmptyDrops 候选通常还需达到 max(500, 背景区最高 UMI + 1)；Cell Ranger v9.0 将 3′/5′分析的 FDR 调整为 0.001。FDR 是统计判定阈值，不是“取 UMI 前 0.1%”。

Cell ranger版本差异要写进分析记录

cell calling 规则会随cell ranger版本变化。v8.0 移除了一个会在高 UMI 样本中过度删除低 RNA 细胞的限制；复现或比较项目时，应记录 Cell Ranger 版本、chemistry 和参数，并以对应版本的算法页与发布说明为准。

版本	cell calling 相关变化	对使用者的影响
v7.0	开始自动估计 expect-cells，也允许用户手动提供	多数样本不再必须先猜回收细胞数
v7.1	NextGEM 的自动搜索范围记录为约 2 到 45000 个细胞	超大规模实验要留意化学版本和搜索范围
v8.0	EmptyDrops 最低候选阈值改为 max(500, 背景区最大 UMI 加 1)，并移除旧的初始细胞 UMI 中位数 1% 限制	更有机会保留与背景不同的低 RNA 细胞
v9.0	3′ 和 5′分析的 EmptyDrops FDR 调为 0.001	调用更严格，结果可能比旧版本少
v10.1 当前	当前算法说明仍采用 OrdMag 加 EmptyDrops；未见上述规则被再次改写	报告中仍应记录准确版本和参数

 

同一数据由不同 Cell Ranger 版本处理，filtered 细胞数可能不同。报告结果时应同时说明软件版本和 cell-calling 参数。

三种从 raw 开始的处理方法 

方法	结果数量	优点	主要风险	建议用途
固定 UMI 阈值	不固定	简单直观	样本间不可直接复用，易漏低 RNA 细胞	快速检查
按 UMI 取前 N	可精确到 N	可得到约 1 万或恰好 1 万	把目标数量当真值，忽略表达谱	审计基线或明确 force-cells
emptyDrops	由数据决定	利用 UMI 与表达组成，可找到低表达细胞	可能检出碎片或受损细胞，仍需 QC	只有 raw 时的首选

 

方法一 按固定 UMI 阈值过滤

固定 UMI 阈值适合快速检查，不适合作为跨样本通用规则。阈值应结合 barcode rank、样本类型和实验信息确定；它不会判断表达谱是否像 ambient RNA。

library(Matrix)

library(DropletUtils)

 

raw_dir <- "raw_feature_bc_matrix"

out_dir <- "filtered_by_umi"

min_umi <- 500

 

sce <- read10xCounts(raw_dir, col.names = TRUE)

is_gex <- rowData(sce)$Type == "Gene Expression"

total_umi <- Matrix::colSums(counts(sce)[is_gex, , drop = FALSE])

keep <- total_umi >= min_umi

 

sce_f <- sce[, keep]

write10xCounts(

  out_dir, counts(sce_f),

  barcodes = colnames(sce_f),

  gene.id = rowData(sce_f)$ID,

  gene.symbol = rowData(sce_f)$Symbol,

  gene.type = rowData(sce_f)$Type,

  version = "3", overwrite = TRUE

)

cat("保留 barcode 数:", sum(keep), "\n")

方法二 按总 UMI 取topN

若任务必须得到topN（例如10,000个barcodes），按 UMI 排序取前 10,000 个最透明，但这属于人为固定数量，不是统计识别。应同时输出 rank plot 和第 10,000 位的 UMI；并列时，下面的代码仍严格保留 10,000 个 barcode。

target_n <- 10000L

ord <- order(total_umi, decreasing = TRUE)

keep_index <- ord[seq_len(min(target_n, length(ord)))]

keep <- rep(FALSE, length(total_umi))

keep[keep_index] <- TRUE

 

cat("第 10000 位 barcode 的 UMI:",

    total_umi[ord[min(target_n, length(ord))]], "\n")

cat("保留 barcode 数:", sum(keep), "\n")

若已知上样量、预期回收率和 web summary，Cell Ranger 的 --expect-cells=10000 通常比下游硬截取更合适。--force-cells=10000 会绕过自动 cell calling，只适用于确需固定数量且已经检查 rank plot 的情况。

方法三 使用 emptyDrops 一次完成过滤与绘图

只有 raw 矩阵时，推荐用下面这段脚本一次完成 cell calling、filtered 矩阵导出、诊断表保存和 barcode rank plot 绘制。蓝色表示 emptyDrops 判定的候选细胞，灰色表示背景；这些候选细胞仍需接受后续细胞级 QC 和 doublet 检查。

#!/usr/bin/env Rscript

 

suppressPackageStartupMessages({

  library(DropletUtils)

  library(SingleCellExperiment)

  library(Matrix)

})

 

输入与输出

raw_dir      <- "raw_feature_bc_matrix"

filtered_dir <- "filtered_feature_bc_matrix_emptydrops"

out_png      <- file.path(filtered_dir, "barcode_rank_plot.png")

fdr_cutoff   <- 0.001

lower_umi    <- 100

niters       <- 30000

seed         <- 1234

 

读取 raw；多模态矩阵只用 Gene Expression 做 cell calling

set.seed(seed)

sce <- read10xCounts(raw_dir, col.names = TRUE)

feature_type <- rowData(sce)$Type

is_gex <- if (is.null(feature_type)) rep(TRUE, nrow(sce)) else feature_type == "Gene Expression"

gex_counts <- counts(sce)[is_gex, , drop = FALSE]

total_umi <- Matrix::colSums(gex_counts)

 

计算 knee、inflection，并运行 emptyDrops

br <- barcodeRanks(gex_counts, lower = lower_umi)

ed <- emptyDrops(

  gex_counts,

  lower = lower_umi,

  niters = niters,

  test.ambient = FALSE

)

is_cell <- !is.na(edFDR) & edFDR <= fdr_cutoff

 

cat("raw barcodes:", ncol(sce), "\n")

cat("called non-empty droplets:", sum(is_cell), "\n")

cat("FDR cutoff:", fdr_cutoff, "\n")

cat("knee UMI:", metadata(br)$knee, "\n")

cat("inflection UMI:", metadata(br)$inflection, "\n")

print(table(Significant = is_cell, Limited = ed$Limited, useNA = "ifany"))

 

1 输出 filtered 10x 三文件；写出时保留所有 feature

sce_filtered <- sce[, which(is_cell)]

write_args <- list(

  path = filtered_dir,

  x = counts(sce_filtered),

  barcodes = colnames(sce_filtered),

  gene.id = rowData(sce_filtered)$ID,

  gene.symbol = rowData(sce_filtered)$Symbol,

  version = "3",

  overwrite = TRUE

)

if (!is.null(rowData(sce_filtered)$Type)) {

  write_args$gene.type <- rowData(sce_filtered)$Type

}

do.call(write10xCounts, write_args)

 

2 保存每个 barcode 的统计结果

diagnostics <- data.frame(

  barcode = colnames(sce),

  total_umi = unname(total_umi),

  p_value = ed$PValue,

  fdr = ed$FDR,

  limited = ed$Limited,

  called = is_cell

)

write.csv(

  diagnostics,

  file.path(filtered_dir, "cell_calling_diagnostics.csv"),

  row.names = FALSE

)

 

3 生成 rank 数据；零 UMI barcode 不参与 log10 绘图

plot_index <- which(total_umi > 0)

ord <- plot_index[order(total_umi[plot_index], decreasing = TRUE)]

rank_df <- data.frame(

  barcode = colnames(sce)[ord],

  rank = seq_along(ord),

  total_umi = unname(total_umi[ord]),

  called = is_cell[ord]

)

rank_df$log10_rank <- log10(rank_df$rank)

rank_df$log10_umi <- log10(rank_df$total_umi)

rank_df$status <- ifelse(rank_df$called, "细胞", "背景")

write.csv(

  rank_df,

  file.path(filtered_dir, "barcode_rank_values.csv"),

  row.names = FALSE

)

 

4 绘制接近 10x 标注风格的正方形 barcode rank plot

内部坐标已经 log10 转换，但刻度显示原始数量级

cell_blue <- "#004799"

background_gray <- "#DDDDDD"

grid_gray <- "#EEEEEE"

 

x_ticks <- seq(0, floor(max(rank_df$log10_rank)), by = 2)

y_data_min <- min(rank_df$log10_umi, na.rm = TRUE)

y_data_max <- max(rank_df$log10_umi, na.rm = TRUE)

y_min_power <- floor(y_data_min)

y_max_power <- ceiling(y_data_max)

y_ticks <- seq(y_min_power, y_max_power, by = 1)

y_plot_range <- extendrange(

  c(y_min_power, y_max_power),

  f = 0.03

)

y_plot_range[1] <- y_min_power

 

 

x_labels <- ifelse(x_ticks >= 6, paste0(10^(x_ticks - 6), "M"),

                   ifelse(x_ticks >= 3, paste0(10^(x_ticks - 3), "k"),

                          format(10^x_ticks, scientific = FALSE, trim = TRUE)))

y_labels <- format(10^y_ticks, scientific = FALSE, trim = TRUE)

 

连续绘制相同类别的区段，避免跨过另一类别错误连线

draw_runs <- function(flag, colour, line_width = 2.7) {

  run_info <- rle(flag)

  run_end <- cumsum(run_info$lengths)

  run_start <- c(1, head(run_end, -1) + 1)

  for (k in which(run_info$values)) {

    idx <- run_start[k]:run_end[k]

    if (length(idx) == 1) {

      points(rank_df$log10_rank[idx], rank_df$log10_umi[idx],

             pch = 16, cex = 0.3, col = colour)

    } else {

      lines(rank_df$log10_rank[idx], rank_df$log10_umi[idx],

            col = colour, lwd = line_width)

    }

  }

}

 

png(out_png, width = 1260, height = 1260, res = 300, pointsize = 9)

par(mar = c(5.2, 7.0, 5.3, 0.8), mgp = c(2.7, 0.75, 0),

    tcl = -0.25, las = 1, xaxs = "i", yaxs = "i",

    fg = "#444444", col.axis = "#444444", col.lab = "#444444",

    cex.axis = 0.72)

plot(rank_df$log10_rank, rank_df$log10_umi, type = "n", axes = FALSE,

     xlab = "", ylab = "",

     xlim = range(rank_df$log10_rank, na.rm = TRUE),

     ylim = y_plot_range)

 

abline(v = x_ticks, h = y_ticks, col = grid_gray, lwd = 1)

 

axis(

  1, at = x_ticks, labels = x_labels,

  col = "#444444", lwd = 0, lwd.ticks = 1

)

axis(

  2, at = y_ticks, labels = y_labels,

  col = "#444444", lwd = 0, lwd.ticks = 1

)

 

左侧和底部边框随绘图范围完整延伸

box(bty = "l", col = "#444444", lwd = 1)

 

 

draw_runs(!rank_df$called, background_gray)

draw_runs(rank_df$called, cell_blue)

mtext("Barcodes", side = 1, line = 3.0, cex = 0.85, col = "#444444")

mtext("UMI Counts", side = 2, line = 4.0, las = 0, cex = 0.85, col = "#444444")

mtext("Barcode Rank Plot", side = 3, line = 3.2, cex = 1.0, col = "#444444")

legend("top", inset = c(0, -0.12), xpd = NA, horiz = TRUE,

       legend = c("Cells", "Background"), lty = 1, lwd = 3,

       col = c(cell_blue, background_gray), bty = "n", cex = 0.75)

dev.off()

 

 

 推荐的处理流程

1.确认输入确实是 raw_feature_bc_matrix，并检查 matrix.mtx.gz、features.tsv.gz、barcodes.tsv.gz 是否齐全。

2.运行方法三脚本：一次完成 emptyDrops、filtered 矩阵导出、诊断表保存和 barcode rank plot 绘制。

3.检查输出中的 FDR、Limited、knee、inflection、候选细胞数及蓝灰曲线分界，再进入细胞级 QC。

4.写出新的 filtered 矩阵和诊断表，再做基因数、UMI、线粒体比例、环境 RNA 与 doublet 检查。

5.对候选细胞再做 nFeature RNA、nCount RNA、线粒体比例、复杂度和双细胞检查。emptyDrops 评判的是 non-empty droplet，不保证每个calling都是完整单细胞。

过滤之后还要做什么

filtered 矩阵只是下游分析的起点。仍需检查每个细胞的基因数、UMI、线粒体比例、环境 RNA 污染和 doublet；应依据本样本分布与组织类型确定阈值。

环境 RNA 校正与 cell calling 是两件事。通常先从 raw 确定候选细胞，再做细胞级 QC；是否使用 SoupX、DecontX 等方法，应由样本背景污染情况决定。

常见误区

· raw 不是“更完整的细胞矩阵”。它是更完整的液滴计数空间。

· CreateSeuratObject 的 min.features 不是 Cell Ranger cell calling 的等价替代。

· FDR 越小不一定越好。过严会先将低 RNA 细胞丢掉。

· emptyDrops 显著代表非空液滴，不自动排除细胞碎片、受损细胞或双细胞。

· 想要 1 万个细胞与数据支持 1 万个细胞是两回事。

参考资料

10x Genomics  Cell Ranger Barcode Rank Plot 

10x Genomics  Cell Ranger Gene Expression Algorithm   

10x Genomics  Cell Ranger Command Line Arguments

10x Genomics  Cell Ranger Release Notes  Cell calling

Lun 等  EmptyDrops  Genome Biology  2019  

Bioconductor  DropletUtils 

DropletUtils Reference Manual  

Salcher 等：肺癌组织中性粒细胞单细胞图谱

  

微生信助力高分文章，谷歌学术10000+