作者：椰子糖
审稿：童蒙
编辑：amethyst

可变剪切能够产生多种类型的mRNA，因此一个基因就可以产生多种不同的蛋白。这个过程极大的增加了mRNA和蛋白质的多样性。可变剪切（alternative splicing）是一种后转录生物学过程，对细胞活动和疾病过程具有重要的且广泛的影响。研究表明人的基因组中有超过90-95%的多外显子基因存在可变剪切。到目前为止，也有很多软件可对其进行检测，今天我们就来了解一下这款常用可变剪切软件rMATS的最新版详情。

1. 软件介绍

rMATS是检测可变剪切事件的常用软件之一，其可以从RNA测序数据中，检测出多种类型的可变剪切事件，并提供了定量和组间差异分析的功能，可对生物学重复的样本进行组间分析。2020年6月更新的4.1版中更是对软件功能进行了完善：

1. 添加参数--task、--tmp等，以在不同的计算机上运行部分计算；
2. 添加参数--variable-read-length，能够允许不同长度的长度的reads进行分析；
3. 添加参数--paired-stats，进行成对统计分析；
4. 添加参数--novelSS, --mil, --mel，以检测新发可变剪切；
5. 输出文件中用fromGTF.novelJunction 和 fromGTF.novelSpliceSite 代替 fromGTF.novelEvents；
6. 版本兼容了python2和python3；
7. 在仅一个样本的组别或仅一个组别时，务必添加参数--statoff；
8. 修改了部分之前版本的bug。

软件网页链接：http://rnaseq-mats.sourceforge.net/

其检测的可变检测的事件类型如下：

2. 软件安装

rMATS turbo是rMATS的C/Cython版本。主要的差别在于速度和存储资源上，相比较rMATS turbo要快100倍，输出文件要小1000倍。具体可以参考文档：https://github.com/Xinglab/rmats-turbo/blob/v4.1.0/README.md，因此我们安装rMATS turbo。

安装依赖：Python (either 2.7 or 3.6),BLAS,LAPACK,GNU Scientific Library,GCC,gfortran,CMake等。保证以上依赖均存在的情况下就可以进行安装了。其实安装好conda，这些基础的包均已包括了。

conda create --name py2 python=2.7
conda activate py2
conda install -c bioconda rmats

安装好以后就可以进行软件测试啦。

3. 软件使用及测试

参数说明:

python rmats.py -h
usage: rmats.py [options]
optional arguments:
 -h, --help            show this help message and exit
 --version             show program's version number and exit
  --gtf GTF             An annotation of genes and transcripts in GTF format
 --b1 B1               A text file containing a comma separated list of the
                        BAM files for sample_1. (Only if using BAM)
 --b2 B2               A text file containing a comma separated list of the
                        BAM files for sample_2. (Only if using BAM)
  --s1 S1               A text file containing a comma separated list of the
                        FASTQ files for sample_1. If using paired reads the
                        format is ":" to separate pairs and "," to separate
                        replicates. (Only if using fastq)
  --s2 S2               A text file containing a comma separated list of the
                        FASTQ files for sample_2. If using paired reads the
                        format is ":" to separate pairs and "," to separate
                        replicates. (Only if using fastq)
  --od OD               The directory for final output
  --tmp TMP             The directory for intermediate output such as ".rmats"
                        files from the prep step
  -t {paired,single}    Type of read used in the analysis: either "paired" for
                        paired-end data or "single" for single-end data.
                        Default: paired
  --libType {fr-unstranded,fr-firststrand,fr-secondstrand}
                        Library type. Use fr-firststrand or fr-secondstrand
                        for strand-specific data. Default: fr-unstranded
  --readLength READLENGTH
                        The length of each read
  --variable-read-length
                        Allow reads with lengths that differ from --readLength
                        to be processed. --readLength will still be used to
                        determine IncFormLen and SkipFormLen
  --anchorLength ANCHORLENGTH
                        The anchor length. Default is 1
  --tophatAnchor TOPHATANCHOR
                        The "anchor length" or "overhang length" used in the
                        aligner. At least "anchor length" NT must be mapped to
                        each end of a given junction. The default is 6. (Only
                        if using fastq)
  --bi BINDEX           The directory name of the STAR binary indices (name of
                        the directory that contains the SA file). (Only if
                        using fastq)
  --nthread NTHREAD     The number of threads. The optimal number of threads
                        should be equal to the number of CPU cores. Default: 1
  --tstat TSTAT         The number of threads for the statistical model.
                        Default: 1
  --cstat CSTAT         The cutoff splicing difference. The cutoff used in the
                        null hypothesis test for differential splicing. The
                        default is 0.0001 for 0.01% difference. Valid: 0 <=
                        cutoff < 1. Does not apply to the paired stats model
  --task {prep,post,both,inte}
                        Specify which step(s) of rMATS to run. Default: both.
                        prep: preprocess BAMs and generate a .rmats file.
                        post: load .rmats file(s) into memory, detect and
                        count alternative splicing events, and calculate P
                        value (if not --statoff). both: prep + post. inte
                        (integrity): check that the BAM filenames recorded by
                        the prep task(s) match the BAM filenames for the
                        current command line
  --statoff             Skip the statistical analysis
  --paired-stats        Use the paired stats model
  --novelSS             Enable detection of novel splice sites (unannotated
                        splice sites). Default is no detection of novel splice
                        sites
  --mil MIL             Minimum Intron Length. Only impacts --novelSS
                        behavior. Default: 50
  --mel MEL             Maximum Exon Length. Only impacts --novelSS behavior.
                        Default: 500     

单个样本运行时

将NA12878的bam文件的具体路径写入到/path/to/b1.txt文件中

condadir/envs/py2/bin/python condadir/envs/py2/rMATS/rmats.py --nthread 4 --b1 /path/to/b1.txt --gtf Homo_sapiens.hg19_ucsc.gtf --od NA12878 -t paired --readLength 101 --libType fr-unstranded --statoff

--b1 为bam文件的路径，若有生物学重复则bam文件路径用逗号隔开，为单比较组时，仅给b1或者给s1即可；
--gtf 为已知的基因及转录本的gtf文件；--od 即为输出路径；-t 测序类型为单端或者双端 ;
--readLength 每条reads的长度，若长度不一致时，可使用--variable-read-length参数与readLength结合使用将reads截取到给定的数值；--libType 文库类型，可选择是否为链特异性；
--statoff 加上该参数则跳过统计部分，单样本或者单比较组时，跳过统计步骤。

比较组运行时

##/path/to/b1.txt
/path/to/1_1.bam,/path/to/1_2.bam
##/path/to/b2.txt
/path/to/2_1.bam,/path/to/2_2.bam
python rmats.py --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 --paired-stats

--b1 为组别1的bam文件的路径，若有生物学重复则bam文件路径用逗号隔开，为单比较组时，仅给b1或者给s1即可；
--b2 为组别2的bam文件的路径，若有生物学重复则bam文件路径用逗号隔开；
--gtf 为已知的基因及转录本的gtf文件；
--od 即为输出路径；
-t 测序类型为单端或者双端 ;
--readLength 若长度不一致时，可使用该参数将reads截取到给定的数值；
--libType 文库类型，可选择是否为链特异性；
--tmp 暂存目录；
--paired-stats 使用成对统计模型。

除了bam文件可作该软件的输入外，还可以使用fq文件做为输入，使用-s1和-s2参数即可，同一样本的双端reads使用冒号分隔，生物学重复间使用逗号分隔。

4. 结果说明

每一种可变剪切事件有相关的一系列的输出文件，每一种事件的相关文件以事件名作为前缀之一，以下文件中以[AS_Event]代替了[SE (skipped exon)，MXE (mutually exclusive exons)，A3SS (alternative 3' splice site)，A5SS (alternative 5' splice site)，RI (retained intron)] 中各事件：

• [AS_Event].MATS.JC.txt：检出的junction区域的reads数（Junction Counts）；
• [AS_Event].MATS.JCEC.txt：检出的junction区域的reads数（Junction Counts）和不跨越的外显子上read数（Exon Counts），考虑已知可变剪切事件时，可重点参考这个文件；
• fromGTF.[AS_Event].txt：从RNA和GTF中检出的所有可变剪切事件；
• fromGTF.novelJunction.[AS_Event].txt：仅使用RNA鉴定的可变剪切事件，与gtf的分析分离，其中并不包含未注释的可变剪切位点；
• fromGTF.novelSpliceSite.[AS_Event].txt：文件中仅包含未知的可变剪切位点的可变剪切事件，仅使用--novelSS参数时产生该文件；
• JC.raw.input.[AS_Event].txt：[AS_Event].MATS.JC.txt文件的input raw文件；
• JCEC.raw.input.[AS_Event].txt：[AS_Event].MATS.JCEC.txt文件的input raw文件。

01. 事件文件中共同的属性列

ID：rMATS 事件的ID；
GeneID：Gene ID；
geneSymbol：Gene 名称；
chr：染色体；
strand：基因的正负链情况；
IJC_SAMPLE_1：sample 1中包含剪切区域的reads数，生物学重复以逗号分隔；
SJC_SAMPLE_1：sample 1中不包含剪切区域的reads数，生物学重复以逗号分隔；
IJC_SAMPLE_2：sample 2中包含剪切区域的reads数，生物学重复以逗号分隔；
SJC_SAMPLE_2：sample 2中不包含剪切区域的reads数，生物学重复以逗号分隔；
IncFormLen：包含区域的长度，用于校正；
SkipFormLen：跳过区域的长度，用于校正；
PValue：两个比较组可变剪切差异的显著性（仅在使用statistical model时存在）；
FDR：由 p-value计算的错误发现率（仅在使用statistical model时存在）；
IncLevel1：由校正后reads数得到的sample 1的区域等级，生物学重复以逗号分隔；
IncLevel2：由校正后reads数得到的sample 2的区域等级，生物学重复以逗号分隔；
IncLevelDifference：average(IncLevel1) - average(IncLevel2)。

02. 事件文件中特异的属性列

• SE：exonStart_0base，exonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE
  包含形式中的目标外显子（该外显子的起始位置, 终止位置）
• MXE：1stExonStart_0base，1stExonEnd，2ndExonStart_0base，2ndExonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE
  +链，包含形式是包含第1个外显子（外显子的起始位置, 终止位置），跳跃第2个外显子
  -链，包含形式是包含第2个外显子（外显子的起始位置, 终止位置），跳跃第1个外显子
• A3SS, A5SS：longExonStart_0base，longExonEnd，shortES，shortEE，flankingES，flankingEE
  包含形式中使用长外显子（长外显子的起始位置, 终止位置）代替短的外显子（短外显子的起始位置 ，终止位置）
• RI：riExonStart_0base，riExonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE
  包含形式中包含内含子区域一般使用（上游外显子的终止位置 , 下有外显子的起始位置）

5. 总结

总体上说目前rMATS4.1版不受限于单双端测序，reads长度不一，是否存在生物学重复，是否有比较组，是否需要检测新转录本，是否链特异性等条件，并且其可以进行分步，分机器计算，功能完善，主要可变剪切事件检测完整的一款软件。在二代测序可变剪切检测的软件中可以算佼佼者，希望小编的介绍能给大家的可变剪切分析带来帮助。

6. 参考文献

• Mehmood A , Laiho A , Venlinen M S , et al. Systematic evaluation of differential splicing tools for RNA-seq studies[J]. Briefings in Bioinformatics, 2019.
• Shen S , Park J W , Lu Z , et al. rMATS: Robust and flexible detection of differential alternative splicing from replicate RNA-Seq data[J]. Proc Natl Acad Sci U S A, 2014, 111(51):5593-601.
• Park J W , Tokheim C , Shen S , et al. Identifying Differential Alternative Splicing Events from RNA Sequencing Data Using RNASeq-MATS[M]// Deep Sequencing Data Analysis. Humana Press, 2013.
• Shihao S , Won P J , Jian H , et al. MATS: a Bayesian framework for flexible detection of differential alternative splicing from RNA-Seq data[J]. Nucleic Acids Research, 2012(8):e61.
• http://rnaseq-mats.sourceforge.net/
• https://github.com/Xinglab/rmats-turbo/blob/v4.1.0