DEXSeq使用说明（1/3）

惠真-市二医院

DEXSeq使用说明（1/3）
最近在学习RNA-seq的分析了，事情最开始是跟着楼主的一个帖子《一个RNA-seq实战-超级简单-2小时搞定！》模仿学习，但是我在学习过程中了不知怎么搞的跑题了，本来是准备学习voom、DEseq2和edgeR包，结果现在来到了DEXSeq包。所以，不管那多，开始吧！
DEXSeq包是用来分析RNA-seq实验数据中exon usage的差异，这里exon usage的差异指的是由于实验条件导致的相对不同的exon usage。DEU（By differential exon usage）定义如下：

基本思想是：针对每个样本的外显子（或部分外显子），我们对比对到该外显子的read数和比对到同一基因（包含多个转录本）其他外显子的read数，然后计算这两个统计数的比值，最后根据不同实验条件下的比值变化情况，推导出相对的exon usage改变。对于基因内的一个外显子，该exon usage同步于该exon被剪接到转录组（可变剪接）的比例，它也包含了发生在转录本5 ‘端和3’端的可导致转录本边界的差异性的exon usage的可变剪接。因此，DEU（By differential exon usage）相比于可变剪接更直观。
最后有我们需要注意的两点：1，DEXSeq实际上并不是使用比值（比率）大小来计算分析，而是使用比值中的分子和分母的数目来计算分析差异的水平；2，DEXSeq并不局限于二个组的比较，它针对复杂的实验设计，使用广义线性模型来带入类似方差分析最终实现差异性的分析。
###实验条件====》转录本不同的剪接===》单个或多个外显子的表达差异
百度一图直观看看可变剪接和外显子的关系：


那么现在，就根据DEXSeq的说明书来翻译解释下：（http://bioconductor.org/packages/3.4/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.pdf ）


数据来源：
这里使用的RNA-seq数据源于文章：Conservation of an RNA regulatory map between Drosophila and mammals. 研究siRNA敲除基因pasilla的在果蝇S2-DRSC细胞的转录组中的影响，因为RNA结合蛋白pasilla蛋白被认为涉及到剪接的调节。这里共7个细胞组，处理组3个，对照组4个，均进行RNA测序，数据可在NCBI下载（GSE18508）。

参考基因组的比对：
分析的第一步是将reads比对到参考序列，这里需要注意的是比对的参考序列是基因组，而不是转录组，同时使用splice-aware的比对软件（例如，TopHat2，GSNAP或STAR）。

Python脚本查找：
这里使用的是DEXSeq包提供的两个Python脚本，这里需要注意的是，用户需要安装Python包HTSeq。安装完成后便可开始做分析前数据的准备了。


准备注释文件:
这里推荐使用的gtf文件是来自Ensembl的，同时确保该gtf文件使用的坐标系统和使用的参考基因组的坐标系统相同（这里最有保证的方式就是同样在Ensembl下载该物种的参考基因组序列）。
在gtf文件中，许多外显子会在多个转录本中出现。下图就是一个gtf文件经过DEXSeq的Python脚本dexseq_prepare_annotation.py转换后的gff文件截图，由图可见，多个转录本常会包含相同的外显子区间：

这里我们是分析差异exon usage，那么我们就将这些外显子所出现的转录本进行收集，定义一个计算外显子统计的bin，也就是一个区间范围，该区间范围针对了两种情况，一个外显子范围，就是该范围和多个转录本的外显子范围吻合；部分外显子范围，由于不同转录本具有不同的剪接范围，当这种情况出现时，就要针对不同的剪接范围设置一个外显子的范围（下图见文章Detecting differential usage of exons from、RNA-seq data）：

由图可见，该基因具有3个转录本，为了对所有外显子的区间范围进行区别统计，这里就设置了4个用于外显子统计的bin。
gtf文件转换命令如下：

在转换gff文件生成统计外显子bin的过程中，可能会遇到重叠的基因。假如两个基因在同一条链上，且外显子区间重叠了，该脚本会默认的将两个基因合为单个”aggregate gene”，这在随后涉及到单个基因的IDs时，会在前面加上一个加号（+）。如果不想使用这种处理办法，那么在转换时加上参数”-r no”，那么来自不同基因重合的外显子就会被跳过。

对reads进行计数：
这里使用的是python_count.py脚本：

该脚本是根据转化为gff格式的文件对给定的外显子区域做reads计数统计，使用方法很简单，但是针对不同的测序数据及应用我们需要选择适当的参数：

双末端测序reads：如果测序数据是采用的双端测序，那就要加上参数”-p yes”，告知采用的测序类型，因为默认是单端reads；此时前面得到的SAM文件需经过排序，可以是根据read名称排序或者是比对位置排序。这里我们可以根据之前的排序方式对该双端测序数据进行对应的排序选择，”-r pos”或”-r name”。双端测序read一定要经过排序，而单端测序read可以忽略。
文库链型：默认条件下”-s yes”，该程序认为你的文库是具有链特异性的，也就是说reads比对到的链和基因所在的链相同，”-s yes”；如果你的文库构建类型不支持链特异性，那么选择相应参数”-s no”；如果你的文库构建方式将reads比对到了基因的反向链上了，那么选择”-s reverse”。针对双端测序，默认是”-s yes”，那么测序得到的第一条序列应比对到基因所在的链，得到的第二条序列应比对到基因所在链的互补链。
这里百度解释下链特异性转录组测序：链特异性转录组测序合成cDNA第二条链时，加入的dNTPs试剂中用dUTP代替dTTP，使cDNA第二链中仅包含A/U/T/G。在PCR扩增前，用UNG酶（UNG酶的作用原理是选择性水解断裂含有dU的双链或单链DNA中的尿嘧啶糖苷键，形成的有缺失碱基的DNA链，在碱性介质以及高温下会进一步水解断裂，从而被消除）将cDNA第二链消化，从而使文库中仅包含cDNA第一链。
比对质量：默认为-a 10，可以根据自己需要选择。
这里打开一个比对后的结果：

该信息为：转录本：基因号：比对reads数目！

惠真-市二医院

感觉有很多知识理解错误，希望大家多多指正！