mats软件只要你运行成功, 结果还是喜人的, 不过目前TCGA数据库的可变剪切都是一个java软件,叫做spliceseq。我们下次再分享spliceseq咯,这次先让学徒带领大家摸索一下mats软件哈!
本次出场的学徒是前面拖后腿学徒的师弟:拖后腿学徒居然也完成作业,理解RNA-seq数据分析结果,目前先不署名吧,等他完成第二个笔记再说!
Part 1
转录本--可变剪切
1、原理
2、查看外显子区域
如果想研究一个基因所有外显子区域,而不是单独一个转录本的外显子区域,因此需要获取该基因的所有转录本信息,这里备选三个数据库(NCBI、Ensembl和UCSC)供使用,以BRCA1为例,得到BRCA1基因区域所有“feature”的物理位置,包括外显子:
使用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
SRR6974562 DCIS KO | SRR6974566 CA1a KO |
---|---|
SRR6974563 DCIS KO | SRR6974567 CA1a KO |
SRR6974564 DCIS NC | SRR6974568 CA1a NC |
SRR6974565 DCIS NC | SRR6974569 CA1a NC |
在PRJNA449427
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的注意事项:
ST
4. 结果解读
放到R里面看会更清晰一点:
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的差异,如果一组内多个样本,那么则是各自组的均值之间差值。