一张融合结果表里,支持数增加通常让人更有信心。但增加的究竟是新的独立片段,还是原有片段的重复副本,必须分清。本次在同一段人工融合序列上增加四十对重复读段,实际运行比对和融合检测后,有效支持数没有增加,断点附近的一侧覆盖度却明显升高。如果只根据覆盖深度排序,就会把这次人为重复当成更强的证据。
这篇文章用一个可以完整重建的小实验,说明如何从原始读段追到比对日志,再追到融合结果中的计数与过滤原因。它解决的是结果审阅问题。实验里的基因、参考序列和读段全部由脚本生成,不涉及真实样本,也不能据此估算真实项目的检出率或假阳性率。

先明确 Arriba 接收的证据
Arriba 从比对结果中寻找异常连接。真正跨越连接位置的分裂比对,和一对读段分别落在两个区域的跨接证据,提供的信息精度不同。前者能帮助定位连接边界,后者提供片段层面的空间约束。审阅时应分别检查,不宜把所有支持压成一个数字后就停止追踪。
尤其要确认传给程序的比对文件确实含有嵌合比对。普通表达定量流程产生的文件,即使文件名相同,也不保证具有融合检测需要的信息。官方输入说明区分了把嵌合记录写入主文件与另存文件的方式,本次使用前一种。相关要求见 Arriba 输入文件说明。
实践中可以先保存比对命令和日志,再保存融合命令。这样,当结果为空时,可以判断问题发生在候选生成之前,还是发生在后续过滤阶段。仅留下最终表格,会让两种完全不同的失败看起来一样。对于需要复查的候选,读段名称还能帮助回到原始数据,但名称列表本身未必等同于过滤后的有效数量。
实验怎样构造
我们安装并运行了 Arriba 2.5.1,配合 STAR 2.7.11b。参考基因组只有两条各一万碱基的人工染色体,每条各放置一个双外显子转录本。融合转录本取甲基因的第一个外显子,接上乙基因的第二个外显子,预设连接坐标为 chr1:1600 与 chr2:4001。
脚本生成一百六十对正常转录本读段和四十一对融合读段。每端长度为一百五十碱基,片段长度为三百五十碱基,不加入测序错误。融合片段沿连接位置错开起点,使它们具有不同的覆盖方式。第二份输入在基础数据上额外加入四十对副本,这些副本来自同一条融合片段,只改名称,不改序列及片段位置。
这种设计把变量限定得很清楚。正常背景、参考文件、比对参数和原来的融合片段全部保留,新增内容只是重复。它不模拟真实实验中所有重复来源,也不讨论重复一定来自扩增还是高表达。我们要验证的只是,在这套明确输入下,程序如何计数以及输出哪些诊断信息。
本次没有使用人类黑名单、已知融合库、蛋白结构域或基因组结构变异。人工参考与那些资源没有对应关系,因此实验关闭黑名单过滤,其余过滤保持默认,另外安排一组专门关闭重复过滤的对照。这一设置只适用于解释本实验,不能照搬为真实人类样本的默认方案。
安装成功要以可执行和版本输出为准
首次检查 STAR 时遇到一个容易忽略的问题。包内名为 STAR 的文件是启动包装器,单独取出它会因找不到同目录的实际程序而失败。本次改为提取包内通用版本二进制,并在记录中注明实际成员。随后在固定镜像提供的依赖环境中运行版本查询、索引生成、比对及融合检测,全部成功。
这一经验适用于手动提取软件包的场景。下载完成不等于安装完成,看到可执行文件也不等于依赖完整。检查应至少包含版本输出和一次能够产生预期结构文件的小任务。直接把一个文件复制到系统目录、再根据文件名判断版本,容易遗漏包装器、动态库和伴随数据之间的关系。
本次的下载地址、压缩包校验值、二进制校验值和运行镜像均保存在结果记录中。实验使用两个计算线程和四吉字节的容器上限。这只是隔离练习的运行条件,不能用于推算真实人类数据的容量需求。小参考上近乎瞬时完成的步骤,也没有作为性能测试发表。
比对日志先回答第一层问题
基础输入共有二百零一对读段,其中人为构造的融合片段为四十一对。STAR 日志报告的嵌合读段数为三十七,Arriba 随后读取到的嵌合记录也是三十七。最终保留的目标融合恰好有三十七个有效支持片段。
这个过程提示我们,知道生成了多少融合片段,并不意味着下游一定接收到同样数量的融合证据。比对对锚定长度、位置和嵌合报告方式有自己的限制。这里没有把少掉的四个片段归咎于某一个未经逐条验证的原因,也没有写成 Arriba 丢掉了四个真阳性,因为日志已经显示差异发生在它接收候选之前。
第二份输入共有二百四十一对读段。Arriba 读取的嵌合记录增加到七十七,执行重复过滤后又回到三十七。这个中间日志把因果链补齐了。若只比较最终两张表,只能看到计数相同。连同比对及过滤日志一起查看,才能确认新增副本确实进入了流程,又被哪个步骤识别。
同一候选的三组结果
三组均保留一条甲基因到乙基因的目标融合,断点与预设位置一致。下面只摘取解释重复效应所需的字段。两类分裂支持分别对应结果中的 split_reads1 和 split_reads2,跨接支持对应 discordant_mates。
| 运行条件 | 分裂支持甲/乙 | 跨接支持 | 有效支持合计 |
|---|---|---|---|
| 基础输入 | 12/17 | 8 | 37 |
| 增加副本并过滤 | 12/17 | 8 | 37 |
| 增副本不去重 | 12/57 | 8 | 77 |
关闭重复过滤后,增加的四十个支持全部进入乙侧分裂支持栏。断点没有因此变得更精确,参考序列也没有增加新信息。这是一个直接的反例,说明支持数变化必须结合预处理和计数规则解释。同一个阈值用于两套过滤规则不同的结果时,含义可能已经改变。
| 运行条件 | 甲侧覆盖度 | 乙侧覆盖度 | 过滤栏 |
|---|---|---|---|
| 基础输入 | 47 | 54 | 无 |
| 增加副本并过滤 | 47 | 94 | 重复副本 40 |
| 增加副本但不去重 | 47 | 94 | 无 |
覆盖度与有效支持数遵循不同的计算口径,官方也提醒覆盖度计算包含重复片段。字段定义可查 Arriba 输出说明。在本次对照中,只有乙侧覆盖度增加,不能进一步推广成每个重复片段都会同时增加两侧覆盖度。应根据实际输出和比对位置解释具体变化。
由此还可以检验一种常见但不严谨的做法。如果将有效支持除以两侧覆盖度之和,基础输入得到约零点三六六,加入副本并过滤后约为零点二六二。融合序列并未改变,比值却下降了。这个比值不能未经定义就命名为肿瘤细胞比例、融合等位基因频率或真实转录本占比,因为分子和分母并不是同一种计数对象。
名称列表也需要核对口径
更容易漏看的现象出现在读段名称栏。基础输入的名称列表有三十七项,增加副本并执行过滤后,列表仍包含七十七个名称,而有效支持数保持三十七。也就是说,本次输出中的名称列表保留了被标记重复的相关读段,并不是一份已经排除全部过滤证据的独立片段清单。
因此,写统计脚本时不要简单对逗号分隔后的名称计数,再把结果覆盖到支持数字段。那样会把程序已经处理过的副本重新算进去。名称更适合作为回查入口。要统计最终证据,应先明确所用字段与过滤状态,并在少量候选上做人工对照,确认自编代码没有悄悄改变口径。
这个核对还可以帮助发现另一个问题。最终候选只有一行,不代表中间只生成过一个候选。本次日志一度记录九个候选,经过相对支持过滤后剩下一条。被丢弃表保留了审查入口。整理结果时同时保留主表与被丢弃表,比只复制通过的那一行更有利于解释为什么某个邻近断点没有进入最终结果。
高置信度和缺失字段分别意味着什么
三组目标融合都标为高置信度,关闭重复过滤并没有让这个标签发生变化。标签是离散分级,不是每增加一个读段就连续变化的数值。也不能把高置信度理解为有特定致病作用。官方明确区分候选可信程度与肿瘤生物学意义,参见 Arriba 结果解释。
本次注释没有提供编码区和结构域信息,因此阅读框、肽序列和结构域相关字段没有给出有效预测。面对这样的空值,正确结论是当前输入与运行没有提供相应判断依据。它不能改写为不编码蛋白、没有功能结构域,或一定不能形成融合蛋白。空值的解释应当回到输入条件,而不是仅凭表格中的一个点号下结论。
同理,两端落在外显子边界只能说明与本次注释中的位置吻合。若要进一步研究功能,需要具体转录本、编码相位及连接后序列等信息。若要验证真实样本中的事件,还需要独立证据。把检测、注释和生物学解释分开,能让结果报告更准确,也便于后续知道下一步应补什么数据。
怎样复现并迁移到自己的流程
完整脚本会生成参考与读段、安装所需可执行程序、建立索引、完成两份输入的比对并执行三组 Arriba 对照。下载 本组实跑脚本与结果记录 后,在有 Docker 的 Linux 环境运行对应目录中的 run.py 驱动脚本即可。压缩包中的记录包含实际输出,而不是只有预期答案。
python3 run.py run156
实验里 Arriba 的关键对照命令如下,输入文件由驱动脚本生成。真实数据应先按自己的物种、参考和注释准备资源,不要直接沿用这里关闭黑名单的教学设置。
./arriba -x duplicates/Aligned.out.bam \
-a genome.fa -g genes.gtf \
-o duplicates_filtered.tsv \
-O duplicates_filtered.discarded.tsv -f blacklist
./arriba -x duplicates/Aligned.out.bam \
-a genome.fa -g genes.gtf \
-o duplicates_keep_duplicates.tsv \
-O duplicates_keep_duplicates.discarded.tsv \
-f blacklist,duplicates
迁移流程时,建议先用一份已知来源的小输入检查字段及计数关系,再扩大数据量。检查目标应是结果能否追溯,而不只是程序是否返回成功。至少要能回答输入含多少片段、上游报告多少嵌合证据、下游保留多少有效支持、哪些过滤解释了差额,以及最终候选对应哪一版参考与注释。
本次还出现了无法用少量嵌合读段估计片段间隔分布的提示,程序使用了默认值。我们保留了这个提示,没有因为任务完成就把它忽略。它提醒读者,小实验可以验证字段和控制变量,却不能覆盖大样本中的全部估计过程。将这类限制写在结果旁边,比把一条成功命令包装成通用参数更有用。
优点、局限与实际选用
Arriba 的实用优点是能够把候选、支持类型和过滤原因放在可追溯的输出中。本次重复对照正是靠这些信息解释清楚的。它适合已有合适比对文件、又需要逐条审阅融合候选的流程。官方提供的完整流程可以作为配置入口,见 固定版本运行脚本。
它的局限也很具体。上游未报告的证据无法凭空恢复,参考与注释会影响定位和解释,人工参考中的漂亮结果不能代表真实数据表现。过滤器提供的是处理规则与优先级,不能替代对实验设计、样本来源和异常模式的理解。对候选有疑问时,应追查证据链,而不是不断关闭过滤直到出现想要的一行。
一次有价值的融合分析,应当能够解释支持为什么增加、为什么减少,以及哪些信息尚未确定。本次最值得带走的习惯,是把有效支持、覆盖度、名称列表和过滤记录同时核对。四个数字或字段都看起来合理,并不代表它们能互相替换,真正可靠的解释来自口径一致和输入可追溯。
当这套证据核对流程需要扩展到多批次样本时,可按实际分析任务了解 云生信生物信息学分析服务器服务。
