●可变剪切(differential splicing)也叫做选择性剪切alternative splicing, 指的是在mRNA前体到成熟mRNA的过程当中,不同的剪切方式使得同一个基因可以产生多个不同的成熟mRNA, 最终产生不同的蛋白质
https://zhuanlan.zhihu.com/p/409865441
rMATs软件可以识别5种可变剪接事件:Skippedexon (SE) 外显子跳跃、Alternative5’ splice site (A5SS) 5’端可变剪切、Alternative3’ splice site (A3SS) 3’端可变剪切、Mutuallyexclusive exons (MXE) 互斥可变外显子、Retainedintron (RI) 内含子保留
1、安装
可以用conda直接安装
conda install -c bioconda rmats
2、使用
1)、计算可变剪接事件
rmats.py --b1 b1.txt \
--b2 b2.txt \
--gtf /home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf \
--od AS \
--tmp tmp \
-t paired \
--readLength 150 \
--cstat 0.0001 \
--nthread 10
--b1 b1.txt 输入sample1的txt格式的文件,文件内以逗号分隔重复样本的bam文件名
--b2 b2.txt 输入sample2的txt格式的文件,文件内以逗号分隔重复样本的bam文件名
-t readType 双端测序则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
--statoff,进行单样本或者是单组的分析,并跳过统计分析
2)可视化
1 整体可视化
rmats2sashimiplot --b1 SRR17709921_sort.bam,SRR17709920_sort.bam,SRR17709917_sort.bam \
--b2 SRR17709910_sort.bam,SRR17709918_sort.bam,SRR17709919_sort.bam \
-t SE \
-e SE.MATS.JC.txt \
--l1 DP_L \
--l2 Han_L \
-o SE_plot &
可以将需要可视化的基因进行筛选,重新做成SE.MATS.JC.txt这种文件,然后可视化就可以了
rmats2sashimiplot --b1 SRR17709911_sort.bam,SRR17709912_sort.bam,SRR17709913_sort.bam \
--b2 SRR17709916_sort.bam,SRR17709915_sort.bam,SRR17709914_sort.bam \
-t SE \
-e SE.MATS.JC.txt \
--l1 DP_M \
--l2 Han_M \
-o M_SE_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 (Onlyif using rMATSformat result as event file).
--l1 L1 The label for first sample.
--l2 L2 The label for second sample.-o OUT_DIR The output directory.
2 基因坐标输入时
报错了
Gene: FGF1 in muscle
rmats2sashimiplot --b1 SRR17709911_sort.bam,SRR17709912_sort.bam,SRR17709913_sort.bam \
--b2 SRR17709910_sort.bam,SRR17709918_sort.bam,SRR17709919_sort.bam \
-c chrNC_040256.1:-:55979601:56069122:/home/sll/genome-sheep/Oar_rambouillet_v1.0-ncbi/GCF_002742125.1_Oar_rambouillet_v1.0_genomic.gtf \
--l1 DP_M \
--l2 Han_M -o ./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 (Onlyif using rMATSformat result as event file).
--l1 L1 The label for first sample.
--l2 L2 The label for second sample.-o OUT_DIR The output directory.
3、结果展示
会输出好几种文件,其中.MATS.JC.txt是我们要用到的文件
以MXE.MATS.JC.txt为例说明每列的意义,引自CSDN博主「次亚硫酸钠」的原创文章https://blog.csdn.net/weixin_42910678/article/details/123587203:
ID GeneID geneSymbol chr strand 1stExonStart_0base 1stExonEnd 2ndExonStart_0base 2ndExonEnd upstreamES upstreamEE downstreamES downstreamEE ID IJC_SAMPLE_1 SJC_SAMPLE_1 IJC_SAMPLE_2 SJC_SAMPLE_2 IncFormLen SkipFormLen PValue FDR IncLevel1 IncLevel2 IncLevelDifference
0 "MS.gene23798" NA chr8.4 - 30758609 30758704 30759122 30759182 30758025 30758095 30760455 30760527 0 1 11 7 9 209 244 0.0120878457309 0.0604392286545 0.096 0.476 -0.38
1 "MS.gene61989" NA chr7.2 - 80697619 80697769 80704270 80704420 80697113 80697232 80705567 80706851 1 1 1 0 3 298 298 0.102057409464 0.34019136488 0.5 0.0 0.5
●ID: 官网描述“rMATS event id”,其实就是序号
●GenelD: 可变剪接事件所在基因编号
●geneSymbol: 可变剪接事件所在基因名称
●chr: 可变剪接事件所在染色体
●strand: 可变剪接事件所在染色体链的方向
●1stExonStart_0base: 第一个可变剪接事件跳跃外显子的起始位置,以0开始计数
●1stExonEnd: 第一个可变剪接事件跳跃外显子的终止位置
●2ndExonStart_0base:第二个可变剪接事件跳跃外显子的起始位置,以0开始计数
●2ndExonEnd: 第二个可变剪接事件跳跃外显子的终止位置
●upstreamES: 可变剪接事件跳跃外显子的上游exon起始位置
●upstreamEE: 可变剪接事件跳跃外显子的上游exon终止位置
●downstreamES: 可变剪接事件跳跃外显子的下游exon起始位置
●downstreamEE: 可变剪接事件跳跃外显子的下游exon终止位置
●ID: 同上
●IJC_SAMPLE_1: 样本一在inclusion junction(IJC)下的count数,重复样本的结果以逗号分隔
●SJC_SAMPLE_1: 样本一在skipping junction(SJC)下的count数,重复样本的结果以逗号分隔
●IJC_SAMPLE_2: 样本二在inclusion junction(IJC)下的count数,重复样本的结果以逗号分隔
●SJC_SAMPLE_2: 样本二在skipping junction(SJC)下的count数,重复样本的结果以逗号分隔
●IncFormLen: 可变剪接事件Exon Inclusion Isoform的有效长度
●SkipFormLen: 可变剪接事件Exon Skipping Isoform的有效长度
●PValue: 两组样本间可变剪接事件表达差异显著性p值
●FDR: 可变剪接事件表达差异显著性FDR值
●IncLevel1: 处理组可变剪接事件Exon Inclusion Isoform在两个Isoform总表达量的比值,也就是PSI
●IncLevel2: 对照组可变剪接事件Exon Inclusion Isoform在两个Isoform总表达量的比值,也就是PSI
●IncLevelDifference: IncLevel1与IncLevel2的差值,和dPSI(different percent spliced in)差不多