CN103984879B

CN103984879B - 一种测定待测基因组区域表达水平的方法及系统

Info

Publication number: CN103984879B
Application number: CN201410096063.1A
Authority: CN
Inventors: 杨力; 朱闪闪; 薛尉
Original assignee: Shanghai Institutes for Biological Sciences SIBS of CAS
Current assignee: Shanghai Institute of Nutrition and Health of CAS
Priority date: 2014-03-14
Filing date: 2014-03-14
Publication date: 2017-03-29
Anticipated expiration: 2034-03-14
Also published as: CN103984879A

Abstract

本发明提供了一种检测基因组区域表达水平(RPKM)的方法和系统，采用本发明，一方面，可以检测出整个基因的表达水平及其所有外显子各自的表达水平；另一个方面可以检测出同一个基因不同的同源异构体的表达水平及其所有外显子各自的表达水平；最后还可以检测出基因组任意指定区间的表达水平。

Description

一种测定待测基因组区域表达水平的方法及系统

技术领域

本发明涉及生物技术和生物信息学领域，具体涉及一种测定基因组区域表达水平的方法及系统。

背景技术

生命遗传信息的表达调控既是生物学研究的重点领域，也是揭示生物学各种生命现象的重要手段，尤其是随着21世纪大量物种基因组序列的测定以及大量测序技术推陈出新，使得基因表达定量方面的研究突飞猛进。测序技术也从传统Sanger测序技术，迅速发展为多种第二代高通量测序技术，如罗氏454、IlluminaHiSeq和AB公司的SOLiD，以及第三代的单分子实时DNA测序技术。其中，Sanger测序技术和罗氏454测序技术的测序读长在700-1000bp，Illumina测序技术的测序读长平均100bp左右，而单分子实时DNA测序技术的读长达到了2500-3000bp。

第二代测序技术也被称为新一代测序技术(NGS,Next Generation Sequencing)，目前主要是Illumina公司出的HiSeq为主，它通过从物种中提取出的RNA转录本中随机进行的短片段测序(通常平均读长50bp、75bp、100bp)获得所测样本的整体表达谱。转录本是通过以连续性基因组为模板进行转录，然后剪切去除内含子，拼接剩余的外显子而形成的。测序过程中，如果一个转录本的丰度高，则测序后定位基因组区域的测序读段也就多，可以通过对定位到基因上的外显子区的测序读段数来估计基因表达水平。测序读段数除了与基因真实表达水平成正比，还与基因长度成正比，同时也与测序深度即测序实验中得到的总读段数正相关。为了保持对不同基因和不同实验间估计的基因表达值的可比性，Mortazavi等人提出了RPKM(Reads Per Kilo-base per Million reads)的概念，并成为RNA-seq应用早期估计基因表达水平和外显子表达水平的主要方法。RPKM是每百万读段中来自于某基因每千碱基长度的读段数，考虑了测序深度对读段计数的影响。

新一代测序技术的广泛普及，使得RNA测序(RNA-seq)已成为基因表达和转录组分析的重要手段。在NGS测序技术出现之前，不同基因表达水平测量的主要手段是基因芯片，利用在基因芯片上高密度集成特点的寡核苷酸，可以对不同组织或者不同发育阶段的特定基因表达差异和模式进行分析。但是与基因芯片数据相比，RNA-seq得到的是全基因组转录水平的数字化信号，具有高灵敏度、高分辨率、无饱和区等优势。

随着新一代测序技术的不断进步，产生的RNA-seq数据通量高、周期短和成本低，越来越多的人选择转录组测序作为科学研究的首选。RPKM在评估基因表达水平上的作用越来越显著，人们通过基因包含的外显子信息，和转录组测序数据在基因组上的定位信息，来计算出RPKM值。FPKM(fragments per kilobase of exon per million fragmentsmapped)也可以用来表示基因表达水平。FPKM与RPKM计算方法基本一致。不同点就是FPKM计算的是片段(fragments)，而RPKM计算的是测序读段(reads)。目前cufflinks软件包中的cufflinks模块和cuffdiff模块及eXpress软件可以计算相关基因表达水平，具体计算过程为，首先统计出映射定位到基因组上的所有测序读段数目，然后统计出定位到各个基因外显子区间上的所有测序读段的数目，再计算出基因包含的外显子的长度，最后计算出基因的FPKM值。

但是，上述软件存在以下问题：

(1)目前大部分计算RPKM的程序，仅支持TopHat、Bowtie、bwa等少数常用的序列比对定位程序，不能支持所有的Illumina/Solexa测序平台的读段定位程序；

(2)在选择注释文件的时候，通常仅支持已知的基因注释文件，不能支持多种文件格式；

(3)在计算基因表达水平的时候，通常计算的是片段的表达水平值，而不是整个基因的表达水平值；

(4)在计算表达水平的时候，没有计算出单个外显子的表达水平；

(5)在计算表达水平的时候，不能够计算出基因组任意指定区间的表达水平；

(6)在计算表达水平的时候，通常仅支持计算一个转录组测序结果，不能够同时支持多个转录测序结果的基因表达水平的计算。

因此，本领域期待一种能够检测基因表达水平和基因组任意指定区间表达水平的方法。

发明内容

本发明的目的是提供一种检测基因组区域表达水平(RPKM)的方法和系统。

本发明的第一方面提供了一种测定待测基因组区域表达水平的方法，包括以下步骤：

(1)对待测样本进行测序，获得包含待测基因组区域转录本的转录组测序数据；

(2)将获得的转录组测序数据与同一物种的基因组序列进行比对；

(3)对定位到基因组的转录组测序读段进行筛选，所述筛选包括去除测序质量≤99.9%的转录组测序读段；

(4)将筛选后的转录组测序读段，按照其定位到基因组上的起始位置进行排序，并对排序结果建立索引；

(5)根据待测基因组区域的位置信息，构建出计算RPKM的基因注释文件；

(6)计算能够映射到基因组上的所有测序读段的总数M；

(7)根据上述步骤(5)构建的基因注释文件计算出定位至待测DNA区间上所有测序读段的总数R；

(8)根据上述步骤(5)构建的基因注释文件，计算出待测DNA区间所有被测序读段定位的序列长度L；和

(9)根据上述步骤(6)-(8)的计算结果，将步骤(7)得到的R除以步骤(6)得到的M与步骤(8)得到的L乘以10⁹，得待测基因组区域的RPKM值，即为待测基因组区域的表达水平，计算公式如下，

在另一优选例中，所述待测基因组区域包含N个同源异构体，且N≥2。如N可以为2、3、4、5、6、7、8、9、10或大于10。

在另一优选例中，所述方法还包括结果验证步骤：提取待测样品的总RNA，经过反转录得到其cDNA，以cDNA作为模板进行PCR检测，验证待测基因组区域的表达水平。

在另一优选例中，所述步骤(5)中所述注释文件整合有已有的基因注释信息、新预测的基因注释信息和/或基因组任意指定区间的注释信息。

在另一优选例中，所述待测基因组区域表达水平，可以为单个基因的表达水平、同一个基因不同的同源异构体的表达水平、所有外显子的表达水平、单个外显子的表达水平以及基因组任意指定区间的表达水平。

在另一优选例中，当所述待测基因组区域中包含两个以上的同源异构体基因序列时，在测定过程中还包括步骤：将各同源异构体的所有外显子进行整合，对于重复的序列区间，仅保留单一序列，从而将同一待测基因组区域中的不同同源异构体的外显子整合成单一序列，将该单一序列的长度作为计算该基因组区域表达水平时的序列长度L。

在另一优选例中，所述步骤(1)中，所述转录组序列数据由罗氏454测序技术、Illumina测序技术、AB公司的SOLiD技术、或者第三代的单分子实时DNA测序技术获得。

在另一优选例中，所述步骤(2)中，序列比对程序为tophat2，以程序默认参数进行比对。

在另一优选例中，所述步骤(2)中，比对结果存储为SAM(Sequence Alignment/Map)格式或其二进制版本BAM格式的定位文件。

在另一优选例中，所述步骤(4)中，所述排序方法为：

a.按照每条测序读段定位到基因组的起始位置进行排序；

b.如果测序读段在基因组位置中的起始位置相同，按照其定位到基因组的先后顺序进行排序，并且保留所有的测序读段；

最后对排序结果建立索引。

在另一优选例中，所述步骤(5)中，基因注释文件存储格式为refFlat或者bed格式。

在另一优选例中，所述步骤(7)中，计算定位至待测DNA区间上所有测序读段的总数R时，如果一个转录组测序读段定位到两个外显子上，每个外显子都会对这个测序读段进行统计，以保证RPKM计算的准确性。

在另一优选例中，所述基因组区域选自如下的组：癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域或其他任意指定基因组区域。

本发明的第二方面提供了一种检测基因组区域表达水平的系统，所述系统包括：

(1)比对单元，用于转录组测序读段与基因组序列进行比对；

(2)筛选单元，用于对定位到基因组的转录组测序读段进行筛选；

(3)排序单元，用于对转录组测序读段，按照其定位到基因组上的起始位置进行排序；

(4)基因注释文件构建单元，用于构建和整合基因注释文件；和，

(5)计算单元，包括：

a.第一模块，用于计算能够映射到基因组上的所有测序读段的总数M；

b.第二模块，用于计算定位至待测DNA区间上所有测序读段的总数R；

c.第三模块，用于计算待测基因组表达区域的序列长度之和L；和，

d.第四模块，用于计算待测基因组区域的RPKM值，计算公式为，

在另一优选例中，所述比对单元中，序列比对程序为tophat2，以程序默认参数进行比对，比对结果存储为SAM(Sequence Alignment/Map)格式或其二进制版本BAM格式的定位文件。

在另一优选例中，所述筛选单元中，所述筛选包括去除测序质量≤99.9%的转录组测序读段。

在另一优选例中，所述排序单元的排序方法为：

a.按照每条测序读段定位到基因组的起始位置进行排序；

最后对排序结果建立索引。

在另一优选例中，所述基因注释文件构建单元中，所述基因注释文件存储格式为refFlat或者bed格式。

在本发明的利用转录组测序结果计算基因表达水平(RPKM)的方法中，所述转录组序列读段可以由罗氏454测序技术、Illumina测序技术和AB公司的SOLiD技术获得的双末端(pair-end)测序读段和单末端(single-end)测序读段；并且可以采用TopHat、TopHat2、Bowtie、Bowtie2、bwa(Burrows-Wheeler Aligner)、SOAP2、SOAP3等多种序列比对定位程序，以程序的最优参数进行比对；在构建基因注释文件时，我们下载已知物种的RefSeq注释文件、KnownGene注释文件，并且也可以通过de novo拼接转录组测序结果构建预测出的新的基因注释文件，因此我们的方法不仅可以对已知的基因、同源异构体、外显子的表达水平进行计算，还可以对预测新的基因、同源异构体、外显子的表达水平；同时也可以计算中用户给定的任意基因组区间表达值；另外，如果同时给定多个的转录组测序结果，我们可以对多个转录组测序结果同时计算，每个转录组测序结果会计算出相应的基因、同源异构体和外显子的RPKM值或者基因组任意指定区间的RPKM值。

应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。

附图说明

图1显示计算基因表达量和外显子表达量示意图。

图2为实施例1中基因GREM1的表达量示意图。

具体实施方式

本发明人通过广泛而深入的研究，获得一种检测基因组区域表达水平(RPKM)的方法和系统，采用本发明，一方面，可以检测出整个基因的表达水平及其所有外显子各自的表达水平；另一个方面可以检测出同一个基因不同的同源异构体的表达水平及其所有外显子各自的表达水平；最后还可以检测出基因组任意指定区间的表达水平。

在进一步描述本发明具体实施方式之前，应理解，本发明的保护范围不局限于下述特定的具体实施方案；还应当理解，本发明实施例中使用的术语是为了描述特定的具体实施方案，而不是为了限制本发明的保护范围。为了对本发明作出清楚的说明，首先针对本说明书中使用的技术术语如下进行定义。

转录组测序(RNA-seq)数据：研究特定组织或细胞在某一功能状态下所能转录出来的RNA的总和，主要包括mRNA和非编码RNA。转录组研究是基因功能及结构研究的基础，通过新一代高通量测序，能够全面快速地获得某一物种特定组织或器官在某一状态下的几乎所有转录本序列信息，在本说明书中主要指通过NGS测到的特定个体的转录组数据。

测序读段(Sequence Reads)，测序技术所产生的单个测序片段，在本说明书中是指转录组测序中的测序片段。

DNA片段(DNA Flagments)，在本说明书中是指最终用于测序的DNA片段。在单末端测序中，DNA片段等价于测序读段；但是在双末端测序中，一个DNA片段，可以得到两个测序读段，而这两个测序读段可能只有一个或全部定位到基因组上，进而形成FPKM和RPKM两种计算方法的差别。

基因表达水平(Gene Expression Level)，也叫基因表达量，是细胞在生命过程中，把储存在DNA中遗传信息经过转录和翻译，转变成具有生物活性的蛋白质分子过程中表达水平的高低丰度，本说明书中是指基因、外显子、同源异构体的转录表达水平。

基因名称(Gene Symbol)，也被称为参考基因，是基因组浏览器中通用基因名称，本说明书中是指用来计算基因表达时的基因名称。以下简称基因。

外显子(Exon)，是真核生物基因中能转录，且在剪接后会被保存下来相应DNA区域。所有的外显子一同组成了遗传信息，其中编码蛋白的信息会体现在蛋白质上，在本说明书中是指用来计算外显子表达水平时的外显子名称。

同源异构体(Isoform)，来自一个基因的mRNA前体因选择性剪接而产生多种mRNA，并翻译出的不同蛋白质，在本说明书中是指由于选择性剪切产生的来源于同一个基因不同的外显子组成的多个mRNA。

基因组任意指定区间(Target Genome Region)，在本说明书中是指用户给定的特定基因组位置信息，包含染色体名称、基因组起始位置和基因组终止位置。

基因组任意指定区间长度(Target Genome Region Length)，在本说明书中是指根据用户给定的任意基因组区间，将基因组起始位置减去基因组终止位置的长度。

注释基因外显子起始位置(Exon Start Position)，在本说明书中是指外显子在基因组中的起始位置。以下简称为起始位置。

注释基因外显子终止位置(Exon End Position)，在本说明书中是指外显子在基因组中的终止位置。以下简称为终止位置。

外显子长度(Exon Length)，在本说明书中是指每个外显子在基因组中的起始位置和终止位置之差，用来表示外显子在基因组中长度。

基因长度(Gene Length)，在本说明书中是指基因在基因组中包含的所有外显子的起始位置和终止位置之差，用来表示基因在基因组中所有外显子的长度之和。

同源异构体长度(Isoform Length)，在本说明书中是指同源异构体在基因组中包含的所有外显子的起始位置和终止位置之差，用来表示同源异构体在基因组中所有外显子的长度之和。

基因组匹配读段数(Total Genome Mapped Reads)，在本说明书中是指比对到基因组区域上的所有转录组测序读段总数。

外显子匹配读段数(Total Exon Mapped Reads)，在本说明书中是指比对到外显子区域上的转录组测序读段的数目。

基因组任意指定区间匹配读段数(Total Genome Region Mapped Reads)，在本说明书中是指比对到基因组任意指定区间的转录组测序读段数目。

每百万读段中来自于某基因每千碱基长度的读段数(RPKM,Reads Per Kilobases per Million reads)，将映射(mapping)定位到基因的读段数除以映射(mapping)定位到基因组的所有读段数(以million为单位)与基因外显子的长度(以KB为单位)。本发明中“每百万读段中来自于某基因每千碱基长度的读段数(RPKM)”是指将映射定位到基因、同源异构体或外显子上的外显子匹配读段数除以基因组匹配读段数与基因、同源异构体或外显子的长度，用来表示基因、同源异构体或外显子表达水平的值，以下简称RPKM。

基因注释文件，在本说明书中是指计算基因表达时输入的基因、同源异构体、外显子以及特定基因组区域的注释文件格式，可以是refFlat格式或者bed格式。

本发明的主要优点在于：

(1)根据步骤(5)-(9)，在计算基因表达量的时候，不仅能够计算单个同源异构体的表达量，而且能够更加准确地计算出包含多个同源异构体的整个基因的表达量以及所有外显子的表达量。每个外显子表达量的准确定量也将有益于基因中不同同源异构体的差异分析。对于同一个基因有多个同源异构体时，如图1所示，假设基因I有三个同源异构体(I、II、III)，各个同源异构体的外显子有差异，我们将该基因的所有同源异构体包含的外显子先整合为一个最完整的转录本，外显子1和外显子2整合为外显子A，外显子3整合为外显子B，外显子4整合为外显子C、外显子5-7整合为外显子D(各外显子的重复序列部分仅保留单一序列)，并计算整合的转录本的整体表达量(作为基因的表达量)；同时对于所有的外显子，可以分别计算出各自的表达量；

(2)根据步骤(5)-(9)，在计算基因表达量的时候，还能够批量计算出用户给定的基因组任意指定区间的表达量，这样可以根据用户需要，进行特定基因组区间的分析；

(3)根据步骤(1)，本发明能够对多种序列定位程序的结果进行计算，包括TopHat、bowtie、bwa和SOAPaligner/soap2等，用户在选择序列比对程序的时候就会有多种选择；

(4)根据步骤(5)，在选择基因注释文件时，我们能够支持多种基因注释文件，包括KnownGene、RefSeq以及de novo拼接转录组测序结果构建预测出的新的基因注释文件；

(5)本发明能够同时计算多个转录组测序结果的基因表达量，不局限于单一转录组测序结果的基因表达量的计算。

实施例1

利用人的胚胎干细胞(H9)的转录组测序结果来计算人的已知参考基因(refseq)的基因、同源异构体和外显子的表达水平。

材料：从美国生物技术信息中心(NCBI，National Center for BiotechonlogyInformation)网站(http://www.ncbi.nlm.nih.gov/sra/)下载人的转录组Illumina测序数据(登录号:SRX243742)，转录组测序数据共有32633419条测序读段，平均读长100bp。从美国加州大学圣克鲁斯分校网站(http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/)下载人的基因组序列(版本号:hg19)，下载(http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/)人的已知参考基因的注释文件(refFlat.txt,版本号:2013-11-24)，已知注释参考基因24082条。

步骤01：从马里兰大学(The University of Maryland)生物信息学和计算生物学中心网站(http://tophat.cbcb.umd.edu/downloads/tophat-2.0.9.Linux_x86_64.tar.gz)下载序列比对定位程序tophat(版本号:2.0.9)。从开源软件开发网站(sourceforge)下载bowtie2(http://netcologne.dl.sourceforge.net/project/bowtie-bio/bowtie2/2.1.0/bowtie2-2.1.0-linux-x86_64.zip，版本号2.1.0)和samtools(http://nchc.dl.sourceforge.net/project/samtools/samtools/0.1.18/samtools-0.1.18.tar.bz2，版本号:0.1.18)。

步骤02：将32633419条人的转录组测序读段与人的基因组进行比对，结果显示有25167366条转录组测序读段能够定位到人的基因组上，占所有测序读段的77.1%，其中1206344条测序读段定位到基因组的多个位置上，占所有测序读段的3.7%。

步骤03：对于仅定位到基因组的一个位置上的测序读段，按照每条测序读段定位到基因组的位置进行排序，排序过程中如下：

(1)按照每条测序读段定位到基因组的起始位置进行排序；

(2)如果测序读段在基因组位置中的起始位置相同，我们按照其定位到基因组的先后顺序进行排序，并且保留所有的测序读段；

最后对排序结果建立索引。

步骤04：将下载的已知注释参考基因文件refFlat.txt作为计算基因表达水平的基因注释文件，根据refFlat.txt文件构建出计算基因同源异构体表达水平的文件refFlat_Isoform.txt，具体方法是在refFlat.txt文件中每个基因按照其同源异构体的顺序分别在基因名称后面标记上同源异构体_1、同源异构体_2、…、同源异构体_N，用来区分同一个基因不同的同源异构体，构建后的基因注释文件注释参考基因24073条，对应的同源异构体47369条。

步骤05：通过统计定位到基因组上的所有转录组测序读段，来确定计算RPKM时的基因组匹配读段数，这一数值在后续计算过程是恒定的，具体统计方法是，先计算定位到每条染色体上的所有转录组测序读段数，然后求所有染色体上的转录组测序读段数之和，最后得到25167366条转录组测序读段作为基因组匹配读段数M(Total Genome MappedReads)。

步骤06：根据上述步骤04构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该位置上所有转录组测序读段的总数，作为计算RPKM时的外显子匹配读段数R(Total exon Mapped Reads)。如果一个转录组测序读段定位到两个外显子上，每个外显子都会对这个测序读段进行统计，以保证RPKM计算的准确性。

步骤07：根据上述步骤04构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该基因所有外显子的长度，作为计算RPKM时的外显子长度，具体计算方法是，对同一个基因，先统计基因包含的外显子数，然后再将每个外显子的终止位置减去起始位置加上1，就得到了每个外显子长度，然后对所有的外显子长度求和，得到的就是基因外显子长度，用来作为计算RPKM时的外显子长度L(exon length)。

步骤08：根据上述步骤05-07的计算结果，将得到的基因的全部外显子匹配读段数除以基因组匹配读段数与基因的全部外显子长度乘以10⁹，这样就得到了基因总的RPKM值，即基因的整个表达水平。

RPKM计算公式：

步骤09：根据上述步骤05-07的计算结果，对于基因的不同同源异构体的表达水平，将得到的同源异构体上的全部外显子匹配读段数除以基因组匹配读段数与同源异构体的全部外显子长度乘以10⁹，这样就得到了同源异构体总的RPKM值，即同源异构体的表达水平。

步骤10：根据上述步骤05-07的计算结果，对于基因的单个外显子的表达水平，将得到的单个外显子的外显子匹配读段数除以基因组匹配读段数与单个外显子的外显子长度乘以10⁹，这样就得到了单个外显子的RPKM值，即外显子的表达水平。

结果：注释参考基因24073个，有12424个基因计算出基因表达水平，对应的外显子一共152317个，其中RPKM大于1的共有12424个基因和129657外显子。同源异构体47369条，有26298个同源异构体计算出基因表达水平，对应的外显子一共296281个，其中RPKM同时大于1的26298个同源异构体和261945外显子。

以基因GREM1为例(图2)，GREM1有NM_001191322(外显子3个，分别为外显子1、2、4)、NM_001191323(外显子3个，分别为外显子1、3、4)和NM_013372(外显子2个，分别为外显子1、5)等3个同源异构体。本发明在计算GREM1基因表达量的时候，对3个同源异构体的所有外显子进行整合，如图2所示，整合后组成的转录本由2个外显子组成，两个外显子长度分别为158bp和3980bp，整合的转录本长度为4138bp，以整合的转录本来计算GREM1基因的表达量，整合后的GREM1基因表达量为2.371。从图2上面部分所示测序读段的分布图也可以得出，GREM1基因的5个外显子中有且仅有共同存在的一个外显子(chr15:33023128-33026870,长度3742bp)有表达量，由于这个外显子都包含在这3个同源异构体中，因此3个同源异构体的表达量差异变化应该不大。通过步骤9计算出同源异构体NM_001191322、NM_001191323和NM_013372的表达量分别为2.498、2.444和2.371；而通过现有cufflinks软件包，3个同源异构体的表达量分别为2.64266、0.00019和3.93723e^-10，基因表达量为2.64285，各异构体的表达量之间存在着较大的差异，与预期结果不符。由于一个基因对应多个转录本的现象大量存在于现有人类基因组中，因此相比现有的cufflinks软件包，本发明能够更加准确计算出基因的表达量，并且准确的计算出所有外显子的表达量，从而有益于各种基因的同源异构体的差异分析。

为了验证基因GREM1的在mRNA水平上存在NM_001191322、NM_001191323和NM_013372这3种同源异构体，并且彼此表达差异变化不大。我们首先提取了人源胚胎干细胞H9的总RNA，经过反转录得到其cDNA。然后根据图2所示GREM1的同源异构体结果，设计两对引物：

引物对1，以跨外显子1和外显子3的序列区域设计上游引物，在外显子4上设计下游引物；

引物对2，在外显子1上设计上游引物，在外显子5的chr15:33023006-33023127区域上设计下游引物。

具体的引物序列如表1所示，由于仅同源异构体NM_013372上含有chr15:33023006-33023127区域，基于以上设计，引物对1(引物1F和引物1R)用于扩增同源异构体NM_001191322和NM_001191323，引物对2(引物2F和引物2R)用于扩增同源异构体NM_013372。

以人源胚胎干细胞H9的cDNA作为模板，以表1所示引物序列进行扩增，结果显示，引物对1成功扩增出了预期的同源异构体NM_001191322和NM_001191323，引物对2成功扩增出了同源异构体NM_013372，电泳结果显示，扩增产物的量没有明显差异。实验结果与本文计算的基因GREM1的3种同源异构体的表达量差异不大的结果相符，而cufflinks的结果显示NM_001191322表达量为2.64266、NM_001191323和NM_013372分别为0.00019和3.93723e^-10，这与事实存在较大误差。

上述结果提示，相比现有的cufflinks软件包，本发明的计算结果更加准确。

表1扩增引物

实施例2

利用人的海拉细胞(Hela)转录组测序结果来计算预测出的新的基因、同源异构体和外显子的表达水平。

材料：从美国生物技术信息中心(NCBI，National Center for BiotechonlogyInformation)网站(http://www.ncbi.nlm.nih.gov/sra/)下载人的转录组Illumina测序读段(登录号:ERX103445)，转录组测序数据共有58076910条测序读段，平均读长72bp。

步骤01：从马里兰大学(The University of Maryland)生物信息学和计算生物学中心网站(http://cufflinks.cbcb.umd.edu/downloads/cufflinks-2.1.1.Linux_x86_64.tar.gz)下载序列拼接程序cufflinks(版本号:2.1.1)。

步骤02：将58076910条人的转录组测序读段进行过滤，去除低质量的数据之后，保留了47569394条转录组测序读段，将这些测序读段与人的基因组进行比对，结果显示有43750301条转录组测序读段能够定位到人的基因组上，占清洗后的所有测序读段的92%，其中2548627条测序读段定位到基因组的多个位置上，占清洗后的所有测序读段的5.4%，有41201674条测序读段仅定位到基因组的一个位置上，占所有测序读段的84%。

(1)按照每条测序读段定位到基因组的起始位置进行排序；

最后对排序结果建立索引。

步骤04：根据步骤02的比对结果，用De novo拼接程序cufflinks进行拼接，结果会得到新的基因注释文件transcript.gtf；结合已知的基因注释文件，对transcript.gtf重新融合为一个新的基因注释文件merge.gtf，该文件包括已知基因的预测基因的注释信息，为了利用转录组测序结果计算预测出的新的基因、同源异构体和外显子的表达水平，我们将已知的基因注释信息删除，仅保留预测出的新的基因注释文件，用于后续的基因注释文件的构建。

步骤05：根据步骤04的结果对预测出的新的基因文件构建出新的基因注释文件CUFF.txt作为计算基因表达水平的基因注释文件，根据CUFF.txt文件构建出计算基因同源异构体表达水平的文件CUFF_Isoform.txt，具体方法是在CUFF.txt文件中每个基因按照其同源异构体的顺序分别在基因名称后面标记上同源异构体_1、同源异构体_2、…、同源异构体_N，用来区分同一个基因不同的同源异构体，构建后的基因注释文件注释参考基因1031条，对应的同源异构体998条。

步骤06：通过统计定位到基因组上的所有转录组测序读段，来确定计算RPKM时的基因组匹配读段数，这一数值在后续计算过程是恒定的，具体统计方法是，先计算定位到每条染色体上的所有转录组测序读段数，然后求所有染色体上的转录组测序读段数之和，最后得到41201674条转录组测序读段作为基因组匹配读段数M(Total Genome MappedReads)。

步骤07：根据上述步骤05构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该位置上所有转录组测序读段的总数，作为计算RPKM时的外显子匹配读段数R(Total exon Mapped Reads)。如果一个转录组测序读段定位到两个外显子上，每个外显子都会对这个测序读段进行统计，以保证RPKM计算的准确性。

步骤08：根据上述步骤05构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该基因所有外显子的长度，作为计算RPKM时的外显子长度，具体计算方法是，对同一个基因，先统计基因包含的外显子数，然后再将每个外显子的终止位置减去起始位置加上1，就得到了每个外显子长度，然后对所有的外显子长度求和，得到的就是基因长度，用来作为计算RPKM时的外显子长度L(exon length)。

步骤09：根据上述步骤06-08的计算结果，将得到的基因的全部外显子匹配读段数除以基因组匹配读段数与基因的全部外显子长度乘以10⁹，这样就得到了基因总的RPKM值，即基因的整个表达水平。

步骤10：根据上述步骤06-08的计算结果，对于基因的不同同源异构体的表达水平，将得到的同源异构体上的全部外显子匹配读段数除以基因组匹配读段数与同源异构体的全部外显子长度乘以10⁹，这样就得到了同源异构体总的RPKM值，即同源异构体的表达水平。

步骤11：根据上述步骤06-08的计算结果，对于基因的单个外显子的表达水平，将得到的单个外显子的外显子匹配读段数除以基因组匹配读段数与单个外显子的外显子长度乘以10⁹，这样就得到了单个外显子的RPKM值，即外显子的表达水平。

结果：注释参考基因1031个，有600个基因计算出基因表达水平，对应的外显子一共887个，其中RPKM大于1的共有600个基因和820外显子。同源异构体998条，有630个同源异构体计算出基因表达水平，对应的外显子一共963个，其中RPKM同时大于1的630个同源异构体和907外显子。本发明的方法可以计算预测出新基因的表达量，这些新基因的表达量，可以用来修正已知基因注释文件中存在的基因注释信息的错误，另一方面可以用来注释新的基因，用于后续的新基因实验验证。

实施例3

利用多组人的类淋巴母细胞(lymphoblastoid)转录组测序结果来计算人的已知参考基因(knownGene)的基因、同源异构体和外显子的表达水平，以及基因组任意指定区间的表达水平。

材料：从欧洲生物信息研究所(EMBL-European Bioinformatics Institute)网站(http://www.ebi.ac.uk/ena/)下载三组人的转录组Illumina测序读段(登录号:ERR188021,ERR188022ERR188023)。从美国加州大学圣克鲁斯分校(http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/)下载人的已知参考基因的注释文件(knownGene.txt,版本号:2013-6-30)，已知注释参考基因28252条。

步骤01：将转录组测序数据ERR188021、ERR188022和ERR188023(测序读段数分别为65015656、44647732和73003930)与人的基因组进行比对，结果显示分别有58128352、41244172和67068647条转录组测序读段能够定位到人的基因组上，占所有测序读段的89.41%、92.38%和91.87%，其中分别有53959738、38274777和61221960条测序读段定位到基因组的一个位置上，占所有测序读段的82.99%、85.73%和83.86%。

步骤02：对于仅定位到基因组的一个位置上的测序读段，按照每条测序读段定位到基因组的位置进行排序，排序过程中如下：

(1)按照每条测序读段定位到基因组的起始位置进行排序；

最后对排序结果建立索引。

步骤03：将下载的已知注释参考基因文件knownGene.txt作为计算基因表达水平的基因注释文件，根据knownGene.txt文件构建出计算基因同源异构体表达水平的文件knownGene_Isoform.txt，具体方法是在knownGene.txt文件中每个基因按照其同源异构体的顺序分别在基因名称后面标记上同源异构体_1、同源异构体_2、…、同源异构体_N，用来区分同一个基因不同的同源异构体，构建后的基因注释文件注释参考基因28252条，对应的同源异构体80343条。

步骤04：根据下载的已知注释参考基因文件knownGene.txt模拟地构建出计算基因组任意指定区间表达水平的文件knownGene_bed.txt，具体方法是选取knownGene.txt文件的第2列染色体名称、第4列基因组起始位置和第5列基因组终止位置，构建出基因组任意指定区间注释文件，基因组任意指定区间共有82960个区间。

步骤05：通过统计定位到基因组上的所有转录组测序读段，来确定计算RPKM时的基因组匹配读段数，这一数值在后续计算过程是恒定的，具体统计方法是，先计算定位到每条染色体上的所有转录组测序读段数，然后求所有染色体上的转录组测序读段数之和，最后分别得到53959738(ERR188021)、38274777(ERR188022)和61221960(ERR188022)条转录组测序读段作为基因组匹配读段数M(Total Genome Mapped Reads)。

步骤06：根据上述步骤03构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该位置上所有转录组测序读段的总数，作为计算RPKM时的外显子匹配读段数R(Total exon Mapped Reads)，如果一个转录组测序读段定位到两个外显子上，每个外显子都会对这个测序读段进行统计，以保证RPKM计算的准确性。

步骤07：根据上述步骤03构建的基因注释文件提供的每一行基因在基因组的位置信息，以及该基因包含的所有外显子在基因组中的位置信息，计算出该基因所有外显子的长度，作为计算RPKM时的外显子长度，具体计算方法是，对同一个基因，先统计基因包含的外显子数，然后再将每个外显子的终止位置减去起始位置加上1，就得到了每个外显子长度，然后对所有的外显子长度求和，得到的就是基因长度，用来作为计算RPKM时的外显子长度L(exon length)。

步骤11：根据上述步骤04构建的基因注释文件提供的基因组任意指定区间，计算出该区间上所有转录组测序读段的总数，作为计算RPKM时的基因组任意指定区间匹配读段数(Total Genome Region Mapped Reads)，如果一个转录组测序读段定位到两个不同的基因组区间上，每个基因组区间都会对这个测序读段进行统计，以保证RPKM计算的准确性。

步骤12：根据上述步骤04构建的基因注释文件提供的基因组任意指定区间，计算出该区间的基因组任意指定区间长度，作为计算RPKM时的基因组任意指定区间长度，具体计算方法是，将每个区间的基因组终止位置减去基因组起始位置加上1，就得到了计算RPKM时每个基因组任意指定区间长度(Target Genome Region Length)。

步骤13：根据上述步骤05和步骤11-12，将得到的基因组任意指定区间的转录组测序读段数除以基因组匹配读段数与基因组任意指定区间长度乘以10⁹，这样就得到了基因组任意指定区间的RPKM值，即基因组任意指定区间的表达水平。

基因组任意指定区间的RPKM计算公式：

结果：注释参考基因28252个，有14349个基因计算出基因表达水平，对应的外显子一共175642个，其中三组转录组数据中RPKM同时大于1的共有11244个基因和119579外显子。同源异构体80343条，有47141个同源异构体计算出基因表达水平，对应的外显子一共470495个，其中三组转录组数据中RPKM同时大于1的38647个同源异构体和362702外显子。基因组任意指定区间82960个，有24465区间计算出表达水平，其中三组转录组数据中RPKM同时大于1的共有14162个基因组任意指定区间。

计算基因组任意区间RPKM的意义：1、为长非编码RNA的研究提供表达量的计算；2、为microRNA的研究提供表达量的计算；3、可以根据用户自身的需要，计算任意区间的表达量，可以计算该区域所有基因的表达量。多个组织、样品同时计算，可以节省单独计算再整合的时间，并且可以为后续的不同样品中同一基因表达量的比较提供前提。

上述实施例中是以人源胚胎干细胞、海拉细胞和类淋巴母细胞为例进行了说明，当然，本发明的利用转录组测序结果计算基因表达水平的方法还可以用于其他生物的基因、同源异构体和外显子的基因表达水平的计算或者基因组任意指定区间的表达水平计算。另外，上述实施例中采用的所述转录组序列读段由Illumina测序技术获得，也可以由罗氏454测序技术、AB公司的SOLiD技术、或者第三代的单分子实时DNA测序技术获得。

在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

Claims

1.一种测定待测基因组区域表达水平的方法，其特征在于，包括以下步骤：

(3)对定位到基因组的转录组测序读段进行筛选，所述筛选包括去除测序质量≤99.9％的转录组测序读段；

(6)计算能够映射到基因组上的所有测序读段的总数M；

R P K M = \frac{R}{M \times L} \times 10^{9};

其中，所述待测基因组区域包含N个同源异构体，且N≥2；并且，在测定过程中还包括步骤：将各同源异构体的所有外显子进行整合，对于重复的序列区间，仅保留单一序列，从而将同一待测基因组区域中的不同同源异构体的外显子整合成单一序列，将该单一序列的长度作为计算该基因组区域表达水平时的序列长度L。

2.如权利要求1所述的方法，其特征在于，N为2、3、4、5、6、7、8、9、10或大于10。

3.如权利要求2所述的方法，其特征在于，所述方法还包括结果验证步骤，所述结果验证步骤包括：提取待测样品的总RNA，经过反转录得到其cDNA，以cDNA作为模板进行PCR检测，验证待测基因组区域的表达水平。

4.如权利要求3所述的方法，其特征在于，所述待测基因组区域表达水平，为单个基因的表达水平、同一个基因不同的同源异构体的表达水平、所有外显子的表达水平、单个外显子的表达水平以及基因组任意指定区间的表达水平，其中所述基因组任意指定区间包含染色体名称、基因组起始位置和基因组终止位置。

5.如权利要求1所述的方法，其特征在于，所述步骤(1)中，所述转录组序列数据由罗氏454测序技术、Illumina测序技术、AB公司的SOLiD技术、或者第三代的单分子实时DNA测序技术获得。

6.如权利要求1所述的方法，其特征在于，所述步骤(4)中，所述排序方法为：

a.按照每条测序读段定位到基因组的起始位置进行排序；

最后对排序结果建立索引。

7.如权利要求1所述的方法，其特征在于，所述基因组区域选自如下的组：癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域。

8.如权利要求1所述的方法，其特征在于，所述基因组区域为基因组任意指定区间，所述基因组任意指定区间指给定的特定基因组位置信息，包含染色体名称、基因组起始位置和基因组终止位置。

9.一种测定待测基因组区域表达水平的系统，其特征在于，所述系统包括：

(1)比对单元，用于转录组测序读段与基因组序列进行比对；

(5)计算单元，包括：

R P K M = \frac{R}{M \times L} \times 10^{9};

其中，

所述筛选单元中，所述筛选包括去除测序质量≤99.9％的转录组测序读段；和/或，

所述排序单元的排序方法为：

a.按照每条测序读段定位到基因组的起始位置进行排序；

最后对排序结果建立索引。

10.如权利要求9所述的系统，其特征在于，所述基因组区域选自如下的组：癌基因基因组区域、遗传疾病基因组区域和/或长非编码基因区域。

11.如权利要求9所述的系统，其特征在于，所述基因组区域为基因组任意指定区间，所述基因组任意指定区间指给定的特定基因组位置信息，包含染色体名称、基因组起始位置和基因组终止位置。