一、 基本介紹

rMATS是一款對(duì)RNA-Seq數(shù)據(jù)進(jìn)行差異可變剪切分析的軟件宰啦。其通過rMATS統(tǒng)計(jì)模型對(duì)不同樣本（有生物學(xué)重復(fù)的）進(jìn)行可變剪切事件的表達(dá)定量，然后以likelihood-ratio test計(jì)算P value來表示兩組樣品在IncLevel（Inclusion Level）水平上的差異，并利用Benjamini Hochberg算法對(duì)p value進(jìn)行校正得FDR值掉房。（lncLevel是Exon Inclusion Isoform在總Isoform（Exon Inclusion Isoform + Exon Skipping Isoform）中所占比例。）

rMATS可識(shí)別的可變剪切事件有5種慰丛，分別是外顯子跳躍（skipped exon, SE）卓囚，選擇性5’剪接位點(diǎn)（alternative 5' splice site, A5SS），選擇性3’剪接位點(diǎn)（alternative 3' splice site, A3SS）诅病，外顯子互斥（mutually exclusive exons, MXE）和 內(nèi)含子保留（retained intron, RI）哪亿。

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

二粥烁、 使用方法

(1) 輸入

rMATS支持兩種格式文件的輸入。

• 第一種是fastq格式蝇棉，那么在安裝的時(shí)候還需要安裝STAR比對(duì)軟件以及提供比對(duì)的索引文件（STAR的索引文件異常的大）讨阻，所以rMATS其實(shí)是建議使用第二種方式；
• 第二種是bam格式篡殷，rMATS支持其他比對(duì)軟件比對(duì)后的結(jié)果bam文件作為輸入钝吮，比如tophat/ hisat2等，這樣也能減少rMATS的運(yùn)行時(shí)間板辽。

(2) 方法

目前rMATS 4.1版不受限于單雙端測(cè)序奇瘦，reads長(zhǎng)度不一，是否存在生物學(xué)重復(fù)劲弦，是否有比較組耳标，是否需要檢測(cè)新轉(zhuǎn)錄本，是否鏈特異性等條件邑跪，并且其可以進(jìn)行分步次坡，分機(jī)器計(jì)算，功能完善画畅，主要可變剪切事件檢測(cè)完整的一款軟件砸琅。

$ 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

構(gòu)建索引：

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

運(yùn)行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) 參數(shù)

--b1 b1.txt：輸入sample1的txt格式的文件，文件內(nèi)以逗號(hào)分隔重復(fù)樣本的bam文件名（/path/to/1_1.bam,/path/to/1_2.bam）
--b2 b2.txt：輸入sample2的txt格式的文件轴踱，文件內(nèi)以逗號(hào)分隔重復(fù)樣本的bam文件名（/path/to/2_1.bam,/path/to/2_2.bam）
-t {paired,single}：雙端測(cè)序則readType為paired明棍，單端測(cè)序則為single
--readLength：測(cè)序reads的長(zhǎng)度
--gtf gtfFile：需要輸入的gtf文件
--od outDir：所有輸出文件的路徑（文件夾）
--nthread：設(shè)置線程數(shù)
--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. 代表的是兩個(gè)樣本中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：使用配對(duì)統(tǒng)計(jì)模型
--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) 結(jié)果

• AS_Event.MATS.JC.txt 僅用跨剪接點(diǎn)的reads來計(jì)算剪接
• AS_Event.MATS.JCEC.txt 使用跨剪接點(diǎn)的reads和在條紋區(qū)域的reads來計(jì)算剪接
• fromGTF.AS_Event.txt 所有來源于GTF和RNA的可能的選擇性剪接(AS)事件
• JC.raw.input.AS_Event.txt 事件計(jì)數(shù)包括僅跨越由rMATS定義的剪接點(diǎn)的reads
• JCEC.raw.input.AS_Event.txt 事件計(jì)數(shù)包括跨越由rMATS定義的剪接點(diǎn)的reads和不跨越外顯子邊界的reads

其中JC和JCEC的區(qū)別在于前者考慮跨越剪切位點(diǎn)的reads，而后者不僅考慮前者的reads還考慮到只比對(duì)到第一張圖中條紋的區(qū)域（也就是說沒有跨越剪切位點(diǎn)的reads）寇僧，但是我們一般使用JC的結(jié)果就夠了（如果只是單純的比較兩組樣品間可變剪切的差異的話）摊腋。

以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數(shù)，重復(fù)樣本的結(jié)果以逗號(hào)分隔嘁傀；列名分別為IJC_SAMPLE_1兴蒸，SJC_SAMPLE_1，IJC_SAMPLE_2细办，SJC_SAMPLE_2橙凳。

下面幾列為：

17）lncFormLen：length of inclusion form, used for normalization（inclusion形式的有效長(zhǎng)度）
18）SkipFormLen：length of skipping form, used for normalization（skipping形式的有效長(zhǎng)度）
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（從歸一化計(jì)數(shù)計(jì)算的SAMPLE_1重復(fù)的inclusion水平，逗號(hào)分隔）
22）IncLevel2：inclusion level for SAMPLE_2 replicates (comma separated) calculated from normalized counts（從歸一化計(jì)數(shù)計(jì)算的SAMPLE_2重復(fù)的inclusion水平笑撞，逗號(hào)分隔）
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

在輸出的結(jié)果文件中還有一個(gè)summary文本文檔岛啸，在該文件中統(tǒng)計(jì)了上述各個(gè)文件中的可變剪接事件數(shù)，其中SignificantEventsJC默認(rèn)為JC文件中滿足FDR值小于等于0.05的事件數(shù)茴肥。

(5) 可視化

rmats2sashimiplot是專門用來對(duì)rMATS分析結(jié)果做可視化的軟件坚踩。可以預(yù)先用samtools對(duì)bam文件建立索引瓤狐，不然rmats2sashimiplot會(huì)先建索引瞬铸，比較慢批幌。rMATS的結(jié)果文件中，如SE.MATS.JC.txt中染色體都是默認(rèn)加上chr的嗓节，而bam文件中的染色體根據(jù)基因組注釋文件來源不同不一定有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/  &

參數(shù)：
--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結(jié)果中產(chǎn)生的可變剪切類型 {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：需要畫圖展示的基因的位置信息拦宣，基因組區(qū)域的坐標(biāo)和注釋截粗，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：縮小外顯子的大小鸵隧。默認(rèn)值為1
--intron_s 1：縮小內(nèi)含子的大小绸罗。例如，如果--intron_s為5掰派，表示內(nèi)含子的大小為5:1（如果內(nèi)含子的實(shí)際大小是5从诲，那么就將縮小為1）左痢。默認(rèn)值為1靡羡。
--group-info：如果用戶想要將樣本分組，可以使用* .gf文件指定此參數(shù)
--color：可以使用自定義繪圖的顏色序列俊性。顏色的數(shù)量應(yīng)該與bam_files對(duì)應(yīng)
--font-size：更改字體大小略步，默認(rèn)等于8
--no-text-background：透明文本背景
--hide-number：隱藏連接數(shù)
--min-counts 0：用于指定展示的最小count數(shù)，如果實(shí)際的counts數(shù)小于該閾值定页，則不會(huì)在圖中顯示趟薄。

結(jié)果：

圖上方標(biāo)題為：可變剪接事件3個(gè)外顯子所在的基因組坐標(biāo)及正負(fù)鏈信息。
跨外顯子比對(duì)的 reads 使用連接外顯子 junction 邊界的弧線表示典徊『技澹弧線的粗細(xì)和比對(duì)到 junction 上的 reads 數(shù)成正比，同時(shí)弧線上的數(shù)字指出了 junction reads 的數(shù)目卒落。
右上方標(biāo)注了各個(gè)樣本可變剪接事件所在的基因及IncLevel（最終會(huì)用這2個(gè)值做組間的差異AS分析）羡铲。
最下方為根據(jù)gtf推斷出的選擇性剪接異構(gòu)體。