转录组----可变剪接分析-rMATs

作者: 千万别加香菜 | 来源:发表于2022-11-08 10:29 被阅读0次

转录组----可变剪接分析-rMATs
生信菜鸟团-专题学习目录
2021-06-23 剪接因子数据集
笔记：rmats绘图结果怎么看
转录组分析(8) - 可变剪接
可变剪切分析（三）rMATS差异可变剪切分析
肿瘤可变剪切相关概念
IsoformSwitchAnalyzeR 可变剪切分析
可变剪接
详细介绍最新版可变剪接软件rMATS

●可变剪切（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) 内含子保留 5种rMATs识别的可变剪接事件.png

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、结果展示

rMATs.png

会输出好几种文件，其中.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）差不多