NGS可变剪切之STAR+rmats软件使用
mats软件只要你运行成功, 结果还是喜人的, 不过目前TCGA数据库的可变剪切都是一个java软件,叫做spliceseq。我们下次再分享spliceseq咯,这次先让学徒带领大家摸索一下mats软件哈!
本次出场的学徒是前面拖后腿学徒的师弟:拖后腿学徒居然也完成作业,理解RNA-seq数据分析结果,目前先不署名吧,等他完成第二个笔记再说!
Part 1
转录本--可变剪切
组成型拼接
外显子跳跃拼接
内含子保留拼接
相互排斥的外显子拼接
替代5’端剪切
替代3'端剪切
如果想研究一个基因所有外显子区域,而不是单独一个转录本的外显子区域,因此需要获取该基因的所有转录本信息,这里备选三个数据库(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 # 测试安装成功
(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
样本名 做成两个分组
下载gtf注释文件:https://www.gencodegenes.org/human/
下载地址: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长度
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的结果就够了(如果只是单纯的比较两组样品间可变剪切的差异的话)
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参数。
放到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的差异,如果一组内多个样本,那么则是各自组的均值之间差值。
友情推广
如果你也对学徒培养或者实习职位感兴趣,想在我们的指导下完成肿瘤外显子等NGS数据分析,可以先看看我是如何培养学徒的:
当然了,学徒培养看缘分!发邮件给我申请:jmzeng1314@163.com