RSEM 实跑，从共享读段到基因汇总，读懂期望计数和可信区间

一条读段同时匹配两个转录本，能否给两个转录本各加一？如果这样做，同一份测序证据就被重复计入。全部丢弃又会损失共享区域的信息。RSEM 处理的核心问题，是在明确模型假设后，把不能唯一定位的证据纳入表达估计，并输出基因和转录本两个层级的结果。

本篇实际准备三条人工转录本、生成一千对读段、完成转录本坐标比对，再运行两组 RSEM 定量。两组只改变转录本与基因的对应关系，观察期望计数、基因内比例和可信区间。运行中还遇到一个真实的构建问题，我们保留诊断记录并修复后重新跑通，没有把异常中间文件当成有效结果。

先明确输出数字的单位

RSEM 的期望计数可以包含小数，因为它表示根据模型分配到某个目标的读段或片段贡献。本篇使用双端数据，所以讨论的是片段。一个共享片段可以按权重影响多个候选，但它不是在每个候选中都被完整复制一次。

转录本丰度、基因丰度和基因内比例也不是同一个对象。前两者分别对应单条转录本与同一基因下的一组转录本，基因内比例则需要明确分母是哪一组成员。若分组文件错误，即使序列比对和转录本计数看起来正常，基因层级解释仍然可能失真。

这一点比背记输出列名更重要。读者在整理结果时，应先写出每一列的对象和单位，再决定怎样比较样本或者汇总基因。方法背景可参阅 RSEM 原始论文，其中多重匹配证据的处理是理解工具定位的关键。

一千个片段里，四百个没有唯一来源

人工转录本 A 和 B 都长两千碱基，前一千七百碱基完全相同，末端三百碱基不同。第三条转录本 C 独立生成，也长两千碱基。双端读长为一百，片段长度在二百二十到二百八十之间，固定随机种子生成并打乱顺序。

我们从共享区域生成四百个片段，从 A 的特有区域生成一百一十个，从 B 的特有区域生成九十个，再从 C 生成四百个。因为 A 和 B 的共享部分相同，那四百个片段不能仅靠序列被唯一分到其中一方。它们需要与特有区域的证据一起解释。

这里故意按区域选择片段，覆盖并非均匀的真实转录组抽样，也没有加入测序错误。因此，本次设计适合说明分配和汇总，不适合拿估计结果评判工具的真实表达准确率。特别是不能把人工指定的特有片段比例直接当成所有实验都满足的表达真值。

同一批转录本，准备两种基因关系

第一组把 A 和 B 放在基因 AB 下，C 单独属于基因 C。第二组不提供分组表，也不提供基因组注释，直接使用转录本序列建立参考。按这一路径的默认规则，每条序列会被看作单独的基因。

RSEM 的对应文件先写基因编号，再写转录本编号，中间用制表符分隔。文件名虽然含有转录本到基因的含义，列顺序却不能凭名称猜。复用其他分析工具的对应表之前，必须核对方向和是否包含表头。官方参考准备说明

printf 'geneAB\ttxA\ngeneAB\ttxB\ngeneC\ttxC\n' > gene_map.tsv

rsem-prepare-reference \
  --transcript-to-gene-map gene_map.tsv \
  transcripts.fa grouped/ref

rsem-prepare-reference transcripts.fa separate/ref

两组生成的转录本序列索引摘要完全相同，差别在分组关系。因此，后面看到比例变化时，可以排除参考序列本身被修改这一解释。真实项目也应采用类似思路，分别保存序列版本和注释对应关系，不能只记录一个含糊的参考名称。

输入要求比“已经有 BAM”更具体

RSEM 接收已有比对时，对记录组织有明确要求。同一读段的多个比对应放在一起，双端每个比对位置的两个 mate 还应相邻，序列和质量字段不能缺失。常见的基因组坐标排序 BAM 并不自动满足这些要求。

本篇使用 HISAT2 对转录本序列进行非剪接比对，禁用软剪切、混合配对和不一致配对，并用很高的缺口惩罚避免插入删除。在这份无错误的设计数据中，实际 CIGAR 全部是完整的一百碱基匹配。两组输入都通过 RSEM 自带验证器。

hisat2-build -p 2 grouped/ref.idx.fa grouped/hisat
hisat2 -p 2 --no-spliced-alignment --no-softclip \
  --no-mixed --no-discordant \
  --rdg 999999,999999 --rfg 999999,999999 \
  -k 20 --seed 144 -x grouped/hisat \
  -1 reads1.fq -2 reads2.fq -S grouped/alignments.sam
rsem-sam-validator grouped/alignments.sam

最终 SAM 有二千八百条记录，但只有一千个片段名称。六百个具有唯一匹配的片段各产生两条记录，四百个共享片段各产生两个位置、每个位置两条记录。记录数增加并不意味着测序量增加。这次核对同时确认了候选比对没有被误当成额外分子。

这些参数用于本次输入约束演示，并不是建议对所有 RNA 数据关闭剪接。这里比对的是已经剪接好的转录本序列，若把相同命令直接用在真实基因组上，输入对象就变了。正式流程应选择与参考构建方式匹配的比对路径，再核对 RSEM 定量输入说明。

安装时遇到的崩溃怎样处理

本次基础容器的软件包元数据为一点三点三，构建号以十二结尾，但命令行版本横幅自报一点三点一。我们同时保存了软件包信息、镜像摘要和实际横幅，没有任选其中一个当成完整的环境记录。

更重要的是，原始镜像在计算可信区间时退出码为一百三十九，中间抽样文件出现非有限数值。换用另一个较早构建也没有解决。此时即使输出目录已经有表达表，也不能把整个分析视为成功，尤其不能引用其中异常的后验丰度。

我们检查固定版本官方源码，发现有些读文件操作位于断言表达式内部。对同一份源码做控制编译，显式关闭断言能够复现崩溃，保留断言则正常完成。原镜像中的相关程序也未检出断言调用符号。这些证据支持断言移除导致运行异常的判断，但并非对所有发布版本的概括，也不等于拿到了原始构建日志。

最终从官方一点三点三源码重新编译 Gibbs 抽样和可信区间两个程序，明确保留断言，使用只读挂载替换本次运行的对应程序。源码没有修改，基础镜像也没有被覆盖。随后从头重跑两组完整定量，均成功退出并得到有限的后验估计和区间。本次核对的固定版本源码

公开脚本包含源码下载摘要检查、这两个程序的编译、输入生成和全部运行步骤，需要 Linux、Docker、Python，以及具备静态链接库的 C++ 编译环境。它会在输出目录中构建所需文件。单独拉取下面的原始镜像并不能复现本文修复后的可信区间结果，完整复现应执行附件中的脚本。

docker pull quay.io/biocontainers/rsem:1.3.3--pl5321h077b44d_12

# 解压复现包后，在对应工具目录执行
# 脚本会按固定源码重建 Gibbs 和可信区间程序
python3 run.py ./practice

共享片段被怎样分配

修复环境后，分组参考下 A 的期望计数为三百二十九点九九，B 为二百七十点零一，C 为四百。A 与 B 的和为六百，基因 AB 的输出也为六百。这说明四百个共享片段被纳入了分配，没有向 A 和 B 各自重复增加四百份。

这组简单设计中，两条相同长度转录本的特有证据分别为一百一十和九十，共享证据的分配方向与它们提供的信息一致。读者可以用五十五与四十五的比例理解本次结果，但不能把这个简化算式当成 RSEM 在所有数据上的完整算法。真实输入还涉及片段长度、有效长度、匹配质量和模型估计。

小数计数不是文件格式错误，也不应该为了看起来像测序条数就直接四舍五入。它保留了模型分配的信息。进入下游分析时，应使用该方法支持的导入方式，并明确当前处理的是期望计数还是相对丰度。把两类数值混在一起，往往比某一个参数选择更难补救。

为什么不提供分组表后全部变成百分之百

第二组中，A、B、C 的期望计数和 TPM 点估计与第一组相同，但三条转录本的基因内比例全部变为百分之百。原因很直接，每条转录本现在都被当成一个独立基因，它在自己的集合里就是唯一成员。

因此，基因内比例从五十五变成一百，不能解释成该异构体表达翻倍，更不能解释成剪接偏好改变。这次读段根本没有变化，改变的只是分母对应的成员集合。比例列需要连同基因编号、转录本编号和注释版本一起阅读。

实际项目合并不同来源的参考时，也容易发生类似问题。某些条目被错误拆成独立基因，另一些条目仍然按正确关系归组，最后会得到看似合理却不可比较的比例。只查看计数是否为正值无法发现这一类错误，应抽查一个明确具有多个异构体的基因，确认它们在结果中确实属于同一组。

基因可信区间不能相加转录本端点得到

分组参考下，基因 AB 的 TPM 点估计为六十万，其默认百分之九十五可信区间约为五十六点九万到六十二点九万。A 的区间约为二十八点六万到三十七点五万，B 约为二十二点八万到三十一点五万。

把 A 与 B 的下限直接相加得到五十一万三千四百五十，上限相加得到六十九万零三百零七，都不是程序输出的基因区间。共享片段在 A 与 B 之间发生的分配波动具有联系，一个目标多分一些，另一个就可能少分一些。汇总时应考虑联合变化，而不是把各自最极端的端点机械拼接。

本次 A 与 B 的后验计数标准差都约为十七点六六，基因 AB 的后验计数标准差为零，但基因 TPM 仍然有非零宽度的区间。这也提醒读者，计数分配的波动与归一化丰度的不确定性并非同一列。不能看到某个计数标准差为零，就推断所有表达参数都确定无疑。

这些区间来自给定模型和样本内证据，不代表生物学重复之间的差异范围。增加抽样次数可以帮助计算更稳定，却不能凭空补上缺失的参考结构或生物学重复。公开结果记录了随机种子、线程和默认抽样设置，多线程抽样及执行环境可能带来末位差异，复现时应先核对主要关系和结果有效性。

适合怎样的问题，有哪些限制

RSEM 的优势是把共享证据纳入概率分配，同时提供基因和转录本层级结果，便于分析异构体和同源序列带来的归属不确定性。明确的参考准备步骤、输入验证器和多种输出，也使流程能够逐层检查。对于需要理解“这份计数是怎样得到的”的项目，它提供了比简单重叠计数更丰富的信息。

限制同样来自模型和参考。候选转录本不完整、错误分组、缺少区分序列或者不合适的比对输入，都可能影响解释。点估计看起来稳定，不等于真实来源已经被唯一确定。复杂度增加时还需要实际评估比对文件、抽样和存储开销，本次小实验没有测量大规模性能，不能据此给出通用配置结论。

这次实跑的经验可以归结为一条可执行的检查顺序。先核对参考成员与分组，再检查候选比对和配对组织，随后确认程序完整退出及输出数值有效，最后才解释计数、比例和区间。这个顺序既能发现把转录本误当成基因的分析错误，也能避免把软件构建故障留下的半成品带入下游结论。

下载本组复现脚本和实际输出，其中包含基础镜像身份、修复构建信息、异常对照证据与两组完整结果，方便读者逐项核对。

当参考体系和定量流程已经核验，需要为多样本项目保留可复现环境和中间结果时，可了解云生信生物信息学分析服务器服务。