ChengyangWang的个人博客分享 http://blog.sciencenet.cn/u/ChengyangWang

博文

同样用Ensembl算TPM,结果还是不一样? | Ensembl的注意事项

已有 4711 次阅读 2017-12-15 14:53 |个人分类:RNA-seq|系统分类:科普集锦| RNA-seq


本文转载自嘉因微信公众号,已获得授权。查看最新文章,敬请关注嘉因,微信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


第三期:heatmapheatmap画不好会得出错误结论 | 数据预处理、聚类分析,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文件的基因组版本




https://wap.sciencenet.cn/blog-3372875-1089848.html

上一篇:Ensembl、UCSC、Refseq,该用哪个
下一篇:TPM、read counts、RPKM/FPKM你选对了吗?
收藏 IP: 124.77.56.*| 热度|

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-3-29 23:39

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部