图 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） 也会进入空液滴，形成低水平计数。

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

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

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

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

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

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

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

固定 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")

若任务必须得到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 的情况。

只有 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(ed$FDR) & ed$FDR <= 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+

支持单位：