搜索
查看: 2974|回复: 0

[mRNA-seq] DEXSeq使用说明(2/3)

[复制链接]

13

主题

66

帖子

276

积分

中级会员

Rank: 3Rank: 3

积分
276
QQ
发表于 2017-4-7 14:58:58 | 显示全部楼层 |阅读模式
本帖最后由 惠真-市二医院 于 2017-4-7 15:04 编辑

DEXSeq使用说明(2/3)

使用R处理比对结果:
剩下的比对就在R里面完成了,这里使用数据是文章介绍的GEO数据,首先看看有哪些文件,我这里将转换后的gff文件重名命了下(Fruit.chr.gff):

分析的第一步是创建分析的DEXSeqDataSet对象,也就是说根据实验来设计比对的方式,该说明使用的是默认的方式:


countfiles:一个包含路径的字符型的向量,指向前面dexseq_count.py分析出的来reads计数结果;
sampleData:包含了说明信息的数据框(例如处理方式、组织类型、表型信息等)。该数据框的行数目要和countData矩阵相同(countData涉及到另一种分析方式DESeqDataSet,对应分配每个样本说明信息)。这里sampleTable表格包含了样本名称(treated vs untreated),样本处理条件(condition)和文库类型(libType),由于Rdefault.stringsAsFactors()TRUE,所以这里的因子为样本处理条件和文库类型分别包含了2个水平。
design:该公式指定了实验的设计方式。它必须包含一个交互项,该交互项需包含sampleData中的一个说明信息(列)和变量”exon”,例如图示的”coditon:exon”,这里就意味着在”condition”条件下的exon usage的改变。
flattenedfile:一个包含路径的flattened注释文件(就是我们之前使用dexseq_prepare_annotaiton.py转换后得到的gff文件)。
这里产生的DEXSeqDataSet对象,我的理解是一个数据集文件(DEXSeqDataSetFromHTSeq),该对象包含了所有输入数据,接下来将进行一系列的分析过程。

标准分析流程:
为了方便展示接下来的分析流程,这里只使用了部分的DEXSeqDataSet对象来分析:

DEXSeqDataSet对象核心的内容就是针对每个外显子的计数,该DEXSeqDataSet对象的每一行包含了给定的外显子的计数数据同时还有该基因其他外显子的计数数据。根据说明顺序依次打开看下:

7个样本,这里前7列是当前外显子的比对数,剩下7列是该基因其他外显子的比对数。
可以看到同一个基因的所有外显子总数相等,将第一行第一列的90加上第一行第八列的10971187,和第二行第一列与第二行第八列的和相等,因为该计数都是来自同一个样本。

这里可以看到DEXSeqDataSet对象的列的注释信息,根据计数方式和样本数量,将外显子分为了thisothers,分别为当前的exon区域的和该基因其他外显子的区域。这也对应了countdxd)的计数信息。

这里我们仅查看了当前基因外显子的reads比对数目情况,可见这里行的样本信息顺序同与之前sampleData列的样本顺序,这个列表其实就是把我们之前使用python_count.py脚本处理得到的信息综合到了一起。
最后我们详细看下用于统计的bin(区间)的信息:

信息一目了然,这里的exonBaseMean就是当前基因的外显子区域在7个样本上比对数目的均值,exonBaseVar就是7个统计数目的方差,后面的transcripts就是该基因1号外显子所在的转录本名称
最后可以通过sampleAnnotation函数查看到生成DEXSeqDataSet对象时使用的sampleData信息:

其实还有很多可以查看DEXSeqDataSet对象的函数,就不一一介绍了。

标准化:
针对不同的实验样本测序,我们会产生不同的测序深度,为了消除这种测序深度的不同导致的表达水平的偏差,我们这里评估实验样本的相对测序深度。这里使用的函数为:estimateSizeFactors

这里可以看到标准化后的DEXSeqDataSet对象会在注释信息列表多了一栏sizeFactor。前后比较下,可以发现到该数值和该样本的测序数据量是约成正比的。具体怎么计算而来的,我还不理解,

离散度评估:
为检测差异性exon usage,就要评估整体数据的变异情况。这里就有必要区分来自技术和生物学(背景噪音)所带来的exon usage差异。这里使用离散度来描述DEXSeqDataSet对象中由于生物学的复制所带来的噪音强度,所说的生物学复制指的就是同样处理条件下的样本重复。但是在RNA-seq实验中,同处理条件下的重复样本数目太小,由此来单个地评估外显子间的离散系数或者变异系数不可靠。因此,这里通过外显子和基因来获得变异信息。
这里仅讨论了单一的实验设计方式,从design公式信息上来看,相同状态(condition)下的样本就是生物学上的复制,那么knockdown条件有3个样本,control条件有4个样本。
DEXSeq使用了DESeq2评估离散度的方式。简而言之,针对每个外显子的离散度,首先Cox-Reid adjusted profile likelihood estimation计算每个外显子的离散度,接着根据每个个体的离散度值来构建离散度均值的关联模型,最后将该关联模型作为先验值来对每个外显子进行收缩评估。总之,和前面的estimateSizeFactors一样,都是统计学知识,我只能根据已有的统计知识进行尝试性理解翻译,肯定有不对的地方,欢迎大家指正。


这里根据绘制的离散度图我们应该就好理解estimateSizeFactorsestimateDispersinos的处理思想了,图示横坐标是标准化后的基因计数均值,纵坐标是离散度值。通俗点就是由于不同样本测序深度不同,那么我们就要根据测序深度的不同来消除不同深度所带来的差异;然后又由于不同实验条件下存在样本的重复例数,重复样本在每个外显子范围内的比对数目就会出现不同,当然单个外显子上的差异太小了,这里就放大到了整个基因水平的差异。我理解的离散度就是方差值,统计上会使用方差来表示离散度。为了消除由于样本重复所导致的差异,那么就要根据整体的离散度来构建一个模型,简单点说就是一个方程式Y=aX+b那种,然后用这个模型来对数据进行校正。

筛选差异性外显子表达量
经过测序深度和离散度的评估校正后,就到了最后检出环节。针对每一个基因,DEXSeq均会使用一个广义的线性模型公式(full formula):

这公式就是前面design的公式内容,那么我们将之与一个小的模型公式(reduced formula)进行比较:

这样我们就得到了conditionexon的交互项,在最开始的生成DEXSeqDataSet对象时,condition是具有两个水平的因子,knockdowncontrolexoncount bin也被分为了thisothers两类别。以上两个公式均针对每个exon的两个类别(thisothers)进行read计数统计。由于相同实验条件下的样本处理条件都是相同的,那么这里condition就是主要影响差异的因素了,针对每个样本使用condition水平对exon的计数统计再进行细分。
这里使用卡方分布来分别检测两个模型的差异性,根据p值,我们可以决定(reduced formula~sample+exon模型足够解释数据的差异性(根据sample不同来比较exon差异),还是存在(conditionexon)相互协同系数(full formula~sample+exon+condition:exon模型能更好的解释数据的差异性(根据condition的不同来比较exon差异),该模型意味着不同实验条件下相同基因外显子的比对数目差异显著。

函数testForDEU对每个基因每个外显子不同实验条件下比对数目执行卡方检验:

有时,我们会评估相对exon usage的倍数改变。计算exon usage倍数的改变是根据广义线性回归公式的系数来计算的:

个人理解,在线性回归中,系数值越大那么该自变量对应变量的影响也就越大,这里的系数就出现在交互项conditionexon这里,该函数计算了exon在不同condition下的usage倍数的改变。
需要注意的是在condition条件下exon usage的差异和gene expression的差异是不同的。例如考虑一个基因在不同condition下差异表达且包含有一个不同condition下差异表达的exon,根据模型公式的系数,它可能区分整体的gene expression效果,该效果的改变是来自该基因所有exon的比对数目;但是来自exon usage的作用,可以通过交互项conditionexon来获得,因为它单独地影响每一个外显子。因此这里使用了参数fitExpToVar=conditon”,通过conditon对外exon进行分类:

至此,所有结果都包含在了DEXSEqDataSet对象(dxd文件)中了,为了总结最终得到的结果,使用函数DEXSeqResults来获得对DEXSEqDataSet对象进行转换,经过DEXSeqResults转换得到的dxr1是一个数据框格式的子类别,这样就方便我们查看!
对于这些结果值的解释可以使用函数mcols来查询得到:



我的理解,其中的[5] "LRT statistic: full vs reduced"[6] "LRT p-value: full vs reduced"是用似然比来评估两个模型中哪个模型更适合当前数据分析的,[7] "BH adjusted p-values"根据后面的内容应该是两相关样本率的卡方检验值!

最后我们可以根据错误率来筛选选具有统计显著性的外显子区域个数:

这里可以得知padj的意思是假阳性率了,也就是错误率。同时还可以查看多少基因受到了影响:

这里可以通过函数plotMA来可视化查看基于比对到外显子数目的read来检出的差异性外显子表达量的情况。该图会将差异性显著的外显子加以红色表示,函数默认的调整后的p值阈值为0.1,该值可在使用函数时根据需要调整。

该图x轴对应的是每个外显子的平均的表达量,y轴是log2为底对表达量的倍数改变值(knockdown/control)求值,此图是将错误率小于0.1的外显子用红色表示。MA图就很简单了,感兴趣的可以查看维基百科(https://en.wikipedia.org/wiki/MA_plot )。





上一篇:读取affymetix的基因表达芯片数据-CEL格式
下一篇:DEXSeq使用说明(3/3)
苛求远离完美
回复

使用道具 举报

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

本版积分规则

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

GMT+8, 2019-10-22 22:58 , Processed in 0.038717 second(s), 27 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.