做完单细胞 ATAC 质控，最容易看到的是“剩下多少细胞”，最容易漏掉的是这个数字的分母。原始片段文件里出现过的条码、写入 Arrow 的细胞、进入最后分析的细胞，是三个不同的集合。把它们混在一起，可能会误以为调低 TSS 阈值就能找回全部细胞。要找到真正影响样本的那一步，可以先把条码名单逐层对齐，再把两条过滤条件拆开计数。
我们用 ArchR 作者提供的一对小文件做这个检查。片段文件里有 150 个不同条码，示例 Arrow 中能读到 127 个条码。对后者再应用一组教学用门槛，留下 86 个。下文既能复算这三个数，也能回答第三个数怎么变来的。这个小文件只含两条染色体上的片段，适合观察数据结构和计数逻辑。ArchR 官方测试文件
先对齐片段文件和 Arrow 的条码
在同一个目录中准备 PBSmall.tsv.gz 和 PBSmall.arrow。下面代码会下载固定版本的官方测试文件，使用 Bioconductor 的 rhdf5 读取已经生成的 Arrow。只需提取几列元数据就能核查这个例子，无须先建立完整的 ArchR 项目。本文使用 R 4.5.2、rhdf5 2.54.1。示例文件内部记录的 ArchR 版本为 1.0.2。rhdf5 使用说明
library(rhdf5)
base <- paste0("https://raw.githubusercontent.com/GreenleafLab/ArchR/",
"6feec354ad6c8052ddbc4626a2ca2d858ed465bf/inst/testdata/")
for (nm in c("PBSmall.tsv.gz", "PBSmall.arrow")) {
if (!file.exists(nm)) download.file(paste0(base, nm), nm, mode = "wb")
}
frag <- read.delim(gzfile("PBSmall.tsv.gz"), header = FALSE,
col.names = c("chr", "start", "end", "barcode", "support"))
get_qc <- function(k) as.vector(h5read("PBSmall.arrow", paste0("Metadata/", k)))
qc <- data.frame(cell = get_qc("CellNames"), nFrags = get_qc("nFrags"),
TSS = get_qc("TSSEnrichment"))
stopifnot(!anyDuplicated(qc$cell), all(is.finite(qc$nFrags)),
all(is.finite(qc$TSS)))
raw_ids <- unique(frag$barcode)
c(fragment_barcodes = length(raw_ids), arrow_cells = nrow(qc),
common = length(intersect(raw_ids, qc$cell)))
输出依次为 150、127、127。这里原始文件有 29,399 条片段记录，同一条码可有许多行，直接数行得到的是片段记录数。Arrow 的 127 个条码全部能在输入名单找到，剩余 23 个条码没有对应的 Arrow 质控值。由这两个文件可以定位差异发生在最终二次筛选之前。要查某一个缺失条码的原因，还需要建档时的参数、允许进入的条码名单和日志。把这 23 个直接登记为“TSS 不合格”会给后续排查指错方向。
在多样本项目里对条码时还要保留样本身份。两个文库完全可能都出现同一串 10x 条码。若项目导出的细胞名带有样本前缀，先按文库建对应关系，再在各文库内部比较条码。不加区分地删除所有前缀再取交集，会让另一个文库中的同名条码看起来像是找回了丢失的细胞。
两个指标回答的是两件不同的事
nFrags 描述每个细胞可用的核片段数量。片段很少时，即使细胞有清晰的启动子信号，许多调控区域仍可能采样不足。TSSEnrichment 反映转录起始位点附近的插入信号相对于两侧背景的富集。一个条码可以积累很多片段，但如果这些片段在起始位点附近没有形成清晰的富集，就需要回看背景信号和注释。两个指标联合使用，才能区分“信息量有限”和“片段不少但信号与背景难以分开”。ArchR 对逐细胞质控指标的定义
本例有意用 nFrags >= 150、TSS >= 5 对现有 127 行做一次二次筛选。150 只是让这个裁剪过的小文件出现四种情况，不能直接用于完整文库。最简单的诊断是保留两列真假值，做一张二维计数表。如果接连筛两次只保留最后的行数，就会丢掉两条条件重叠的信息。
frag_ok <- qc$nFrags >= 150
tss_ok <- qc$TSS >= 5
table(nFrags_pass = factor(frag_ok, levels = c(FALSE, TRUE)),
TSS_pass = factor(tss_ok, levels = c(FALSE, TRUE)))
c(total = nrow(qc), kept = sum(frag_ok & tss_ok),
fails_either = sum(!frag_ok | !tss_ok))
| 在这 127 个条码中 | TSS 未达 5 | TSS 达到 5 | 合计 |
|---|---|---|---|
| 核片段数未达 150 | 23 | 13 | 36 |
| 核片段数达到 150 | 5 | 86 | 91 |
| 合计 | 28 | 99 | 127 |
如果先按片段数筛，127 个条码里剩 91 个，这时提高 TSS 要求又排除 5 个。如果先按 TSS 筛，会先剩 99 个，再排除 13 个。最后都是 86 个。两种顺序的中间数字不同，是因为 23 个条码同时未达两项要求。写质控记录时，可以用“仅片段数不足”“仅 TSS 不足”“两项都不足”三个互斥原因，分别填 13、5、23。读者既能复算合计剔除的 41 个，也能看懂究竟是哪条边界值得仔细检查。这里两项都不足的 23 个仍位于 Arrow 内，与前面完全没有进入 Arrow 的 23 个是不同的条码。
把这张表移到自己的 ArchR 项目
已有项目可以通过公开接口 getCellColData() 提取当前细胞的质控列。下面假设你的项目对象名为 proj。先统计每个样本当前可见的细胞数，再看中位数和四分位数。同一个合并散点图里，如果一个大文库占了绝大多数点，小文库整体的位移很容易被遮住。ArchR 项目数据提取接口
md <- as.data.frame(ArchR::getCellColData(
proj, select = c("Sample", "nFrags", "TSSEnrichment")))
print(table(md$Sample))
by(md, md$Sample, function(d) {
rbind(nFrags = quantile(d$nFrags, c(.25, .5, .75), na.rm = TRUE),
TSS = quantile(d$TSSEnrichment, c(.25, .5, .75), na.rm = TRUE))
})
如果某个样本整体核片段偏少，先看输入的有效片段、测序覆盖和片段长度分布。如果核片段不少而整批的 TSS 分布都偏低，先核对比对参考组装、染色体命名以及计算富集使用的起始位点注释，再看该文库的聚合 TSS 曲线。若只有少数点处在筛选边缘，可单独查看它们在后续降维中的位置，以及阈值轻微变化后目标群的成员是否明显改变。单靠一个低分无法区分文库质量、注释错配与细胞本身的差异。
还要记住这张表的起点。项目对象只含当前纳入的细胞，之后放宽筛选只能找回仍有记录的条码。若需要重新评估建档阶段排除的条码，应回到每个样本对应的 fragments 文件，核对排序、bgzip 压缩、参考基因组与 createArrowFiles() 当时使用的门槛。在独立目录重建一份可比较的 Arrow，保存原始结果和参数，再逐条码做差集。最后交付每个文库的输入条码数、入档数、最终保留数及两项指标的交叉表，同事就能沿原路线核查一次。createArrowFiles 参数说明
把门槛挪一小步，实际增加的是哪些条码？
对同一份 Arrow，可以固定一个指标，只移动另一个。以片段数一百五十、富集分五为基线，保留八十六个条码。将片段数要求降到一百而保持富集分不变，保留九十九个。把片段数升到两百，保留四十九个。再把片段门槛恢复到一百五十，分别把富集分设为四和六，保留数为九十一和七十二。以下代码直接从上文的数据框计算。把两个门槛一次只动一个，才能看清哪些变化来自片段量，哪些来自起始位点信号。
grid <- data.frame(frag = c(100, 150, 200, 150, 150), tss = c(5, 5, 5, 4, 6))
grid$kept <- mapply(function(f, t) sum(qc$nFrags >= f & qc$TSS >= t),
grid$frag, grid$tss)
print(grid)
rescued <- subset(qc, nFrags >= 100 & nFrags < 150 & TSS >= 5)
stopifnot(nrow(rescued) == 13)
这些增加的条码值得单独保存。例如放宽片段要求时找回的十三个，全都跨过了富集分门槛，可以逐条查看片段长度、检测到的开放区域，以及这些条码在同一个文库中的分布位置。反过来，从一百五十提高到两百所减少的三十七个条码都有达标的起始位点信号。若某类稀少细胞恰好集中在这个深度区间，门槛会连带改变它的观察频率。要判断是否确有这种情况，需要独立标记或后续注释来核对。本例只展示哪些条码受边界影响，不能从条码名称猜测组织成分。
实际项目最好把每个文库的敏感性表做成相同的行列，然后同时报告绝对数与比例。甲文库从一万个细胞增加五十个，乙文库从一百个增加五十个，绝对数看起来一样，影响却完全不同。把两者合并再计算比例会让深度较大的甲文库掩盖乙文库的变化。如果不同条件由不同批次完成，还可在样本表中列出每个条件对应的建库日期和处理方式。某一条件的保留率始终偏低时，先查条件和批次有没有分开安排，再讨论细胞本身的染色质状态。
从二维计数回到片段信号
当一个文库整体富集分偏低时，可再画聚合起始位点曲线。选择与比对参考一致的起始位点注释，把插入信号按相对位置汇总，观察中心是否有峰、两侧背景是否均匀。不同文库并排展示时沿用同一套窗口与归一化方式。假如一份文库曲线看起来几乎平坦，先核对注释和片段能否在相同染色体上对齐，再查细胞筛选的先后顺序。若在合并图上画的是全体条码，而单细胞统计来自严格筛选后的子集，两者会显示不同的信号强度。
交付给合作者时，给每个文库一张条码级结果表会比只存最终通过名单更方便。用一列保存输入名，一列保存项目中的带样本名前缀的名字，同时保存原始片段数、富集分和两个独立的通过标记。上游已筛掉、无法取得富集分的条码，在这一列留缺失并注明来源。后续若要重定门槛，同事可以从表中直接重得四格计数，并挑出要再检查的条码。真正需要重新运行上游建档的情况也会一目了然。
下载代码和示例数据，在 Linux 上复现
下载本篇 R 代码（ZIP，4.0 KiB） · 下载配套示例数据（ZIP，836.5 KiB）
将两个 ZIP 解压到同一位置，进入 027-archr 目录，运行 Rscript --vanilla run.R。打开 results/metrics.tsv 查看结果，代码包中的 expected/metrics.tsv 可供逐行核对。
本次从下载包复现的步骤是“用 rhdf5 从作者的 Arrow 与片段核对过滤结果”。Linux x86_64、R 4.3.3、rhdf5 2.40.0 的输出为片段文件条码 150，Arrow 细胞 127，两项质控通过 86。
下载版另存 arrow-qc.csv。把 TSS 富集和片段数分别设为通过标记，再取两个标记的交集。只有 Arrow 中的条码有这两项可供统计。原片段文件中未写入 Arrow 的条码应保留在另一张清单中。
实跑后的经验总结
在这份真实小型 Arrow 文件里，两项质控的先后顺序会改变中间淘汰数，最后同时达标的条码始终是八十六个。复核时先保留每个样本的完整条码名单，再给两项阈值各留一列真假值。对在建档前已消失的条码另列来源，后来改变任意一个门槛都能直接找回实际受影响的细胞。
当门槛已经在代表文库上核对清楚，后面有整批 fragments 需要按相同流程重建 Arrow，可以按项目周期租用环境继续跑。云生信生物信息学分析服务器租赁
