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

  • A+
所属分类:Bioinformatics

对于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的表达情况。

原文来自:http://pgfe.umassmed.edu/ou/archives/3690

avatar

发表评论

:?: :razz: :sad: :evil: :!: :smile: :oops: :grin: :eek: :shock: :???: :cool: :lol: :mad: :twisted: :roll: :wink: :idea: :arrow: :neutral: :cry: :mrgreen:

目前评论:1   其中:访客  1   博主  0

    • avatar liuqs1990 0

      您好!看了您的文章感觉受益匪浅。我现在正在学习用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 或者我可能什么地方出现了错误。
      不胜感激。