rMATS这款差异可变剪切分析软件的使用体验

今天在全国第8届生物信息学大会有幸听到了rMATS软件开发实验室PI的演讲，正好推送一波顾兄关于rMATS这款差异可变剪切分析软件的使用体验

PPT镇楼：

rMATS是一款对RNA-Seq数据进行差异可变剪切分析的软件。其通过rMATS统计模型对不同样本（有生物学重复的）进行可变剪切事件的表达定量，然后以likelihood-ratio test计算P value来表示两组样品在IncLevel（Inclusion Level）水平上的差异(从公式上来看，IncLevel跟PSI的定义也是类似的)，lncLevel并利用Benjamini Hochberg算法对p value进行校正得FDR值。rMATS可识别的可变剪切事件有5种，分别是skipped exon (SE)外显子跳跃，alternative 5' splice site (A5SS)第一个外显子可变剪切，alternative 3' splice site (A3SS)最后一个外显子可变剪切，mutually exclusive exons (MXE)外显子选择性跳跃和 retained intron (RI)内含子滞留，展现形式如下图（来自官网http://rnaseq-mats.sourceforge.net/index.html）

rMATS_splicing

软件下载及安装

rMATS最近刚现在出了rMATS 4.0.1版，相比之间的rMATS 3.2.5版，其用C，Python，Cython重写了该软件，运算速度提升了100倍，并且可支持多线程执行（明显感觉到计算速度的提升），并且新版的安装也简便好多了。PS.老版的rMATS我那时都是用bioconda安装的，不然太折腾了。。

进官网或者下述网站进行下载

https://sourceforge.net/projects/rnaseq-mats/files/MATS/rMATS.4.0.1.tgz/download

然后按照官网说明安装一些库文件以及python库（以ubuntu为例）

pip install numpy
sudo apt-get install libblas-dev liblapack-dev
sudo apt-get install gfortran

如果python的numpy装的有问题，可以使用bioconda来装下旧版的rMATS，其会顺便把numpy也装好，然后将其放置在环境变量中就行了（一般也不用这样）

如果运行时报错：error while loading shared libraries: libgsl.so.0: cannot open shared object file: No such file or directory，则说明缺少libgsl.so.0库文件，按照下述安装下就好

sudo apt-get install libgsl0-dev

软件的使用

软件使用也很简单，rMATS支持两种格式文件的输入。

第一种是fastq格式，那么在安装的时候还需要安装STAR比对软件以及提供比对的索引文件（STAR的索引文件异常的大），所以rMATS其实是建议使用第二种方式；

第二种是bam格式，rMATS支持其他比对软件比对后的结果bam文件作为输入，比如tophat等（那么hisat2也没啥问题，我试过），这样也能减少rMATS的运行时间。

新版rMATS下载解压后，你会发现有两个rmats.py执行脚本，这是由于rMATS v4.0.1 (turbo) was built with two different settings of Python interpreter，所以我们需要先测试下自己的系统支持那种，进入python，输入下述命令

>>> import sys
>>> print sys.maxunicode

如果出现1114111则说明需要使用rMATS-turbo-Linux-UCS4文件夹下rmats.py；如果出现65535则说明使用rMATS-turbo-Linux-UCS2文件夹下rmats.py

rMATS的参数设置不多，我一般使用以下几个，其他具体可参考官网

--b1 b1.txt 输入sample1的txt格式的文件，文件内以逗号分隔重复样本的bam文件名 --b2 b2.txt 输入sample2的txt格式的文件，文件内以逗号分隔重复样本的bam文件名 -t readType 双端测序则readType为paired，单端测序则为single --readLength 测序reads的长度 --gtf gtfFile 需要输入的gtf文件 --od outDir 所有输出文件的路径（文件夹） --nthread 设置线程数 --cstat The cutoff splicing difference. The cutoff used in the null hypothesis test for differential splicing（这个我一直不太理解是怎么卡的阈值，是在算法识别一些新的可变剪切的时候的差异性吗）

python rMATS-turbo-Linux-UCS4/rmats.py --b1 b1.txt --b2 b2.txt --gtf Homo_sapiens.Ensembl.GRCh37.72.gtf --od AS -t paired --readLength 151 --cstat 0.0001 --nthread 10

结果文件

rMATS的结果文件是以各个可变剪切事件的分布的，主要由AS_Event.MATS.JC.txt，AS_Event.MATS.JCEC.txt，fromGTF.AS_Event.txt，JC.raw.input.AS_Event.txt，JCEC.raw.input.AS_Event.txt这几类；其中JC和JCEC的区别在于前者考虑跨越剪切位点的reads，而后者不仅考虑前者的reads还考虑到只比对到第一张图中条纹的区域（也就是说没有跨越剪切位点的reads），但是我们一般使用JC的结果就够了（如果只是单纯的比较两组样品间可变剪切的差异的话）。

这几类文件中比较重要的要数S_Event.MATS.JC.txt，因为其他文件的信息差不多最终都整合在这个文件里面，以SE.MATS.JC.txt为例：

1-5列看列名就能懂其意思的，分别ID，GeneID，geneSymbol，chr，strand

6-11列分别为外显子的位置信息，分别为exonStart_0base，exonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE；网上有张图能很好的解释其含义，如下所示；其他可变剪切文件的这几列有点略微不同，但都可以类似的理解

skipped_exon.png

12列又为ID，不知道为啥重复再来一次，可能为了布局美观吧。。。

13-16列展示两组样品在inclusion junction（IJC）和skipping junction counts（SJC）下的count数，重复样本的结果以逗号分隔；列名分别为IJC_SAMPLE_1，SJC_SAMPLE_1，IJC_SAMPLE_2，SJC_SAMPLE_2，可以从下图来理解下

SE.png

下面几列信息rMATS认为是极为重要的信息：

• lncFormLen ：length of inclusion form, used for normalization
• SkipFormLen : length of skipping form, used for normalization
• P-Value : Significance of splicing difference between two sample groups（两组样品可变剪切的统计学显著差异指标）
• FDR ： False Discovery Rate calculated from p-value（对p-value的FDR校正）
• lncLevel1 ： inclusion level for SAMPLE_1 replicates (comma separated) calculated from normalized counts
• IncLevel2 : inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized counts
• IncLevelDifference ： average(IncLevel1) - average(IncLevel2)

lncFormLen和SkipFormLen分别是inclusion form和skipping form的有效长度值，虽然有计算公式但是还是要根据reads跨越时的具体情况来定的，具体解释可见https://groups.google.com/forum/#!topic/rmats-user-group/d7rzUBKXF1U(需访问外国网站。。。才能看哦)

IncLevel1可被认作是exon inclusion level（ψ），是Exon Inclusion Isoform在总（Exon Inclusion Isoform + Exon Skipping Isoform）所占比例

For a skipped exon, the exon inclusion level (denoted as percent spliced in or ψ) is calculated by the number of reads uniquely mapped to the exon inclusion isoform or the exon skipping isoform

因此计算公式如下：

ψ = (I/lI)/(S/lS)

• I 是指mapping到exon inclusion isoform上的reads数，对应结果表格中的IC_SAMPLE_1
• S 是指mapping到exon skipping isoform上的reads数，对应结果表格中的SC_SAMPLE_1
• lI 是指effective length of the exon inclusion isoform，对应结果表格中的IncFormLen，主要用于对I标准化，用于计算ψ
• lS 是指effective length of the exon skipping isoform，对应结果表格中的SkipFormLen，主要用于对S标准化，用于计算ψ

对于effective length的计算公式则是：

For a segment whose length is l, and the length of the read is r, the effective length of the segment is defined by the number of unique read intervals in this region, which is l − r + 1

这里的计算也可以参照https://www.biostars.org/p/256949/

IncLevelDifference则是指两组样本IncLevel的差异，如果一组内多个样本，那么则是各自组的均值之间差值

rMATS采用的是Likelihood-ratio test来比较两组样本可变剪切是否有差异

If the maximum-likelihood estimations (MLEs) of ψi1,ψi2 have a difference smaller than or equal to the user-defined cutoff (i.e., |ψi1 −ψi2|≤ c), we set the P value to be 1

所以跟我们参数--cstat设置有关（默认是0.0001），基于marginal distribution（在零假设前提下，近似χ2 distribution），备择假设是|ψi1 −ψi2|> c，因此粗略的理解就是ψi1和ψi2差异比c越大P值就越小

rMATS可视化-rmats2sashimiplot

rmats2sashimiplot这款专门用来对rMATS分析结果做可视化的软件，最近尝试用了下，操作比较简洁

rmats2sashimiplot，网址：rmats2sashimiplot

下载

git clone https://github.com/Xinglab/rmats2sashimiplot

or

https://pypi.org/project/rmats2sashimiplot/2.0.2/#files

安装

python setup.py install

ubuntu 16.04 报错：error while loading shared libraries: libgsl.so.0: cannot open shared object file: No such file or directory

这是因为ubuntu 16.04不支持libgsl0ldbl，而是替换为libgsl2，所以先安装libgsl2，然后在将其软链接/复制到/usr/lib目录下

sudo apt-get install libgsl-dev
dpkg -L libgsl2
sudo cp /usr/lib/x86_64-linux-gnu/libgsl.so.19 /usr/lib/libgsl.so.0

解决方法参照自：https://jingyan.baidu.com/article/63acb44a13cc7d61fdc17e6d.html

如果直接使用rMATS从fastq输入文件到结果输出文件，那么rMATS是不保留中间的bam文件的

后来看了下rmats.py中的代码，发现只要将shutil.rmtree(args.tmp)注释掉就会保留STAR比对的bam结果文件了，但是其是放在服务器上/tmp目录中，也不太方便

而rmats2sashimiplot可视化则需要bam文件作为输入，所以需要我们先用STAR比对得到bam文件再用rMATS做差异可变剪切分析，如下所示：

1. 用STAR对各个样本做比对生成bam文件，比对参数参考rMATS软件调用STAR时所用参数，其实就是比默认比对参数多了： --chimSegmentMin 2 --outFilterMismatchNmax 3 --outSAMstrandField intronMotif --alignEndsType EndToEnd --alignSJDBoverhangMin 6 --alignIntronMax 299999 STAR --readFilesIn A1_1.clean.fq A1_2.clean.fq --chimSegmentMin 2 --outFilterMismatchNmax 3 --alignEndsType EndToEnd --runThreadN 12 --outSAMstrandField intronMotif --outSAMtype BAM SortedByCoordinate --alignSJDBoverhangMin 6 --alignIntronMax 299999 --genomeDir ~/index/STAR/Ensembl/human/GRCh38.92/ --sjdbGTFfile ~/annotation/Ensembl/human/Homo_sapiens.GRCh38.92.chr.gtf --outFileNamePrefix A1
2. rMATS做差异可变剪切分析（例如A.txt包含了A组bam的文件名） python rmats.py --b1 A.txt --b2 B.txt --gtf Homo_sapiens.GRCh38.92.chr.gtf --od ./test -t paired --readLength 150 --cstat 0.0001 --nthread 10

rmats2sashimiplot作图

rmats2sashimiplot --b1 A1.bam,A2.bam,A3.bam --b2 B1.bam,B2.bam,B3.bam -t SE -e ./SE.MATS.JC_top20.txt --l1 A --l2 B --exon_s 1 --intron_s 1 -o SE_plot

rmats2sashimiplot的注意事项：

• 可以预先用samtools对bam文件建索引，不然rmats2sashimiplot也会先建索引（比较慢）
• rMATS的结果文件中，如SE.MATS.JC.txt中染色体都是默认加上chr的，而bam文件中的染色体根据基因组注释文件来源不同，有些的没chr，而是直接以数字表示，这时需要预先将两者保持一致
• rmats2sashimiplot默认是出PDF图片的，如果需要出png，则需要自己修改下源码（由于我也不懂python作图，所以用了最无脑的方式就是将所有跟出图相关的脚本中pdf都替换为png。。。）

rMATS可视化的图片如下：

rmats2sashimiplot

参考资料： https://www.biostars.org/p/256949/#274012 http://rnaseq-mats.sourceforge.net/user_guide.htm#as_events