CN104657628A - 基于Proton的转录组测序数据的比较分析方法和系统 - Google Patents

Abstract

本发明提供了一种基于Proton的转录组测序数据的比较分析方法和系统。其中，该方法包括通过Proton测序平台获取某一物种至少两个转录组的原始测序数据；过滤掉不合格数据，获得clean reads；然后进行第一步分析和第二步分析，其中，所述第一步分析包括：将clean reads分别比对到所述物种的参考基因组；进行转录本定量分析；筛选出显著差异表达的基因；进行显著差异表达基因功能分析；所述第二步分析包括：将clean reads分别比对到所述物种的参考基因组；进行可变剪接分析；进行基因结构优化分析；进行新转录本预测；进行SNP分析；进行基因融合检测。本发明的方法和系统可以提高转录组测序数据比较分析的准确度和可靠性。

Description

基于Proton的转录组测序数据的比较分析方法和系统

技术领域

本发明涉及生物信息技术领域，尤其涉及一种基于Proton的转录组测序数据的比较分析方法和系统。

背景技术

转录组(transcriptome)广义上指某一生理条件下，细胞内所有转录产物的集合，包括信使RNA(mRNA)、核糖体RNA、转运RNA及非编码RNA；狭义上指所有信使RNA的集合。由于转录组代表了生物在某一时刻的基因表达状态，因此，对转录组的研究具有极大的生物学意义。

基因可以具有产生mRNA的多个不同拼接的转录本，且mRNA也可以由基因组上的其他区域产生。测序技术是用于调查和理解mRNA表达谱的有力工具。相对于传统的芯片杂交平台，转录组测序无需预先针对已知序列设计探针，即可对任意物种的整体转录活动进行检测，提供更精确的数字化信号，更高的检测通量以及更广泛的检测范围，且还允许发现新的序列变异和转录本的可能性，是目前深入研究转录组复杂性的强大工具。基于Ion Proton平台的RNA-Seq产品可以实现在2～4个小时完成转录组测序，产出高达10Gb的产量，读长为100-200bp。测序成本方面，Ion Proton摆脱了昂贵的光学测序的束缚，测序成本大大降低。较低的测序成本，更快的测序速度、更长的数据片段以及良好的一致准确性使得Ion Proton平台的转录组测序在生命科学研究领域，尤其是临床诊断方面的优势十分明显。

由于基于Proton的mRNA测序技术是高吞吐量的并产生相当大量的基因数据，mRNA测序数据可能需要使用各种计算技术以根据所获得的的大量数据对测序的mRNA进行准确和有意义的推断。此外，来自不同试验(其可以是来自不同组织或不同实验条件下的样本)的mRNA测序数据可能还需要被比较，分析相当大量的mRNA的数据是相当复杂的过程。对不同试 验之间的mRNA测序数据的分析检测会受到测序平台特征的影响，也受到用于分析数据的计算技术的影响，需要大量的人工检查和验证后续工作。

发明内容

针对现有技术中存在的问题，本发明的主要目的在于提供一种基于Proton的转录组测序数据的比较分析方法，可以提高转录组测序数据比较分析的准确度和可靠性。

一方面，本发明提供了一种基于Proton的转录组测序数据的比较分析方法，包括：

通过Proton测序平台获取某一物种至少两个转录组的原始测序数据；

过滤所述各组原始测序数据中的不合格数据，获得所述各个转录组的待分析数据(即clean reads，它是原始测序数据经过过滤处理后留下的数据)；

对所述各个转录组的待分析数据进行第一步分析和第二步分析，其中，所述第一步分析包括：

(1)将所述各个转录组的待分析数据分别比对到所述物种的参考基因组；

(2)利用RPKM标准化方法对所述各个转录组的待分析数据进行转录本定量分析；

(3)对转录本定量分析后的所述各个转录组的待分析数据进行两两差异比较，并根据设定的阀值筛选出显著差异表达的基因；

(4)对所述显著差异表达的基因进行功能分析；

所述第二步分析包括：

(1)选取合适的参数，将所述各个转录组的待分析数据分别比对到所述物种的参考基因组；

(2)可变剪接分析：得到各个基因的可信的剪接位点，对剪接位点信息进行计算，获得含有可变剪接形式的基因注释结果；和/或

基因结构优化分析：根据比对结果进行基因结构的优化分析；和/或

新转录本预测：根据比对结果进行新转录本的预测；和/或

SNP分析：根据比对结果进行SNP预测；和/或

基因融合检测：根据比对结果进行融合基因检测。

在本发明的一个实施例中，所述基因结构优化分析的步骤为：将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与所述物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的重叠部分(overlap)超过阀值，且与其对应位置的已知转录本覆盖度大于阀值，则将该转录本作为最终结果，往其已知的注释结果的两端延伸。在本发明的又一个实施例中，基因结构优化分析中过滤掉组装结果小于阀值的转录本中的阀值为180bp；如果某个转录本与其对应的基因上游或下游的重叠部分(overlap)超过阀值中的阀值为200bp；且与其对应位置的已知转录本覆盖度大于阀值中的阀值为0.9。

在本发明的又一个实施例中，所述新转录本预测的步骤为：将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与所述物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的重叠部分(overlap)没有达到阀值，且满足一定的测序深度，则称为新转录本。在本发明的又一个实施例中，新转录本预测中过滤掉组装结果小于阀值的转录本中的阀值为180bp；如果某个转录本与其对应的基因上游或下游的重叠部分(overlap)超过阀值中的阀值为200bp。

在本发明的又一个实施例中，不合格数据包括序列长度低于30的reads；修剪掉测序接头后序列长度低于30的reads；修剪掉低质量的碱基序列后序列长度小于30的reads，其中所述低质量碱基序列是指从3’端开始，以15为窗口统计reads平均质量，若平均质量值小于9，则继续向5’端滑动，直到窗口平均质量大于9时该窗口3’端到该reads 3’端的所有碱基序列。

在本发明的又一个实施例中，第一步分析中的比对是使用Tmap软件进行的，所述第一步分析中的所述两两差异比对是使用泊松分布模型进行的。

在本发明的又一个实施例中，利用所述Tmap软件比对时，参数的具体设置如下：设置比对结果输出类型a为2，即为所有最好的位点；设置程序运行时的线程数为1～16；设置输出文件类型为压缩的bam格式；设置比对策 略为map4，即为bwa快速比对。

在本发明的又一个实施例中，第一步分析中的功能分析包括聚类分析、GO富集分析和/或Pathway分析。

在本发明的又一个实施例中，第二步分析中的所有比对是使用tophat软件进行的。

在本发明的又一个实施例中，利用所述tophat软件比对时，参数的具体设置如下：设置比对reads的错配数为0～5；设置reads开gap最大的缺口数为0～20；设置reads的最大的编辑距离为5～20；设置reads剪接比对时最大错配数为0～4；设置reads比对最大插入长度为1～10；设置reads比对最大缺失长度为1～10；设置片段的比对最大错配数为0～5；设置Bowtie2比对内设参数为1～3；设置线程数为1～16。

在本发明的又一个实施例中，利用所述tophat软件比对时，参数的具体设置如下：设置比对reads的错配数为3～5；设置reads开gap最大的缺口数为6～15，设置reads的最大的编辑距离为10～20，设置reads剪接比对时最大错配数为2～4；设置reads比对最大插入长度为5～10；设置reads比对最大缺失长度为5～10；设置片段的比对最大错配数为2～5；设置Bowtie2比对内设参数为1～3；设置线程数为6～16。

在本发明的又一个实施例中，可变剪接分析中的对剪接位点信息进行计算包括以下步骤：利用自己编写的perl脚本对剪接位点信息进行计算。

在本发明的又一个实施例中，基因结构优化分析中和新转录本预测中reads组装是使用cufflink软件进行的；所述SNP分析中SNP检测是使用pyrohmmvar软件进行的；所述基因融合检测是使用tophat-fusion软件进行的。

在本发明的又一个实施例中，使用所述cufflink软件组装时，参数的具体设置如下：设置参数u，对多reads进行矫正；设置线程数p为1～16。

在本发明的又一个实施例中，使用所述pyrohmmvar软件进行SNP检测时，参数的具体设置如下：滑动窗口的长度设为3～30；最小的变异质量值分数设为20～50；最小碱基质量值设为5～20；最小的reads长度设为15～30；设置忽略indels I；设置最小比对打分为0；设置最小的比对一致性为 0.50～0.90；设置用于计算SNP最小支持数reads的系数为0.1～0.3；设置支持SNP的最大覆盖度为400～1000。

在本发明的又一个实施例中，使用所述tophat-fusion软件进行融合基因检测时，参数的具体设置如下：设置线程数为1～16；设置跳过blast比对，即skip-blast；设置支持断点最小的reads数目为1～10；设置支持断点的最小的成对reads数目为0。

另一方面，本发明还提供了一种基于Proton测序平台转录组测序数据的比较分析的系统，其特征在于，包括：

序列接收模块，用于获取某一物种至少两个转录组的基于Proton的测序平台原始测序数据；

序列过滤模块，与序列接收模块相连，用于过滤所述各组原始测序数据中的不合格数据，获得所述各个转录组的clean reads；

第一序列比对模块，与序列过滤模块相连，用于将所述各个转录组的clean reads分别比对到所述物种的参考基因组；

定量分析模块，与第一序列比对模块相连，用于利用RPKM标准化方法对所述各个转录组的clean reads进行转录本定量分析；

差异表达筛选模块，与定量分析模块相连，用于对转录本定量分析后的所述各个转录组的clean reads进行两两差异比较，并根据设定的阀值筛选出显著差异表达的基因；

功能分析模块，与差异表达筛选模块相连，用于对所述显著差异表达的基因进行功能分析；

第二序列比对模块，与序列过滤模块相连，用于将所述各个转录组的clean reads分别比对到所述物种的参考基因组；

可变剪接基因注释模块，与第二序列比对模块相连，用于得到各个基因的可信的剪接位点，对剪接位点信息进行计算，获得含有可变剪接形式的基因注释结果；和/或

基因结构优化分析模块，与所述第二序列比对模块相连，用于根据比对结果进行基因结构的优化分析；和/或

新转录本预测模块，与所述第二序列比对模块相连，用于根据比对结 果进行新转录本的预测；和/或

SNP分析模块，与第二序列比对模块相连，用于根据比对结果进行SNP预测；和/或

基因融合检测模块，与第二序列比对模块相连，用于根据比对结果进行融合基因检测。

本发明的基于Proton的转录组比较分析方法以不同组织或者不同实验条件下的个体为研究对象，从整体水平研究基因功能以及基因结构，揭示特定生物学过程或者疾病发生过程中的分子机理，可广泛应用于基础研究、临床诊断和药物研发等领域。相对于其他的测序平台和数据比较分析方法，本发明的方法可以提高转录组测序数据比较分析的准确度和可靠性。

附图说明

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

图2是本发明一种基于Proton的转录组测序数据的比较分析方法的一个实施例的流程示意图。

图3是本发明可变剪接分析中Exon Skipping(外显子跳读)算法示意图。

图4是本发明可变剪接分析中Intron Retentio(内含子保留)算法示意图。

图5是本发明可变剪接分析中Alternative 5’Splice Site(可变5’剪接)算法示意图。

图6是本发明可变剪接分析中Alternative 3’Splice Site(可变3’剪接)算法示意图。

图7是本发明基因结构优化分析流程示意图。

图8是本发明新转录本预测流程示意图。

图9是本发明一种基于Proton的转录组测序数据的比较分析系统的一个实施例的框图。

具体实施方式

以下结合附图和优选实施例对本发明进行更全面地描述，应当理解的是，本文所描述的优选实施例仅用于说明和解释本发明，并不构成对本发明的限制。

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

1.1测序

采集不同组织或不同实验条件下的样品，对其进行测序得到比较分析所需的数据。在一个实施例中，采用Ion Proton平台进行测序，测序基本过程如图1所示：

1)总RNA提取：提取样品总RNA，使用DNaseⅠ酶消化总RNA中的DNA；

2)mRNA富集：对于真核生物，用Oligo(dT)的磁珠富集mRNA；

3)mRNA片段化：将mRNA与适量打断试剂相混合，高温作用一定时间进行打断；

4)cDNA合成：以片断后的mRNA为模板，用六碱基随机引物进行反转录，合成双链cDNA；

5)末端修复：用末端修复复合酶对双链cDNA进行末端修复，修复产物用磁珠纯化，连接接头，连接产物片段进行胶纯化回收；

6)PCR扩增：配置PCR反应体系对回收的连接产物进行扩增，并用磁珠纯化，文库构建完成；

7)文库质控：使用Agilent 2100对文库的大小和浓度进行检测；

8)上机测序：使用Ion Proton平台进行测序。

1.2转录组比较分析

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

在本发明一个实施例中，比较分析流程如图2所示。

S101过滤，由于某些原始测序序列带有接头(adaptor)序列或含有少量低质量序列，首先需要经过一系列数据过滤以去除杂质数据，原始序列数据经过去除杂质后得到的数据称为clean reads，后续分析都基于clean reads。在本发明的一个实施例中，过滤按按以下方式进行：去除某些含adaptor的reads的adaptor，如果去除adaptor后reads的序列长度低于30，那就过滤掉该reads；从3’端开始，以15为窗口统计reads的平均质量，若平均质量值小于9，则继续向5’端滑动，直到窗口平均质量大于9为止，修剪掉该窗口3’端到该reads的3’端的所有碱基序列，若修剪后的reads序列长度低于30，那就过滤掉该reads。

接下来，进行第一步分析：

S102Tmap比对，Tmap软件主要用作短序列的比对，最大的特点是支持长短不一的reads,并支持开gap的快速比对，这些特点是其它很多软件所没有的。我们使用Tmap比对软件(https://github.com/iontorrent:/TMAP)分别将clean reads比对到参考基因组上，通过这一步骤，可以将测序得到的reads对应到基因及基因组上，第一步分析中的后续分析都是基于上述比对结果。为了获得良好的比对效果，加快比对速度，节省硬盘内存，在利用Tmap比对软件比对到基因或基因组时，参数的具体设置如下：设定比对结果输出类型a为2，即为所有最好的位点；设定程序运行时设置的线程数为1～16，优选为6；设定输出文件类型为1，以便输出bam压缩格式文件，一方面能节省硬盘存储，另一方面能跟下游的比对软件直接对接；设定比对策略为map4，设置map4比对策略是因为该参数内嵌BWA比对软件，比对算法为巴罗斯比对算法，较为快速，更适合短序列比对。也即，优选地，比对的详细参数为：a 2-n 6-v-Y-u-o 1 stage1 map4。Proton平台测序作为一个新的测序平台，支持其数据分析的软件很少甚至没有，目前在分析这点上选择Tmap，是因为它的参数设置较其他软件更加灵活，选择余地更大。前期的测试表明用它进行比对及定量分析其结果是可靠和 科学的。

S103表达定量，使用RPKM法(Reads Per Kb per Million reads)计算基因的表达量，其计算公式为：

$\mathrm{RPKM}=\frac{{10}^{6}C}{\mathrm{NL}/{10}^3}$

设RPKM(A)为基因A的表达量，则C为唯一比对到基因A的reads数，N为唯一比对到基因组的总reads数，L为基因A编码区的碱基数。RPKM法能消除基因长度和测序量差异对计算基因表达的影响，计算得到的基因表达量可直接用于比较不同样品间的基因表达差异。

如果一个基因存在多个转录本，则用该基因的最长转录本计算其测序覆盖度和表达量。

S104差异基因筛选，差异基因筛选旨在找出不同样本间存在的差异表达的基因，在一个实施例中，使用了泊松分布模型对样本进行两两差异比较。假设观测到基因A对应的reads数为x，已知在一个大文库中，每个基因的表达量只占所有基因表达量的一小部分，在这种情况下，p(x)的分布服从泊松分布：

(λ为基因A的真实转录数)

已知，样本1中唯一比对到基因组的总reads数为N1，样本2中唯一比对到基因组的总reads数为N2，样本1中唯一比对到基因A的总reads数为x，样本2中唯一比对到基因A的总reads数为y，则基因A在两个样本中表达量相等的概率可由以下公式计算：

$2{\mathrm{\Σ}}_{i=0}^{i-y}p\left(i\right|x)$

或 $2\×(1-{\mathrm{\Σ}}_{i=0}^{i-y}p\left(i\right|x))$ (如果 ${\mathrm{\Σ}}_{i=0}^{i-y}p\left(i\right|x)>0.5$ ) 

$p\left(y\right|x)={\left(\frac{N_2}{N_1}\right)}^y\frac{(x+y)!}{x!y!{(1+\frac{N_2}{N_1})}^{(x+y+1)}}$

然后，对差异检验的p value作多重假设检验校正，通过控制FDR(False  Discovery Rate)来决定p value的域值。假设挑选了R个差异表达基因，其中S个是真正有差异表达的基因，另外V个是其实没有差异表达的基因，为假阳性结果。希望错误比例Q＝V/R平均而言不能超过某个可以容忍的值，比如1％，则在统计时预先设定FDR不能超过0.01。在得到差异检验的FDR值同时，我们根据基因的表达量(RPKM值)计算该基因在不同样本间的差异表达倍数。FDR值越小，差异倍数越大，则表明表达差异越显著。在一个实施例的分析中，差异表达基因定义为FDR≤0.001且倍数差异在2倍以上的基因。

S105差异表达基因功能分析，在一个实施例中，得到差异表达基因之后，对差异表达基因做聚类分析、GO富集分析和KEGG Pathway分析。

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

$P=1-\underset{i=0}{\overset{m-1}{\mathrm{\Σ}}}\frac{\left(\begin{array}{c}M\\ i\end{array}\right)\left(\begin{array}{c}N-M\\ n-i\end{array}\right)}{\left(\begin{array}{c}N\\ n\end{array}\right)}$

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

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

接下来，进行第二步分析：

S106tophat比对，为了适应Proton测序平台的测序数据，以进行更好的比对，利用tophat比对软件与基因组比对时，参数的具体设置如下：设置比对reads的错配数为3～5；设置reads开gap最大的缺口数为6～15，设置reads的最大的编辑距离为10～20，设置reads剪接比对时最大错配数为2～4；设置reads比对最大插入长度为5～10；设置reads比对最大缺失长度为5～10；设置片段的比对最大错配数为2～5；设置Bowtie2比对内设参数为1～3；设置线程数为6～16。以上参数与比对的准确性密切相关，数值偏大或过小均会影响比对的准确性。尤其是过大的reads开gap最大的缺口数和reads的最大的编辑距离对比对的准确性影响较大，过大的reads开gap最大的缺口数或reads的最大的编辑距离可能会产生算法上假的比对结果，降低比对的准确性；由于proton本身测序数据有一定的indel，过小的reads开gap最大的缺口数或reads的最大的编辑距离可能会造成许多reads比对不上的情况，也会影响比对的准确性；各个参数值在此优选范围内时，可实现更佳的比对效果。

更优选地，参数的具体设置如下：设置比对reads的错配数，即N，为5；设置reads开gap最大的缺口数，即read-gap-length，为12；设置reads的最大的编辑距离，即read-edit-dist，为16；设置reads剪接比对时最大错配数，即splice-mismatches，为3；设置reads比对最大插入长度，即 max-insertion-length，为8；设置reads比对最大缺失长度，即max-deletion-length，为10；设置片段的比对最大错配数，即segment-mismatches，为3；设置Bowtie2比对内设参数，即b2-N，为2；设置线程数，即p，为15。考虑到proton测序数据本身存在一定的Indel，这里我们设置read-gap-length为12，read-edit-dist设置为16都比默认值要大，是为了开较大的gap，放松比对环境，提高reads的利用率。也即，比对的参数为：-N 5--read-gap-length 12--read-edit-dist 16--splice-mismatches 3--max-insertion-length 8--max-deletion-length 10--segment-mismatches 3--b2-N 2-p 15--b2-very-sensitive。关于参数的选取，这里我们使用hiseq测序的数据和proton测序的数据，并且设定梯度参数进行了评估，以保证可变剪切分析结果的准确性。该套参数为proton数据分析的宽松参数，能满足我们对数据分析的预期要求。关于tophat比对软件的使用具体可参考tophat的使用方法。

使用tophat比对软件分别将clean reads比对到参考基因组上，后续只要涉及到tophat与参考基因组比对，均使用该比对参数。通过这一步骤，可以将测序得到的reads对应到基因及基因组上，第二步分析中的后续分析都是基于上述比对结果。

S1071可变剪接分析，可变剪接使一个基因产生多个mRNA转录本，不同mRNA可能翻译成不同蛋白。因此，通过可变剪接一个基因可能产生多个蛋白，极大地增加了蛋白多样性。在将测序的reads比对到参考基因组上，得到可信的剪接(junction)位点后，利用自己编写的proton_AS_v1.3.pl程序根据junction位点信息进行计算，找出与已知的注释位点不同的四种可变剪接方式，其具体过程如下：

如图3所示，以转录本1为例，根据tophat比对找出的junction(剪接)结果，如果存在junction1(外显子1和外显子2拼接过)和junction2(外显子2和外显子3拼接过)，同时又满足存在junction3(外显子1和外显子3拼接过)的情况，则可判断发生外显子跳读剪接。

如图4所示，根据tophat比对找出的junction结果，如果存junction1(外显子1和外显子2拼接过)，并且根据比对统计得出的基因组深度统计文件，当 这段内含子满足：1)90％的区域存在覆盖；2)覆盖深度大于或等于两侧较小外显子的15％；3)该内含子与两侧外显子相邻5bp内必须有覆盖；4)该junction必须唯一存在，则可判断发生内含子保留的可变剪接。

如图5所示，根据tophat比对找出的junction结果，如果存junction1的类型(外显子1和外显子2拼接过)，同时存在junction2或者junction3的类型，这里只要外显子1的3’端终止位置不一致，且外显子2的起始和终止位置一致，则可判断发生可变5’端剪接。

如图6所示，根据tophat比对找出的junction结果，如果存junction1的类型(外显子2和外显子3拼接过)，同时存在junction2或者junction3的类型，这里只要外显子3的5’端起始位置不一致，且外显子2的起始和终止位置一致，则可判断发生可变3’端剪接。

使用软件“tophat”鉴定转录本的剪接位点(使用软件默认参数)，剪接位点给出了转录本不同外显子的边界及组合关系。下面，概述检测可变剪接的算法：

(1)外显子跳跃(Exon Skipping)剪接，如图3所示。转录本1和转录本2分别同时检测到如图3所示三个剪接位点，可认为转录本1的Exon1、Exon2和Exon3存在Exon Skipping剪接方式；转录本2的Exon1、Exon3和Exon4也存在Exon Skipping剪接方式。

(2)内含子保留(Intron Retention)剪接，如图4所示。1)检测到Junction1的存在，表明在某个成熟mRNA中Exon1和Exon2之间的Intron被剪接下来；2)Exon1和Exon2之间的Intron有90％以上的区域均有唯一比对reads覆盖，说明在某个成熟mRNA中该intron被保留下来了(考虑到转录的exon通常也不是100％被reads覆盖到，所以在这里以90％为阈值)。若同时满足以上两个条件，则认为该基因Exon1和Exon2之间存在Intron Retention的可变剪接方式。

(3)可选5’端剪接位点(Alternative 5’Splice Site)，如图5所示。如图5，一个转录本的Junction 1位点被检测到，并且Junction 2和Junction 3中有一个被检测到(它们共同点是3’剪接位点和Junction 1相同，但5’剪接位点和Junction 1不同)，那么就认为Exon1和Exon2存在Alternative 5’ Splice Site的剪接方式。

(4)可选择性3’端剪接位点(Alternative 3’Splice Site)，如图6所示如图6，一个转录本的Junction 1位点被检测到，并且Junction 2和Junction 3中有一个被检测到(它们共同点是5’剪接位点和junction 1相同，但3’剪接位点和junction 1不同)，那么就认为Exon1和Exon2存在Alternative 3’Splice Site的剪接方式。

S1072基因结构优化，通过比较转录本测序结果和现有基因注释结果，对基因的5’端或3'端进行延长。基因结构优化的步骤如图7所示，将reads比对到参考基因组后，采用cufflink软件组装比对上reads。过滤组装结果小于180bp的转录本，并与已知的转录本进行比对，过滤同源性较高(即相似性＞0.8)的转录本。将符合的转录本与基因组进行比对，如果与基因上下游200bp有overlap，并且与相应位置已知的转录本覆盖度大于0.9，则作为最终的结果对基因的5’端和3’端进行延长。基于cufflink的组装结果，选取合适长度的转录本及覆盖度，并根据其在基因组的位置，过滤假阳性结果，提高了准确性，并最终提供给客户可方便阅读的报告。使用cufflink软件进行组装时，参数的具体设置如下：设置参数u，以对多reads进行矫正；设置线程数p为1～16，优选为8。也即，优选地，具体参数为：-u-p 8。Cufflink软件能很好的与tophat比对结果相兼容，考虑到同一位置的reads可能存在较高的冗余，这里设置了参数对其进行矫正。

S1073新转录本预测，现有数据库中对转录本的注释可能还不全面，通过高通量测序还可能检测到新的转录本。新转录本预测的步骤如图8所示，采用cufflink软件组装比对上reads，过滤组装结果小于180bp的转录本，并与已知的转录本进行比对，过滤同源性较高的转录本。将符合的转录本与基因组进行比对，如果与基因上下200bp都没有overlap，且测序深度不小于2，则认为其为新的转录本。最后，用CPC软件预测该新转录本的编码能力，按设定的阈值将新转录本分成编码RNA和长的非编码RNA。

S1074SNP分析，在tophat比对软件将测序的reads比对到参考基因组后，根据比对结果，使用pyrohmmvar软件根据HMM模型对比对结果进行训练，推断出单体型，再根据贝叶斯模型，计算出潜在基因型的后验概率 值，值最高的便被推断为此位点的基因型。

使用所述pyrohmmvar软件进行SNP检测时，参数的具体设置如下：滑动窗口的长度设为3～30；最小的变异质量值分数设为20～50；最小碱基质量值设为5～20；最小的reads长度设为15～30；设置忽略indels I；设置最小比对打分为0；设置最小的比对一致性为0.50～0.90；设置用于计算SNP最小支持数reads的系数为0.1～0.3；设置支持SNP的最大覆盖度为400～1000。以上参数与检测的准确性密切相关，尤其是最小碱基质量值和支持SNP的最大覆盖度对检测准确性的影响较大，碱基质量设置为5-20是考虑了测序错误，过滤由于测序错误而造成的假阳性结果。覆盖度设置为400-1000是为了考虑测序中可能引起的冗余，如果过高，可能是测序的重复引起。各个参数值在此优选范围内时，可实现更佳的检测效果。

优选地，使用pyrohmmvar软件call SNP的具体参数为：滑动窗口的长度(w)设为25；最小的变异质量值分数(t)设为35.0；最小碱基质量值(B)设为13；最小的reads长度(L)设为25；设置忽略indels I；设置最小比对打分(M)为0；设置最小的比对一致性(E)为0.75；设置用于计算SNP最小支持数reads的系数(e)为0.2；设置支持SNP的最大覆盖度(c)为700。也即，具体参数为：-w 25-t 35.0-B 13-L 25-I-M 0-E 0.75-e 0.2-c 700。这里我们使用GATK软件与pyrohmmvar软件进行过评估测试，结果显示GATK虽然在检测出SNP位点的数量上具有优势，但存在较高的假阳性。而pyrohmmvar能很好的基于存在假indel的数据进行重矫正，最终能得到较为准确的结果。

S1075基因融合检测,基因融合是指两个或两个以上基因的部分或全部的序列构成一个新的杂合基因的过程。由于支持SE数据(注：SE数据即为单链测序的数据)基因融合检测的软件不多，在一个实施例中，使用tophat系列软件进行检测。在tophat比对软件将测序的reads比对到参考基因组后，根据比对结果，使用tophat-fusion(http://tophat-fusion.sourceforge.net/tutorial.html)来进行融合基因检测。使用tophat-fusion软件进行融合基因检测时，参数的具体设置如下：设置线程数p为1～16，优选为10；设置跳过blast比对，即skip-blast；设置支持断点最小的reads数目，即num-fusion-reads， 为1～10，优选为2；设置支持断点的最小的成对reads数目，即num-fusion-pairs，为0。也即，具体参数为：-p 10--skip-blast--num-fusion-reads 2--num-fusion-pairs 0。该软件能在linux上很好的移植与兼容，并能很好的分析proton的下游数据。tophat-fusion的使用可参考：

Edgren,H.,et al.,Identification of fusion genes in breast cancer by paired-end RNA-sequencing.Genome biology,2011.12(1):p.R6.。

在实际具体操作中，上述的第一步分析和第二部分析无先后顺序之分，也可以同时进行。

图9是本发明一种基于Proton的转录组测序数据的比较分析系统的一个实施例的框图。如图9所示，该实施例的比较分析系统包括序列接收模块11，用于获取某一物种至少两个转录组的基于Proton的原始测序数据；序列过滤模块12，与序列接收模块11相连，用于过滤各组原始测序数据中的不合格数据，获得各个转录组的clean reads；第一序列比对模块13，与序列过滤模块12相连，用于将各个转录组的clean reads分别比对到所述物种的参考基因组；定量分析模块14，与第一序列比对模块13相连，用于利用RPKM标准化方法对各个转录组的clean reads进行转录本定量分析；差异表达筛选模块15，与定量分析模块14相连，用于对转录本定量分析后的各个转录组的clean reads进行两两差异比较，并根据设定的阀值筛选出显著差异表达的基因；功能分析模块16，与差异表达筛选模块15相连，用于对显著差异表达的基因进行功能分析；第二序列比对模块17，与序列过滤模块12相连，用于将各个转录组的clean reads分别比对到所述物种的参考基因组；可变剪接基因注释模块181，与第二序列比对模块17相连，用于得到各个基因的可信的剪接位点，对剪接位点信息进行计算，获得含有可变剪接形式的基因注释结果；基因结构优化分析模块182，与第二序列比对模块17相连，用于将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的overlap超过阀值，且与其对应位置的已知转录本覆盖度 大于阀值，则将该转录本作为最终结果，往其已知的注释结果的两端延伸；新转录本预测模块183，与第二序列比对模块17相连，用于将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与所述物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的overlap没有达到阀值，且满足一定的测序深度，则预测其为新转录本；SNP分析模块184，与第二序列比对模块17相连，用于根据比对结果进行SNP预测；基因融合检测模块185，与第二序列比对模块17相连，用于根据比对结果进行融合基因检测。

相对于其他的测序平台和数据比较分析系统，该实施例的系统可以提高转录组测序数据比较分析的准确度和可靠性，该系统的主要优点有以下几方面：

1、国内首创基于proton测序平台一整套的数据处理方法。从软件的选择与测试，到数据处理程序的编写，包过结题报告的产生，相比其它平台无明显缺点或者劣势。

2、基于proton平台的数据从实验到测序再到数据分析，这系统相比hiseq平台，整个项目周期可在减少30％左右的时间。

3、尤其对于医院以及一些药企，该系统的应用相比芯片或者hiseq平台，项目的质量无劣势，但价格却价格却要降低不少。

最后应说明的是：以上所述仅为本发明的优选实施例而已，并不用于限制本发明，尽管参照前述实施例对本发明进行了详细的说明，对于本领域的技术人员来说，其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换。凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (10)

1.一种基于Proton的转录组测序数据的比较分析方法，其特征在于，该方法包括：

通过Proton测序平台获取某一物种至少两个转录组的原始测序数据；

过滤所述各组原始测序数据中的不合格数据，获得所述各个转录组的待分析数据；

对所述各个转录组的待分析数据进行第一步分析和第二步分析，其中，所述第一步分析包括：

(1)将所述各个转录组的待分析数据分别比对到所述物种的参考基因组；

(2)利用RPKM标准化方法对所述各个转录组的待分析数据进行转录本定量分析；

(3)对转录本定量分析后的所述各个转录组的待分析数据进行两两差异比较，并根据设定的阀值筛选出显著差异表达的基因；

(4)对所述显著差异表达的基因进行功能分析；

所述第二步分析包括：

(1)将所述各个转录组的待分析数据分别比对到所述物种的参考基因组；

(2)可变剪接分析：得到各个基因的可信的剪接位点，对剪接位点信息进行计算，获得含有可变剪接形式的基因注释结果；和/或

基因结构优化分析：根据比对结果进行基因结构的优化分析；和/或

新转录本预测：根据比对结果进行新转录本的预测；和/或

SNP分析：根据比对结果进行SNP预测；和/或

基因融合检测：根据比对结果进行融合基因检测。

2.如权利要求1所述的方法，其特征在于：所述不合格数据包括序列长度低于30的reads；修剪掉测序接头后序列长度低于30的reads；修剪掉低质量的碱基序列后序列长度低于30的reads，其中所述低质量碱基序列是指从3’端开始，以15为窗口统计reads平均质量，若平均质量值小于9，则继续向5’端滑动，直到窗口平均质量大于9时该窗口3’端到该reads 3’端的所有碱基序列。

3.如权利要求1所述的方法，其特征在于：所述第一步分析中的比对是使用Tmap软件进行的，所述第一步分析中的所述两两差异比对是使用泊松分布模型进行的；

优选地，利用所述Tmap软件比对时，参数的具体设置如下：设置比对结果输出类型a为2，即为所有最好的位点；设置程序运行时的线程数为1～16；设置输出文件类型为压缩的bam格式；设置比对策略为map4，即为bwa快速比对。

4.如权利要求1所述的方法，其特征在于：所述第一步分析中的功能分析包括聚类分析、GO富集分析和/或Pathway分析。

5.如权利要求1所述的方法，其特征在于：所述第二步分析中的所有比对是使用tophat软件进行的；

优选地，利用所述tophat软件比对时，参数的具体设置如下：设置比对reads的错配数为0～5；设置reads开gap最大的缺口数为0～20；设置reads的最大的编辑距离为5～20；设置reads剪接比对时最大错配数为0～4；设置reads比对最大插入长度为1～10；设置reads比对最大缺失长度为1～10；设置片段的比对最大错配数为0～5；设置Bowtie2比对内设参数为1～3；设置线程数为1～16；

更优选地，设置比对reads的错配数为3～5；设置reads开gap最大的缺口数为6～15，设置reads的最大的编辑距离为10～20，设置reads剪接比对时最大错配数为2～4；设置reads比对最大插入长度为5～10；设置reads比对最大缺失长度为5～10；设置片段的比对最大错配数为2～5；设置Bowtie2比对内设参数为1～3；设置线程数为6～16。

6.如权利要求1所述的方法，其特征在于：所述可变剪接分析中的对剪接位点信息进行计算包括以下步骤：利用perl脚本对剪接位点信息进行计算。

7.如权利要求1所述的方法，其特征在于：所述基因结构优化分析的步骤为：将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与所述物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的重叠部分超过阀值，且与其对应位置的已知转录本覆盖度大于阀值，则将该转录本作为最终结果，往其已知的注释结果的两端延伸；

优选地，所述基因结构优化分析中过滤掉组装结果小于阀值的转录本中的阀值为180bp；如果某个转录本与其对应的基因上游或下游的重叠部分超过阀值中的阀值为200bp；且与其对应位置的已知转录本覆盖度大于阀值中的阀值为0.9；

优选地，所述基因结构优化分析中的reads组装是使用cufflink软件进行的；

更优选地，使用所述cufflink软件组装时，参数的具体设置如下：设置参数u，对多reads进行矫正；设置线程数p为1～16。

8.如权利要求1所述的方法，其特征在于：所述新转录本预测的步骤为：将比对上的reads进行组装，过滤掉组装结果小于阀值的转录本，且过滤掉与所述物种已知的转录本同源性较高的转录本，再将剩下的转录本与所述物种的参考基因组进行比对，如果某个转录本与其对应的基因上游或下游的重叠部分没有达到阀值，且满足一定的测序深度，则称为新转录本；

优选地，所述新转录本预测中过滤掉组装结果小于阀值的转录本中的阀值为180bp；如果某个转录本与其对应的基因上游或下游的重叠部分超过阀值中的阀值为200bp；

优选地，所述新转录本预测中reads组装是使用cufflink软件进行的；

更优选地，使用所述cufflink软件组装时，参数的具体设置如下：设置参数u，对多reads进行矫正；设置线程数p为1～16。

9.如权利要求1所述的方法，其特征在于：所述SNP分析中SNP检测是使用pyrohmmvar软件进行的；所述基因融合检测是使用tophat-fusion软件进行的；

优选地，使用所述pyrohmmvar软件进行SNP检测时，参数的具体设置如下：滑动窗口的长度设为3～30；最小的变异质量值分数设为20～50；最小碱基质量值设为5～20；最小的reads长度设为15～30；设置忽略indels I；设置最小比对打分为0；设置最小的比对一致性为0.50～0.90；设置用于计算SNP最小支持数reads的系数为0.1～0.3；设置支持SNP的最大覆盖度为400～1000；

优选地，使用所述tophat-fusion软件进行融合基因检测时，参数的具体设置如下：设置线程数为1～16；设置跳过blast比对，即skip-blast；设置支持断点最小的reads数目为1～10；设置支持断点的最小的成对reads数目为0。

10.一种基于Proton的转录组测序数据的比较分析的系统，其特征在于，包括：

序列接收模块，用于获取某一物种至少两个转录组的基于Proton的原始测序数据；

序列过滤模块，与所述序列接收模块相连，用于过滤所述各组原始测序数据中的不合格数据，获得所述各个转录组的clean reads；

第一序列比对模块，与所述序列过滤模块相连，用于将所述各个转录组的clean reads分别比对到所述物种的参考基因组；

定量分析模块，与所述第一序列比对模块相连，用于利用RPKM标准化方法对所述各个转录组的clean reads进行转录本定量分析；

差异表达筛选模块，与所述定量分析模块相连，用于对转录本定量分析后的所述各个转录组的clean reads进行两两差异比较，并根据设定的阀值筛选出显著差异表达的基因；

功能分析模块，与所述差异表达筛选模块相连，用于对所述显著差异表达的基因进行功能分析；

第二序列比对模块，与所述序列过滤模块相连，用于将所述各个转录组的clean reads分别比对到所述物种的参考基因组；

可变剪接基因注释模块，与所述第二序列比对模块相连，用于得到各个基因的可信的剪接位点，对剪接位点信息进行计算，获得含有可变剪接形式的基因注释结果；和/或

基因结构优化分析模块，与所述第二序列比对模块相连，用于根据比对结果进行基因结构的优化分析；和/或

新转录本预测模块，与所述第二序列比对模块相连，用于根据比对结果进行新转录本的预测；和/或

SNP分析模块，与所述第二序列比对模块相连，用于根据比对结果进行SNP预测；和/或

基因融合检测模块，与所述第二序列比对模块相连，用于根据比对结果进行融合基因检测。