alevin-fry 的三个计数层怎么选？把同一基因的 6 和 9 算清楚

同事把一个单细胞定量目录读进 Python，你把同一个目录读进另一份笔记本。细胞条形码相同，基因名也相同，有一个基因却分别显示六个和九个分子。先不要重新比对原始读段。用 alevin-fry 的 USA 模式定量时，同一基因可以有三种剪接状态的计数。读入函数按下游用途组合它们，两个矩阵都叫主表达矩阵，里面装的却可以是不同的分子集合。

这个问题特别容易出现在整细胞和单核测序合并分析时。单核样本通常能捕获较多尚未完全剪接的转录本，沿用整细胞项目的读入方式，会改变这部分证据是否进入表达量。若研究的是转录动态，则可能需要分别保留不同剪接状态。同一个计数目录留有三层信息，正是为了让这些问题都有机会回到原始证据上。alevin-fry 对定量结果和导入方式的说明

先辨认三种计数各自表达什么

可以把三个状态想成一个基因下的三只抽屉。已剪接证据进入一只，未剪接证据进入一只，还有一只用来保存无法在这两种状态之间可靠区分的计数。第三种的“模糊”描述的是剪接状态。一个分子可以有证据归到某个基因，却缺乏足够信息判断它来自哪种剪接形式。看到模糊这一栏，不能直接把它统统解释成别的基因的错配，更不能默认应该丢掉。

这里说的是聚合后的计数，来源于所选参考、比对或映射、分子去重和歧义判定过程。换一个扩展参考，三个抽屉的分配可能跟着变化。比较两次统计前先核对是否用了同一个定量目录和同一批条形码，再检查读入函数的输出格式。若这些条件已经不同，加减分子数通常无法完全解释差别。作者关于扩展转录组参考的说明

亲手算一次，为什么主矩阵会差三个分子

设一个细胞中的某个基因有四个已剪接分子、三个未剪接分子、两个状态模糊的分子。下面的计算专为读懂文件而设，数字不代表哪个组织的典型剪接比例。先在纸上把三类单独写清楚，再决定上游计数保持不动时，怎样组合成你的观察量。

六和九相差的三个，恰好是未剪接计数。如果软件界面同时显示主矩阵和一个已剪接层，也不要想当然地把这两份相加。有的输出预设会让主矩阵和某一层保存同一组合，它们是为了方便不同函数读取而提供的两个入口。把两个入口当成独立观测相加，会把部分分子重复计算。保存对象时，把每个层的中文定义也写进分析记录，常常比只留一个简称更有用。

从实际目录中抽一条记录核账

下面用 pyroe 同一版的两种预设读入一个 USA 定量结果。路径应指向 alevin-fry quant 的输出目录。例子先对齐两个对象的细胞和基因名称，再查同一个坐标。如果想查看一个特定基因，可把下标改成该基因的唯一编号，并选择一个确实有未剪接计数的细胞。若选中的坐标刚好全为零，这个例子验证不了两种定义之间的差别。

from pyroe import load_fry

one = load_fry("quant_dir", output_format="scRNA")
two = load_fry("quant_dir", output_format="snRNA")
assert one.obs_names.equals(two.obs_names)
assert one.var_names.equals(two.var_names)
i, j = 0, 0
print(one.X[i, j], one.layers["unspliced"][i, j], two.X[i, j])

按这里所用预设，打印出来的第一个值与第二个值相加，应等于第三个值。若等式失败，先检查目录是否确实来自 USA 模式，然后确认有没有在其中一份对象上提前做归一化、取对数或过滤基因。原始分子计数可以相加，归一化或经过非线性变换后的两个数一般不能这样核账。若需要保留三个抽屉供他人使用，pyroe 还允许通过字典指定主矩阵和附加层分别取哪几类，读入后应随手查看真实存在的层名。pyroe 对预设、计数层及自定义字典的逐项说明

还有一种容易误判的情况。同一个基因符号可能对应多个稳定编号，而不同参考版本又可能更新符号。检查单个坐标时，把细胞条形码、基因稳定编号和参考版本三者一起抄下，不要只截取屏幕上的基因简称。若两个对象的基因集合相同但排序不同，应先重排到共同顺序再做矩阵相减。直接按位置相减会制造出一大片本来不存在的剪接差异。

实跑经验｜六列输入怎样变成两个基因

在 Linux 的独立环境安装 pyroe 0.9.3，固定 NumPy 1.26.4、pandas 1.5.3，用真实的 load_fry 读取自行构造的最小 USA 格式目录。目录含两枚条形码、两个基因和六列计数。第一枚细胞第一个基因设置为已剪接四、未剪接三、模糊二。用常规单细胞预设得到六，切换细胞核预设得到九。第二个基因给出了更醒目的反例，它只含七个未剪接分子，在常规主矩阵恰好为零，切到细胞核预设立刻变为七。两份对象的基因名和条形码均一一相同。

再打开速度预设，主矩阵与 spliced 层逐格一致，第一枚细胞的 unspliced 层为三和七。在四个细胞与基因组合上，我们都执行了“常规主矩阵加未剪接层等于细胞核主矩阵”的断言。初次导入时旧依赖还报缺少 pkg_resources，补装 setuptools 80.9.0 后成功。这个试验实际运行的是定量输出的读取与聚合，没有从原始读段重新比对。排查不同报告里的基因检出率前，可以取同一目录做这四格小测试，先确认差异是否仅由读入定义产生。

从一个基因的差别，追到整批样本的偏移

拿到多样本项目时，建议从每个样本挑几类表达特点不同的基因，再多查几个细胞。一个未剪接计数很多的基因，在两种预设下会出现明显差别。另一个几乎全落在已剪接层的基因，数值则可能完全相同。若你只检查了后者，容易误以为两种导入方式可以互换。把各基因在两种口径下的差额与其未剪接计数并排，就能判断眼前的变化主要来自读入选择，还是还叠加了其他处理。

对照组来自整细胞、处理组来自单核时，这一点尤其重要。即使把两组强行读成相同的主矩阵定义，样本制备造成的差异仍然存在，而且可能与处理条件完全重合。此时应先写明数据收集方式，单独查看各批次和各细胞类型，再讨论哪些对照具有解释空间。如果目的只是制作表达图谱，可以给出统一口径的主要结果，并在代表性样本上另算一版，看看细胞类型的标志是否稳定。如果目的是推断动态，则要把相应状态层按模型需要交付，并另行验证速度方向与实际生物过程能否对应。

下次交付时，除了结果矩阵，把定量模式、参考版本、读入格式、各层定义和一个能复算的细胞—基因坐标放进说明。接手人看到六和九，就知道该查哪三个数。若这些设置要跨批次重复使用，而且课题组成员希望共用一套可复现的 Python 环境，可了解云生信的生物信息分析服务器。