使用DEXSeq分析NGS数据中的exon表达差异

~ admin

对于RNA-seq，除了gene水平的差异分析外，还可以进行exon水平的差异分析。这时可以使用Bioconductor的DEXSeq软件包。

使用DEXSeq软件包，其输入为Count table。要生成这个Count table，与DESeq, edgeR类似的，都需要使用到htseq。DEXSeq提供了两个python角本，来调用htseq生成计数文件。

关于HTSeq的安装，请参考http://www-huber.embl.de/users/anders/HTSeq/doc/install.html

对于DEXSeq的安装，直接使用source(“http://bioconductor.org/biocLite.R”);biocLite(“DEXSeq”)就可以了。

第一步：计数
安装结束以后，使用system.file(“python_scripts”, package=”DEXSeq”)来找到DEXSeq提供的python角本的路径。

python ~/projects/RLib.3.01/DEXSeq/python_scripts/dexseq_prepare_annotation.py Homo_sapiens.GRCh37.75.fixed.gtf DEXSeq.hg19.gene.gff

这里需要注意的是，gtf文件必须与mapping的基因组一致，尤其是染色体的名字要一致，比如说如果mapping时有chr，gtf文件中的染色体一定也需要有chr。这里运行python是在terminal中，而不是R中。
而后使用dexseq_count.py来计数。

python ~/projectGREEN/RLib.3.01/DEXSeq/python_scripts/dexseq_count.py -p yes -s no -a 10 -f bam ~/DEXSeq/DEXSeq.hg19.gene.gff bam.file out.counts

在这里，参数-p指出mapping文件是否是pair end文件。参数-s表示是否是stranded，默认为yes。-f指输入文件的格式，默认为sam。
在运行计数结束之后，需要检查一下最后四行，看看empty的多不多，如果超过20%，可能需要检查一下mapping的结果，当然也可能是计数文件准备错误，比如mapping结果没有index等等。如果以上都不是，那可能是polyA太多了。

第二步：读入数据
在R中读入计数数据，需要准备好计数文件，实验设计，以及前面用到的gff文件。在这里，我们使用Bioconductor中已有的pasilla数据来示例。

> suppressPackageStartupMessages(library("DEXSeq"))
> inDir <- system.file("extdata", package="pasilla")
> countFiles <- list.files(inDir, pattern="fb.txt$", full.names=TRUE)
> gffFile <- list.files(inDir, pattern="gff$", full.names=TRUE) ##注意，如果是自己的数据的话，比如之前示例使用的是DEXSeq.hg19.gene.gff，这里就用DEXSeq.hg19.gene.gff
> ##实验设计
> sampleTable <- data.frame(row.names=c(paste("treated", 1:3, sep=""), paste("untreated", 1:4, sep="")),
+                           condition=rep(c("knockdown", "control"), c(3, 4)))
> 
> dxd <- DEXSeqDataSetFromHTSeq(
+         countFiles,
+         sampleData=sampleTable,
+         design= ˜ sample + exon + condition:exon,
+         flattenedfile=gffFile)

第三步：获得差异表达数据
只需要一步

> dxr <- DEXSeq(dxd)
> dxr

但是，我们似乎并不太明白它在背后做了些什么。于是我们一步一步的来查看它是如何做的。
同DESeq一样，它分为三个步骤：normalization, Dispersion estimation 以及 testing for differential exon usage。

> dxd <- estimateSizeFactors(dxd) #第一步
> dxd <- estimateDispersions(dxd) #第二步，此时可以使用plotDispEsts(dxd)来观察离散情况
> dxd <- testForDEU(dxd) #第三步
> dxd <- estimateExonFoldChanges(dxd, fitExptoVar="condition")
> dxr1 <- DEXSeqResults(dxd) #可以使用plotMA(dxr1)来查看结果

在第二步中，我们可以重新设计formula以符合实验要求，当然第三步也要随之改变。
对于结果dxr，可以直接视为data.frame来操作。也可以使用as.data.frame来转换它。结合使用plotDEXSeq就可以查看自己感兴趣的目标基因中的exon的表达情况。

发布者：admin

12 thoughts on “使用DEXSeq分析NGS数据中的exon表达差异”

  1. 谢谢分享教程，请教您一下，为何我在运行DEXSeqDataSetFromHTSeq时会报错呢。
    运行的结果如下：

    > suppressPackageStartupMessages( library( "DEXSeq" ) )
    > inDir = system.file("extdata",package="pasilla")
    > countFiles = list.files(inDir,pattern="count.txt$",full.name=TRUE)
    > countFiles
    [1] "/home/.../R/x86_64-pc-linux-gnu-library/3.0/pasilla/extdata/b2_dexseq_count.txt"
    [2] "/home/.../R/x86_64-pc-linux-gnu-library/3.0/pasilla/extdata/upf1a_dexseq_count.txt"
    > flattenedFile = list.files(inDir,pattern="mod.gff$",full.names=TRUE)
    > flattenedFile = list.files(inDir,pattern="mod.gff$",full.names=TRUE)
    > flattenedFile
    [1] "/home/tianmiao/R/x86_64-pc-linux-gnu-library/3.0/pasilla/extdata/DEXSeq_mod.gff"
    > sampleTable = data.frame(
    + row.names = c("upf1a","wtb2"),
    + condition = c("knockout","control"),
    + condition = c("knockout","control"),
    + libType = c("paired-end","paired-end"))
    > sampleTable
    condition condition.1 libType
    upf1a knockout knockout paired-end
    wtb2 control control paired-end
    > sampleTable = data.frame(
    + row.names = c("upf1a","wtb2"),
    + condition = c("knockout","control"),
    + libType = c("paired-end","paired-end"))
    > sampleTable
    condition libType
    upf1a knockout paired-end
    wtb2 control paired-end
    > dxd <- DEXSeqDataSetFromHTSeq(
    + countFiles,
    + sampleData=sampleTable,
    + design= ~ sample + exon + condition:exon,
    + flattenedfile= flattenedFile)
    Error: could not find function "DEXSeqDataSetFromHTSeq"

    回复

    1. 请检查你的DEXSeq的版本。如果是低于1.10.0，那请使用read.HTSeqCounts或者参与http://bioconductor.org/packages/2.13/bioc/vignettes/DEXSeq/inst/doc/DEXSeq.pdf。但是我建议还是升级至最新版为好。

      回复

      1. 你说的非常正确！
        我没有注意到我的R版本太低（3.0）,所以下载了对应的DEXSeq1.6版本。但是却是参考的是1.10版本的说明进行操作的，所以导致了上述错误。
        现已更新新的R，下载了最新的DEXSeq。
        非常感谢您的回复！

        回复

  2. 博主您好。我想问一下：
    我如果有不止2个组织的话，比如10个组织，想分析组织间是否存在差异表达。
    我能想的方法是：10个分别两两比较，若10个组织两两比较中有1个（或几个）存在DE，则说明这些组织在这个基因/exon存在差异。
    但总觉得这个方法不是很好，不知道有没有方法可以看多组织间的一次性比较？
    谢谢！

    回复

      1. 谢谢博主！ 昨天看了下DEXSeq，发现里面有可以一次针对多个样本。但最后只给出一个pvalue，无法体现这一组中两两比较的情况。所以我决定还是用最简单直接的两两比较，再判断是否有>=1对样本出现差异。

        不过还有个问题就是，比如5个样本，两两比较，会有10次pairwise，这样每个基因/exon会在10次比较中都得到一个pval或padj。这样得到的10个pvalue还需要再较正，以得到整体是否有差异的结论吗？ 谢谢！

        回复

          1. 好的，谢谢！我打算先按照最后得到的个数看是否要校正，个数多就校正一下，少的话就不校正了。

  3. 您好！看了您的文章感觉受益匪浅。我现在正在学习用DEXseq这个包。。但是我用了里面的那两个py script后跑出来 得出来的reads都在ambiguous_readpair_position 和 empty categories里面。。。我是使用的是蜜蜂（Apis mellifera）的基因组 并且都是在ensemble上面下载的。 但是由于似乎ensemble上面没有gtf文件 所以第一步的prepare annotation的那个script我就使用了从NCBI上下载的文件。然后得到的就是reads都在ambiguous_readpair_position 和 empty categories里面。
    求大神帮忙解答 如果没有gtf文件的话 怎么来用DEXseq 或者我可能什么地方出现了错误。
    不胜感激。

    回复

    1. 这种情况我遇到的比较少。我一般会检查一下assembly的版本是不是一样，染色体的名字的风格是不是一致。然后就是你有没有正确的使用DEXseq的python程序准备好完全没有overlaps的注释信息。

      回复

  4. 您好！根据您的步骤已经开始了DEXSeq分析，万分感谢！
    很不幸，在dxd <- estimateDispersions(dxd)遇到问题了，报错如下：
    using supplied model matrix
    Error in estimateDispersionsFit(object, fitType = fitType, quiet = quiet) :
    all gene-wise dispersion estimates are within 2 orders of magnitude
    from the minimum value, and so the standard curve fitting techniques will not work.
    One can instead use the gene-wise estimates as final estimates:
    dds <- estimateDispersionsGeneEst(dds)
    dispersions(dds) <- mcols(dds)$dispGeneEst
    …then continue with testing using nbinomWaldTest or nbinomLRT
    请问这步可以用DESeq代替吗？希望您能帮忙解决。万分感谢！！！

    回复

发表评论

电子邮件地址不会被公开。 必填项已用*标注