利用球衣对转录本水平进行差异分析

欢迎关注”生信修炼手册”!

常见的转录组差异分析有两种策略，一种是基于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种水平的差异分析

1. exon
2. intron
3. gene
4. 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·

—如果喜欢，快分享给你的朋友们吧—