本文轉載自嘉因微信公眾號，已獲得授權。查看最新文章，敬請關注嘉因，微信ID：rainbow-genome
作者：小哈 來源：嘉因
大家都會做方便面，有人做辛拉面，有人做三鮮伊面，工藝有何不同？
大家都會做RNA-seq，有人能篩出有意義的基因，有人能找出有價值的線索，有人。。。差別在哪？
前五期介紹了回帖、均一化處理、差異基因篩選、畫heatmap和富集分析的合理方法：
第五期：gene model：同樣算read counts，為什么公司跟師兄算的結果不一樣？| Ensembl、UCSC、Refseq，該用哪個
第一期：數據預處理：同一套RNA-seq，為什么公司做的跟師兄跑的結果不一樣？ | TPM、read counts、RPKM/FPKM你選對了嗎？
第二期：差異基因篩選：同一套RNA-seq，公司篩出的差異基因跟師兄篩出的為什么不一樣？| Pvalue, FDR, cutoff
第三期：heatmap：heatmap畫不好會得出錯誤結論 | 數據預處理、聚類分析，HCL、 K means里的講究
第四期：富集分析：富集分析，倆人做的結果差5歲 | 你用的注釋文件有多老？
本文先解決去除rRNA的RNA-seq問題，后面談ensembl基因組使用的注意事項。
實驗和分析去除rRNA
通常做lncRNA-seq時會用到特殊的建庫方式，去除rRNA，而不是抓polyA尾巴，這樣可以看到不帶polyA尾巴的lncRNA。
第五期說Ensembl的基因注釋最准確，基因最全；
好！我用Ensembl的GTF做mapping；
第一期說衡量基因的豐度要用TPM；
好！我用TPM排序，看看哪個基因表達豐度最高。
Q：暈！怎么這么多rRNA ！？哥花了大價錢買kit去除rRNA，怎么測序還是測到這么多rRNA？！
A：淡定！通常實驗去除rRNA，還是會剩0.1-10%的rRNA。
sample間rRNA去除效率不同，剩下那些rRNA的read count會影響TPM的計算。因此，在mapping前，先刪掉GTF文件里的rRNA。ensembl的GTF有一列專門標注了gene biotype，用“grep -v rRNA”一條命令搞定。
用刪掉rRNA的GTF做mapping，不影響rRNA mapping到genome，也就是不影響回帖率；只是在算TPM的時候，不把rRNA的read counts算在內。這樣就不會因為rRNA去除效率的差異而影響其他基因豐度的計算。
看看用刪掉rRNA的GTF文件做mapping后，Top 10的基因是誰吧！
暈！怎么最多的還是無聊的snRNA？
Q：我想把snRNA、snoRNA也刪掉！Ensembl的GTF gene biotype有這么多種ncRNA。只留下lincRNA，其余的都刪掉多清凈？
A：不能刪！
tRNA（76-90 nt）、scRNA（100 nt）、miRNA（22 nt）因為長度太短撈不到，所以read counts很低，影響不到TPM的計算。剩下的snoRNA、snRNA，你確定它們在group間沒有差異表達？也許他們也有着重要的功能，在你的處理條件或組織類型里表達量提高或被抑制，從而發揮了調控作用，或被調控。
所以，如果rRNA去除效率差異太大，刪掉GTF文件里的rRNA就好。其實，通常也不需要這樣處理，rRNA read count的差異對DESeq2和edgeR結果不會影響很大。
關於GTF里gene biotype的處理，沒有找到文獻支持，只有comment，例如：
https://www.biostars.org/p/130420/
https://www.biostars.org/p/106590/
Ensembl基因組使用注意事項
Q：我又想到一個辦法：rRNA通常是repeat，ensembl有三種fasta文件，其中dna_rm的fasta文件已經將repeat變成N，所以，rRNA的read就無法mapping回去，這是個好主意嗎？
A：推薦dna或dna_sm版本的fasta文件，不要用dna_rm。點擊左下角“閱讀原文”直達原貼
Q：同樣是dna或dna_sm，又有Primary assembly和toplevel，該選哪個？
A：通常推薦Primary assembly。因為Primary assembly去掉了toplevel當中的單倍體和patch，它們會干擾分析結果。
Q：我只想要chr1-22和X、Y，那些單個chromosome文件是primary assembly？還是toplevel？
A：是Toplevel。所以，還是用primary assembly吧。
如果你不想要scaffolds，可以提取primary assembly文件里的chromosome。不過，不建議這樣做，因為那些scaffolds上面還有300多個基因。
Q：ensembl的版本號更新太頻繁，我需要每次mapping都下載最新版本嗎？
A：對於基因組成熟的物種，例如人和小鼠等模式生物，最近版本間差異不大，hg19、hg38和mm9、mm10對應的版本都可以。這取決於你想要一起分析的數據用的是哪個基因組版本。ENCODE項目最新版本都用hg38和mm10了。不過geo上的大多數數據都是hg19和mm9的。用liftover轉換基因組版本，簡單快捷，但會損失信息。實在不行，就把公共數據用最新版本重新跑一下。
其他該避免的陷阱：
-
確保文件下載完整
-
fasta文件和gtf文件版本對應
-
注意其他來源的index文件的基因組版本
