搜索
查看: 12136|回复: 3

[mRNA-seq] 小白生信学习记(三)转录组分析流程

[复制链接]

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
发表于 2017-2-21 19:53:15 | 显示全部楼层 |阅读模式
注:本文是生信媛微信公众号原创文章
作者:生物女博士
原文链接:
http://mp.weixin.qq.com/s?__biz=MzI1MjU5MjMzNA==&mid=2247483788&idx=1&sn=a8afed11acc6123ad9e5d210d85332c5&chksm=e9e0282dde97a13b656daee65bd5d7aa62d0c29f078275f15422da180f58faec0dcbfe75198c#rd


这周开始,正式进入RNAseq分析部分。
      在分析数据的过程中,其实会有很多需要了解的细节,这些细节需要大家通过自己查找文献等方式进行补充。这个系列的写作目的主要是带着大家实际上机操作,快速走一个分析的框架,学会软件安装、基本linux命令及R语言等等。而分析流程中的其他许多重要的东西,会逐渐在后续的写作中去详细讲解、逐步完善。


      转录组测序分析可以分为有参转录组分析和无参转录组分析。有参无参的意思是,有/无参考基因组。通常模式动植物基本都会有参考基因组。而这里仅以有参转录组分析为例。
      有参转录组的分析流程大体是这样的:
        我们从下机数据开始说起。
        在测序之前,通常合同上会规定交付的数据量,以及测序质量符合什么样的要求:比如平均Q20≥85% ,平均Q30≥80% 等(Q20、Q30是什么意思,我们等下解释)。

        所以在拿到下机数据后,通常测序公司会进行数据过滤(比如去除接头序列、舍去质量太差的reads等质控),最后交付符合合同要求测序质量的数据。
        我们先来看看到我们手里的Illumina的数据长什么样子,我们可以用more命令查看:

        more xiaobai.fq


       上图显示的就是fastq(fq)格式的测序文件。这个格式曾在《小白生信学习记(一)》里简单说过。我们来复习一下,第1行是序列的名字。第2行是测的DNA序列。第4行则是测序质量值。1到4行组成了一条reads的完整信息。
        我们重点看一下这些质量值是什么意思。
        这些质量值是用下列的ASCII码表示的,ASCII码与数值的对应关系如下表(参考资料1):


        ASCII码与测序质量值的对应关系,每个公司不一样,而Illumina公司本身不同时期也使用过不同的换算方式。如下图(参考资料2):


       目前使用的,都是Illumina1.8+对应的换算方式。
       也就是L表示的那行:!是最低的质量,为0,而J是最高的质量,为41。以@为例,计算:查表可以知道 @ 对应的数值为64,但对于测序质量值来说,多了33(因为测序质量是从!开始用的),所以@ 对应的测序质量值应该为64-33=31。
       那为什么搞得这么复杂,而不用直接用数字表示?因为当质量比较高的时候,比如38,这个数值首先和原来的DNA序列的对应关系会变得复杂,且占用的存储也更多。
       这个质量值是什么意思?最初Sanger中心用Phred Quality Score来衡量该read中每个碱基的质量,Q=-10logP ,其中P代表该碱基被测序错误的概率,如果该碱基测序出错的概率为0.001,则Q应该为30,那么30+33=63,那么63对应的ASCii码为“?”,则在该碱基对应的质量值即“?”(参考资料1)。对应的,Q20表示碱基测序错误率为0.01。如下图(参考资料3):


        我们继续往下走流程。
        在拿到这些质控后的数据后,我们也可以自己亲自用软件(如:Fastqc)查看一下数据是否符合要求。如果感觉还不满意,也可以自己使用软件进行质控(比如:trimmomatic,Trim Galore等软件)
       这些经过质量控制的数据,我们称其为clean data 。对应的,那些未经过质控的数据,我们称为raw data。
       这些clean data ,我们会使用一些软件(比如:HISAT2,tophat2,STAR等)把它们和参考基因组序列进行比对,寻找每个reads的最佳匹配位置。我们把这个过程称为Mapping。
       接着,统计每个基因或者转录本的表达量(软件:HTSeq,RSEM等)后,进行差异表达分析(edgeR, DEseq2, EBseq等)。
        为什么我们需要寻找差异表达的基因?举例来说,在正常条件和某种试剂处理的时候,出现了某些我们关心的特性,但是我们可能不知道这些特性的发生是什么机制。那我们可以通过寻找这些有变化的基因,分析它们大致是参与了什么过程的基因起了变化,就可以为进一步深入研究提供线索。
        那如何知道这些基因参与了什么过程?我们可以通过GO分析和KEGG通路分析来提供这方面信息。这部分内容今天暂不介绍。

       今天我们先简单介绍到这里。上面提到的内容,我之后都会更为详细地进行一一讲解,并且给出命令以及数据供大家练习。
       生物信息学的发展是非常迅速,软件也在快速迭代的。以GO分析为例:2002年-2008年之间的软件至少有68个之多(参考资料4)!所以本系列文章不在于教会大家某个软件非常具体的操作,而是希望通过某些软件的学习,大家逐渐了解生物信息以及linux系统的使用和编程入门。


参考资料4:   Da Wei Huang等,Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists



欢迎到微信公众号订阅我们
生信媛
bio_sxy


本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x



上一篇:find和xargs连用虽好,但用起来要小心哦~
下一篇:R语言ggplot2画图系列——Pathway富集分析气泡图
回复

使用道具 举报

0

主题

15

帖子

215

积分

中级会员

Rank: 3Rank: 3

积分
215
QQ
发表于 2017-2-24 17:31:48 | 显示全部楼层
每天学一点儿,每天练一点儿...
回复 支持 反对

使用道具 举报

33

主题

46

帖子

230

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
230
 楼主| 发表于 2017-2-27 23:23:24 | 显示全部楼层
Martin 发表于 2017-2-24 17:31
每天学一点儿,每天练一点儿...

这孩子不多时必成大器!
回复 支持 反对

使用道具 举报

5

主题

26

帖子

291

积分

中级会员

Rank: 3Rank: 3

积分
291
发表于 2017-3-20 17:38:36 | 显示全部楼层
小白想问一下, We then aligned the remaining transcripts to the reference grape genome in the Swiss-Prot database using BLAST with the parameters “-e 1e-5 -F F -a 5”. The transcripts that could be aligned to reference sequences were selected and scanned to define coding regions (CDS).这一段里,说和Swiss-prot数据库做比对,用blast,可是,怎么用啊。。。
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 00:17 , Processed in 0.036605 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.