精選推文 | 基于三代轉錄組的基因注釋踩坑經歷以及GSAman使用
邀請并收到一位「GSAman」用戶的稿件，非常詳盡且實在。相信這份推文可以為一些做功能基因組方面工作的朋友，提供實用參考。?-- CJ-陳程杰
前言
隨著測序技術的進步和普及，現如今已經步入到“功能基因組學”（functional genomics)的時代。基因組注釋是各種基因功能研究的基礎，在之前的《生信石頭》推文中，CJ大佬已經多次強調一個準確可靠的基因注釋的重要性。
基因注釋可以分為結構注釋和功能注釋，結構注釋描述了一個基因在基因組上具體的位置信息，是功能注釋的前提。
目前做結構注釋大約有三種途徑：
從頭注釋(de novo prediction)：通過已有的概率模型來預測基因結構，在預測剪切位點和UTR區(qū)準確性較低
同源預測(homology-based prediction)：有一些基因蛋白在相近物種間的保守性高，所以可以使用已有的高質量近緣物種注釋信息通過序列聯配的方式確定外顯子邊界和剪切位點
基于轉錄組預測(transcriptome-based prediction)：
通過物種的RNA-seq數據輔助注釋，能夠較為準確的確定剪切位點和外顯子區(qū)域。
----?引徐洲更，https://blog.csdn.net/u012110870/article/details/82500684
相較于前兩種方法，基于轉錄組預測用的是實實在在檢測到的RNA，結果最為真實可靠。比如EST（現在好像比較少見），各種類型的RNA-seq，最常見的是普通mRNA-seq。RNA-seq最大的不足是，由于基因轉錄表達受內外條件影響，單樣本測序只能獲得一部分基因表達信息。此外，二代測序短讀長的特點會導致不能獲得轉錄本的完整序列，引起信息的錯誤或丟失。好在三代測序的出現給注釋提供了極大便利，三代測序技術無需打斷mRNA序列，可以完整測通，在理想情況下，只需要將測得序列比對回基因組，直接就可以得到完整的基因結構注釋。
言歸正傳，我的一個課題就是對某個物種進行基因組注釋，最開始測了以下數據：
多個樣品分別提取mRNA等量混合后建一個庫測全長轉錄本，pacbio平臺。
不同樣品分別建庫，鏈特異性建庫，illumina平臺。
本來計劃測基因組的，但是國外有實驗室已經基于三代測序公開了一個品種A的基因組，還順便把我們實驗室的品種B用二代測了。所以就不測基因組，計劃用公開的數據做參考基因組進行注釋。（重點，此處挖一大坑?。?/span>
上文已經說了基于轉錄組預測的缺點，所以從頭注釋(de novo prediction)仍是必不可少的，同源預測(homology-based prediction)就見仁見智了，由于與本物種相關的公開注釋質量一般，且自己已經測了足量的轉錄組數據，所以同源預測就略過?；诙D錄組數據的基因注釋流程非常多，常見的有maker，braker等，網上一搜一大把，此處不贅述。
三代轉錄組注釋流程
目前三代轉錄組數據分析流程中，pacbio官方發(fā)布的工具包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序列（坑?。?/h4>
這一步看上去沒什么問題，因為flnc序列中確實會存在重復或冗余。官網上cluster的邏輯好像也沒什么毛病。然而實際比對時發(fā)現，該步驟合成所謂的transcript可能會產生錯誤。
看圖：
紅色箭頭指向的轉錄本明顯有一處比對不正常的，根據cluster report找到對應的transcript和flnc序列，分別比對回基因組，如下圖：
第一條藍色的是cluster合并得到的transcript，下面兩條紅色的是flnc序列，白色空缺是內含子區(qū)域。很明顯，兩條可變剪切產生的轉錄本很不幸的被合并成了一個錯誤的轉錄本。這種問題非常隱秘，如果不是后期分析時仔細檢查比對結果，是很難發(fā)現的。
目前官方軟件并沒有可供調整的參數選項。我想到的第一個解決方法是使用別的聚類軟件來去冗余，比如cd-hit。cd-hit是找出一個最具代表性的序列，不會變更序列信息，但是實際使用后還是存在問題，可能會合并掉可變剪切差異比較小，或者3’端長度有差異的的轉錄本。第二個解決方法是比對之后再合并。
4. 序列矯正（polish）
三代測序都有錯誤率高的問題，如果是hifi reads，可能不需要糾錯，但是按以上流程走出來的flnc reads，包含了許多低質量ccs序列，需要做糾錯。我從二代測序數據中各抽出一部分，用lordec軟件進行糾錯。
可以看到，糾錯的效果還是很明顯的，錯誤比較多的一般都是低質量的flnc序列，另外pacbio測序錯誤貌似堿基插入居多。要確保二代轉錄組能覆蓋三代測序的內容，最好就是把三代測序的RNA分一部分去做二代。
5. 比對回基因組（坑?。?/h4>
如果做的是無參轉錄組注釋，把flnc序列折疊一下就可以做后續(xù)篩選和功能注釋了。但是對于有參注釋，一定是要比對回基因組，確定轉錄本在基因組上的位置信息。官方給的軟件是pbmm2（官方對許多三代轉錄組分析工具進行了整合，例如cDNAcupcake，SQANTI3）。Pbmm2是minimap2的封裝，比對獲得SAM文件后，用isoseq3 collapse對冗余轉錄本折疊， 因為三代測序雖然能確保3’端信息完整，但是mRNA的5’端降解是不可避免的（解決的方法是做CAGE-seq，把轉錄起始位點測出來），collapse這一步會把能完全匹配上更長轉錄本的序列折疊掉，順手解決了之前聚類的問題。如果這個基因組有公開的注釋，用Pigeon基于公開注釋，對新注釋進行分類和統(tǒng)計。
看上去很順利，很簡單對不對，然而比對這一步就充滿了問題。
首先，由于參考基因組來自別的實驗室，由二代數據組裝，幾百個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規(guī)律，有的地方為了符合GTAG規(guī)律，又搞出deletion和snp變異。
6.基因注釋結構矯正-apollo（坑?。?/h4>
怎么辦呢？既然軟件不行，那就只能上人工了，在GSAman開發(fā)之前，人工做基因結構注釋的軟件相對較好的大概就是apollo了。安裝復雜，好在有docker，勉強也用上了。然而及其難用，首先它對gff文件的兼容性極差，也沒有報錯信息，只能來來回回的調整gff，好不容易正確顯示出基因結構了，想改也不是那么容易，首先要把目的轉錄本一個個調整到注釋欄，可供矯正的操作非常少，遇到復雜的情況，修一個轉錄本要花好幾分鐘。修改過的轉錄本與沒修改的是隔離的，只能單獨導出修改后的轉錄本，再合并。我猜開發(fā)者可能沒有想過有人會大量修改基因注釋，apollo改一兩個注釋還算輕松，改整個基因組就是折磨了（也可能是我沒怎么看軟件文檔，或許有更快捷的手段）。我花了一星期，才檢查和矯正了約一千個基因。而且當時用的還是錯誤的參考基因組，矯正起來非常頭疼。因為我不能把所有精力放在用apollo做矯正上，于是就擱置了。
GSAman登場
之前曾搜索到的GSAme的文章，很可惜CJ大佬暫停開發(fā)了一段時間，沒了消息。然而幸運的是，今年10月初，就在我重新測完基因組，準備再次死磕apollo的時候，意外刷到的公眾號內測的信息，第一時間報名參加（當時群里才20~30人，現在已經300多人），這一個月的時間里，CJ大佬高強度迭代了幾十個版本，我也是用著不斷升級的GSAman，花了約二十天完成了矯正，這當中還因為上文提到的各種坑完全重做了一次。效率可見一斑！
關于如何具體使用GSAman矯正注釋，公眾號里已經發(fā)了好些推文，我這里就不再演示。
對于基因組注釋矯正，一個難點是定位錯誤，肉眼一個個看注釋太慢了。一般情況下，出錯的是少數（如果大多數都有錯建議重頭找原因?。?，很多時間會浪費正確的注釋上。
我的方法是基于比對得到的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，甚至兩個基因混為一談的。
另外，如果讀者中有大佬對三代轉錄組注釋有更好的方法和見解，懇請在留言區(qū)留言指點一二，感謝！
最后非常感謝CJ大佬的開發(fā)工作，造福了我們這些搞基因組注釋的，10月以來幾乎每天都會發(fā)布新的版本，非常nb！點贊?。?！
