外顯子定量 可視化
DEXSeq包是用來分析RNA-seq實驗數(shù)據(jù)中exon的差異钉凌；
這里exon的差異指的是由于實驗條件導致的相對不同的exon娜扇，DEU（By differential exon usage）定義如下：
- 針對每個樣本的外顯子（或部分外顯子）跑筝，我們對比對到該外顯子的read數(shù)和比對到同一基因（包含多個轉(zhuǎn)錄本）其他外顯子的read數(shù)，
- 計算這兩個統(tǒng)計數(shù)的比值滚局，最后根據(jù)不同實驗條件下的比值變化情況工猜，推導出相對的exon usage改變摊趾。
- 對于基因內(nèi)的一個外顯子币狠，該exon同步于該exon被剪接到轉(zhuǎn)錄組（可變剪接）的比例，它也包含了發(fā)生在轉(zhuǎn)錄本5 ‘端和3’端的可導致轉(zhuǎn)錄本邊界的差異性的exon usage的可變剪接砾层。因此漩绵，DEU（By differential exon usage）相比于可變剪接更直觀。
Attention：
- DEXSeq實際上并不是使用比值（比率）大小來計算分析肛炮，而是使用比值中的分子和分母數(shù)目來分析差異的水平止吐；
- DEXSeq并不局限于二個組的比較，它針對復雜的實驗設(shè)計侨糟，使用廣義線性模型來帶入類似方差分析最終實現(xiàn)差異性的分析碍扔。
- 在轉(zhuǎn)換gene.gtf生成gff過程中，可能會遇到重疊的基因粟害。假如兩個基因在同一條鏈上蕴忆，且外顯子區(qū)間重疊了颤芬，該腳本會默認的將兩個基因合為單個”aggregate_gene”悲幅。
如果不想使用這種處理辦法，那么在轉(zhuǎn)換時加上參數(shù) -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.分析
整理數(shù)據(jù),準備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
## 簡單看一下每個樣本的數(shù)據(jù)要求聚蝶，其實就是每個基因的每個外顯子的reads數(shù)量
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"
##注意，如果是自己的數(shù)據(jù)的話碘勉，比如之前示例使用的是bean.gene.gtf巷挥，這里就是bean.gff
- 最后是分組文件
##實驗設(shè)計
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
構(gòu)造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
對自己的數(shù)據(jù)，完全模仿這個形式即可
- 做差異分析验靡，并解讀結(jié)果,這一步較慢
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
## 一步法做差異分析倍宾，分步法做差異分析比較繁瑣復雜。
###結(jié)果表頭
## [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 )
##可視化全部是顯著差異基因
##所有的顯著性差異的基因的相關(guān)圖片生成一個網(wǎng)頁胜嗓。
DEXSeqHTML( dxr, FDR=0.1, color=c("#FF000080", "#0000FF80") )
*# 并行分析高职，多差異分組
DEXSeq 分析的計算量可能很大，尤其是對于包含大量樣本的數(shù)據(jù)集或包含具有大量外顯子基因的基因組而言辞州。
我們使用包BiocParallel怔锌，調(diào)用BPPARAM=函數(shù)的參數(shù)estimateDispersions()，testForDEU()以及estimateExonFoldChanges()：
BPPARAM = MultiCoreParam(4)
dxd = estimateSizeFactors( dxd )
dxd = estimateDispersions( dxd, BPPARAM=BPPARAM)
dxd = testForDEU( dxd, BPPARAM=BPPARAM)
dxd = estimateExonFoldChanges(dxd, BPPARAM=BPPARAM)
- 表格數(shù)據(jù)導出
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)來查看結(jié)果
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") )
