软件13 —— rMATS

一、 基本介绍

rMATS是一款对RNA-Seq数据进行差异可变剪切分析的软件。其通过rMATS统计模型对不同样本（有生物学重复的）进行可变剪切事件的表达定量，然后以likelihood-ratio test计算P value来表示两组样品在IncLevel（Inclusion Level）水平上的差异，并利用Benjamini Hochberg算法对p value进行校正得FDR值。（lncLevel是Exon Inclusion Isoform在总Isoform（Exon Inclusion Isoform + Exon Skipping Isoform）中所占比例。）

rMATS可识别的可变剪切事件有5种，分别是外显子跳跃（skipped exon, SE），选择性5’剪接位点（alternative 5' splice site, A5SS），选择性3’剪接位点（alternative 3' splice site, A3SS），外显子互斥（mutually exclusive exons, MXE）和 内含子保留（retained intron, RI）。

官方文档：
https://github.com/Xinglab/rmats-turbo/blob/v4.1.2/README.md

二、 使用方法

(1) 输入

rMATS支持两种格式文件的输入。

• 第一种是fastq格式，那么在安装的时候还需要安装STAR比对软件以及提供比对的索引文件（STAR的索引文件异常的大），所以rMATS其实是建议使用第二种方式；
• 第二种是bam格式，rMATS支持其他比对软件比对后的结果bam文件作为输入，比如tophat/ hisat2等，这样也能减少rMATS的运行时间。

(2) 方法

目前rMATS 4.1版不受限于单双端测序，reads长度不一，是否存在生物学重复，是否有比较组，是否需要检测新转录本，是否链特异性等条件，并且其可以进行分步，分机器计算，功能完善，主要可变剪切事件检测完整的一款软件。

$ conda create -n bioinf3
$ conda activate bioinf3  #python 3.9
$ conda install -c bioconda rmats=4.1.2
$ conda install -c bioconda rmats2sashimiplot
$ rmats.py –version
$ rmats2sashimiplot --help

构建索引：

$ for i in `ls ./04_bam/*.bam` ; do (samtools index $i -@ 4) ; done  &

运行rmats.py：

$ rmats.py  --b1 ./06_rmats/LF.txt --b2 ./06_rmats/CF.txt  --gtf ./annotation/genome/Drosophila_melanogaster.BDGP6.32.105.gtf  --od ./06_rmats/  -t paired  --readLength 150  --cstat 0.0001  --nthread 4  --libType fr-firststrand  --tmp ./06_rmats/tmp/  &

(3) 参数

--b1 b1.txt：输入sample1的txt格式的文件，文件内以逗号分隔重复样本的bam文件名（/path/to/1_1.bam,/path/to/1_2.bam）
--b2 b2.txt：输入sample2的txt格式的文件，文件内以逗号分隔重复样本的bam文件名（/path/to/2_1.bam,/path/to/2_2.bam）
-t {paired,single}：双端测序则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. The default is 0.0001 for 0.01% difference. Valid: 0 ≤ cutoff < 1. Does not apply to the paired stats model. 代表的是两个样本中inclusion level的差值
--libType {fr-unstranded,fr-firststrand,fr-secondstrand}：Library type. Default is unstranded (fr-unstranded). Use fr-firststrand or fr-secondstrand for strand-specific data.
--bi STARIndexFolder：The folder name of the STAR binary indexes (i.e., the name of the folder that contains SA file). For example, use ~/STAR_index/hg19 for hg19. (Only if using fastq)
--paired-stats：使用配对统计模型
--variable-read-length：Allow reads with lengths that differ from –readLength to be processed. --readLength will still be used to determine IncFormLen and SkipFormLen
--tmp TMP：The directory for intermediate output such as ".rmats" files from the prep step

(4) 结果

• AS_Event.MATS.JC.txt 仅用跨剪接点的reads来计算剪接
• AS_Event.MATS.JCEC.txt 使用跨剪接点的reads和在条纹区域的reads来计算剪接
• fromGTF.AS_Event.txt 所有来源于GTF和RNA的可能的选择性剪接(AS)事件
• JC.raw.input.AS_Event.txt 事件计数包括仅跨越由rMATS定义的剪接点的reads
• JCEC.raw.input.AS_Event.txt 事件计数包括跨越由rMATS定义的剪接点的reads和不跨越外显子边界的reads

其中JC和JCEC的区别在于前者考虑跨越剪切位点的reads，而后者不仅考虑前者的reads还考虑到只比对到第一张图中条纹的区域（也就是说没有跨越剪切位点的reads），但是我们一般使用JC的结果就够了（如果只是单纯的比较两组样品间可变剪切的差异的话）。

以SE.MATS.JC.txt为例：

1）ID
2）GeneID
3）geneSymbol
4）chr
5）strand

外显子的位置信息：

6）exonStart_0base
7）exonEnd
8）upstreamES
9）upstreamEE
10）downstreamES
11）downstreamEE
12）ID

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

下面几列为：

17）lncFormLen：length of inclusion form, used for normalization（inclusion形式的有效长度）
18）SkipFormLen：length of skipping form, used for normalization（skipping形式的有效长度）
19）P-Value：Significance of splicing difference between two sample groups（两组样品选择性剪接差异的显著性）
20）FDR：False Discovery Rate calculated from p-value
21）lncLevel1：inclusion level for SAMPLE_1 replicates (comma separated) calculated from normalized counts（从归一化计数计算的SAMPLE_1重复的inclusion水平，逗号分隔）
22）IncLevel2：inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized counts（从归一化计数计算的SAMPLE_2重复的inclusion水平，逗号分隔）
23）IncLevelDifference：average(IncLevel1) - average(IncLevel2)（两组样本IncLevel均值的差异）

其中：

lncLevel = Exon Inclusion Isoform / (Exon Inclusion Isoform + Exon Skipping Isoform)
lncLevel = (IJC_SAMPLE / lncFormLen) / (SJC_SAMPLE / SkipFormLen)
effective length = segment length / read length + 1

在输出的结果文件中还有一个summary文本文档，在该文件中统计了上述各个文件中的可变剪接事件数，其中SignificantEventsJC默认为JC文件中满足FDR值小于等于0.05的事件数。

(5) 可视化

rmats2sashimiplot是专门用来对rMATS分析结果做可视化的软件。可以预先用samtools对bam文件建立索引，不然rmats2sashimiplot会先建索引，比较慢。rMATS的结果文件中，如SE.MATS.JC.txt中染色体都是默认加上chr的，而bam文件中的染色体根据基因组注释文件来源不同不一定有chr，需要使两者一致。

方法：

$ rmats2sashimiplot  --b1  ./04_bam/LF1.bam,./04_bam/LF2.bam,./04_bam/LF3.bam  --b2 ./04_bam/CF1.bam,./04_bam/CF2.bam,./04_bam/CF3.bam  -t SE  -e ./06_rmats/plot/SE.MATS.JC_top5.txt  --l1 LF  --l2 CF  --exon_s 1  --intron_s 1  -o ./06_rmats/plot/  &

参数：
--b1 B1：sample_1 in bam format (s1_rep1.bam[,s1_rep2.bam])
--b2 B2：sample_2 in bam format (s2_rep1.bam[,s2_rep2.bam])
-t：rMATS结果中产生的可变剪切类型 {SE,A5SS,A3SS,MXE,RI}
-e EVENTS_FILE：The rMATS output event file (Only if using rMATS format result as event file) 在*.MATS.JC.txt文件里把需要画图的基因信息提取出来
-c：需要画图展示的基因的位置信息，基因组区域的坐标和注释，gff3文件，-c {chromosome}:{strand}:{start}:{end}:{/path/to/gff3}
--l1 L1：The label for first sample
--l2 L2：The label for second sample
-o OUT_DIR：The output directory
--exon_s 1：缩小外显子的大小。默认值为1
--intron_s 1：缩小内含子的大小。例如，如果--intron_s为5，表示内含子的大小为5:1（如果内含子的实际大小是5，那么就将缩小为1）。默认值为1。
--group-info：如果用户想要将样本分组，可以使用* .gf文件指定此参数
--color：可以使用自定义绘图的颜色序列。颜色的数量应该与bam_files对应
--font-size：更改字体大小，默认等于8
--no-text-background：透明文本背景
--hide-number：隐藏连接数
--min-counts 0：用于指定展示的最小count数，如果实际的counts数小于该阈值，则不会在图中显示。

结果：

图上方标题为：可变剪接事件3个外显子所在的基因组坐标及正负链信息。
跨外显子比对的 reads 使用连接外显子 junction 边界的弧线表示。弧线的粗细和比对到 junction 上的 reads 数成正比，同时弧线上的数字指出了 junction reads 的数目。
右上方标注了各个样本可变剪接事件所在的基因及IncLevel（最终会用这2个值做组间的差异AS分析）。
最下方为根据gtf推断出的选择性剪接异构体。