CN107766696A

CN107766696A - 基于RNA‑seq数据的真核生物可变剪接分析方法和系统

Info

Publication number: CN107766696A
Application number: CN201610707885.8A
Authority: CN
Inventors: 张翼; 程超
Original assignee: ABLIFE (WUHAN) Inc
Current assignee: ABLIFE (WUHAN) Inc
Priority date: 2016-08-23
Filing date: 2016-08-23
Publication date: 2018-03-06

Abstract

本发明提供一种基于RNA‑seq数据的真核生物可变剪接分析方法和系统。包括通过illumina二代测序平台获取某一具有参考基因组和注释的真核生物的一个或多个样品的转录组原始测序数据；将质量不合格的数据过滤掉，留下的数据作为待分析的数据；接着进行基础分析：将各个转录组样本待分析数据分别比对到所述物种的参考基因组，筛选出唯一比对的结果；计算各样本基因的表达量；筛选出显著差异表达的基因；对差异基因进行功能注释和分析；然后进行可变剪接分析：已知可变剪接事件的鉴定；新的可变剪接事件的鉴定；样品（组）间可变剪接事件差异分析；可变剪接与基因表达关联分析；可变剪接分析结果统计和报表生成；可变剪接可视化图生成。

Description

基于RNA-seq数据的真核生物可变剪接分析方法和系统

技术领域

本发明涉及生物信息技术领域，尤其涉及一种基于RNA-seq数据的真核生物可变剪接分析方法和系统。

背景技术

在真核生物中，有些基因的一个mRNA前体通过不同的剪接方式(选择不同的剪接位点)产生不同的mRNA剪接异构体，这一过程称为可变剪接(或选择性剪接，alternativesplicing)。可变剪接是调节基因表达和产生蛋白质组多样性的重要机制，是导致真核生物基因和蛋白质数量较大差异的重要原因。

高通量测序技术(High-throughput sequencing)又称“下一代”测序技术("Next-generation"sequencing technology)，以能一次并行对几十万到几百万条DNA分子进行序列测定和一般读长较短等为标志。而RNA-seq测序技术利用高通量测序平台已经成为非常广泛被使用的研究RNA的技术。而使用RNA-seq技术来研究可变剪接比其他技术拥有更多的好处。

到目前为止,RNA-seq分析方法主要集中在基因表达水平的评估和发现新的外显子,以及转录本的组装和注释,对外显子选择从表达水平进行量化和研究。可变的剪接事件的研究和分析依然具有挑战性。

发明内容

针对现有技术中存在的问题，本发明的主要目的是提供一种基于RNA-seq数据的可变剪接分析方法，对可变剪接提供更全面的分析，并分析了可变剪接对生物体功能产生的影响。另外，本发明基于以上的方法建立了一个自动化的分析系统。

一方面，本发明提供了一种基于RNA-seq数据的可变剪接分析方法，包括：

1)通过illumina二代测序平台获取某一具有参考基因组和注释的真核生物的一个或多个样品的转录组原始测序数据；

2)对上述各组原始测序数据进行过滤，将质量不合格的数据过滤掉，留下的数据作为待分析的数据，过滤的条件是：截掉adapter接头及之后的序列；截掉序列末尾质量低于20的碱基；丢掉序列长度小于16的序列；去掉50％碱基质量低于20的序列；

3)对各个转录组的待分析数据进行基础分析和可变剪接分析，其中，所述的基础分析包括：

(1)将所述各个转录组样本待分析数据分别比对到所述物种的参考基因组，获取发生剪接的比对结果，并筛选出唯一比对的结果；

(2)计算各样本基因的表达量：基于RPKM标准化方法使用python编写程序，计算基因表达量信息；

(3)将各样品按照样品间或样品组间进行差异分析，筛选出显著差异表达的基因：样本(组)间差异分析使用R软件包edgeR进行，显著差异基因的筛选标准为：p value小于等于0.01，fold change大于等于2；

(4)对差异基因进行功能注释和分析：包括样品间相关性分析，差异基因聚类分析，差异基因GO富集分析，差异基因KEGG Pathway分析。

所述的可变剪接分析包括：

(1)参考基因组注释文件中已知可变剪接事件的鉴定；

(2)新的可变剪接事件的鉴定；

(3)样品(组)间可变剪接事件差异分析；

(4)可变剪接与基因表达关联分析；

(5)可变剪接分析结果统计和报表生成；

(6)可变剪接可视化图生成：使用perl编写程序，绘制可变剪接事件的可视化图。

在本发明的一个实施例中，基础分析的比对使用tophat2软件进行，软件的参数具体设置如下：设置比对reads的错配数为4；设置Bowtie2片段比对最大错配数为1；设置reads最大的多位置比对结果输出个数为2；设置线程数为16；其他均使用软件默认设置。

在本发明的又一个实施例中，各个转录组样本待分析数据分别比对到所述物种的参考基因组得到结果后，筛选出唯一比对结果的方法如下：检查bam文件中每条比对结果TAG的NH，如果匹配“NH:i:1”，则表示该reads是唯一比对的结果，保留下来，否则就丢掉。最后筛选留下的结果使用samtools工具转换为bam文件，并建立index；该bam用于后续分析；Tophat2可以提取出发生剪接的reads比对结果，并生成bed文件：junctions.bed，该文件是后序可变剪接分析的输入文件。

在本发明的又一个实施例中，可变剪接分析中研究的可变剪接事件包括的类别如下：外显子跳跃事件(ES/cassetteExon)，互斥外显子事件，可变3’剪接事件，可变5’剪接事件，可变的第一个外显子事件，可变的最后一个外显子事件，同时外显子跳跃和可变3’剪接的事件，同时外显子跳跃和可变5’剪接的事件，内含子保留事件。

在本发明的又一个实施例中，参考基因组注释文件中已知可变剪接事件的鉴定步骤为：首先为每个基因定义一个基因模型，也就是gene model，默认选择注释文件中该基因第一个转录本作为gene model；将gene model的剪接位置(splice junction)提取出来，及exon与exon的连接位点信息，再将所有基因的剪接位置提取出来，再使用自己编写的perl程序分析这些剪接位置，鉴定出相对于gene model存在的可变剪接事件并进行筛选和分类，按照规范格式进行输出；已知可变剪接事件的筛选条件如下：对于已知的非内含子保留事件，相对于model可变的剪接型支持reads要大于等于2，model的支持reads要大于等于2，两者加起来要大于等于10。

在本发明的又一个实施例中，新的可变剪接事件的鉴定的步骤如下：以genemodel作为参考，从tophat2检测出的发生可变剪接的比对结果中筛选出新的splicejunction结果，使用自己编写的perl程序鉴定相对于gene model发生的可变剪接事件，并进行筛选和分类，按照规范格式进行输出；新的可变剪接事件的筛选条件如下：相对于model可变的剪接型支持reads要大于等于2，该剪接型比例需大于等于0.15。

在本发明的又一个实施例中，样品(组)间可变剪接事件差异分析的步骤如下：分别计算每个样品每个可变剪接事件发生的比例，计算样品间该比例的差值，使用费雪尔正确性检定(Fisher’s Exact Test)或t检验(T-Test)计算差异的可靠性，即得到p value，差异可变剪接事件的筛选条件如下：样品间比例差值需大于等于0.2，p value需小于等于0.05。

另一方面，本发明还提供了一种基于RNA-seq数据的可变剪接分析的系统，包括：

数据处理模块，用于过滤各样本原始测序数据中质量不合格的序列，获得各个转录组的clean reads；

比对模块，用于将各个转录组数据比对到参考基因组上，并提取唯一比对序列；

表达量分析模块，用于对各样本的基因进行定量分析，计算各基因的RPKM；

差异基因分析模块，用于样本(组)间差异分析，筛选出显著差异表达的基因，并按上下调进行分类；

功能注释分析模块，用于对选定基因进行GO富集分析和KEGG Pathway分析；

已知可变剪接事件检测和注释模块，用于gene model的提取和定义，以及已知可变剪接事件的检测和注释；

新的可变剪接事件检测和注释模块，用于鉴定新的可变剪接事件，并进行定量和筛选；

可变剪接事件差异分析模块，用于样本(组)间可变剪接事件差异分析，筛选出发生显著差异的可变剪接事件，并与发生表达显著差异的基因做关联分析；

可变剪接事件统计模块，用于统计和生成报表；

可变剪接事件全基因组可视化模块，用于绘制可变剪接事件的可视化图。

本发明提供了一整套可变剪接分析方案，从数据处理到RNA-seq数据基础分析，以及可变剪接的各种分析，相比其他平台无明显缺点或劣势；本发明能够快速准确的鉴定已知和新的可变剪接事件，可变剪接事件类别细分为10种，在鉴定和筛选新的可变剪接事件时使用的参数是根据多年可变剪接研究经验得到的优化值；对可变剪接事件进行了定量，计算获取可变型发生的比率。针对两个以上样品，可以进行可变剪接差异分析，以获取有可变剪接进行调控的基因，进而研究可变剪接调控导致的功能变化；可视化展示可变剪接事件。可广泛应用于基础研究、临床诊断和药物研究等领域。

附图说明

图1是本发明利用illumina Nextseq500测序平台进行转录组测序的一个实施例的流程示意图。

图2是本发明基于illumina二代测序数据的真核生物可变剪接分析方法的一个实施例的流程示意图。

图3是不同可变剪接事件鉴定算法的示意图。

图4是可变剪接事件差异分析的流程示意图。

图5是本发明一个实施例中产生的外显子跳跃事件可视化图。

图6是本发明一个实施例中产生的内含子保留事件可视化图。

图7是本发明一个实施例中产生的可变3’剪接位点事件和互斥外显子事件可视化图。

图8是可变剪接分析系统的一个实施例的框图。

图9是外显子跳跃事件的示意图。

具体实施方式

通过以下详细说明结合附图可以进一步理解本发明的特点和优点。所提供的实施例仅是对本发明方法的说明，而不以任何方式限制本发明揭示的其余内容。

除非另有说明，否则在这些实施例中阐述的部件和步骤的相对布置、数字表达式和数值不构成对本发明的限制。对于本领域普通技术人员已知的技术、方法和设备可能不作详细讨论，但在适当情况下，技术、方法和设备应当被视为本说明的一部分。

【实施例1】基础分析

本发明利用illumina Nextseq500测序平台进行转录组测序的流程示意如图1。得到样品基于NextSeq500平台转录组测序数据后，寻找样品的参考数据库及相应的注释文件(物种本身的基因、基因组)，然后利用如图2所示的比较分析流程对数据进行详细的分析。由于以下所有的流程都是基于参考序列进行的，所以选择合适的参考数据库(如NCBI、UCSC等公共数据库的基因组序列和cDNA序列)十分重要。

数据过滤，由于某些原始测序序列带有接头(adaptor)序列或含有少量低质量序列，首先需要经过一系列数据过滤以去除杂质数据，原始序列数据经过去除杂质后得到的数据称为clean reads，后续分析都基于clean reads。过滤按以下方式进行：截掉adapter接头及之后的序列；截掉序列末尾质量低于20的碱基；丢掉序列长度小于16的序列；去掉50％碱基质量低于20的序列。

接下来对各个转录组的待分析数据进行基础比对：

1.1序列比对

本发明使用tophat2进行序列比对，tophat可以将剪掉过内含子的序列比对到基因组上，即进行可变剪切的比对。对于nextseq500的数据，由于数据读长一般为150，所以在参数的选择上需要read的mismatch要设置稍大些，本发明中read的错配数设为4，设置Bowtie2片段比对最大错配数为1；设置reads最大的多位置比对结果输出个数为2；设置线程数为16；其他均使用软件默认设置。也即，比对的参数为：-b2-N 1-read-mismatches 2-a8-m 0-g 2-p 16。对每个样品按照设定的参数进行比对，比对的结果生成对应的文件夹下，比对结束后对比对结果进行筛选，提取出uniq map的reads，同时对筛选结果进行统计，包括total raw reads，比对效率，uniq map reads数，multiple map reads数。提取的方法如下：检查bam文件中每条比对结果TAG的NH，如果匹配“NH:i:1”，则表示该reads是唯一比对的结果，保留下来，否则就丢掉。最后筛选留下的结果使用samtools工具转换为bam文件，并建立index。该bam用于后续分析。Tophat2可以提取出发生剪接的reads比对结果，并生成bed文件：junctions.bed，该文件是后序分析的输入文件。

1.2表达量计算

首先读取基因组注释文件(gff文件)，按照设定的window大小对处于同一window的gene进行cluster；然后依次读取uniq map的比对结果，将read比对位置和所在cluster的gene做比较，如果落在gene的exon上，则将gene的reads数按比例增加，同时统计gene的exon长度，intron长度，coverage，depth，uniq map reads数等信息，最后根据RPKM的公式计算gene的RPKM表达量，RPKM的公式为：

RPKM为基因的表达量，total exon reads为唯一比对到基因A的总reads数，mapped reads为唯一比对到基因组上的所有的reads数，exon length为基因A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

1.3差异基因分析

该分析的目的是找到两个样本间或两个样本组间发生显著表达差异的基因。本发明使用R的edgeR包对样本进行差异比较。EdgeR主要基于负二项分布模型，被主要用做RNA-seq数据的表达量分析。首先将表达量分析中获取的每个样本落在每个基因exon区域(与gene的任一exon有overlap)的reads整理成一个总表，行为基因，列为样本。然后使用edgeR包的exactTest(精确检验)函数计算表达量变化，然后通过结果中的fold-change和pvalue对结果进行筛选，获得显著差异的gene及其差异信息。fold-change的阈值为2，即两个样本间的表达量倍数要大于2或小于0.5，p value的阈值选择为0.01。p value的值越小表示差异越显著。

1.4功能注释分析

得到差异表达基因之后，对差异表达基因做聚类分析、GO富集分析和KEGGPathway分析。

聚类分析给出差异表达基因的功能分类注释；GO富集分析给出差异表达基因的GO功能显著性富集分析。聚类分析给出具有某个功能的基因列表及基因数目统计。GO富集分析给出与基因组背景相比，在差异表达基因中显著富集的GO功能条目，从而给出差异表达基因与哪些生物学功能显著相关。在一个实施例中，聚类分析和GO富集分析也可以整合到GO功能分析中，以方便地分析具有某一功能的所有差异基因的表达模式。GO功能分析首先把所有差异表达基因向Gene Ontology数据库(http://www.geneontology.org/)的各个term映射，计算每个term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的GO条目，其计算公式为：

其中，N为所有基因中具有GO注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定GO term的基因数目；m为注释为某特定GO term的差异表达基因数目。计算得到的pvalue通过Bonferroni校正之后，以corrected p value≤0.05为阈值，满足此条件的GO term定义为在差异表达基因中显著富集的GO term。通过GO功能显著性富集分析能确定差异表达基因行使的主要生物学功能。

不同基因相互协调行使其生物学功能，Pathway分析有助于更进一步了解基因的生物学功能。KEGG是有关Pathway的主要公共数据库，Pathway显著性富集分析以KEGGPathway为单位，应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著性富集的Pathway。该分析的计算公式同GO功能显著性富集分析，在这里N为所有基因中具有Pathway注释的基因数目；n为N中差异表达基因的数目；M为所有基因中注释为某特定Pathway的基因数目；m为注释为某特定Pathway的差异表达基因数目。FDR≤0.05的Pathway定义为在差异表达基因中显著富集的Pathway。通过Pathway显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。

【实施例2】可变剪接分析

各种可变剪接事件的鉴定算法(如图3)：

1、基因模型：

可变剪接是一种相对的事件，一个可变剪接事件中至少含有两种剪接型，一种剪接型相对于另一种是可变的。由于许多注释文件中一个基因不只一个转录本，所以为了方便研究可变剪接这种相对性，我们从每个基因中都会选取一个转录本作为参考模型，即我们可变剪接研究的基因模型，我们认为这种模型是发生的，如果检测出支持其他剪接型的证据，如新的剪接位点，则认为这里可能发生了一个可变剪接事件。

2、外显子跳跃事件(ES)：

在基因模型中含有N个外显子，N>＝3，若存在剪接点跨越在外显子[0]和外显子[N-1]之间，外显子[1]—外显子[N-2]被跳跃，则可以预测这里发生了外显子跳跃事件，存在一个跨越了N-2个外显子的剪接型。

3、外显子跳跃事件(CE):

也是一种外显子跳跃事件，与ES的区别在于:此处基因模型含有两个外显子：A和B。存在剪接点j 1以外显子A的结尾为起始，剪接点j2以外显子B的起始为结尾，j1和j2没有重叠，且有序列落在j1的结尾和j2的起始之间。则我们认为在外显子A和B之间发现了一个外显子是在基因模型中被跳跃掉的，存在一个cassette外显子事件。

4、互斥外显子：

即外显子互斥事件，基因模型中含有三个外显子：A,B和C，使用AC应用上面CE中判断新外显子的方法，如果发现AC间存在新的外显子，且该外显子在AB之间或者在BC之间，则可以判断这里存在一个互斥外显子事件，外显子B和新发现的外显子就是互斥的两个外显子。

5、可变3’剪接位点:

基因模型含有两个外显子，存在sj 1以外显子A的结尾为开始，sj1的末端落在内含子区域或外显子B的内部，则可能发生了一个可变3’剪接位点事件，这里要求sj 1在寻找其他的可变剪接事件(除可变5’剪接位点)时没有用到过。

6、可变5’剪接位点:

基因模型含有两个外显子，存在sj1以外显子B的开始为结尾，sj1的起始落在内含子区域或外显子A的内部，则可能发生了一个可变5’剪接位点事件，这里要求sj1在寻找其他的可变剪接事件(除可变3’剪接位点)时没有用到过。

7、可变的第一个外显子:

在基因模型中外显子A是转录本5’端第一个外显子，如果存在剪接点j 1以外显子B的开始为结尾，j1的开始在外显子A的左侧，则认为存在一个可变的第一个外显子事件。

8、可变的最后一个外显子:

在基因模型中外显子A是转录本3’端第一个外显子，如果存在剪接点j 1以外显子B的结尾为开始，j1的结尾在外显子A的右侧，则认为存在一个可变的最后一个外显子事件。

9、内含子保留事件:

以整个转录本为基因模型，考虑其中的一个内含子，如果满足以下条件，则认为此处发生了一个内含子保留事件：

a、内含子区碱基平均深度要大于总内含子区碱基平均深度的2倍；

b、内含子左右两边边界序列都要大于1；

c、内含子区碱基平均深度要大于左右两边外显子区碱基平均深度的15％；

基于上面提到的可变剪接检测算法模型，整个可变剪接的具体流程描述如下：

1.序列比对和剪接点检测：

现在有许多已经很成熟的工具用来做序列比对和剪接点检测，例如tophat/tophat2,mapsplice等等，这些软件找到的剪接点具有较高的正确率和可信度。例如在我们的流程所用到的tophat2，它检测到的Correct junction序列比例可达到95％以上。这里我们使用tophat2将clean序列比对到基因组，获取到整个转录组比对的情况以及有tophat2检测到的剪接点s。根据基因组注释信息我们使用一个自己编写的perl脚本将这些剪接连接分类为已知的剪接连接，即在注释信息中出现的外显子-外显子连接；新的剪接连接。这些新的剪接连接则代表可能有新的可变剪接事件发生。

2.已知可变剪接事件检测：

首先将基因组注释信息中的剪接点提取出来和已经注释的gene(转录本)信息作为输入，这里假设每个剪接点的序列数为5，将已知基因组注释信息中包含的可变剪接事件检测出来并分类，已知的类别包括SE，可变3’剪接位点,可变5’剪接位点,可变的第一个外显子,可变的最后一个外显子,内含子保留,互斥外显子。ABLas1算法对于已知AS事件的检测有些特别之处：(1)因为注释文件中的外显子是确定的，所以算法中关于推测某区域是否为外显子直接查看注释文件即可。(2)互斥外显子、可变3’剪接位点、可变5’剪接位点事件需要同时找到互为模型的事件。在确定所有的已知可变剪接事件后，计算每个AS事件支持两种剪接型的剪接点是否在该实验数据比对得到的剪接点中包含以及检测出的剪接点序列数。如果支持两种剪接型的剪接点满足阈值条件，则认为这个可变剪接事件在该数据中被检出。阈值条件的设置为：可变型支持的junction reads至少为2，model型支持的junction reads至少为2，两者相加至少为10。最后可以得到数据中检出的所有的在注释信息中已经包含的可变剪接事件库。利用这个库可以分析注释文件中已有的可变剪接事件在数据中表达情况。

3.新的可变剪接事件检测：

新的可变剪接事件的检测分为两个部分。首先将转录本中新的剪接点以及注释文件作为输入可以得到非内含子保留的AS事件。另外将转录本的比对结果作为输入，可检测内含子保留事件。将两个算法得到的结果合并则是原始的新的可变剪接事件。再经过筛选和去重复则得到最后的新的可变剪接事件。筛选的标准如下：可变型支持的junctionreads至少为2，可变型占可变剪接事件总比例的15％。

4.可变剪接事件的可视化展示：

本发明开发了一个针对可变剪接的基因组可视化软件，软件读取基因组注释文件和基础分析中tophat2唯一比对的bam结果文件，以及检测出的可变剪接事件，将reads分布及可变剪接以可视化的形式进行展示。图5是本发明一个实施例中产生的外显子跳跃事件可视化图。图6是本发明一个实施例中产生的内含子保留事件可视化图。图7是本发明一个实施例中产生的可变3’剪接位点事件和互斥外显子事件可视化图。

5.可变剪接调控分析：

在检出所有的已知和新的可变剪接事件后，通过观察可变剪接事件的两种可变型在不同样本中的表达比例，可以发现样本间可变剪接模式的变化，及其所介导的基因表达的变化。图4是可变剪接调控分析的流程图。在非内含子保留事件中，我们使用支持每种剪接型的junction序列数计算剪接型表达的密度，通过两种剪接型密度的比较得到每种剪接型表达的比例。剪接型密度的计算通过用junction序列数除以junction序列起始位点可能落在的所有位点总数，例如在外显子跳跃事件中，如图9所示，序列长度为rl,junction的overhang为ol,则每个junction start可能落下的位点长度为rl-ol，type1的密度ψ1为type2的密度ψ2为Type1的密度比例为两个样本中type1密度比例的差值Ψ_diff＝Ψ_sample1-Ψ_sample2即表示剪接型tpye1在两样本中表达比例的差异大小。Ψ_diff越大(绝对值)说明在两个样本中该可变剪接事件发生的调控越明显。

另外如何衡量Ψ_diff的可信度，因为如果j1,j2,j3的值偏小，则观察值就没有统计意义，从而得到的Ψ_diff存在较大的错误率，因此我们使用fisher exact test方法来计算Ψ_diff的p value值。

fisher exact test的四格表如下：

【实施例3】

图8是本发明一种基于RNA-seq数据的可变剪接分析系统的一个实施例的框图。如图8所示，分析系统包括数据处理模块11，用于过滤各样本原始测序数据中质量不合格的序列，获得各个转录组的clean reads；比对模块12，与数据处理模块11相连，用于将各个转录组数据比对到参考基因组上，并提取唯一比对序列；表达量分析模块13，与比对模块12相连，用于计算基因的表达量RPKM值；差异基因分析模块14，与表达量分析模块13相连，用于筛选发生显著表达差异的基因；功能注释分析模块15，与差异基因分析模块14和可变剪接事件差异分析模块18相连，用于对筛选出的基因做功能注释分析，包括GO分析和KEGG分析；已知可变剪接事件检测和注释模块16，与比对模块12相连，用于将注释文件中已知的可变剪接事件进行提取，并且进行定量；新的可变剪接事件检测和注释模块17，与已知可变剪接事件检测和注释模块16相连，用于鉴定注释文件中没有的新的可变剪接事件，并进行定量和筛选；可变剪接事件差异分析模块18，与新的可变剪接事件检测和注释模块17及差异基因分析模块14相连，用于计算和筛选在样本间剪接发生显著变化的差异可变剪接事件，并与表达差异基因做overlap分析；可变剪接事件统计模块19，与可变剪接事件差异分析模块18相连，用于对上面所检出的已知可变剪接事件，新的可变剪接事件，差异事件等做统计，生成总表；可变剪接事件全基因组可视化模块20，与可变剪接事件统计模块19相连，用于展示reads分布及各样本可变剪接情况。

与现有技术相比，本发明具有的有益效果有：

(1)能够准确的获取splice junction信息。

(2)能够快速准确的鉴定已知和新的可变剪接事件，在鉴定和筛选新的可变剪接事件时使用的参数是根据多年可变剪接研究经验得到的优化值。

(3)对可变剪接事件进行了定量，计算获取可变型发生的比率。

(4)针对两个以上样品，可以进行可变剪接差异分析，以获取有可变剪接进行调控的基因，进而研究可变剪接调控导致的功能变化。

(5)可视化展示可变剪接事件。

Claims

1.一种基于RNA-seq数据的可变剪接分析方法，其特征在于，包括：

1）通过illumina二代测序平台获取某一具有参考基因组和注释的真核生物的一个或多个样品的转录组原始测序数据；

2）对上述各组原始测序数据进行过滤，将质量不合格的数据过滤掉，留下的数据作为待分析的数据，过滤的条件是：截掉adapter接头及之后的序列；截掉序列末尾质量低于20的碱基；丢掉序列长度小于16的序列；去掉50%碱基质量低于20的序列；

3）对各个转录组的待分析数据进行基础分析和可变剪接分析，其中，所述的基础分析包括：

（1）将所述各个转录组样本待分析数据分别比对到所述物种的参考基因组，获取发生剪接的比对结果，并筛选出唯一比对的结果；

（2）计算各样本基因的表达量：基于RPKM标准化方法使用python编写程序，计算基因表达量信息；

（3）将各样品按照样品间或样品组间进行差异分析，筛选出显著差异表达的基因：样本（组）间差异分析使用R软件包edgeR进行，显著差异基因的筛选标准为：p value小于等于0.01，fold change大于等于2；

（4）对差异基因进行功能注释和分析：包括样品间相关性分析，差异基因聚类分析，差异基因GO富集分析，差异基因KEGG Pathway分析；

所述的可变剪接分析包括：

（1）参考基因组注释文件中已知可变剪接事件的鉴定；

（2）新的可变剪接事件的鉴定；

（3）样品（组）间可变剪接事件差异分析；

（4）可变剪接与基因表达关联分析；

（5）可变剪接分析结果统计和报表生成；

（6）可变剪接可视化图生成：使用perl编写程序，绘制可变剪接事件的可视化图。

2.根据权利要求1所述的方法，其特征在于，所述的基础分析的比对使用tophat2软件进行，软件的参数具体设置如下：设置比对reads 的错配数为4 ；设置Bowtie2片段比对最大错配数为1；设置reads 最大的多位置比对结果输出个数为2 ；设置线程数为16；其他均使用软件默认设置。

3.根据权利要求1所述的方法，其特征在于，各个转录组样本待分析数据分别比对到所述物种的参考基因组得到结果后，筛选出唯一比对结果的方法如下：检查bam文件中每条比对结果TAG的NH，如果匹配“NH:i:1”，则表示该reads是唯一比对的结果，保留下来，否则就丢掉；

最后筛选留下的结果使用samtools工具转换为bam文件，并建立index；该bam用于后续分析；Tophat2可以提取出发生剪接的reads比对结果，并生成bed文件：junctions.bed，该文件是后序可变剪接分析的输入文件。

4.根据权利要求1所述的方法，其特征在于，可变剪接分析中研究的可变剪接事件包括的类别如下：外显子跳跃事件（ES/cassetteExon），互斥外显子事件，可变3’剪接事件，可变5’剪接事件，可变的第一个外显子事件，可变的最后一个外显子事件，同时外显子跳跃和可变3’剪接的事件，同时外显子跳跃和可变5’剪接的事件，内含子保留事件。

5.根据权利要求1所述的方法，其特征在于，参考基因组注释文件中已知可变剪接事件的鉴定步骤为：首先为每个基因定义一个基因模型，也就是gene model，默认选择注释文件中该基因第一个转录本作为gene model；将gene model的剪接位置（splice junction）提取出来，及exon与exon的连接位点信息，再将所有基因的剪接位置提取出来，再使用自己编写的perl程序分析这些剪接位置，鉴定出相对于gene model存在的可变剪接事件并进行筛选和分类，按照规范格式进行输出；已知可变剪接事件的筛选条件如下：对于已知的非内含子保留事件，相对于model可变的剪接型支持reads要大于等于2，model的支持reads要大于等于2，两者加起来要大于等于10。

6.根据权利要求1所述的方法，其特征在于，新的可变剪接事件的鉴定的步骤如下：以gene model作为参考，从tophat2检测出的发生可变剪接的比对结果中筛选出新的splicejunction结果，使用自己编写的perl程序鉴定相对于gene model发生的可变剪接事件，并进行筛选和分类，按照规范格式进行输出；新的可变剪接事件的筛选条件如下：相对于model可变的剪接型支持reads要大于等于2，该剪接型比例需大于等于0.15。

7.根据权利要求1所述的方法，其特征在于，样品（组）间可变剪接事件差异分析的步骤如下：分别计算每个样品每个可变剪接事件发生的比例，计算样品间该比例的差值，使用费雪尔正确性检定（Fisher’s Exact Test）或t检验（T-Test）计算差异的可靠性，即得到pvalue，差异可变剪接事件的筛选条件如下：样品间比例差值需大于等于0.2，p value需小于等于0.05。

8.一种基于RNA-seq数据的可变剪接分析的系统，其特征在于，包括：

分析系统包括数据处理模块11，用于过滤各样本原始测序数据中质量不合格的序列，获得各个转录组的clean reads；比对模块12，与数据处理模块11相连，用于将各个转录组数据比对到参考基因组上，并提取唯一比对序列；表达量分析模块13，与比对模块12相连，用于计算基因的表达量RPKM值；差异基因分析模块14，与表达量分析模块13相连，用于筛选发生显著表达差异的基因；功能注释分析模块15，与差异基因分析模块14和可变剪接事件差异分析模块18相连，用于对筛选出的基因做功能注释分析，包括GO分析和KEGG分析；已知可变剪接事件检测和注释模块16，与比对模块12相连，用于将注释文件中已知的可变剪接事件进行提取，并且进行定量；新的可变剪接事件检测和注释模块17，与已知可变剪接事件检测和注释模块16相连，用于鉴定注释文件中没有的新的可变剪接事件，并进行定量和筛选；可变剪接事件差异分析模块18，与新的可变剪接事件检测和注释模块17及差异基因分析模块14相连，用于计算和筛选在样本间剪接发生显著变化的差异可变剪接事件，并与表达差异基因做overlap分析；可变剪接事件统计模块19，与可变剪接事件差异分析模块18相连，用于对上面所检出的已知可变剪接事件，新的可变剪接事件，差异事件等做统计，生成总表；可变剪接事件全基因组可视化模块20，与可变剪接事件统计模块19相连，用于展示reads分布及各样本可变剪接情况。