rMATS差异可变剪切分析

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.jpg

软件下载及安装

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(需翻墙。。。才能看哦)

lncLevel1和IncLevel2分别为sample1和sample2的inclusion level(粗略可理解为PSI),从公式上看比较好理解

ψ = (I/LI) / (I/LI + S/LS)

ψ = Inclusion Level(sample1)I = IJC_SAMPLE_1S = SJC_SAMPLE_1LI = lncFormLenLS = SkipFormLen

从公式上理解为lncLevel是在最终的成熟mRNA中,这个respective exon(SE事件中则是对应那个被跳跃的exon)出现的平均频率是多少,或者说所占的比例?也可以认为是read counts(标准化后的)在可变剪切的事件中的各自exon上所占的比例?简单的说就是lncLevel越小说明出现外显子的跳跃的比例越高(这个是以SE文件为例的)???可能我还理解的不是很准确。。。。

最后还有一个软件能绘制可变剪切的图片的软件rmats2sashimiplot

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

本文出自于http://www.bioinfo-scrounger.com转载请注明出处