【原】精选推文 | 基于三代转录组的基因注释踩坑经历以及GSAman使用

邀请并收到一位「GSAman」用户的稿件，非常详尽且实在。相信这份推文可以为一些做功能基因组方面工作的朋友，提供实用参考。 -- CJ-陈程杰

前言

随着测序技术的进步和普及，现如今已经步入到“功能基因组学”（functional genomics)的时代。基因组注释是各种基因功能研究的基础，在之前的《生信石头》推文中，CJ大佬已经多次强调一个准确可靠的基因注释的重要性。
基因注释可以分为结构注释和功能注释，结构注释描述了一个基因在基因组上具体的位置信息，是功能注释的前提。

目前做结构注释大约有三种途径：

1. 从头注释(de novo prediction)：通过已有的概率模型来预测基因结构，在预测剪切位点和UTR区准确性较低
   同源预测(homology-based prediction)：

2. 有一些基因蛋白在相近物种间的保守性高，所以可以使用已有的高质量近缘物种注释信息通过序列联配的方式确定外显子边界和剪切位点

3. 基于转录组预测(transcriptome-based prediction)：
   通过物种的RNA-seq数据辅助注释，能够较为准确的确定剪切位点和外显子区域。

---- 引徐洲更，https://blog.csdn.net/u012110870/article/details/82500684

相较于前两种方法，基于转录组预测用的是实实在在检测到的RNA，结果最为真实可靠。比如EST（现在好像比较少见），各种类型的RNA-seq，最常见的是普通mRNA-seq。RNA-seq最大的不足是，由于基因转录表达受内外条件影响，单样本测序只能获得一部分基因表达信息。此外，二代测序短读长的特点会导致不能获得转录本的完整序列，引起信息的错误或丢失。好在三代测序的出现给注释提供了极大便利，三代测序技术无需打断mRNA序列，可以完整测通，在理想情况下，只需要将测得序列比对回基因组，直接就可以得到完整的基因结构注释。

言归正传，我的一个课题就是对某个物种进行基因组注释，最开始测了以下数据：

1. 多个样品分别提取mRNA等量混合后建一个库测全长转录本，pacbio平台。

2. 不同样品分别建库，链特异性建库，illumina平台。
   本来计划测基因组的，但是国外有实验室已经基于三代测序公开了一个品种A的基因组，还顺便把我们实验室的品种B用二代测了。所以就不测基因组，计划用公开的数据做参考基因组进行注释。（重点，此处挖一大坑！）

上文已经说了基于转录组预测的缺点，所以从头注释(de novo prediction)仍是必不可少的，同源预测(homology-based prediction)就见仁见智了，由于与本物种相关的公开注释质量一般，且自己已经测了足量的转录组数据，所以同源预测就略过。基于二代转录组数据的基因注释流程非常多，常见的有maker，braker等，网上一搜一大把，此处不赘述。

三代转录组注释流程

目前三代转录组数据分析流程中，pacbio官方发布的工具包PacificBiosciences/pbbioconda: PacBio Secondary Analysis Tools on Bioconda (github.com)最为全面。

1.从sureads生成ccs序列

这里就遇到了第一个小坑，要质量还是要数量？官方把高质量的ccs序列称为HiFi序列，生成HiFi序列要求多次pass的subreads，按照早期版本软件的默认的处理步骤，会过滤掉很多只有一两个pass的subreads，当测序量不足的时候，会丢失许多低频转录本，加大测序量可以缓解这个问题（￥得加钱￥）。好在目前最新版本的ccs软件reads.bam | CCS Docs考虑到这个问题，--all参数生成reads.bam，包含所有质量的CCS。低质量的CCS可以留到后面步骤处理。

2.从ccs序列生成FLNC序列

这步没什么好说的，软件limahttps://lima.how去接头，isoseq3Schematic Workflow | Iso-Seq Docs (isoseq.how)refine步骤去嵌合体，polyA尾巴。

3. isoseq3 cluster 合并flnc序列（坑！）

这一步看上去没什么问题，因为flnc序列中确实会存在重复或冗余。官网上cluster的逻辑好像也没什么毛病。然而实际比对时发现，该步骤合成所谓的transcript可能会产生错误。
看图：

红色箭头指向的转录本明显有一处比对不正常的，根据cluster report找到对应的transcript和flnc序列，分别比对回基因组，如下图：

第一条蓝色的是cluster合并得到的transcript，下面两条红色的是flnc序列，白色空缺是内含子区域。很明显，两条可变剪切产生的转录本很不幸的被合并成了一个错误的转录本。这种问题非常隐秘，如果不是后期分析时仔细检查比对结果，是很难发现的。
目前官方软件并没有可供调整的参数选项。我想到的第一个解决方法是使用别的聚类软件来去冗余，比如cd-hit。cd-hit是找出一个最具代表性的序列，不会变更序列信息，但是实际使用后还是存在问题，可能会合并掉可变剪切差异比较小，或者3’端长度有差异的的转录本。第二个解决方法是比对之后再合并。

4. 序列矫正（polish）

三代测序都有错误率高的问题，如果是hifi reads，可能不需要纠错，但是按以上流程走出来的flnc reads，包含了许多低质量ccs序列，需要做纠错。我从二代测序数据中各抽出一部分，用lordec软件进行纠错。

可以看到，纠错的效果还是很明显的，错误比较多的一般都是低质量的flnc序列，另外pacbio测序错误貌似碱基插入居多。要确保二代转录组能覆盖三代测序的内容，最好就是把三代测序的RNA分一部分去做二代。

5. 比对回基因组（坑！）

如果做的是无参转录组注释，把flnc序列折叠一下就可以做后续筛选和功能注释了。但是对于有参注释，一定是要比对回基因组，确定转录本在基因组上的位置信息。官方给的软件是pbmm2（官方对许多三代转录组分析工具进行了整合，例如cDNAcupcake，SQANTI3）。Pbmm2是minimap2的封装，比对获得SAM文件后，用isoseq3 collapse对冗余转录本折叠， 因为三代测序虽然能确保3’端信息完整，但是mRNA的5’端降解是不可避免的（解决的方法是做CAGE-seq，把转录起始位点测出来），collapse这一步会把能完全匹配上更长转录本的序列折叠掉，顺手解决了之前聚类的问题。如果这个基因组有公开的注释，用Pigeon基于公开注释，对新注释进行分类和统计。

看上去很顺利，很简单对不对，然而比对这一步就充满了问题。
首先，由于参考基因组来自别的实验室，由二代数据组装，几百个scaffold，我就基于参考基因组用ragtag拼接了一下。等到转录组比结果导入igv一看，满屏的indel和snp变异，
好像这个材料突变了几万年一样，人都看傻了。后来用sanger测序检查了备份材料的部分序列，确定是做转录组测序的材料和别的实验室做基因组测序的材料是不一样的。只能赶紧把手头材料做三代基因组重测序，前前后后浪费了大半年时间。
拿到三代基因组数据，马不停蹄的做组装，再来比对，总算没有满屏的变异了。但是仍然有问题。
看图：

默认参数下，minimap2将一些微小外显子判定为上下游外显子边界的变异，最常见的就是外显子5’端的“插入突变”。这种情况下，由SAM文件生成的gff文件就丢失了外显子信息，并生成错误的外显子边界。这个问题github上也有人询问过，李恒大神表示20bp以下的微外显子很难判断Documentation: Limitations "small exons" · Issue #948 · lh3/minimap2 (github.com)。既然minimap2不行，那其他软件行不行？能做三代转录组比对的软件除了minimap2，还有GMAP，deSALT，STAR等，搜一搜相关文章也能找到不少软件，但是就目前我尝试的结果来说，deSALT相对好一点，10bp的外显子没有问题，再小的外显子就无能为力了，此外该软件会犯一些其他错误，比如有的地方剪切位点不符合GTAG规律，有的地方为了符合GTAG规律，又搞出deletion和snp变异。

6.基因注释结构矫正-apollo（坑！）

怎么办呢？既然软件不行，那就只能上人工了，在GSAman开发之前，人工做基因结构注释的软件相对较好的大概就是apollo了。安装复杂，好在有docker，勉强也用上了。然而及其难用，首先它对gff文件的兼容性极差，也没有报错信息，只能来来回回的调整gff，好不容易正确显示出基因结构了，想改也不是那么容易，首先要把目的转录本一个个调整到注释栏，可供矫正的操作非常少，遇到复杂的情况，修一个转录本要花好几分钟。修改过的转录本与没修改的是隔离的，只能单独导出修改后的转录本，再合并。我猜开发者可能没有想过有人会大量修改基因注释，apollo改一两个注释还算轻松，改整个基因组就是折磨了（也可能是我没怎么看软件文档，或许有更快捷的手段）。我花了一星期，才检查和矫正了约一千个基因。而且当时用的还是错误的参考基因组，矫正起来非常头疼。因为我不能把所有精力放在用apollo做矫正上，于是就搁置了。

GSAman登场

之前曾搜索到的GSAme的文章，很可惜CJ大佬暂停开发了一段时间，没了消息。然而幸运的是，今年10月初，就在我重新测完基因组，准备再次死磕apollo的时候，意外刷到的公众号内测的信息，第一时间报名参加（当时群里才20~30人，现在已经300多人），这一个月的时间里，CJ大佬高强度迭代了几十个版本，我也是用着不断升级的GSAman，花了约二十天完成了矫正，这当中还因为上文提到的各种坑完全重做了一次。效率可见一斑！
关于如何具体使用GSAman矫正注释，公众号里已经发了好些推文，我这里就不再演示。
对于基因组注释矫正，一个难点是定位错误，肉眼一个个看注释太慢了。一般情况下，出错的是少数（如果大多数都有错建议重头找原因！），很多时间会浪费正确的注释上。
我的方法是基于比对得到的SAM文件，通过CIGAR 来找错误注释。例如：用grep从sam中筛选出符合M[0-9]+N[0-9]I[0-9]+M的reads。
比如这些reads的比对信息：

这是从sam提取的前6列。其中外显子5’端插入被我标注了出来，应该是一个3bp的微外显子被软件当成了插入突变。
实际如何呢，看图：

这个就可以确定是一个3bp的微外显子了。用GSAman右键直接插入一个外显子，然后调整到合适大小。在筛选时如果遇到同一位置上很多reads一样的变异，那大概率是需要矫正的地方。
其它外显子边界的indel类错误大都可以这样来定位，难的是错配，因为错配在CIGAR中看不出来，deSALT输出的sam没有错配的位置信息。不过Minimap2可以在比对时增加MD参数获得错配位置从而进行定位。
把定位出来的错误位置做成表，结合GSAman的便捷操作，我差不多一星期就能把整个基因组粗略的矫正完，更精细的矫正就不太需要了，具体基因具体分析。GSAman还提供了BLAT等手段进行精细预测，非常实用。

---

强烈建议做单（少量）基因功能研究老师同学也拿GSAman好好检查一下手头的基因，特别是有转录组数据的情况下。如果基因结构一开始就有问题，那后面的工作就要大打折扣。就我手头的结果来看，很多参考基因都有问题，有的微外显子丢失，有的UTR边界错误或者干脆没有UTR，只有CDS，甚至两个基因混为一谈的。

另外，如果读者中有大佬对三代转录组注释有更好的方法和见解，恳请在留言区留言指点一二，感谢！

最后非常感谢CJ大佬的开发工作，造福了我们这些搞基因组注释的，10月以来几乎每天都会发布新的版本，非常nb！点赞！！！