NGS可变剪切之STAR+rmats软件使用

mats软件只要你运行成功， 结果还是喜人的， 不过目前TCGA数据库的可变剪切都是一个java软件，叫做spliceseq。我们下次再分享spliceseq咯，这次先让学徒带领大家摸索一下mats软件哈！

本次出场的学徒是前面拖后腿学徒的师弟：拖后腿学徒居然也完成作业，理解RNA-seq数据分析结果，目前先不署名吧，等他完成第二个笔记再说！

Part 1

转录本--可变剪切

1、原理

• 组成型拼接
• 外显子跳跃拼接
• 内含子保留拼接
• 相互排斥的外显子拼接
• 替代5’端剪切
• 替代3'端剪切

2、查看外显子区域

如果想研究一个基因所有外显子区域，而不是单独一个转录本的外显子区域，因此需要获取该基因的所有转录本信息，这里备选三个数据库（NCBI、Ensembl和UCSC）供使用，以BRCA1为例,得到BRCA1基因区域所有“feature”的物理位置，包括外显子:

• NCBI

使用Ensembl数据库获取BRCA1基因的所有外显子区域: https://www.cnblogs.com/yahengwang/p/9361101.html

Part 2

STAR 软件使用

参考：

jimmy的教程：https://mp.weixin.qq.com/s/3JV2ykRqJi_sBHXpa2ONXg

官网帮助：http://rnaseq-mats.sourceforge.net/user_guide.htm

# 以防万一，先创建小环境：
create -n rmats python=2
ca rmats
# 安装rMATS：
conda install -y rmats # 4.02版本
conda install -y rmats2sashimiplot #rmats的可视化软件
RNASeq-MATS.py -h # 测试安装成功

两种输入方式，建议第2种！！

（1）输入文件为fastq文件时，需安装STAR比对软件并提供对应的基因组索引文件，命令如下：

python rmats.py --s1 s1.txt --s2 s2.txt --gtf gtfFile --bi STARindexFolder -od outDir -t readType -readLength readLength [options]* 

（2）输入文件为bam文件时，命令格式如下：

python rmats.py \
    --b1 b1.txt \
    --b2 b2.txt \
    --gtf gtfFile \
    --od outDir \
    -t readType \
    --nthread nthread \
    --readLength readLength \
    --tstat tstat [options]* #可选
    ## group1 即sample1；group2 即sample2

在PRJNA449427

• 样本名 做成两个分组

• 下载gtf注释文件：https://www.gencodegenes.org/human/

1. 安装rmats4.0.1

下载地址：https://sourceforge.net/projects/rnaseq-mats/files/MATS/rMATS.4.0.1.tgz/download

# 下载完后解压：
cd /home/kof/rna/8.alternative-splicing/data
gunzip rMATS.4.0.1.tgz
tar -vxf rMATS.4.0.1.tar
cd rMATS.4.0.1/

# 进入python 看系统版本
python
>>> import sys
>>> print(sys.maxunicode)
# 返回 1114111
#如果出现1114111则说明需要使用rMATS-turbo-Linux-UCS4文件夹下rmats.py；如果出现65535则说明使用rMATS-turbo-Linux-UCS2文件夹下rmats.py
>>> exit()
cd rMATS-turbo-Linux-UCS4
cp rmats.py /home/kof/rna/mapping
cd /home/kof/rna/mapping
cp /home/kof/rna/8.alternative-splicing/data/gencode.v32lift37.annotation.gtf ./
cp /home/kof/rna/8.alternative-splicing/b1.txt ./
cp /home/kof/rna/8.alternative-splicing/b2.txt ./

建立bam文件名的txt文件，分布为b1.txt 、b2.txt —readLength 要看质控报告的read长度

2. 下面直接用 conda 安装的rmats 做:

RNASeq-MATS.py --b1 b1.txt --b2 b2.txt --gtf gencode.v32lift37.annotation.gtf    --od /home/kof/rna/8.alternative-splicing/outDir \
    -t paired --nthread 6 --cstat 0.0001 --readLength 150
# ^Z 
# bg 1  #挂后台

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.raw.input.AS_Event.txt的结果就够了（如果只是单纯的比较两组样品间可变剪切的差异的话）

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

cd /home/kof/rna/8.alternative-splicing/outDir
cd /home/kof/rna/mapping
cp /home/kof/rna/8.alternative-splicing/outDir/ ./

cat > plot.command<<EOF
rmats2sashimiplot \
--b1 SRR6974562.sort.bam,SRR6974563.sort.bam,SRR6974566.sort.bam,SRR6974567.sort.bam \
--b2 SRR6974564.sort.bam,SRR6974565.sort.bam,SRR6974568.sort.bam,SRR6974569.sort.bam \
-t SE \
-e /home/kof/rna/8.alternative-splicing/outDir/SE.MATS.JC.txt \
--l1 SRR697456_1 \
--l2 SRR697456_2 \
--exon_s 1 \
--intron_s 5 \
-o SE_plot
EOF

• 可以改的: --b1 #样本1的全部bam --b2 #样本2的全部bam -t #查看的剪切类型: SE \ A5SS \ A3SS \ MXE \ RI -e # 对应剪切类型的txt文件 --l1 #样本1 bam文件共有的文件名前几个字母 --l2 #样本2 bam文件共有的文件名前几个字母 --exon_s 1 # 显示范围 --intron_s 1 -o # 输出位置,建立一个新的文件夹

rmats2sashimiplot的注意事项：

• 可以预先用samtools对bam文件建索引，不然rmats2sashimiplot也会先建索引（比较慢）
• rmats2sashimiplot默认是出PDF图片

• -e参数后的结果文件里有时会有超级多的基因，如果你想都画图的话，就不用对文件做别的操作，慢慢画就好了（运行时间较长）；如果不想都画，就把需要画图的基因信息提出来就好。
• -e和-c参数不能同时使用；
• 其实直接给定坐标信息画图很直接，很好；但是rMATS的结果文件中的染色体都是默认加上chr的，而bam文件中的染色体信息来源于基因组注释文件，有的没chr，而是直接以数字表示，这时和rmats2sashimiplot预设的不一样，就会报错。可以预先修改基因组文件（加chr）将两者保持一致或者修改代码更改这一块内容，小编觉得这样就比较麻烦了，所以 不建议使用-c参数。

ST

4. 结果解读

放到R里面看会更清晰一点：

以SE.MATS.JC.txt为例：

• 1-5列分别: ID，GeneID，geneSymbol，chr，strand
• 6-11列分别为外显子的位置信息，分别为exonStart_0base，exonEnd，upstreamES，upstreamEE，downstreamES，downstreamEE；网上有张图能很好的解释其含义:

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

• 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 根据归一化计数计算得出的SAMPLE_1重复（以逗号分隔）
• 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）所占比例；IncLevelDifference则是指两组样本IncLevel的差异，如果一组内多个样本，那么则是各自组的均值之间差值。