可变剪接(Alternative splicing):一个基因的外显子以不同的组合方式剪接形成不同的成熟RNA,由此产生的不同的 mRNA 可能被翻译成不同的蛋白质构体,因此,一个基因可能编码多种蛋白质。常见的可变剪接软件包括rMATS,Asprofile以及miso等。
rMATS介绍
rMATS是一个从RNA-Seq数据中检测差异选择性剪接事件的计算工具,根据RNA-Seq数据,rMATS可以自动检测和分析与所有主要类型的可变剪接模式相对应的可变剪接事件。rMATS可识别的可变剪切事件有5种。
skipped exon (SE),外显子跳跃,指一个或多个外显子连同其两端的内含子一起被剪切,在成熟mRNA中不存在。
alternative 5' splice site (A5SS),5’端可变剪接,它们的3’端剪接位点一致但5’端剪接位点不同,产生不同长度的5’端外显子。
alternative 5' splice site (A3SS),3’端可变剪接,它们的5’端剪接位点一致但3’端剪接位点不同,产生不同长度的3’端外显子。
mutually exclusive exons (MXE),外显子互斥,成熟的mRNA变体中,彼此特有的外显子,这些外显子不能同时出现在同一成熟mRNA中。
retained intron (RI),内含子保留,在一些转录本中内含子不会被剪切掉,保留在成熟的mRNA。
定量
rMATS采用exon inclusion level 来定义样本中可变剪切事件的表达量,以外显子跳跃(Skipped Exon)为例,正常的转录本称之为Exon Inclusion Isofrom, 发生了外显子跳跃的转录本则称之为Exon Skipping Isofrom。
用 I 表示比对到Exon Inclusion Isofrom上的reads,S表示比对到Exon Skipping Isofrom上的reads, 则该外显子跳跃的可变剪切事件比例可以表示为:
可以看到,exon inclusion level实际上是inclusion isofrom所占的比例,计算时,用长度校正了原始的reads数。其他类型的可变剪切事件也可以划分成上述两种isoform, 示意图如下
可以看到,rmats在计算isofrom的长度时,提供了两种方式,二者的区别就在于是否考虑跳过的exon的长度。
差异分析与统计检验
rmats 在差异分析时,比较的就是两组样本中inclusion level的差异,给定阈值c, 判断两个样本中对应inclusion level 的是否发生了变化,公式如下
c这个阈值通过--cstat参数自定义,取值范围为0-1,代表的是两个样本中inclusion level的差值,0.1表示两个样本中该可变剪切事件的inclusion level相差10%。当然,实际计算过程是非常繁琐的,需要考虑数据的分布,对应的统计模型等各种因素,最终会给出每个可变剪切事件的p值和多重假设检验校正后的FDR值。
rMATS安装
非root权限下安装rMATS总是出现各种问题,使用常规的conda安装貌似只能安装低版本的rMATS,运行的时候也是各种不兼容。幸得高手指点,在最新版rmats-turbo-4.1.1,有一个快速安装的方式:./build_rmats --conda(再次提醒阅读使用说明的重要性:))。如无意外,应该可以安装成功。
rMATS使用
!!!重要
使用./build_rmats --conda方式安装的rmats似乎必须在安装目录下使用./run_rmats --b1 /path/to/b1.txt --b2 /path/to/b2.txt --gtf /path/to/the.gtf -t paired --readLength 50 --nthread 4 --od /path/to/output --tmp /path/to/tmp_output进行运行,反而在b1.txt目录下是用绝对路径引用run_rmats(例如/path/to/run_rmats)进行运行总是提示:FileNotFoundError: [Errno 2] No such file or directory: 'b1.txt'
具体程序执行代码参考原文。
rMATS结果解读
在输出目录下,有很多的文件,我们重点关注其中两种文件即可:AS_Event.MATS.JC.txt, AS_Event.MATS.JCEC.txt。这里的AS_Event对应五种不同类型的可变剪切事件,每种类型是一个单独的文件,而JC和JCEC对应的是isoform effective length的两种计算方式。由于两种计算方式没有绝对的孰优孰劣的区分,根据需要进行选择。在这些文件中,包含了定量和差异的结果,其中InclevelDifference就是两组样本中表达量的差值,通过Pvalue和FDR可以对结果进行过滤和筛选。
针对exonStart_0base,exonEnd,upstreamES,upstreamEE,downstreamES,downstreamEE如下所示。