Cell Ranger 给了两套矩阵，怎么确认少掉的是哪些条形码？

同一个文库跑完 Cell Ranger，原始矩阵里可能有几十万个条形码，过滤矩阵里只剩几千个。第一次接手数据的人常问，是测序丢了细胞，还是软件删掉了基因？这两个目录回答的是同一批分子属于哪些条形码，只是交付给下游的条形码范围不同。选一个熟悉的基因，逐个核对它在保留列里的计数，可以把这个问题从猜测变成一项文件检查。

以下以常规基因表达文库的 cellranger count 输出为例。多样本拆分、探针型文库以及额外做过环境核糖核酸校正的矩阵，文件层级和筛选含义可能不同。拿到外部交付时，先找分析命令、软件版本和文库编号。目录名都是“过滤后”，也有可能一个来自原始软件，另一个已经被分析人员再次质控。

分清行、列和矩阵里存的值

在 outs 下面找到 raw_feature_bc_matrix 和 filtered_feature_bc_matrix。各目录通常含压缩的特征表、条形码表和稀疏矩阵。特征表给行的身份，条形码表给列的身份，矩阵里的整数是归到某特征和某条形码的去重分子计数。两个目录应成套读取。拿原始矩阵配过滤版条形码，即使程序没有立刻报错，列号也已经失去原来的含义。10x 对特征与条形码矩阵的说明

原始矩阵保留的是满足该流程条件、有读段支持的白名单条形码，其中可能有细胞，也可能只有环境中的游离分子。过滤矩阵保留的是软件判为与细胞关联的条形码。因而不能把两者相减后直接称作“死掉的细胞数”。对同一次常规计数输出，应该先看到列集合收缩，再验证保留列上的计数能够一一对上。特征表有时还含抗体等其他类型，计算“检测到多少基因”时应先限定基因表达特征。

用四列数据看懂变化发生在哪里

设原始矩阵只有三个基因和四个条形码。为了演示，假设软件保留了第一、第三个条形码。这个设定只演示如何查文件，不模拟软件如何做细胞判定。表内的零也要读进去，一列中某个基因为零，只表示该组合没有计到分子。

甲列在过滤前是五、一、零，过滤后仍应是这三个数。第一行基因的全矩阵总数却会从七降到五，因为乙和丁各贡献的一个分子随列退出。基因没有从特征表消失，保留细胞的表达值也没有被重新归一化。若有人说“过滤使这个基因的表达下降”，先追问他比的是该基因跨所有条形码的总和，还是同一个细胞中的数。这个例子里只有前者发生了变化。

在真实文件中验证列身份

进入同一文库的 outs 目录，在装有 Python 的环境运行下面这段代码。它只读两个压缩条形码文件，打印两组列数、过滤版独有的数量，以及原始版独有的数量。如果第三项大于零，先把样本、导出的流程和条形码后缀核对清楚，再继续读矩阵。对跨文库合并数据，末尾的后缀可能承担区分文库的作用，不要为了对上数量直接删除。

import gzip
from pathlib import Path

def barcodes(folder):
    with gzip.open(Path(folder) / "barcodes.tsv.gz", "rt") as fh:
        values = [line.strip() for line in fh]
    assert len(values) == len(set(values))
    return set(values)

raw = barcodes("raw_feature_bc_matrix")
kept = barcodes("filtered_feature_bc_matrix")
print(len(raw), len(kept), len(kept - raw), len(raw - kept))

集合检查通过后，再抽几个保留条形码。先通过条形码表找到它在两个目录中各自的列号，再按稳定的特征编号对齐行、逐项比较数值。不要拿基因简称当作唯一键，多个特征可能拥有相同的显示名称。稀疏矩阵还容易造成一个视觉错觉，文件只写非零项，找不到某个“行号、列号”组合通常代表该处为零。这和特征压根没有进入参考注释，是两种需要分别记录的情况。

实跑经验｜公开 PBMC 数据的原矩阵有多少列真的有分子

我们下载了 10x 公开的一千个外周血单个核细胞示例的两份 HDF5 文件，使用 h5py 3.16.0 和稀疏矩阵逐列核对。厂商的 Cell Ranger 3.0.0 原始输出有 33,538 个特征、6,794,880 个条形码列，过滤矩阵保留相同的特征和 1,222 个条形码。在原表中按实际字符串重新定位这 1,222 列，所有基因与条形码组合的计数差异均为零。这个实测覆盖读取和核对厂商发布的计数结果，未重新从 FASTQ 进行比对。10x 公开 PBMC 样本、实验条件及许可

数量最容易看错的是过滤掉的部分。原表中有 6,793,658 个条形码未进过滤表，其中只有 233,378 个至少有一个分子。仅原表中出现的条形码总 UMI 中位数是零，最高却达到 2,372。未保留列既包括大量零列，也含值得进一步复核的非零条形码。统计原表时把“所有列”“非零列”和“判定为细胞的列”分别写开，比把六百多万列直接叫作六百多万个液滴或细胞更准确。较高的未保留计数本身依然不能证明它是一枚应该救回的细胞。

细胞数偏低，回到排名图判断下一步

在配套网页报告中打开条形码排名图。横轴按分子数为条形码排序，纵轴表示每个条形码的分子数。先看高分子数区域和长尾之间是否有清晰落差，再看边界处哪些条形码被认作细胞。曲线上同一小段可混有不同判定，图上的颜色与图例要一起读。把肉眼看到的拐点当作手工通用阈值，会漏掉本来分子含量较低的细胞类型，也可能纳入含环境分子的液滴。10x 对条形码排名图的解读

如果预期回收量和图上差距很大，挑三个小集合比较，分别取明显保留的条形码、边界附近未保留的条形码、长尾中的条形码。查看总分子数、检出基因数，并检查一组有组织背景依据的标记是否共同出现。某个高丰度基因在三组都零星出现，可能来自环境背景。边界条形码若同时出现一组连贯的谱系标记，才有理由进一步调查。预期上样数可以帮助发现问题，不能单独决定具体哪些列应补回来。若改动判定规则，应从同一份原始矩阵生成另一版结果，并保留最初的清单。

交给下一位分析人员时，最有用的是让他知道每一次“减少”发生在哪一步。写下原始条形码数、软件判定的细胞数、下游质控后剩余细胞数，并说明每个数字的文件来源。若后续只收到最后一套矩阵，上游已经排除的条形码就无从逐列检查。把两个目录和网页报告一起归档，将来有人质疑一类低分子数细胞是否被漏掉，仍然可以回到原始计数复核。

一份文库核对清楚后，可以沿用同样的脚本逐批检查。需要让多人长期复核原始矩阵并重跑同一流程时，可了解云生信的生物信息分析服务器。