StringTie 实跑，参考引导组装与只按注释定量究竟有什么不同

同一个 RNA-seq 比对文件，一次分析报告两个转录本，另一次只报告一个，是否说明其中一次丢失了数据？对于 StringTie，首先应检查分析任务本身。使用参考注释引导组装，与限制在给定注释中估计表达，是两种不同的问题。后者没有输出新结构，并不能证明样本不存在这种结构。

本篇实际安装 StringTie，构造一个包含外显子跳跃的人工转录区域，生成读段并完成剪接比对、坐标排序和三组分析。我们只改变是否限制参考模型，以及参考中包含哪些转录本，观察结构和丰度怎样变化。所有结果来自实际运行，人工设计用于解释参数，不用于宣称真实样本中的发现率或定量准确率。

先把研究问题写成两句话

如果问题是样本里可能存在哪些转录结构，应允许程序结合比对证据和参考注释组装候选模型。参考提供已有结构信息，读段仍然可能支持参考之外的连接方式。如果问题是这些既定转录本在各样本中的表达如何，则需要固定模型集合，在同一套坐标与标识下估计。

StringTie 的 -G 提供参考注释。单独使用时，它可以帮助组装已知和新结构。加入 -e 后，程序进入针对给定参考的表达估计模式，需要同时提供参考注释。这里的差别影响候选模型范围，不能仅把它理解成提高速度的开关。StringTie 官方运行说明

这种区分在多样本项目中很重要。各样本先分别组装，可能得到不一致的边界和临时编号。如果直接把同名临时转录本放在同一行，就可能把不同结构拼成一张矩阵。比较表达之前，应先建立适合研究目的的统一模型集合，再让各样本在相同集合上定量。

设计一个能逐段核对的外显子跳跃

我们生成一条六千碱基的人工染色体，在正链设置三个外显子。已知转录本包含三个外显子，另一条转录本跳过中间外显子，只连接第一段与第三段。两个转录本的首尾坐标相同，内部结构不同，正好说明为什么不能只按起止位置判断转录本是否一致。

两者的剪接后序列分别长一千和八百碱基。

初始参考只写入三外显子的已知模型。另准备一份完整参考，同时写入已知模型和跳跃模型。注意，完整参考之所以知道跳跃结构，是因为这是我们设计的对照，不代表真实项目可以提前知道所有新转录本。

根据两条剪接后的序列，脚本以固定间隔取一百五十碱基的单端读段，每个位置生成三份，然后打乱顺序。已知模型贡献一千二百七十八条，跳跃模型贡献九百七十八条，总计二千二百五十六条。序列无测序错误，质量值一致，覆盖也受人为取样方式影响。这些条件便于核对结构，不能代替有批次差异和真实噪声的生物学重复。

安装后先确认输入坐标体系

本次 StringTie 实际版本为三点零点三，HISAT2 为二点二点三，排序使用 Samtools 一点二十二点一。复现脚本锁定三个容器摘要，并记录参考序列摘要、读段来源和实际命令。读者可以从公开附件生成同样的输入，再完整执行三组对照。

docker pull quay.io/biocontainers/stringtie:3.0.3--h29c0135_0
docker run --rm quay.io/biocontainers/stringtie:3.0.3--h29c0135_0 \
  stringtie --version

# 解压复现包后，在对应工具目录执行
python3 run.py ./practice

StringTie 使用按基因组坐标排列的 RNA 比对结果。本次先用 HISAT2 将读段比对到人工染色体，再用 Samtools 排序。不能拿到一个文件扩展名为 BAM 的结果就直接认为适用，转录本坐标上的比对与基因组坐标上的剪接比对回答的是不同问题。参考序列名和注释坐标也必须匹配。

hisat2_extract_splice_sites.py known.gtf > known.ss
hisat2-build -p 2 genome.fa genome
hisat2 -p 2 --rna-strandness F \
  --known-splicesite-infile known.ss --seed 146 \
  -x genome -U reads.fq -S reads.sam
samtools sort -@ 1 -o reads.bam reads.sam
samtools quickcheck -v reads.bam

上面命令中的链方向来自人工读段的生成方式，不是所有实验都应使用的默认答案。真实项目应核对文库方向，以及剪接比对中有关转录方向的标记。排序检查能够排除一部分文件问题，但无法证明读段已经正确支持了每一个候选结构。

三组分析只改变模型范围

三组 StringTie 命令使用同一份排序文件，均保留端部位置、不按覆盖衰减修剪，以便比较设计边界。第一组用旧参考引导组装，第二组在旧参考上限制定量，第三组在包含两条转录本的完整参考上限制定量。

stringtie reads.bam -p 2 --fr -t \
  -G known.gtf -o guided.gtf -A guided.genes.tsv

stringtie reads.bam -p 2 --fr -t -e \
  -G known.gtf -o known_only.gtf -A known_only.genes.tsv

stringtie reads.bam -p 2 --fr -t -e \
  -G union.gtf -o union_only.gtf -A union_only.genes.tsv

这里关闭端部修剪是为了减少本次对照中的额外变量，不应推广成所有样本都必须照搬的参数。对于覆盖渐降、降解明显或者端部证据不足的数据，边界本来就需要结合读段分布和其他证据判断。本篇只检验外显子连接方式和参考集合对输出的影响。

第一组输出的两条模型，外显子坐标分别与设计的三外显子和两外显子结构一致。第二组只留下已知模型。第三组恢复两条模型，但它们作为参考中的条目被估计。三个结果不是简单的正确与错误之分，而是不同任务边界的直接表现。

判断一个新结构，要看连接证据

此次两千多条读段全部成功且唯一比对，但这还不够。解析实际比对中的跳跃操作后，已知结构的两个内含子分别得到二百三十五次和二百三十六次支持，跳跃结构跨越的长内含子得到一百九十五次支持。这些数字说明，跳跃连接有独立的剪接比对证据，而不是只因为一个输出文件里多了一行。

对应的内含子闭区间分别为一千四百零一到二千五百、二千七百零一到四千，以及一千四百零一到四千。最后一段跨过了原来的中间外显子。检查这种关系时，最好同时写出前一外显子末端、后一外显子起点和中间间隔，避免把内含子长度与坐标端点混在一起。

实际比对还出现二十七个软剪切操作。即使总体比对率达到百分之百，也不能据此说每条读段都完整支持了预期剪接。本文保留这一现象，并没有将比对率包装成结构正确率。真实项目应检查跨连接读段的两侧锚定长度、错配、重复区域及异常覆盖，而不仅是浏览器中看起来存在一条连线。

人工对照中我们知道生成真值，可以确认输出结构与设计一致。真实样本中的新结构则需要更多证据，尤其是低覆盖连接、单一样本独有结构和复杂同源区域。候选结构可以成为后续验证对象，但不能只凭组装输出就宣布发现了稳定表达的新异构体。

为什么旧参考定量得到一百万 TPM

旧参考限制定量组只输出一条转录本，其 TPM 为一百万。这个数字很容易误导读者，以为所有读段都来自该转录本，或者它在真实样本中具有极端高表达。实际上，这次输出集合只有一个可报告的模型，相对丰度归一化的分母也受模型集合约束。

允许两个模型后，已知转录本的 TPM 约为六十点八万，跳跃模型约为三十九点二万。引导组装组分别为六十万七千九百五十三点一八和三十九万二千零四十六点八二，完整参考定量组分别为六十万七千八百五十一点八一和三十九万二千一百四十八点一九。两组接近，但并非逐位相同。

这些差异不应该被直接解释为表达发生了变化，因为输入比对文件完全相同。改变的是模型范围和估计过程。旧参考组将可解释的共享区域证据放在更受限的模型中，结果也不会告诉你参考之外的结构不存在。若一个样本用旧参考、另一个样本用扩充参考，丰度表面上的差异可能首先反映分母和分配方式改变。

还应避免拿人工生成的读段份数直接检验 TPM 比例。两条转录本长度不同，取样位置的数量不同，边缘覆盖和共享区域也不同。读段来源比例、平均覆盖和按有效信息估计的相对丰度，各有自己的统计含义。本次设计适合观察参数造成的变化，并不是经过校准的表达真值基准。

临时编号相同，不代表跨样本结构相同

引导组装组为两条输出生成临时编号，并给已知结构附上参考对应信息。限制定量组则使用参考中的转录本编号。整理结果时，应保存参考对应属性和完整外显子链，而不是只保留一个看起来简洁的名称。

基因丰度表也需要仔细看。此次引导组装输出同时包含参考基因编号与组装基因编号，它们位于同一条染色体、同一链和同一范围。若仅把不同编号当成不同基因相加，就会把一个区域重复解释成两个独立对象。这个现象提醒我们先核对记录之间的关系，再汇总数量或绘图。

批量分析时，统一参考不仅是为了减少表格中的空值，更是为了让同一行在不同样本中代表同一种结构。需要区分已知结构、新连接、端部变化和无法明确对应的模型。过度合并可能掩盖真实差异，完全不合并又可能把边界噪声放大成大量新条目。保留原始组装结果和统一模型之间的对应表，能够让后续解释可追溯。

组装、合并、再定量分别解决什么

一个常见的多样本思路，是先在各样本中获得候选结构，经过质量检查后建立统一集合，再对各样本重新定量。StringTie 的合并模式处理的是转录本注释结果，而不是把多个 BAM 文件拼在一起重新比对。理解输入对象，才能避免把合并参数放到错误阶段。

统一集合应包含哪些候选，取决于研究目标和证据要求。参考充分、主要关注已知基因表达时，可以采用更受限的模型范围。关注特定剪接变化时，需要追查支持该连接的证据和样本间重复性。研究新转录区域时，还需要特别关注注释不完整、读穿转录以及技术噪声带来的假候选。

本次没有运行多样本合并，也没有执行差异表达分析，因此不提供这些环节的实测性能或发现数量。读者迁移流程时，应按官方说明核对相应输出与下游工具的输入要求。GTF 中的相对丰度不能直接冒充原始计数，单个样本的两条模型也不构成有统计重复的差异剪接结论。

优点、限制与实跑经验

StringTie 的优势是把基因组剪接比对、参考信息和转录本结构联系起来，可以同时观察模型形态和表达估计。参考引导与限制定量两种模式，使它既能用于探索候选结构，也能在固定集合上进行比较。对于需要解释为什么某个转录本被报告、另一个没有被报告的项目，这种区分非常有用。

限制是结果依赖上游比对和可用证据。覆盖不足、复杂剪接、同源区域和不完整注释都会影响候选结构与定量。软件更新增加的其他模式也有各自适用条件，本次只运行普通短读段分析，没有测试长读段、混合读段或新生转录模式，因此不会把它们写成已经验证的能力。

这次实跑最值得带走的做法，是把结构检查放在丰度解释前面。先确认参考允许哪些模型，再看实际读段支持哪些连接，随后核对输出编号及其对应关系，最后解释相对丰度。这样面对一条模型的一百万 TPM，或者一次运行多出一个转录本时，读者能够提出可检查的问题，而不是急着给出表达升高或新转录本发现的结论。

下载本组复现脚本和实际输出。组装方法背景可参阅 StringTie 原始论文，实验中的链方向和剪接比对设置可核对 HISAT2 官方手册。

当统一注释和分析步骤已经确定，需要保存中间结果并稳定处理多样本转录组时，可了解云生信生物信息学分析服务器服务。