ballgown是一個差異表達分析RNA-Seq數(shù)據(jù)的R包
對數(shù)據(jù)的要求：
1. RNA-Seqreads應(yīng)已比對到參考基因組上。
2.?轉(zhuǎn)錄組應(yīng)已經(jīng)組裝或下載參考轉(zhuǎn)錄組。
3.轉(zhuǎn)錄組中特征(轉(zhuǎn)錄本、外顯子和內(nèi)含子連接)的表達應(yīng)該處理成ballgown可讀格式。
兩個流程能生成ballgown所需的格式數(shù)據(jù)
1 TopHat2+Stringtie
2?pHat2+Cufflinks+Tablemaker
由Stringtie或?Tablemaker生成的Ballgown可讀的表達文件如下：
e_data.ctab:?外顯子水平表達值
i_data.ctab:內(nèi)顯子水平表達值
t_data.ctab:轉(zhuǎn)錄組水平表達值
e2t.ctab:表中有兩列，e_id和t_id，表示哪些外顯子屬于哪些轉(zhuǎn)錄本。這些id與e_data和t_data表中的id匹配。
i2t.ctab:表中有兩列，i_id和t_id，表示哪些內(nèi)含子屬于哪些轉(zhuǎn)錄本。這些id與i_data和t_data表中的id匹配。
ballgown的安裝：
在R的面板下執(zhí)行以下命令：
source("http://bioconductor.org/biocLite.R")
biocLite(
"ballgown")
導(dǎo)入數(shù)據(jù)到R中：
導(dǎo)入ballgown包
library(ballgown)
載入數(shù)據(jù)，并創(chuàng)建一個ballgown項目
儲存數(shù)據(jù)的文件夾名為：extdata
data_directory=system.file('extdata',?package='ballgown')
bg=ballgown(dataDir=data_directory,samplePattern='sample',meas='all')
提取外顯子，內(nèi)含子，轉(zhuǎn)錄本
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
例如：
提取轉(zhuǎn)錄本的表達，用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)
創(chuàng)建表型表格：
在差異表達分析之前，需要一個表格儲存樣本的表型信息，需要自己手動創(chuàng)建，一行一個樣本。
例如：
指定分組，及重復(fù)樣本數(shù)目：
pData(bg)?=data.frame(id=sampleNames(bg),?group=rep(c(1,0),?each=10))
phenotype_table=pData(bg)
差異表達分析：
stattest?能自動處理兩組比較(例如，病例/對照)、多組比較和“時間過程”比較。對于兩組和多組的比較，顯著的結(jié)果表明，該特征在至少一組中有差異表達。對于時間的比較，顯著的結(jié)果意味著特征的表達隨時間而顯著變化(即(連續(xù)協(xié)變量的值）。
1示例數(shù)據(jù)集bg包含兩個組標簽，0和1。我們可以用stattest?檢驗每個轉(zhuǎn)錄本在不同組之間的差異表達:
stat_results=stattest(bg,feature='transcript',?meas='FPKM',covariate='group')
結(jié)果如下：
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?檢驗每個轉(zhuǎn)錄本在時間刻度上的差異表達
pData(bg)=data.frame(pData(bg),time=rep(1:10,2))，timecourse_?results=stattest(bg,?feature='transcript',?meas='FPKM',?covariate='time',?timecourse=TRUE)
最近寫報告寫到炸~
