做过1000遍RNA-seq的老司机告诉你如何翻车

熟悉我的人都知道RNA-seq是我的拿手好戏(如果你不熟悉我,今天过后请记住)。

但是我今天处理了一个公共数据,比对率低的惊人。

究竟为什么会发生这种小概率事情呢?

  • 是测序数据质量不好?
  • 难道grcm38与mm10有差别?
  • 比对工具的默认参数不行?

这篇文章告诉你,老司机如我,也有翻车的时候。

数据比较新,所以理所当然的认为测序数据肯定是OK的。

数据下载地址

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE81916

下载sra数据,转换为fastq的过程我就不讲解了!(请翻看我以前的文章)

(17)一个ChIP-seq实战-生信菜鸟团博客2周年精选文章集

(18)一个mRNA-seq实战-生信菜鸟团博客2周年精选文章集

(19)一个affymetrix表达芯片实战-生信菜鸟团博客2周年精选文章集

(20)一个WES实战-生信菜鸟团博客2周年精选文章集

总之,这次处理了GEO里面的4个公共数据,数据量如下:

Written 30468155 spots for SRR3589959.sraWritten 52972617 spots for SRR3589960.sraWritten 36763726 spots for SRR3589961.sraWritten 43802631 spots for SRR3589962.sra

我用的是hisat2工具来比对,一般情况下就用默认参数。

reference=/home/jianmingzeng/reference/index/hisat/grcm38/genome~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589959.fastq -S control1.sam 2>control1.log~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589960.fastq -S control2.sam 2>control2.log~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589961.fastq -S Akap951.sam 2>Akap951.log~/biosoft/HISAT/current/hisat2 -p 5 -x $reference -U SRR3589962.fastq -S Akap952.sam 2>Akap952.logls sam |while read id;do (nohup samtools sort -n -@ 5 -o ${id%%.}.Nsort.bam $id &);done

但是让我意外的是比对率出奇的低!请看我结果。。。

0.48% overall alignment rate0.62% overall alignment rate0.48% overall alignment rate0.49% overall alignment rate

起初我怀疑是参考基因组用错了,但是我查看了GEO里面的介绍,的确是mouse的ESC,所以我用grcm38没有问题!然后我怀疑是测序数据质量的问题,但是质量再差也不会导致如此低的比对率啊!所以我还是用fastqc检查了一下:

果然,质量值好到爆!

而且我抽取了几条序列去blat一下,发现也可以比对,而且明显是跨越intron的比对,超级经典的RNA-seq数据呀,这完全没问题啊。( 其实我这个blat结果也没有看仔细,正常的reads不应该被截成比对到基因组的正负链的,这其实预示着我把PE序列拼接了。)

那么就是我hisat2这个步骤的问题,我首先怀疑是不是我下载hisat的index搞错了,虽然看起来我命名是grcm38,但是有可能是我下载错误!我打开了sam文件看了看开头:

貌似的确是mouse基因组的染色体长度呀!很诡异,而且我清楚的记得,我下载的就是mouse的基因的索引。这里也没有问题。

https://ccb.jhu.edu/software/hisat2/index.shtml

难道grcm38与mm10有差别? 我就先用bowtie2测试一下mm10吧,毕竟我还没有hisat2的mm10的index。

head -1000 SRR3589959.fastq >tmp.fq~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 6 -x ~/reference/index/bowtie/mm10 -U tmp.fq -S tmp.sam

结果我挑出来的这1000条序列,全军覆没了,0.00% overall alignment rate,我傻眼了! 没办法呀,逼着我换hg19参考基因组看看:

~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -p 6 -x ~/reference/index/bowtie/hg19 -U tmp.fq -S tmp.hg19.sam

仍然是全军覆没了,0.00% overall alignment rate,心好累,继续傻眼!

回过头看了看fastqc的报告,发现前面10个碱基的确有问题的!如果只是对RNA-seq进行定量,可能需要trim掉,但是我以前从来不trim,照样不影响比对。(如果你认真研究过测序流程,你应该知道前面这十个左右的碱基其实是bias)

不过,暂时看到这个问题,我就试着解决一下吧,先从这个思路来,而且比对工具里面本来就有这个选项,没必要自己来trim的!

具体参数网站:https://ccb.jhu.edu/software/hisat2/manual.shtml

-5/--trim5 <int> trim <int> bases from 5'/left end of reads (0)-3/--trim3 <int> trim <int> bases from 3'/right end of reads (0)

所以我加上了-p 6 -5 10 -3 10 --local 参数,比对人,可以拿到 35.60% overall alignment rate,比对mouse,可以拿到 98.80% overall alignment rate ,我勒个去,问题出来了,看起来好像是应该trim掉呀。以前的万能默认参数不行了?

但是有个问题,虽然我用local模式都比对上了,但是首先100bp的reads我切成了80,而且都是40M,40S,说明只有reads的一半成功比对到了参考基因组序列!

我然后用同样的参数,测试了hisat2工具,但是hisat2里面压根就没有local的选项,仅仅是trim一下,对比对的改善毫无意义,所以重点在于--local这个参数,但它只是表象,本质还是这个测序数据出问题了!

数据为什么会出问题呢?

终于,我想起了最开始的fastqc,决定再回过头看看测序数据的fastqc报告,请看下图,这么重要的图我居然忽略掉了,忽略掉了,略掉了,掉了,了。。。

再联想到前面的40M,40S我瞬间明白了,这肯定是一个双端测序被我搞成了单端测序数据!

而且我再去GEO介绍上面看,上面赫然写着PAIRED! 我死也想不明白,我明明是加了--split-3 参数呀,为什么sra转换成fastq会出这么明显的错误呢?

然后我检查我的脚本,我自己从我的博客里面复制了我的代码。

唯一值得你看的就是上面这个图

是--不—是全角半角害死人,而且这个参数不识别它居然不停止,而是忽略我的参数! 是--不—是全角半角害死人,而且这个参数不识别它居然不停止,而是忽略我的参数! 是--不—是全角半角害死人,而且这个参数不识别它居然不停止,而是忽略我的参数!

更要命的是我把wget跟fastq-dump一起运行,而wget会给出一大堆的log日志把fastq-dump的报错日志给掩盖了。

此处请脑补我内心的呐喊!

因为前面一直处理的是单端的数据,所以这个错误没有被发现。

我痛恨我自己直接拷贝了博客的脚本!我痛恨 --这样的参数设置!

下面是我修改后的代码!

cut -f 3 config.txt |while read id ; do wget $id 2>/dev/null ;donels *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump --gzip --split-3 $id;done

这就是老司机翻车以及我debug的心路历程,希望你们引以为戒!

  1. 因为太熟练了,所以我已经很久没有关心过完整的fastqc报告。
  2. 因为太熟练了,所以我把wget跟fastq-dump一起运行,很久没关心过log文件。
  3. 因为太熟练了,所以我对用了一千遍fastq-dump的命令产生了一万个放心。
  4. 因为太熟练了,所以我先想了各种复杂的情况而没有一开始就从最简单的问题入手。

所以,一切看起来都是因为我是一个老司机!

对了,顺便说一句,我已经把下好的sra数据删除了!(你想笑就笑吧)

如果你认为自己还是一个菜鸟就请记住老司机的这句话:

行走江湖,靠的是一个

文:Jimmy

编辑校对:一只思考问题的熊

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2017-01-13

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏FreeBuf

小心,Android木马工具SpyNote免费啦!远程监听就是这么简单

只要是做生意,都得讲究价值规律,黑市也不例外。某款产品要是搞打折促销,群众们必然蜂拥而至——要是免费大派送,那一传十十传百的速度又怎是门庭若市可形容! 最近Pa...

3617
来自专栏pangguoming

Android平台相机接口的应用

第一部分、前述: Android作为Google移动互联网战略的重要组成部分,将进一步推进“随时随地为每个人提供信息”这一企业目标的实现。Google的目标是让...

3775
来自专栏linux驱动个人学习

MIPI协议-DSI

对于现代的智能手机来说,其内部要塞入太多各种不同接口的设备,给手机的设计 和元器件选择带来很大的难度。下图是一个智能手机的例子,我们可以看到其内部存储、显示、摄...

4096
来自专栏华章科技

7大笔记应用,让你的代码效率翻7倍

但是大多数笔记应用的设计并不是以程序员作为目标受众,这些程序可能会让使用者用起来很难受,甚至完全放弃这些工具。这就是为什么我们为你找来了这些最好的笔记工具。快来...

2102
来自专栏机器之心

资源 | 挑战谷歌,Facebook 发布交互数据可视化工具 Visdom

选自GitHub 机器之心编译 参与:微胖、吴攀 FAIR 发布了 Visdom,一款可在 Torch、PyTorch 以及 NumPy 上实现交互式数据可视化...

3018
来自专栏FreeBuf

浅谈拒绝服务攻击的原理与防御(5) | NTP反射攻击复现

0×01 故事起因 前两天以为freebuf上的网友stream(年龄、性别不详)私信我说他在阿里云上的服务器被NTP攻击了,流量超过10G,希望我帮忙一起分析...

4645
来自专栏Coding01

我也来打造一个个人阅读追踪系统

国庆放假期间,偶然发现这篇文章《Serverless实战:打造个人阅读追踪系统》http://insights.thoughtworks.cn/serverle...

1192
来自专栏安富莱嵌入式技术分享

【原创开源】网络版二代双通道示波器开源发布,支持电脑,手机和Pad等各种OS平台访问

一代示波器发布于3年前,去年年底的时候发布了二代示波器,软件性能已经比较强劲,但依然有值得升级改进的地方,经过今年这半年多努力,在二代示波器的基础上再推出网络版...

2491
来自专栏酷玩时刻

QQ轻游戏入门到精通OR放弃?

注册很简单,使用已有Q号登录「厘米游戏」开放平台按照流程提交资料审核即可 。开发者接入官方说明文档

7794
来自专栏程序人生

分布式系统中的监工:Overseer

最近从无趣的工作中发现了有趣的事情,工作和业余时间都扑了些精力上去,本待上周末最终的成果出来后再写文章的,无奈事情太多,代码还没写完,二月上旬已过,再不写文章春...

3317

扫码关注云+社区

领取腾讯云代金券