欢迎关注”生信修炼手册”!
常见的转录组差异分析有两种策略，一种是基于raw count的定量方式，比如DESeq2, edgeR等；另外一种是基于FPKM/RPKM的定量方式，比如cuffdiff等。
在之前的文章中，我们也提到过基于FPKM值的pipeline由tophat+cufflinks+cuffdiff 升级更新为hisat + stringTie + ballgown。ballgown这个R包也是针对FPKM值的表达量进行差异分析，有两种方式可以得到转录本水平的FPKM值。
1. stringTie
为了方便下游的ballgown分析，在stringTie软件中直接添加-b参数就可以生成ballgown的输入文件，基本用法如下
stringtie -p 10 \
-G hg19.gtf \
-o output.gtf \
-b ballgown_out_dir -e \
align.sorted.bam
2. tablemaker
tablemaker软件通过调用cufflinks软件，也可以生成ballgown的输入文件，该软件可以从以下链接下载
https://figshare.com/articles/Tablemaker_Linux_Binary/1053137
基本用法如下
tablemaker \
-p 4 \
-q -W \
-G hg19.gtf \
-o out_dir \
align.sorted.bam
对于每个样本，都会生成一个文件夹，包含如下5个文件
e_data.ctab
e2t.ctab
i2t.ctab
i_data.ctab
t_data.ctab
e代表exon, i代表intron, t代表transcript，_data的文件为不同水平的表达量值。i2t表示intron和transcript之间的对应关系，e2t表示exon和transcript的对应关系。
输入文件准备好之后，就可以进行差异分析了。现在的R包都是高度封装的，几个函数就可以完成整套分析了。首先是读取所有样本的输入文件，代码如下
library(ballgown)
bg = ballgown(
samples = c("sampleA.dir", "sampleB.dir")，
meas='all')
samples 指定所有样本的ballgown的输入文件夹。导入成功之后，可以通过*expr函数在R中查看样本在不同水平的表达量信息， *的取值范围为i, e, t， g，代表不同水平。
查看转录本水平的表达量的代码示例如下
transcript_fpkm = texpr(bg, 'FPKM')
需要注意的是，intron, exon, transcript 这些水平的表达量信息在原本的ctab文件中都有，而gene水平的表达量信息，需要根据基因对应的转录本的表达量来计算，所以比较费时。
读取之后，需要设置样本分组, 代码如下
pData(bg) <- data.frame(
id=sampleNames(bg),
group=rep(c(1,0), each=3)
)
其实就是一个数据框，第一列为样本名称，第二列为样本对应的分组。
ballgown会自动根据group的种类进行不同类型的差异分析，如果样本分为两组，则进行两组间的差异分析，如果样本为多组，则进行多组间的差异分析。
ballgown通过stattest函数进行差异分析，支持以下4种水平的差异分析
- exon
- intron
- gene
- transcript
通过feature参数指定差异分析的水平。常规用法如下
# 转录本水平的差异分析
stat_results = stattest(bg,
feature='transcript',
meas='FPKM',
covariate='group')
# 基因水平的差异分析
stat_results = stattest(bg,
feature='gene',
meas='FPKM',
covariate='group')
ballgown还支持时间序列的差异分析，用法如下
pData(bg) <- data.frame(
pData(bg),
time=rep(1:10, 2)
)
results <- stattest(bg,
feature='transcript',
meas='FPKM',
covariate='time',
timecourse=TRUE
)
只需要添加timecourse=TRUE即可。ballgown还支持自定义差异分析的模型，更多用法可以参考官方文档。
·end·
—如果喜欢，快分享给你的朋友们吧—
扫描关注微信号，更多精彩内容等着你！
