rmats是一款比较好的分析可变剪切的分析软件，将rmats分析的结果，利用rmats2sashimiplot绘画出的图也是比较漂亮的，

可变剪切事件：

可变剪切事件可以分为：

image.png

我们就讲一下外显子跳跃事件，

image.png
如下图所示：
image.png

上面的那个属于Exon lnclusion isoform的情况。
下面那个属于Exon Skipping isoform的情况
外显子跳跃后，中间的外显子就没有了，两边的连到一起

那么，测序的reads是如何捕捉到的转录本怎么和可变剪切事件的呢？

上图中左边为inclusion isofrom, 右边为skipping isoform, 比对上inclusion isoform的reads 有两种，分别为红色和紫色，其中红色reads只比对到了1号和3号exons上，无法区分这部分reads是来自inclusion isofrom还是skipping reads, 而紫色的reads 比对到了exon2以及exon2与其他exon的连接处，这部分reads可以明确来自inclusion isoform. 所以在统计inclusion isofrom count时，只会计数这部分reads（紫色reads）, 同理，统计skipping isoform的reads数时，只会统计绿色部分的reads。

结果解读：

那么结果该如何解读呢？

image.png

其中RPKM衡量表达量情况，lnclevel是什么呢？

文献中给的定义是：
1.lncLevel1 ： inclusion level for SAMPLE_1 replicates (comma separated) calculated from normalized counts（对于sample1）
2.IncLevel2 : inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized counts（对于sample2）
3.IncLevelDifference ： average(IncLevel1) - average(IncLevel2)

回到刚才，那么
lnclevel是exon inclusion level，是Exon Inclusion Isoform在总（Exon Inclusion Isoform + Exon Skipping Isoform）所占比例

lnclevel的计算公式如下：
ψ = (I/lI)/[(S/lS) + (I/II)]

其中：

1. I是指mapping到exon inclusion isoform的reads数
2. S指mapping到exon skipping isoform的reads数
3. II指exon inclusion isoform的有效长度
4. IS 指exon skipping isoform的有效长度

那么有效长度的定义又是什么呢？？？

image.png

刚才我们介绍了rMATS只会记录连接处（紫色和绿色）的reads
以exon skipping为例，每一个连接处都会有8bp左右的anchor，作者默认为8bp，因此有效比对到两边exon的长度为 r-2a+1（加 1 是为了防止0值）

因此可以定义出 junction length ：the overall junction region covered by reads across junctions, affected by read length, anchor length (by default 8 bps in both upstream and downstream exons) and exon length
junction length is calcualted as (read length - anchor)*2
假设 read length = 150np；anchor = 8bp，则 junction length = 284bp，skipping的有效长度为135bp

rMATs是根据Likelihood-ratio test来比较两组样本可变剪切是否有差异

image.png

所以，我们根据这个就可以看两个样品是否有可便剪切上的差异了
那么，当Exon inclusion isoform数量较高时，lnclevel1较高；当Exon skipping isoform数量较高时，lnclevel2较高。
如果Exon inclusion isoform数量低的离谱，Exon skipping isoform有一定数量且保持不变那么lnclevel1会接近于0；反之，Exon skipping isoform也一样。

但是，我们在利用miso或rmats分析可变剪切的时候，往往用sashimiplot进行可视化，其中基于两款软件画图的reads可能会不一样，那么：

对于rmats2sashimiplot产生的图片，我们不需要太过于纠结reads数目和rmats对应不上，只需要看reads分布的大体趋势即可，重点是关注不同样本中IncLevel的差异。

这是因为sashimiplot的input是miso的输出文件，而rmats的output会被sashimiplot转换成miso的output，才可以继续作图
在转换的过程中由于rmats和miso判断可变剪切的原理不一样，所以会造成一定误差

判断差异可变剪切的统计学模型

由上面的叙述我们知道判断某个基因是否存在差异可变剪切是通过 | Ψi1-Ψi2 | ≤ c，其中 i 表示第 i 个基因；Ψ 即为 lnclevel。
首先rmats先计算全部基因的 | Ψi1-Ψi2 | ，并以此拟合出一个分布：

上面这张图的横坐标表示 | Ψi1-Ψi2 | ，纵坐标表示频率（频数）分布直方图，也就是表示当 | Ψi1-Ψi2 | 为某值是，对应一共有多少个基因满足这个条件
那么rmats里面有一个参数是 --cstat，默认为0.01%，也就是说在全部的基因中，存在差异可变剪切的只占0.01%。因此统计全部基因的 | Ψi1-Ψi2 | 值，选择差异最大的那0.01%，并定义为存在差异可变剪切的基因，那么p_val即为上图阴影部分的面积

如何理解rMAT结果

前几列代表的是基因的基本信息：

1. chr: 代表染色体
2. strand: 代表链的信息

而 exonStart_0base ,exonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE代表的是可变剪切的isoform的位置信息：

exonStart_0base ,exonEnd代表skip的外显子

而：

前面两列代表显著性，而lncLevels1和lncLevels2分别对应这不同处理中各个生物学重复的lncLevels值，而lncLevelDifference代表lncLevels1与lncLevels2的差值，若差值很大并且显著则判断为差异可变剪切

参考：https://cloud.tencent.com/developer/article/1366294
参考：https://www.biostars.org/p/256949/
参考：https://www.jianshu.com/p/d2f00c1c9067