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

mats软件只要你运行成功, 结果还是喜人的, 不过目前TCGA数据库的可变剪切都是一个java软件,叫做spliceseq。我们下次再分享spliceseq咯,这次先让学徒带领大家摸索一下mats软件哈!

本次出场的学徒是前面拖后腿学徒的师弟拖后腿学徒居然也完成作业,理解RNA-seq数据分析结果,目前先不署名吧,等他完成第二个笔记再说!

Part 1

转录本--可变剪切

1、原理
  • 组成型拼接

  • 外显子跳跃拼接

  • 内含子保留拼接

  • 相互排斥的外显子拼接

  • 替代5’端剪切

  • 替代3'端剪切

2、查看外显子区域

如果想研究一个基因所有外显子区域,而不是单独一个转录本的外显子区域,因此需要获取该基因的所有转录本信息,这里备选三个数据库(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 # 测试安装成功

两种输入方式,建议第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

  • 样本名 做成两个分组

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

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的注意事项:

  • 可以预先用samtools对bam文件建索引,不然rmats2sashimiplot也会先建索引(比较慢)

  • rmats2sashimiplot默认是出PDF图片

  • -e参数后的结果文件里有时会有超级多的基因,如果你想都画图的话,就不用对文件做别的操作,慢慢画就好了(运行时间较长);如果不想都画,就把需要画图的基因信息提出来就好。

  • -e和-c参数不能同时使用;

  • 其实直接给定坐标信息画图很直接,很好;但是rMATS的结果文件中的染色体都是默认加上chr的,而bam文件中的染色体信息来源于基因组注释文件,有的没chr,而是直接以数字表示,这时和rmats2sashimiplot预设的不一样,就会报错。可以预先修改基因组文件(加chr)将两者保持一致或者修改代码更改这一块内容,小编觉得这样就比较麻烦了,所以 不建议使用-c参数

ST
4. 结果解读

放到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

(0)

相关推荐