Understanding rMATS：基于 Junction Reads 的剪接Inclusion水平（PSI）估计

Understanding rMATS: 基于 Junction Reads 的剪接Inclusion水平（PSI/ψ）估计

.

rMATS（multivariate analysis of transcript splicing）是我们使用最多的一个alternative splicing events计算软件，其计算原理主要是基于RNA-Seq的junction reads计数，通过分层统计模型对每个选择性剪接事件的两类isoforms (inclusion vs exclusion)比例进行建模，从而在重复样本中准确地检测组间差异剪接事件。

.

1). junction reads对剪接事件的量化

rMATS通过统计跨越剪接位点的junction reads以及落在isoform特异性区域内的reads，来计算剪接包含水平（Percent Spliced In，PSI，常用ψ表示），这个值本质上反映的是特定条件下，某个外显子被“使用”的比例。

rMATS 首先从 RNA-Seq 数据中识别五类基本可变剪接事件：外显子跳跃（SE）、可变 5‘ 剪接位点（A5SS）、可变 3’ 剪接位点（A3SS）、互斥外显子（MXE）和内含子保留（RI）。

对于每一个事件，它关注两类互相竞争的isoforms:包含型（exon inclusion isoform，包含目标外显子或使用特定剪接位点）和排除型（exon skipping isoform，不包含）。

.

以最经典的exon skipping （SE）事件为例，一个SE事件涉及3个外显子：up-stream exon，skipping exon以及down-stream exon。

支持包含型的reads证据是两类junction reads，即一类是跨越 ‘up-stream exon，skipping exon’的连接点，另一类是跨越 ‘skipping exon，down-stream exon‘的连接点，跨这两类连接点的reads直接证明目标外显子被剪接进了成熟mRNA中。

而支持排除型的证据则是跨越‘up-stream exon，down-stream exon’连接点的junction reads，这证明生成的isoform跳过了中间的外显子。

除了junction reads，rMATS还可以统计exon body reads，即完全落在目标外显子(skipped exon)内部的reads，不跨越任何剪接位点。在JCEC模式下，这些reads也被计入包含型的计数中，但它们的特异性不如junction reads，因为一个落在外显子内部的reads，并不直接证明它来自于成熟mRNA。

有效长度校正

IJC（包含型 junction reads 数）和 SJC（排除型 junction reads 数）并不能直接计算**IJC/(IJC+SJC)**，因为exon inclusion isoform比exon skipping isoform多了目标外显子那一段序列，能容纳mapped reads的位置更多。即使两种剪接类型isoform的摩尔浓度完全相同，exon inclusion isoform也会比对到更多的reads，从而产生系统性偏差。

所以，rMATS 引入有效长度（Effective Length）来校正，对一个长度为L的序列片段，read长度为r，这样能容纳一个完整read的起始位置数为L-r+1，这样就是有效长度。

关于有效长度，如果前文尚未完全理解，在此详细解释：一个read要能完整地落在一段序列上，它的起始位置加上reads长度不能超过序列的末端，所以，如果read的起始位置为x，则该read覆盖的范围为[x，x+r-1]，要保证read完全落在序列内部，则需要满足x+r-1 <= L，转换一下，也就是x+r-1 <= L-r+1。也就是说，如果x是从1开始计算（第一个碱基位置），那么满足条件的起始位置就有x=1，2，3，……，L-r+1，总共有L-r+1个可能得起始位置。

.

现在设：

I：支持包含型的读段数

S：支持排除型的读段数

L_I：包含型的有效长度

L_S ：排除型的有效长度

c_I：包含型的摩尔浓度

c_S：包含型的摩尔浓度

read数与摩尔浓度、有效长度的关系为（正相关）：

我们真正想估计的剪接inclusion水平是摩尔浓度的比例：

但从read计数出发，直接得到的是：

这个值和 ψ 并不相等。为了还原 ψ，需要先把reads计数除以各自的有效长度，得到单位有效长度上的read密度：

然后用这两个密度来估计摩尔比例：

以上即为 rMATS 中 ψ 估计公式的来源。其本质是：先把长度带来的系统性偏倚除掉，再计算比例。

.

rMATS 将用有效长度校正的结果嵌入到二项分布的概率参数中，使得后续的统计检验是在校正后的摩尔比例尺度上进行的。这意味着组间比较的是真实的剪接差异，而不是长度差异带来的假像。

这里需要注意的一点是，rMATS程序运行中有一个参数–readLength需要设置，这个参数反映的是测序reads的长度，它是一个全局技术参数，与这里的inclusion和exclusion两种异构体isoform长度不是一回事，容易混淆。rMATS 会根据基因组注释文件自动读取每个外显子、内含子的坐标，然后根据剪接事件的具体结构，自动计算出两种isoform各自的有效长度。