一种测定待测基因组区域表达水平的方法及系统的制作方法
【专利摘要】本发明提供了一种检测基因组区域表达水平(RPKM)的方法和系统，采用本发明，一方面，可以检测出整个基因的表达水平及其所有外显子各自的表达水平；另一个方面可以检测出同一个基因不同的同源异构体的表达水平及其所有外显子各自的表达水平；最后还可以检测出基因组任意指定区间的表达水平。
【专利说明】一种测定待测基因组区域表达水平的方法及系统
【技术领域】
[0001]本发明涉及生物技术和生物信息学领域，具体涉及一种测定基因组区域表达水平的方法及系统。
【背景技术】
[0002]生命遗传信息的表达调控既是生物学研究的重点领域，也是揭示生物学各种生命现象的重要手段，尤其是随着21世纪大量物种基因组序列的测定以及大量测序技术推陈出新，使得基因表达定量方面的研究突飞猛进。测序技术也从传统Sanger测序技术，迅速发展为多种第二代高通量测序技术，如罗氏454、IlluminaHiSeq和AB公司的SOLiD，以及第三代的单分子实时DNA测序技术。其中，Sanger测序技术和罗氏454测序技术的测序读长在700-1000bp, Illumina测序技术的测序读长平均IOObp左右,而单分子实时DNA测序技术的读长达到了 2500-3000bp。
[0003]第二代测序技术也被称为新一代测序技术(NGS, Next Generation Sequencing),目前主要是Illumina公司出的HiSeq为主，它通过从物种中提取出的RNA转录本中随机进行的短片段测序(通常平均读长50bp、75bp、100bp)获得所测样本的整体表达谱。转录本是通过以连续性基因组为模板进行转录，然后剪切去除内含子，拼接剩余的外显子而形成的。测序过程中，如果一个转录本的丰度高，则测序后定位基因组区域的测序读段也就多，可以通过对定位到基因上的外显子区的测序读段数来估计基因表达水平。测序读段数除了与基因真实表达水平成正比，还与基因长度成正比，同时也与测序深度即测序实验中得到的总读段数正相关。为了保持对不同基因和不同实验间估计的基因表达值的可比性，Mortazavi等人提出了 RPKM(Reads Per Kilo-base per Million reads)的概念，并成为 RNA-seq 应用早期估计基因表达水平和外显子表达水平的主要方法。RPKM是每百万读段中来自于某基因每千碱基长度的读段数，考虑了测序深度对读段计数的影响。
[0004]新一代测序技术的广泛普及，使得RNA测序(RNA-seq)已成为基因表达和转录组分析的重要手段。在NGS测序技术出现之前，不同基因表达水平测量的主要手段是基因芯片，利用在基因芯片上高密度集成特点的寡核苷酸，可以对不同组织或者不同发育阶段的特定基因表达差异和模式进行分析。但是与基因芯片数据相比，RNA-seq得到的是全基因组转录水平的数字化信号，具有高灵敏度、高分辨率、无饱和区等优势。
[0005]随着新一代测序技术的不断进步，产生的RNA-seq数据通量高、周期短和成本低，越来越多的人选择转录组测序作为科学研究的首选。RPKM在评估基因表达水平上的作用越来越显著，人们通过基因包含的外显子信息，和转录组测序数据在基因组上的定位信息，来计算出 RPKM值。FPKM(fragments per kilobase of exon per million fragments mapped)也可以用来表示基因表达水平。FPKM与RPKM计算方法基本一致。不同点就是FPKM计算的是片段(fragments),而RPKM计算的是测序读段(reads)。目前cuff links软件包中的cufflinks模块和cuffdiff模块及eXpress软件可以计算相关基因表达水平,具体计算过程为，首先统计出映射定位到基因组上的所有测序读段数目，然后统计出定位到各个基因外显子区间上的所有测序读段的数目，再计算出基因包含的外显子的长度，最后计算出基因的FPKM值。
[0006]但是，上述软件存在以下问题:
[0007](I)目前大部分计算RPKM的程序,仅支持TopHat、Bowtie、bwa等少数常用的序列比对定位程序，不能支持所有的Illumina/Solexa测序平台的读段定位程序；
[0008](2)在选择注释文件的时候，通常仅支持已知的基因注释文件，不能支持多种文件格式；
[0009](3)在计算基因表达水平的时候，通常计算的是片段的表达水平值，而不是整个基因的表达水平值；
[0010](4)在计算表达水平的时候，没有计算出单个外显子的表达水平；
[0011](5)在计算表达水平的时候，不能够计算出基因组任意指定区间的表达水平；
[0012](6)在计算表达水平的时候，通常仅支持计算一个转录组测序结果，不能够同时支持多个转录测序结果的基因表达水平的计算。
[0013]因此，本领域期 待一种能够检测基因表达水平和基因组任意指定区间表达水平的方法。
【发明内容】
[0014]本发明的目的是提供一种检测基因组区域表达水平(RPKM)的方法和系统。
[0015]本发明的第一方面提供了一种测定待测基因组区域表达水平的方法，包括以下步骤:
[0016](I)对待测样本进行测序，获得包含待测基因组区域转录本的转录组测序数据；
[0017](2)将获得的转录组测序数据与同一物种的基因组序列进行比对；
[0018](3)对定位到基因组的转录组测序读段进行筛选，所述筛选包括去除测序质量(99.9%的转录组测序读段；
[0019](4)将筛选后的转录组测序读段，按照其定位到基因组上的起始位置进行排序，并对排序结果建立索引；
[0020](5)根据待测基因组区域的位置信息，构建出计算RPKM的基因注释文件；
[0021](6)计算能够映射到基因组上的所有测序读段的总数M ;
[0022](7)根据上述步骤(5)构建的基因注释文件计算出定位至待测DNA区间上所有测序读段的总数R ；
[0023](8)根据上述步骤(5)构建的基因注释文件，计算出待测DNA区间所有被测序读段定位的序列长度L ;和
[0024](9)根据上述步骤(6)-(8)的计算结果，将步骤(7)得到的R除以步骤(6)得到的M与步骤(8)得到的L乘以109，得待测基因组区域的RPKM值，即为待测基因组区域的表达水平，计算公式如下，
D
[0025]RPKli=^IjX IO9
O
[0026]在另一优选例中，所述待测基因组区域包含N个同源异构体，且N≥2。如N可以为 2、3、4、5、6、7、8、9、10 或大于 10。
[0027]在另一优选例中，所述方法还包括结果验证步骤:提取待测样品的总RNA，经过反转录得到其cDNA，以cDNA作为模板进行PCR检测，验证待测基因组区域的表达水平。
[0028]在另一优选例中，所述步骤(5)中所述注释文件整合有已有的基因注释信息、新预测的基因注释信息和/或基因组任意指定区间的注释信息。
[0029]在另一优选例中，所述待测基因组区域表达水平，可以为单个基因的表达水平、同一个基因不同的同源异构体的表达水平、所有外显子的表达水平、单个外显子的表达水平以及基因组任意指定区间的表达水平。
[0030]在另一优选例中，当所述待测基因组区域中包含两个以上的同源异构体基因序列时，在测定过程中还包括步骤:将各同源异构体的所有外显子进行整合，对于重复的序列区间，仅保留单一序列，从而将同一待测基因组区域中的不同同源异构体的外显子整合成单一序列，将该单一序列的长度作为计算该基因组区域表达水平时的序列长度L。
[0031]在另一优选例中，所述步骤(I)中，所述转录组序列数据由罗氏454测序技术、Illumina测序技术、AB公司的SOLiD技术、或者第三代的单分子实时DNA测序技术获得。
[0032]在另一优选例中，所述步骤(2)中，序列比对程序为tophat2，以程序默认参数进行比对。
[0033]在另一优选例中，所述步骤(2)中，比对结果存储为SAM (Sequence Alignment/Map)格式或其二进制版本BAM格式的定位文件。
[0034]在另一优选例中，所述步骤⑷中，所述排序方法为:
[0035]a.按照每条测序读段定位到基因组的起始位置进行排序；
[0036]b.如果测序读段在基因组位置中的起始位置相同，按照其定位到基因组的先后顺序进行排序，并且保留所有的测序读段；
[0037]最后对排序结果建立索引。
[0038]在另一优选例中，所述步骤(5)中，基因注释文件存储格式为refFlat或者bed格式。
[0039]在另一优选例中，所述步骤(7)中，计算定位至待测DNA区间上所有测序读段的总数R时，如果一个转录组测序读段定位到两个外显子上，每个外显子都会对这个测序读段进行统计，以保证RPKM计算的准确性。
[0040]在另一优选例中，所述基因组区域选自如下的组:癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域或其他任意指定基因组区域。
[0041]本发明的第二方面提供了一种检测基因组区域表达水平的系统，所述系统包括:
[0042](I)比对单元，用于转录组测序读段与基因组序列进行比对；
[0043](2)筛选单元，用于对定位到基因组的转录组测序读段进行筛选；
[0044](3)排序单元，用于对转录组测序读段，按照其定位到基因组上的起始位置进行排序;
[0045](4)基因注释文件构建单元，用于构建和整合基因注释文件；和，
[0046](5)计算单元，包括:
[0047]a.第一模块，用于计算能够映射到基因组上的所有测序读段的总数M ;
[0048]b.第二模块，用于计算定位至待测DNA区间上所有测序读段的总数R ；[0049]c.第三模块，用于计算待测基因组表达区域的序列长度之和L ;和，
[0050]d.第四模块，用于计算待测基因组区域的RPKM值，计算公式为，
[0051]
【权利要求】
1.一种测定待测基因组区域表达水平的方法，其特征在于，包括以下步骤: (1)对待测样本进行测序，获得包含待测基因组区域转录本的转录组测序数据； (2)将获得的转录组测序数据与同一物种的基因组序列进行比对； (3)对定位到基因组的转录组测序读段进行筛选，所述筛选包括去除测序质量(99.9%的转录组测序读段； (4)将筛选后的转录组测序读段，按照其定位到基因组上的起始位置进行排序，并对排序结果建立索引； (5)根据待测基因组区域的位置信息，构建出计算RPKM的基因注释文件； (6)计算能够映射到基因组上的所有测序读段的总数M; (7)根据上述步骤(5)构建的基因注释文件计算出定位至待测DNA区间上所有测序读段的总数R ； (8)根据上述步骤(5)构建的基因注释文件，计算出待测DNA区间所有被测序读段定位的序列长度L ;和 (9)根据上述步骤(6)-⑶的计算结果，将步骤(7)得到的R除以步骤(6)得到的M与步骤(8)得到的L乘以109，得待测基因组区域的RPKM值，即为待测基因组区域的表达水平，计算公式如下，
2.如权利要求1所述的方法，其特征在于，所述方法还包括结果验证步骤，优选地，所述结果验证步骤包括:提取待测样品的总RNA，经过反转录得到其cDNA，以cDNA作为模板进行PCR检测，验证待测基因组区域的表达水平。
3.如权利要求1所述的方法，其特征在于，所述待测基因组区域包含N个同源异构体，且N≥2。
4.如权利要求3所述的方法，其特征在于，在测定过程中还包括步骤:将各同源异构体的所有外显子进行整合，对于重复的序列区间，仅保留单一序列，从而将同一待测基因组区域中的不同同源异构体的外显子整合成单一序列，将该单一序列的长度作为计算该基因组区域表达水平时的序列长度L。
5.如权利要求1所述的方法，其特征在于，所述步骤(1)中，所述转录组序列数据由罗氏454测序技术、Illumina测序技术、AB公司的SOLiD技术、或者第三代的单分子实时DNA测序技术获得。
6.如权利要求1所述的方法，其特征在于，所述步骤(4)中，所述排序方法为: a.按照每条测序读段定位到基因组的起始位置进行排序； b.如果测序读段在基因组位置中的起始位置相同，按照其定位到基因组的先后顺序进行排序，并且保留所 有的测序读段； 最后对排序结果建立索引。
7.如权利要求1所述的方法，其特征在于，所述基因组区域选自如下的组:癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域或任意指定的基因组区域。
8.一种测定待测基因组区域表达水平的系统，其特征在于，所述系统包括:(1)比对单元，用于转录组测序读段与基因组序列进行比对； (2)筛选单元，用于对定位到基因组的转录组测序读段进行筛选； (3)排序单元，用于对转录组测序读段，按照其定位到基因组上的起始位置进行排序； (4)基因注释文件构建单元，用于构建和整合基因注释文件；和， (5)计算单元,包括: a.第一模块，用于计算能够映射到基因组上的所有测序读段的总数M; b.第二模块，用于计算定位至待测DNA区间上所有测序读段的总数R； c.第三模块，用于计算待测基因组表达区域的序列长度之和L;和， d.第四模块，用于计算待测基因组区域的RPKM值，计算公式为，
9.如权利要求8所述的系统，其特征在于， 所述筛选单元中,所述筛选包括去除测序质量< 99.9%的转录组测序读段；和/或， 所述排序单元的排序方法为: a.按照每条测序读段定位到基因组的起始位置进行排序； b.如果测序读段在基因组位置中的起始位置相同，按照其定位到基因组的先后顺序进行排序，并且保留所有的测序读段； 最后对排序结果建立索引。
10.如权利要求8所述的系统，其特征在于，所述基因组区域选自如下的组:癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域或任意指定的基因组区域。
【文档编号】G06F19/22GK103984879SQ201410096063
【公开日】2014年8月13日 申请日期:2014年3月14日 优先权日:2014年3月14日
【发明者】杨力, 朱闪闪, 薛尉 申请人:中国科学院上海生命科学研究院