1. sci666首页
  2. 实用技巧
  3. 生物信息学

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

转录本–可变剪切

1、原理

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

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

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

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

  • NCBI

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

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

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

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

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

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

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

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

  • 样本名 做成两个分组

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

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

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

 
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/

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

# 进入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 ./

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

建立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图片

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

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

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

ST
 
4. 结果解读

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

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

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

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

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

以SE.MATS.JC.txt为例:

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

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

  • 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的差异,如果一组内多个样本,那么则是各自组的均值之间差值。

这些可能会帮助到你: 问答社区 | 共享百度SVIP | 留言建议

欢迎入群交流:生信分析群: 732179952 · Meta分析群: 797345521

发表评论

登录后才能评论