复现-jimmy教程之 外显子可视化
2019-11-22 本文已影响0人
小鹏_哒哒哒
外显子定量 可视化
DEXSeq包是用来分析RNA-seq实验数据中exon的差异；
这里exon的差异指的是由于实验条件导致的相对不同的exon，DEU（By differential exon usage）定义如下：
- 针对每个样本的外显子（或部分外显子），我们对比对到该外显子的read数和比对到同一基因（包含多个转录本）其他外显子的read数，
- 计算这两个统计数的比值，最后根据不同实验条件下的比值变化情况，推导出相对的exon usage改变。
- 对于基因内的一个外显子，该exon同步于该exon被剪接到转录组（可变剪接）的比例，它也包含了发生在转录本5 ‘端和3’端的可导致转录本边界的差异性的exon usage的可变剪接。因此，DEU（By differential exon usage）相比于可变剪接更直观。
Attention：
- DEXSeq实际上并不是使用比值（比率）大小来计算分析，而是使用比值中的分子和分母数目来分析差异的水平；
- DEXSeq并不局限于二个组的比较，它针对复杂的实验设计，使用广义线性模型来带入类似方差分析最终实现差异性的分析。
- 在转换gene.gtf生成gff过程中，可能会遇到重叠的基因。假如两个基因在同一条链上，且外显子区间重叠了，该脚本会默认的将两个基因合为单个”aggregate_gene”。
如果不想使用这种处理办法，那么在转换时加上参数 -r no，那么来自不同基因重合的外显子就会被跳过。
分析流程
1.文件准备
- 输入文件：
gene.gtf，sample.bam，sample.txt
- 输出文件：
sample.exon.counts.txt，diffgene&allgene_html- 软件：
python，HTSeq，DESeq
##Python
pip install HTSeq
---
##R
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DEXSeq")
BiocManager::install("pasilla")
2.分析
整理数据,准备3类输入文件
- 首先是counts矩阵：
inDir <- system.file("extdata", package="pasilla")
countFiles <- list.files(inDir, pattern="fb.txt$", full.names=TRUE)
##countFiles
Phvul.001G000400:001 88
Phvul.001G000400:002 2
Phvul.001G000400:005 15
## 简单看一下每个样本的数据要求，其实就是每个基因的每个外显子的reads数量
counts.txt
末尾冗余行太多的话建议删除
包自带python脚本可以定量，HTseq-count软件也可以
- gtf格式的基因组注释文件
gffFile <- list.files(inDir, pattern="gff$", full.names=TRUE)
##gffFile
Chr01 dexseq_prepare_annotation.py aggregate_gene 706 6715 . - . gene_id "Phvul.001G000400"
Chr01 dexseq_prepare_annotation.py exonic_part 706 1266 . - . transcripts "Phvul.001G000400"; exonic_part_number "001"; gene_id "Phvul.001G000400"
Chr01 dexseq_prepare_annotation.py exonic_part 1705 1945 . - . transcripts "Phvul.001G000400"; exonic_part_number "002"; gene_id "Phvul.001G000400"
##注意，如果是自己的数据的话，比如之前示例使用的是bean.gene.gtf，这里就是bean.gff
- 最后是分组文件
##实验设计
sampleTable <- read.delim("sample.group.txt",sep="\t",header=T,row.names=1)
sampleTable
## condition
## treated1 knockdown
## treated2 knockdown
## treated3 knockdown
## untreated1 control
## untreated2 control
## untreated3 control
## untreated4 control
构造DEXSeqDataSet对象
就是利用上面准备好的3类文件即可
dxd <- DEXSeqDataSetFromHTSeq(
countFiles,
sampleData=sampleTable,
design= ~sample + exon + condition:exon,
flattenedfile=gffFile
)
## converting counts to integer mode
dxd
## class: DEXSeqDataSet
## dim: 70463 14
## metadata(1): version
## assays(1): counts
## rownames(70463): FBgn0000003:E001 FBgn0000008:E001 ...
## FBgn0261575:E001 FBgn0261575:E002
## rowData names(5): featureID groupID exonBaseMean exonBaseVar
## transcripts
## colnames: NULL
## colData names(3): sample condition exon
对自己的数据，完全模仿这个形式即可
- 做差异分析，并解读结果,这一步较慢
dxr <- DEXSeq(dxd)
## Warning in MulticoreParam(workers = 1): MulticoreParam() not supported on
## Windows, use SnowParam()
## using supplied model matrix
## using supplied model matrix
## using supplied model matrix
## 一步法做差异分析，分步法做差异分析比较繁琐复杂。
###结果表头
## [1] "group/gene identifier"
## [2] "feature/exon identifier"
## [3] "mean of the counts across samples in each feature/exon"
## [4] "exon dispersion estimate"
## [5] "LRT statistic: full vs reduced"
## [6] "LRT p-value: full vs reduced"
## [7] "BH adjusted p-values"
## [8] "exon usage coefficient"
## [9] "exon usage coefficient"
## [10] "relative exon usage fold change"
## [11] "GRanges object of the coordinates of the exon/feature"
## [12] "matrix of integer counts, of each column containing a sample"
## [13] "list of transcripts overlapping with the exon"
- 指定基因进行可视化
plotDEXSeq( dxr, "FBgn0000210", legend=TRUE, cex.axis=1.2, cex=1.3,
lwd=2 )
## visualize the transcript models
plotDEXSeq( dxr, "FBgn0000210", displayTranscripts=TRUE, legend=TRUE,
cex.axis=1.2, cex=1.3, lwd=2 )
plotDEXSeq( dxr, "FBgn0000210", expression=FALSE, norCounts=TRUE,
legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )
plotDEXSeq( dxr, "FBgn0000210", expression=FALSE, splicing=TRUE,
legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )
##可视化全部是显著差异基因
##所有的显著性差异的基因的相关图片生成一个网页。
DEXSeqHTML( dxr, FDR=0.1, color=c("#FF000080", "#0000FF80") )
*# 并行分析，多差异分组
DEXSeq 分析的计算量可能很大，尤其是对于包含大量样本的数据集或包含具有大量外显子基因的基因组而言。
我们使用包BiocParallel，调用BPPARAM=函数的参数estimateDispersions()，testForDEU()以及estimateExonFoldChanges()：
BPPARAM = MultiCoreParam(4)
dxd = estimateSizeFactors( dxd )
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM)
- 表格数据导出
write.table(dxr, "DEXSeq_results.xls" , sep="\t", col.names=T, row.names = T)
个人示例
/public/cluster2/works/lipeng/RNA_seq/test/RNAedit/DEX/
image.png
python /path/to/library/DEXSeq/python_scripts/dexseq_prepare_annotation.py gene.gtf gene.gff
htseq-count -i exon_id -f bam -s reverse -r name sample.bam genome/gene.gff >sample.counts.txt
#R
suppressPackageStartupMessages(library("DEXSeq"))
suppressPackageStartupMessages(library("pasilla"))
countFiles <- list.files("./", pattern="counts.txt$", full.names=TRUE)
gffFile <- list.files("./", pattern="gff$", full.names=TRUE)
sampleTable <- read.table("samplegroup",sep="\t",header=T,row.names=1)
dxd <- DEXSeqDataSetFromHTSeq(
countFiles,
sampleData=sampleTable,
design= ~sample + exon + condition:exon,
flattenedfile=gffFile
)
###dxd <- estimateSizeFactors(dxd) #第一步
###dxd <- estimateDispersions(dxd) #第二步，此时可以使用plotDispEsts(dxd)来观察离散情况
###dxd <- testForDEU(dxd) #第三步,The function testForDEU performs these tests for each exon in each gene.
###dxd <- estimateExonFoldChanges(dxd, fitExptoVar="condition")
###dxr1 <- DEXSeqResults(dxd) #可以使用plotMA(dxr1)来查看结果
dxr <- DEXSeq(dxd)
plotDEXSeq( dxr, "FBgn0000210", expression=FALSE, norCounts=TRUE,legend=TRUE, cex.axis=1.2, cex=1.3, lwd=2 )
write.table(dxr, "DEXSeq_results.xls" , sep="\t", col.names=T, row.names = T)
DEXSeqHTML( dxr,FDR = 0.1, color=c("#FF000080", "#0000FF80") )
