ballgown是一個(gè)差異表達(dá)分析RNA-Seq數(shù)據(jù)的R包
對(duì)數(shù)據(jù)的要求：
1. RNA-Seqreads應(yīng)已比對(duì)到參考基因組上好爬。
2.?轉(zhuǎn)錄組應(yīng)已經(jīng)組裝或下載參考轉(zhuǎn)錄組。
3.轉(zhuǎn)錄組中特征(轉(zhuǎn)錄本甥啄、外顯子和內(nèi)含子連接)的表達(dá)應(yīng)該處理成ballgown可讀格式存炮。
兩個(gè)流程能生成ballgown所需的格式數(shù)據(jù)
1 TopHat2+Stringtie
2?pHat2+Cufflinks+Tablemaker
由Stringtie或?Tablemaker生成的Ballgown可讀的表達(dá)文件如下：
e_data.ctab:?外顯子水平表達(dá)值
i_data.ctab:內(nèi)顯子水平表達(dá)值
t_data.ctab:轉(zhuǎn)錄組水平表達(dá)值
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)建一個(gè)ballgown項(xiàng)目
儲(chǔ)存數(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
提取表達(dá)值：
*expr(ballgown_object_name,<EXPRESSION_MEASUREMENT>)
*?is either e for exon, ifor intron, t for transcript, or g for gene
例如：
提取轉(zhuǎn)錄本的表達(dá)茴迁，用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)建表型表格：
在差異表達(dá)分析之前苛预，需要一個(gè)表格儲(chǔ)存樣本的表型信息，需要自己手動(dòng)創(chuàng)建笋熬，一行一個(gè)樣本热某。
例如：
指定分組，及重復(fù)樣本數(shù)目：
pData(bg)?=data.frame(id=sampleNames(bg),?group=rep(c(1,0),?each=10))
phenotype_table=pData(bg)
差異表達(dá)分析：
stattest?能自動(dòng)處理兩組比較(例如胳螟，病例/對(duì)照)昔馋、多組比較和“時(shí)間過程”比較。對(duì)于兩組和多組的比較糖耸，顯著的結(jié)果表明秘遏，該特征在至少一組中有差異表達(dá)。對(duì)于時(shí)間的比較嘉竟，顯著的結(jié)果意味著特征的表達(dá)隨時(shí)間而顯著變化(即(連續(xù)協(xié)變量的值）邦危。
1示例數(shù)據(jù)集bg包含兩個(gè)組標(biāo)簽，0和1舍扰。我們可以用stattest?檢驗(yàn)每個(gè)轉(zhuǎn)錄本在不同組之間的差異表達(dá):
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?檢驗(yàn)每個(gè)轉(zhuǎn)錄本在時(shí)間刻度上的差異表達(dá)
pData(bg)=data.frame(pData(bg),time=rep(1:10,2))倦蚪，timecourse_?results=stattest(bg,?feature='transcript',?meas='FPKM',?covariate='time',?timecourse=TRUE)
最近寫報(bào)告寫到炸~
