ballgown是一个差异表达分析RNA-Seq数据的R包
对数据的要求：
1. RNA-Seqreads应已比对到参考基因组上。
2. 转录组应已经组装或下载参考转录组。
3.转录组中特征(转录本、外显子和内含子连接)的表达应该处理成ballgown可读格式。
两个流程能生成ballgown所需的格式数据
1 TopHat2+Stringtie
2 pHat2+Cufflinks+Tablemaker
由Stringtie或 Tablemaker生成的Ballgown可读的表达文件如下：
e_data.ctab: 外显子水平表达值
i_data.ctab:内显子水平表达值
t_data.ctab:转录组水平表达值
e2t.ctab:表中有两列，e_id和t_id，表示哪些外显子属于哪些转录本。这些id与e_data和t_data表中的id匹配。
i2t.ctab:表中有两列，i_id和t_id，表示哪些内含子属于哪些转录本。这些id与i_data和t_data表中的id匹配。
ballgown的安装：
在R的面板下执行以下命令：
source("http://bioconductor.org/biocLite.R")
biocLite(
"ballgown")
导入数据到R中：
导入ballgown包
library(ballgown)
载入数据，并创建一个ballgown项目
储存数据的文件夹名为：extdata
data_directory=system.file('extdata', package='ballgown')
bg=ballgown(dataDir=data_directory,samplePattern='sample',meas='all')
提取外显子，内含子，转录本
structure(bg)$exon
structure(bg)$intron
structure(bg)$trans
提取表达值：
*expr(ballgown_object_name,<EXPRESSION_MEASUREMENT>)
* is either e for exon, ifor intron, t for transcript, or g for gene
例如：
提取转录本的表达，用FPKM值表示
transcript_fpkm=texpr(bg, 'FPKM')
transcript_cov=texpr(bg, 'cov')
whole_tx_table=texpr(bg, 'all')
exon_mcov=eexpr(bg, 'mcov')
junction_rcount=iexpr(bg)
whole_intron_table=iexpr(bg, 'all')
gene_expression=gexpr(bg)
创建表型表格：
在差异表达分析之前，需要一个表格储存样本的表型信息，需要自己手动创建，一行一个样本。
例如：
指定分组，及重复样本数目：
pData(bg) =data.frame(id=sampleNames(bg), group=rep(c(1,0), each=10))
phenotype_table=pData(bg)
差异表达分析：
stattest 能自动处理两组比较(例如，病例/对照)、多组比较和“时间过程”比较。对于两组和多组的比较，显著的结果表明，该特征在至少一组中有差异表达。对于时间的比较，显著的结果意味着特征的表达随时间而显著变化(即(连续协变量的值）。
1示例数据集bg包含两个组标签，0和1。我们可以用stattest 检验每个转录本在不同组之间的差异表达:
stat_results=stattest(bg,feature='transcript', meas='FPKM',covariate='group')
结果如下：
head(stat_results)
## feature id pval qval
## 1 transcript 10 0.01381576 0.10521233
## 2 transcript 25 0.26773622 0.79114975
## 3 transcript 35 0.01085070 0.08951825
## 4 transcript 41 0.47108019 0.90253747
## 5 transcript 45 0.08402948 0.48934813
## 6 transcript 67 0.27317385 0.79114975
2 用stattest 检验每个转录本在时间刻度上的差异表达
pData(bg)=data.frame(pData(bg),time=rep(1:10,2))，timecourse_ results=stattest(bg, feature='transcript', meas='FPKM', covariate='time', timecourse=TRUE)
最近写报告写到炸~
网友评论