rMATS进行差异可变剪切分析并可视化

可变剪接(Alternative splicing;又称“选择性剪接”)是一种在真核生物中非常普遍的基因表达方式,具体表现为一个基因的外显子以不同的组合方式剪接形成不同的成熟RNA,从而在不同的时空环境和状态下形成不同的蛋白质,执行不同的生物学功能。常见的可变剪接软件包括rMATS,Asprofile以及miso等。本文主要介绍rMATS软件的使用,并对结果利用rmats2sashimiplot可视化。

rMATS是一个从RNA-Seq数据中检测差异选择性剪接事件的计算工具,根据RNA-Seq数据,rMATS可以自动检测和分析与所有主要类型的可变剪接模式相对应的可变剪接事件。rMATS可识别的可变剪切事件有5种(Figure1,http://rnaseq-mats.sourceforge.net/):

1)Skippedexon (SE) 外显子跳跃

2)Alternative5′ splice site (A5SS) 5’端可变剪切

3)Alternative3′ splice site (A3SS) 3’端可变剪切

4)Mutuallyexclusive exons (MXE) 互斥可变外显子

5)Retainedintron (RI) 内含子保留

Figure 1:rMATS可识别的可变剪切类型

1.分析所需软件

1)Linux操作系统

2)rMATS

3)rmats2sashimiplot

安装使用rMATS需预先安装如下软件:

4)InstallPython 2.7.x and corresponding versions of NumPy and SciPy

5)Downloadand install pysam (rMATS was tested with v0.9.1.4)

6)Downloadand install samtools (version 1.2 or later)

7)Downloadand install STAR (version 2.5 or later)

2.软件下载及安装

1)Python 安装

下载:https://www.python.org/下载后得到文件:Python-2.7.12.tgz,安装过程如下:tar-zxvf Python-2.7.12.tgzcd /software/python/Python-2.7.12./configure--prefix=/software/python/python2712installmakemakeinstallvim ~/.bashrcexportPATH="$PATH:/software/python/python2712install/bin";source ~/.bashrc

Python安装完成后,需要安装如下python模块,python模块一般均可从https://pypi.org/网站下载得到。

NumPy,SciPy,pysam

下面以python setup.py install安装的方式介绍Python模块的安装过程(需要预先安装setuptools模块)。NumPy和pysam模块可以选择该安装方式。

首先下载相关模块文件,以NumPy为例:

下载:https://pypi.org/project/scipy/#files.下载后得到文件numpy-1.16.4.zip,安装过程如下:
unzip numpy-1.16.4.zipcd numpy-1.16.4/python2.7 setup.py install

在安装SciPy时,由于其依赖模块较多,建议下载.whl格式文件,并通过pip install的方式安装(需预先安装pip模块)。

下载:scipy-0.19.0-cp27-cp27m-manylinux1_x86_64.whlpip2.7 install scipy-0.19.0-cp27-cp27m-manylinux1_x86_64.whl

2)samtools

下载:https://sourceforge.net/projects/samtools/files/

安装:

tar -xjf samtools-1.3.1.tar.bz2cd samtools-1.3.1/./configure--prefix=/software/samtools/samtools-1.3.1-installmakemakeinstallvim ~/.bashrcexportPATH="$PATH:/software/samtools/samtools-1.3.1-install/bin";source ~/.bashrc

3)STAR

下载:https://github.com/alexdobin/STAR

安装:

unzip STAR-master.zipvim ~/.bashrcexportPATH="$PATH:/software/STAR/STAR-master/bin/Linux_x86_64";source ~/.bashrc

4)rMATS

下载:http://rnaseq-mats.sourceforge.net/

下载后得到rMATS.4.0.2.tgz文件

安装:

cd /software/rMATStar -zxcf /software/rMATScd rMATS.4.0.2

可以发现该文件夹下出现了两个文件(rMATS-turbo-Linux-UCS2 和rMATS-turbo-Linux-UCS4),如何确定使用哪个文件夹下的文件呢?首先进入python控制台,然后输入以下内容:

>>>import sys>>>print sys.maxunicode1114111Thisoutput indicates that your python is built with --enable-unicode=ucs4, and youshould use rMATS-turbo-xxx-UCS4.>>>import sys>>>print sys.maxunicode65535

This output indicates that your python is built with –enable-unicode=ucs2, and youshould use rMATS-turbo-xxx-UCS2.

笔者在操作过程中,操作系统出现的结果是1114111,因此使用rMATS-turbo-Linux-UCS4下的文件。测试一下是否可以运行:

python/software/rMATS/rMATS.4.0.2/rMATS-turbo-Linux-UCS4/rmats.py --help

5)rmats2sashimiplot

下载:https://github.com/Xinglab/rmats2sashimiplot/

安装:

unzip rmats2sashimiplot-master.zipcd rmats2sashimiplot-master/python2.7 setup.py install

在运行过程中发现还需要安装matplotlib模块,该模块的安装可参照上述.whl格式文件安装过程,在此不在进一步描述。

3.运行rMATS

一般情况下,新一代测序数据文件都非常大,本次以官网测试数据为例。

http://rnaseq-mats.sourceforge.net/rmats4.0.2/下载测试数据testData.tgz(Figure 2)

Figure 2:testData.tgz数据文件

Fastq和bam文件均为二代测序过程中常见的文件,此处不做过多描述,主要看一下txt格式文件中的信息。

cat  s1.txt231ESRP.25K.rep-1.R1.fastq:231ESRP.25K.rep-1.R2.fastq,231ESRP.25K.rep-2.R1.fastq:231ESRP.25K.rep-2.R2.fastq
cat  s2.txt231EV.25K.rep-1.R1.fastq:231EV.25K.rep-1.R2.fastq,231EV.25K.rep-2.R1.fastq:231EV.25K.rep-2.R2.fastq

 

1)通过fastq文件开始运行:

python2.7 /software/rMATS/rMATS.4.0.2/rMATS-turbo-Linux-UCS4/rmats.py --s1 s1.txt --s2 s2.txt --gtf Homo_sapiens.GRCh37.75.gtf --bi STARindexFolder -od outDir -tpaired --readLength 50 --libType fr-unstranded --nthread 10

2)通过bam文件开始运行:

python2.7 /software/rMATS/rMATS.4.0.2/rMATS-turbo-Linux-UCS4/rmats.py --b1 b1.txt --b2 b2.txt --gtf Homo_sapiens.GRCh37.75.gtf --od outDir -t paired --readLength 50--libType fr-unstranded --nthread 10

参数详解:

--s1 s1.txt  A text file contains FASTQfile(s) for the sample_1.文件内以逗号分隔重复样本的FASTQ文件名(Onlyif using fastq)--s2 s2.txt  A text file contains FASTQfile(s) for the sample_2.文件内以逗号分隔重复样本的FASTQ文件名(Onlyif using fastq)--b1 b1.txt  A text file records mappingresults for the sample_1 in bam format. 文件内以逗号分隔重复样本的bam文件名(Only if using bam)--b2 b2.txt  A text file records mappingresults for the sample_2 in bam format. 文件内以逗号分隔重复样本的bam文件名(Only if using bam)-treadType  Type of read used in theanalysis. readType is either 'paired' or 'single'. 'paired' is for paired-enddata and 'single' is for single-end data--readLength  测序reads的长度--gtf  gtfFile 需要输入的gtf格式文件--bi  STARIndexFolderThe folder name of the STARbinary indexes (i.e., the name of the folder that contains SA file). Forexample, use ~/STARindex/hg19 for hg19. (Only if using fastq)--odoutDir  所有输出文件的路径--tophatAnchorThe "anchor length" or "overhang length" used in thealigner. At least “anchor length” NT must be mapped to each end of a givenjunction. The default is 6. (This parameter applies only if using fastq)--nthreadThe number of thread.--cstat  The cutoff splicing difference. The cutoffused in the null hypothesis test for differential splicing. The default is0.0001 for 0.01% difference. Valid: 0 ≤ cutoff < 1--tstat  The number of thread for statistical model.--statoff  Turn statistics part off.-libTypelibraryTypeLibrary type.  Default isunstranded (fr-unstranded). Use fr-firststrand or fr-secondstrand forstrand-specific data.

软件运行后主要产生如下结果文件(Figure 3):rMATS的结果文件是以各个可变剪切事件的分布的,主要由以下几类构成,详细可参考http://rnaseq-mats.sourceforge.net/rmats4.0.2/user_guide.htm.

AS_Event.MATS.JC.txtevaluates splicing with only reads that span splicing junctions.

AS_Event.MATS.JCEC.txtevaluates splicing with reads that span splicing junctions and reads on target(striped regions on Figure 1).

fromGTF.AS_Event.txtall possible alternative splicing (AS) events derived from GTF and RNA.

JC.raw.input.AS_Event.txtevaluates splicing with only reads that span splicing junctions.

JCEC.raw.input.AS_Event.txtevaluates splicing with reads that span splicing junctions and reads on target(striped regions on Figure 1).

其中JC和JCEC的区别在于JC考虑跨越剪切位点的reads,而JCEC不仅考虑前者的reads还考虑到比对到Figure 1中条纹的区域(也就是说没有跨越剪切位点的reads),一般情况使用JC的结果就够了(如果只是单纯的比较两组样品间可变剪切的差异的话)

Figure 3:rMATS结果文件

这几类文件中我们需要重点关注的是AS_Event.MATS.JC.txt,因为其他文件的信息大多数都整合在该文件中,以SE.MATS.JC.txt为例,该文件主要包含如下各列信息:

1)ID2)GeneID3)geneSymbol4)chr5)strand6)exonStart_0base7)exonEnd8)upstreamES9)upstreamEE10)downstreamES11)downstreamEE12)ID13)IJC_SAMPLE_1:样品在inclusion junction下的count数,重复样本的结果以逗号分隔14)SJC_SAMPLE_1:样品在skipping junction counts(SJC)下的count数,重复样本的结果以逗号分隔15)IJC_SAMPLE_2:样品在inclusion junction下的count数,重复样本的结果以逗号分隔16)SJC_SAMPLE_2:样品在skipping junction counts(SJC)下的count数,重复样本的结果以逗号分隔17)IncFormLen:lengthof inclusion form, used for normalization18)SkipFormLen:lengthof skipping form, used for normalization19)PValue:Significanceof splicing difference between two sample groups(两组样品可变剪切的统计学显著差异指标)20)FDR:FalseDiscovery Rate calculated from p-value(校正后的PValue)21)IncLevel1:inclusionlevel for SAMPLE_1 replicates (comma separated) calculated from normalizedcounts22)IncLevel2:inclusionlevel for SAMPLE_2 replicates (comma separated) calculated from normalizedcounts23)IncLevelDifference:average(IncLevel1)– average(IncLevel2)

4.利用rmats2sashimiplot可视化

rmats2sashimiplot --b1 231ESRP.25K.rep-1.bam,231ESRP.25K.rep-2.bam --b2 231EV.25K.rep-1.bam,231EV.25K.rep-2.bam -t SE -e /software/rMATS/test/outDir/SE.MATS.JC.txt --l1 231ESRP.25K --l2 231EV.25K -o SE_plot

注意:各个文件中染色体个格式均需保持一致,要么使用1,2,3……21,22,X,Y,要么使用chr1,chr2,chr3……chr21,chr22,chrX,chrY

参数详解:

--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.

程序产生的图片如下所示(Figure 4,Figure5,Figure6):

 

 

Figure 4:HRAS基因可变剪切

 

 

Figure 5:RPLP2基因可变剪切

 

 

Figure 6:RNH1基因可变剪切

好了,可变剪切分析和可视化就这样完成了,貌似有些复杂,不过只要你不抛弃不放弃,还是可以实现的。当然谁都不可能一次就成功,生信分析的坑还是不少的,这次不行,就再试一次,一步步的跨过这些坑才会成为生信大神!

生物信息学

tcR包:T细胞受体和免疫球蛋白数据进行高级分析和可视化(二)

2020-8-28 0:49:53

生物信息学

R包survminer画生存曲线的实用技能

2020-8-28 0:53:30

加入Q群
0 条回复 A文章作者 M管理员
    暂无讨论,说说你的看法吧
个人中心
今日签到
有新私信 私信列表
搜索