使用rmats进行可变剪切的分析

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

rmats是目前使用的最广泛的可变剪切分析软件，该软件不仅可以识别可变剪切事件，还提供了定量和组间差异分析的功能，功能强大，网站链接如下

http://rnaseq-mats.sourceforge.net/rmats4.0.2/index.html

该软件前后经历了多个版本，目前最新版本为v4.0.2, 相比之前的版本，v4.0之后的版本在运行速度，内存消耗，磁盘占用等方面进行了优化，最明显的就是运行速度，比之前快了100多倍。

安装也很简单，直接下载解压缩即可使用，这里不赘述。rmats可以识别以下五种类型的可变剪切事件

软件的基本用法如下

python rmats.py  \
--b1 b1.txt --b2 b2.txt \
--gtf ref.transcript.gtf \
--od out_dir \
-t paired \
--readLength 101 \
--cstat 0.1 \
--libType fr-unstranded

b1.txt中保存的是每个样本比对参考基因组的bam文件的路径，示例如下

/bams/rep1.bam,/bams/rep2.bam

这种用法从bam文件开始，实用性更强，除此之外，也支持从fastq文件开始，用法如下

python rmats.py \
--s1 s1.txt --s2 s2.txt \
--gtf ref.transcript.gtf \
--bi /STARindex/hg19 \
--od out_dir \
-t paired \
--nthread 6 \
--readLength 151

S1.txt中保存的是每个样本fastq文件的路径，rmats会自动调用STAR进行比对，bi参数指定参考基因组STAR的索引，更多参数和细节请参考官方文档。

rmats中核心功能就是定量和差异分析，解释如下

1. 定量

rmats采用exon inclusion level 来定义样本中可变剪切事件的表达量，以外显子跳跃为例，正常的isoform称之为Exon Inclusion Isofrom, 发生了外显子跳跃的转录本称之为Exon Skipping Isofrom, 示意如下

比对到inclusion isoform上的reads用I表示，比对到skipping isoform上的reads用S表示, 则该外显子跳跃的可变剪切事件的表达量如下

可以看到，exon inclusion level实际上是inclusion isofrom所占的比例，计算时，用长度校正了原始的reads数。其他类型的可变剪切事件也可以划分成上述两种isoform, 示意图如下

可以看到，rmats在计算isofrom的长度时，提供了两种方式，二者的区别就在于是否考虑跳过的exon的长度，详细的公式在上图中共也有给出。

2. 差异分析

rmats 在差异分析时，比较的就是两组样本中inclusion level的差异，给定阈值c, 判断两个样本中对应inclusion level 的是否发生了变化，公式如下

c这个阈值通过--cstat参数自定义，取值范围为0-1，代表的是两个样本中inclusion level的差值，0.1表示两个样本中该可变剪切事件的inclusion level相差10%。当然，实际计算过程是非常繁琐的，需要考虑数据的分布，对应的统计模型等各种因素，最终会给出每个可变剪切事件的p值和多重假设检验校正后的FDR值。

在输出目录下，有很多的文件，我们重点关注其中两种文件即可。

1. AS_Event.MATS.JC.txt
2. AS_Event.MATS.JCEC.txt

这里的AS_Event对应五种不同类型的可变剪切事件，每种类型是一个单独的文件，而JC和JCEC对应的是isoform effective length的两种计算方式。由于两种计算方式没有绝对的孰优孰劣的区分，所以同时给出两种结果，在这些文件中，包含了定量和差异的结果

IJC表示inclusion isoform counts, SJC表示是skipping isoform counts, 生物学重复样本用逗号分隔；IncFormLen代表effective inclusion isoform length, SkipFormLen代表effective inclusion isoform length；lencLevel代表定量的结果，InclevelDifference就是两组样本中表达量的差值，通过Pvalue和FDR可以对结果进行过滤和筛选。

除了定量和差异的结果，还给出了每种可变剪切事件对应的exon的区间信息，示意如下

以上示例是外显子跳跃中的区间信息，其他类型的表头会有一定差异，但是意思是一样的。

rmats中，识别可变剪切是以exon为单位的，只需要比较邻近的3到4个exon的表达情况，就可以确定一个可变剪切事件是否发生，这个思路从可变剪切最核心的地方，即exon的变化出发进行分析，直接有效，但是由于其对问题的高度抽象和简化，导致对应的结果看起来不够直观。

·end·

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