长读长转录组里一条读段从一个外显子跨到另一个外显子,浏览器有时画成跨内含子的长线,有时只给出两个单独的比对块。两块都严丝合缝地落在正确的位置,比对质量也很高。这时是出现了新的融合,还是把转录本按基因组读段比了一遍?可以先设计一个答案已知的剪接位点,把同一条序列只改比对参数送入 minimap2,看看两种显示各需要什么条件。这个练习也能说明,做融合或新异构体筛选之前,为什么应该先还原一条分子的全部比对。minimap2 作者的用法说明
造一个跨六千碱基内含子的分子
我们固定随机种子生成一条人工染色体,两个外显子分别长四百八十和五百二十碱基。中间插入六千碱基的内含子,正向参考上以 GT 开头、以 AG 结尾,两端再放入不参与转录本的随机侧翼。直接拼接两个外显子得到一千碱基的成熟分子,另外生成它的反向互补,分别赋予不同的读名。两个输入代表同一条设计序列的两种测序方向,可以用来检查方向翻转会改变哪些列。数据没有人为添加测序错误,也没有第二处相似基因,便于把差异归因到本轮比较的参数。
一基闭区间下,第一段外显子是第七百零一位至第一千一百八十位,第二段从第七千一百八十一位至第七千七百位。将一千碱基的查询比回基因组后,我们预期能看见两个外显子边界,也能查出一条记录到底用掉多少查询碱基。在三个运行目录保留相同的参考、读段名称和输入顺序,用同一版 minimap2 生成 SAM 并写出碱基差异标签。从文件里逐列读取起点、标志位、比对质量、操作串和编辑距离,避免只凭浏览器缩放后的画法判断。
先看同一读名有几行,再决定如何计算覆盖
普通基因组长读段预设下,正向分子产生了两条记录。主比对覆盖较长的第二个外显子,操作串前面软剪切四百八十个碱基,剩下五百二十个碱基全部匹配。补充比对落在第一个外显子,覆盖另外四百八十个碱基,后面硬剪切五百二十个。两行编辑距离都为零,比对质量均为六十。分别核对两段在原查询上的位置,恰好互补。因此按读名还原时整条查询都有去处,只保留主比对再计算每条分子的比对比例,会得到百分之五十二。两者回答的问题不同,尤其不能把两条记录当成两个细胞中的两条转录本。
换成剪接预设,同一个正向读名只剩一条主比对。操作串先匹配四百八十,接着在参考上跳过六千,再匹配五百二十。该行消耗一千个查询碱基,参考跨度有七千个位置。查询序列中没有那六千个碱基,测序深度也不能从跨接线的长度推算。一张表达报告如果要展示外显子覆盖,就要按实际有碱基对齐的区块算。一张剪接报告如果要展示内含子支持,则可以从跨接操作提取供体与受体边界,再以独立分子为单位汇总。
| 运行方式 | 每个读名的记录数 | 主记录 | 补充记录 |
|---|---|---|---|
基因组 map-ont |
2 | 480S520M |
480M520H |
剪接 splice |
1 | 480M6000N520M |
无 |
剪接但限定 -G 1000 |
2 | 480S520M |
480M520H |

两个外显子的支持也不能仅用基因总读数代替。只落在第二个外显子里的另一条读段会增加它的覆盖,却没有穿过接头。做结果表时,在接头两侧分别设置最小锚定长度,并注明同一分子在同一接头只计一次,这样有效支持数才可以与覆盖度图对照。
实跑经验|一个最大间隔就能让剪接连接消失
minimap2 -t 1 -a --cs=short -x map-ont ref.fa reads.fa > genomic.sam
minimap2 -t 1 -a --cs=short -x splice ref.fa reads.fa > splice.sam
minimap2 -t 1 -a --cs=short -x splice -G 1000 ref.fa reads.fa > limited.sam
samtools view splice.sam | cut -f1-6
第三轮仍然启用剪接模式,随后把允许的最大参考间隔改为一千。真实设计的间隔是六千,结果又拆成两段,正反两个读名得到相同的记录数变化。三组运行中,出现在结果里的比对质量始终是六十。这个分数高,可以和“两个外显子是否被同一条比对串连接起来”同时检查,但它本身没有给出生物学结构的答案。真核物种的内含子长度分布不同,选上限时可以先从可靠注释提取长度分布,再挑短、中、长几个已知接头做测试。把参数设置保存在运行日志里,比事后只拿一份已排序的文件解释漏检容易得多。minimap2 参数手册
反向互补读段还给了一个独立的检查。它比回同一条染色体,跨度与外显子边界都一致,记录里标记了反向比对。剪接差异串在两条记录中都显示相同的内含子长度和参考上的供受体序列。注释只在同一个链方向上定义这两个外显子,不能因为查询以反向互补形式输入就另外发明一个反义转录本。实际实验中还要查建库方向、读段是否经过定向处理,再按软件定义解读相应标签。对于没有方向信息的文库,可以把参考上的接头和文库链特异性分别保存,避免后续合并表时丢失判断条件。
如果候选融合只剩两段分开的高分比对,我会先按读名提取主比对与补充比对,检查查询上的片段是否重叠、参考上的顺序是否符合普通剪接、间隔是否刚好超过运行参数,再寻找其他独立读段的支持。人工数据里的每个边界都精确命中,真实样本却可能遇到短外显子、重复序列或非典型接头。把已知阳性位点与这些难点分开记录,能判断调整参数究竟找回了原有连接,还是额外引入了一批需要复核的候选。
需要在多个组织或样本中沿用这批已验证的剪接检查时,可以用云生信的生物信息分析服务器集中保存比对版本及原始读段。
