CN103902852A

CN103902852A - 基因表达的定量方法及装置

Info

Publication number: CN103902852A
Application number: CN201410108121.8A
Authority: CN
Inventors: 黄文潘; 卢志远; 龚梅花; 章文蔚; 席凤; 韩鸿雁
Original assignee: BGI Shenzhen Co Ltd
Current assignee: Wuhan Hua made Dazhi Technology Co. Ltd.
Priority date: 2014-03-21
Filing date: 2014-03-21
Publication date: 2014-07-02
Anticipated expiration: 2034-03-21
Also published as: CN103902852B

Abstract

本发明公开了一种基因表达的定量方法和装置，包括：获取含有核酸序列信息的读段序列；将读段序列与所有参考基因进行比对，获取比对上的读段序列；对比对上的读段序列进行过滤，舍去软剪切比例超过第一预设值，序列长度小于第二预设值，以及比对得分小于第三预设值的读段序列，软剪切比例是指没有比对上的碱基数目占该读段序列总碱基数目的比例；比对得分是按照每个读段序列与参考基因的匹配程度以及读段序列的长度而确定的数值；对于已过滤的读段序列，使用RPKM对目标基因表达进行定量。通过将读段序列与参考基因进行比对，而不是现有的与参考基因组进行比对，可以简化比对过程，提高比对效率。

Description

基因表达的定量方法及装置

技术领域

本发明涉及基因组学及生物信息学技术领域，具体涉及一种基因表达的定量方法及装置。

背景技术

转录组测序技术（RNA-seq，RNA sequencing）是把小RNA（Ribonucleic Acid，核糖核酸）、mRNA和非编码RNA等或者其中一些用高通量测序技术把它们的序列测出来。目前RNA-seq测序平台有多种，包括Hiseq、Roche FLX、IlluminaSolexa、ABI solid等。不同测序平台的测序原理有所不同，但测序步骤基本包括文库制备，聚合酶链式反应（PCR，Polymerase Chain Reaction）扩增等。通过RNA-seq，科研工作者能够获得生物中基因表达的情况，研究不同个体、不同时期、不同形态的组织的基因表达水平的差异。

中国专利申请（申请号：201110283718.2，名称：一种分析基因表达定量的方法）基于Illumina平台公开一种分析基因表达定量的方法，可以克服数字基因表达谱（DGE，Digital Gene Expression）技术对CATG位点和参考基因完整性依赖性强的缺点。但是，该方法测序分析需时较长，劳动效率有待提高。

发明内容

本发明提供一种基因表达的定量方法及装置，可以快速地完成基因表达的定量。

依据本发明的一方面提供一种基因表达的定量方法,包括：获取含有核酸序列信息的读段序列；将读段序列与所有参考基因进行比对，获取比对上的读段序列；对比对上的读段序列进行过滤，舍去软剪切比例超过第一预设值，序列长度小于第二预设值，以及比对得分小于第三预设值的读段序列，软剪切比例是指没有比对上的碱基数目占该读段序列总碱基数目的比例；比对得分是按照每个读段序列与参考基因的匹配程度以及读段序列的长度而确定的数值；对于已过滤的读段序列，使用每百万读段序列中来自目标基因每千碱基长度的读段序列数目RPKM对所述目标基因表达进行定量，定义为RPKM=（比对到目标基因对应的参考基因的读段序列的数目）*10⁹/（比对到所有参考基因的读段序列的数目*目标基因的长度）。

优选地，比对到目标基因对应的参考基因的读段序列的数目是指只能比对到目标基因对应的参考基因上，而且能够比对到所述参考基因的至少一个转录本的读段序列的数目；目标基因的长度是指目标基因的所有转录本中最长的转录本的长度。

依据本发明的另一方面提供一种基因表达的定量装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与数据输入单元、数据输出单元及存储单元数据连接，用于执行存储单元中存储的可执行的程序，该程序的执行包括完成上述基因表达的定量方法。

本发明的有益效果是：通过将读段序列与参考基因进行比对，而不是现有的与参考基因组进行比对，可以简化比对过程，提高比对效率。特别地，比对到目标基因对应的参考基因的读段序列的数目是指只能比对到目标基因对应的参考基因上，而且能够比对到所述参考基因的至少一个转录本的读段序列的数目，则不会认为这部分读段序列是重复比对而需要被过滤，从而提高RPKM和QPCR的相关性，即提高基因表达定量的准确性。

附图说明

图1为现有技术中RNA-seq的流程图；

图2为本发明实施例一的流程图（A）；

图3为本发明实施例一的流程图（B）；

图4为本发明实施例一的读段序列选择示意图；

图5是本发明实施例一的HBRR标准品和QPCR标准的相关性结果图；

图6是本发明实施例一的HBRR标准品的重复性结果图。

具体实施方式

下面通过具体实施方式结合附图对本发明作进一步详细说明。

现有的高通量测序平台有多种，包括Roche454，Ion PGM和Ion Proton等。本发明中的实施例以Ion Proton测序平台作说明，其他测序平台亦同样适用本发明所提供的方法，测序平台并不构成本发明的限制。RNA样本的文库构建一般包括将RNA反转录为DNA来进行文库构建，RNA的提取、构建文库等均可利用现有技术进行，测序文库构建步骤一般包括打断、末端修复、加proton接头、扩增等，请参考图1，测序步骤及参数可以根据不同测序平台的建议操作说明、测试样本种类进行调整，不构成对本发明的限制。实施例中未注明具体条件的，按照常规条件或制造商建议的条件进行；所用试剂或仪器未注明生产厂商的，均为可以通过市面购买获得的常规产品。

实施例一：

本实施例采用RNA样本构建文库。RNA样本使用人组织混合液RNA的微阵列质量控制标准品（UHRR-MAQC，Universal Human Reference RNA-MicroArray Quality Control）和人脑混合液RNA微阵列质量控制标准品（HBRR-MAQC，Human Brain Reference RNA-MicroArray Quality Control），其中UHRR-MAQC标准品采购自安捷伦公司（Agilent Technologies,Inc.），HBRR-MAQC购自Ambion公司。在其他具体实施方式中，亦可以使用其他种类的RNA标准品，或是采购自其他公司所生产的RNA标准品，对本发明不构成限制。

本实施例构建文库的过程如下：取总RNA样品，用DEPC（diethylpyrocarbonate，焦碳酸二乙酯）水稀释，混匀，65℃变性，使用dT（DynalbeadsOligo）₂₅磁珠将总RNA中的信使RNA（mRNA）调取出来并纯化；将所得mRNA与打断试剂混合得到打断的mRNA，再与试剂I混合进行一链合成反应；将一链合成反应后的体系与试剂II混合，进行二链合成反应，反应完成后，用AmpureXP磁珠纯化二链产物；所得二链产物与试剂III混合进行末端修复，并用AmpureXP磁珠纯化末端修复产物；所得末端修复产物与试剂IV混合进行加接头，并用Ampure XP磁珠纯化加接头产物；采用PCR仪扩增，并用Ampure XP磁珠纯化PCR产物，获得测序文库。构建转录本文库或其它RNA文库亦可利用现有方法，文库构建并不构成本发明的限制。

试剂I：0.5μl的100mM二硫苏糖(DTT，DL-Dithiothreitol)、0.5μl的10mM脱氧核糖核苷三磷酸(dNTP Mix，deoxy-ribonucleoside triphosphate)、0.5μl的RNases抑制剂（RNase Inhibitor)。

试剂II：10μl GEX Second Strand Buffer、2μl10mM dNTP Mix，0.2μl逆转录酶RNaseH、2.5μl DNA聚合酶I（DNA Pol I）。

试剂III：5μl10X末端修复缓冲液（End Repair Buffer）、0.4μl25mM dNTPMix、1.2μl T4DNA聚合酶（T4DNA Polymerase）、0.2μl Klenow DNA聚合酶（Klenow DNA Polymerase）、1.2μl T4多聚核苷酸激酶（T4PNK）。

试剂IV：2μl T4DNA连接酶（T4DNA Ligase）、2μl proton Adapter OligoMix(12um)、25μl2X Rapid T4DNA Ligase Buffer。

利用Agilent2100质检构建得的文库，上机测序，获得测序序列，即获得读段序列（reads）。

请参考图2至图6，本实施例提供一种基因定量表达方法，可以快速地完成定量表达。其中在先步骤如文库制备、PCR扩增等采用前述步骤与参数。本实施例具体包括：

S100：获取含有核酸序列信息的读段序列reads

S101：对读段序列进行修剪（trimming）

Trimming可以减少碱基序列在拼接之后产生的错误。在其他具体实施方式中，亦可以不对读段序列进行修剪，直接进行后续步骤；或者使用校正（correct），或修剪与校正结合的方式，以进一步提高测序分析的准确率。Trimming针对读段序列的开头和末尾的的3到4bp，这几个bp通常带有测序接头。包括低质量reads，接头（adapter），基因组3’端位置相同的reads。在高通量测序中，每测一个碱基会给出一个相应的质量值（Q-Value)，可以参考公开号为CN102653784A，名称为《用于多重核酸测序的标签及其使用方法》的中国专利申请。质量值可以反映测序质量的好坏，数值越高表示测序质量越好。因此，低质量reads是指质量值低于y1的碱基的数目超过该reads总碱基数目的y2%，y1的取值范围为15＜y1≤20，y2的取值范围为15＜y2≤25，本实施例取y1为17，y2为20。本领域人员知道，譬如Q20是指质量值大于20的碱基在所有碱基中所占的比例，取值范围为[0,1]，Q20数值越接近1，质量值大于20的碱基在所有碱基中所占的比例越大。因此，低质量reads可以描述为Q(y1)小于（100-y2）%的reads，或其他等同描述方式。譬如本实施例的低质量reads，亦可以描述为Q17小于80%的reads，其中80来源于100-y2=100-20。譬如对于Hiseq测序平台，y1优选设置为20，y2优选设置为20，则低质量reads可以描述为Q20小于80%的reads。y1和y2的取值之间没有必然的数值联系，可以相同或是不同的数值。

在其他具体实施方式中，y1及y2的取值可以根据样品、测试平台等有所调整，y1、y2越高，被筛选的reads越多，即留下的reads越少；y1、y2越低，则被筛选的reads越少，处理效率越慢。

S102：将读段序列与参考基因进行比对，获取比对上的读段序列

基因组作图(genome mapping)是应用界标或遗传学标记对基因组进行精细的划分，进而标示出碱基序列或基因排列。本实施例中利用reads与参考基因进行比对，而不是现有的reads和参考基因组比对，从而提高比对准确性及比对效率。对于真核生物，基因是由基因组中的外显子拼接而成，而测序平台测出来的是拼接之后的序列，直接和参考基因进行比对可以较为直接、准确。另外，在输出比对结果时，本实施例是输出所有的比对匹配结果，即如果有两条以上的读段序列都与参考基因比对匹配，则这两条以上的读段序列都会输出，而不是只输出唯一匹配的reads。一个基因包括多个转录本，很多转录本是来自外显子的不同组合方式，所以有些转录本会有许多同源序列，所以有许多序列会比对到多个转录本上，因此保留所有这些碱基序列，用来判断这些序列是否来自同一个基因。在本实施例中，应用tmap比对工具。tmap是一款适用proton测序平台的商业比对软件，由Life Tech.公司开发。比对的过程主要通过比对得分进行，利用设置基础比对分值，比如本实施例设置基础分为0，reads上的一个碱基位置匹配上参考基因加一分，一个位置错配减一分，该位置缺失计0分等，由此对该read的比对情况进行打分，一般地，一条reads越长，与参考基因匹配程度越高，则其得分越高。在其他具体实施方式中，计分的规则可以根据实现的程序进行调整，譬如基础分为100，每匹配上一个参考基因加0.1分，具体的计分规则不构成本发明的限制。在其他实施方式中，亦可以根据测序平台的不同使用合适的商用比对软件，比如Bowtie、SOAP2、BWA-SW等，或者是自编程序，只要该程序可以达到reads与参考基因进行比对并输出所有的比对匹配结果的目的即可，因此具体的设置参数及比对工具并不构成本发明的限制。

S103：对比对上的读段序列进行过滤

对步骤S102得出的比对读段序列过滤，去掉含软剪切比例超过第一预设值x1的reads，序列长度小于第二预设值x2的reads，以及比对得分小于第三预设值x3的reads。软剪切是指没有比对匹配的reads段，例如一条100bp的reads，共有90bp与参考序列比对匹配，但剩下的10bp没有比对匹配，则这10bp称为软剪切，该reads的软剪切比例为10%。在本实施例中，第一预设值x1为自然数，取值范围是[10%，30%]，优选为20%；x1越大，被过滤的reads数目越多，可能导致后面检测到的基因数目偏少，x1如果过小，则可能导致部分错误的reads没有被过滤掉。第二预设值x2为正整数，取值范围是[15，25]，优选为20，对于过短的序列，如10bp的reads，由于长度较短，可能会比对到参考基因的多个区域。第三预设值x3为正整数，取值范围为[20，50]，x3过低则说明比对匹配的程度过差，易引入错误，x3过高则会导致reads过多的被去掉。值得注意的是，x3的取值范围必然根据步骤S102的比对得分规则而调整，对于本实施例的proton测序平台以及比对得分规则而适用于为[20，50]的取值范围。在其他具体实施方式中，x1、x2、x3的具体数值可以根据测试平台、测试样品进行调整。x1、x2、x3之间没有必然的数值联系，可以相同或是不同的数值。

S104：对基因表达进行定量

本实施例用RPKM来定量，RPKM（reads per kilobase of exon model permillion mapped reads）是目前通用的定量归一化的方法，定义为：

RPKM=（比对到目标基因对应的参考基因的读段序列的数目）*10⁹/（比对到所有参考基因的读段序列的数目*目标基因的长度）。

选取唯一比对到参考基因上的read作为比对到目标基因对应的参考基因的read。对于比对到多个参考基因的read，无法区分来自哪个参考基因，因此将比对到多个参考基因的read都去掉。对于一条read比对到一个参考基因的多个同源转录本，或者一个参考基因的多个位置的情况，则认为只比对到该参考基因一次。当一条read比对到多个转录本时判断所有比对上的转录本是否来自同一个基因，即是否所有比对上的转录本同源，如果判断结果为是，即所有比对上的转录本都是来自同一个基因，则这条read并不是重复比对（multiple map）而不需要去除；如果判断为否，则该条read是multiple map而需要去除，不能作为唯一比对到参考基因上的read。在本实施例中，步骤S102的显示结果可以包括reads比对上哪些转录本，可以有multiple map的显示提示，因此可以利用基因和转录本对应的数据库，来对multiple map的reads进行过滤。然后，统计比对到该参考基因的总reads数目，一个基因可以存在多个转录本或者多个位置，但是这些read都来自同一个参考基因，不会干扰基因的定量，选取该基因的最长转录本代表该基因的长度。基因的长度越长，在同等表达水平下产生的read会比长度短的基因要多。因此在计算RPKM的时候除以基因的长度，能够尽量避免基因长度对定量的影响。请参考图4，以基因A（Gene A）为例进行说明。图4是基因A三个转录本（transcript）的覆盖（coverage）情况，分别是transcript1、transcript2、transcript3。在计算RPKM的时候，覆盖到基因A的read数目为3，包括read1，read2，read3，其中基因的长度我们用最长的转录本3（transcript3）的长度来当做该基因的总长度。对于本实施例中的RPKM计算公式，由于前述步骤中比对、过滤的设置，以及本步骤中对参数的限制选择，使得基因表达的定量变得快速、简单。

本实施例提供的基因表达的定量准确性用QPCR的相关性作评价。这里以皮尔森相关性系数（pearson correlation）作说明。皮尔森相关系数是用来反映两个变量线性相关程度的统计量，皮尔森相关性系数越高，QPCR的相关性越强，基因表达的定量越准确。相关系数用r表示，其中n为样本量，分别为两个变量的观测值和均值。r描述的是两个变量间线性相关强弱的程度，绝对值越大表明相关性越强，具体公式为

在其他具体实施方式中，亦可以与其他相关性系数联合评价，如斯皮尔曼相关性系数(spearman relativity)等。

图5是HBRR标准品和QPCR标准的相关性结果图，其中横坐标是HBRR标准品的proton测序结果计算出来的RPKM值的以10为底的对数值，纵坐标是QPCR值的以10为底的对数值，一个黑点代表一个基因。该标准品的QPCR基因为1000个，即genenum为1000。经计算，pearson correlation可达到0.917，spearman relativity亦可达到0.868。

图6是HBRR标准品的重复性结果图，分别使用了两个HBRR标准品，分别命名为proton_A和proton_B以作说明上的区分，实质并无区别。横坐标是proton_A用proton测序得到的RPKM值的以10为底的对数值，纵坐标是重复proton_B用proton测序得到的RPKM值的以10为底的对数值。基因数目genenum为17463表示，在proton_A和proton_B中均能够检测到的基因个数为17463。图5中的genenum数目与图6的genenum数目不同是因为图5中的genenum中的QPCR结果是标准品RNA提供方Agilent公司提供的经过验证的1000个，而图6中的genenum是proton_A和proton_B都能测出来的基因，但是其中很大部分基因仍未有经过验证的QPCR结果。可以看出，图6的pearson correlation可达到0.997，用spearman relativity亦能达到0.985，说明对于不同的样品的定量结果具有很好的重复性。

对于UHRR的标准品，QPCR的相关性亦达到0.86以上，详细的结果请见表1。以8个样本为例，其中UHRR为4个，HBRR为4个，其中样本的名称不具有实质意义，只是作为不同样本的区分之用。

表1不同样本的基因定量表达评价

然后，可以根据国际标准化的基因功能分类体系Gene Ontology全面描述基因的属性，其中包括基因的分子功能molecular function、所处的细胞位置cellularcomponent、参与的生物过程biological process。亦可以通过比较不同样本间的数据从而筛选出差异表达的基因，后续分析中的差异基因表达模式聚类分析，Gene Ontology功能显著性富集分析，Pathway显著性富集分析，蛋白互作网络分析均是基于差异表达基因。

本领域技术人员可以理解，上述实施方式中各种方法的全部或部分步骤可以通过程序来指令相关硬件完成，该程序可以存储于一计算机可读存储介质中，存储介质可以包括：只读存储器、随机存储器、磁盘或光盘等。

依据本发明的另一方面还提供一种基因定量表达的装置，包括：数据输入单元，用于输入数据；数据输出单元，用于输出数据；存储单元，用于存储数据，其中包括可执行的程序；处理器，与上述数据输入单元、数据输出单元及存储单元数据连接，用于执行存储单元中存储的可执行的程序，该程序的执行包括完成上述实施方式中各种方法的全部或部分步骤。

以上内容是结合具体的实施方式对本发明所作的进一步详细说明，不能认定本发明的具体实施只局限于这些说明。对于本发明所属技术领域的普通技术人员来说，在不脱离本发明构思的前提下，还可以做出若干简单推演或替换。

Claims

1.一种基因表达的定量方法，其特征在于，包括：

获取含有核酸序列信息的读段序列；

将所述读段序列与所有参考基因进行比对，获取比对上的读段序列；

对所述比对上的读段序列进行过滤，舍去软剪切比例超过第一预设值，序列长度小于第二预设值，以及比对得分小于第三预设值的读段序列，所述软剪切比例是指没有比对上的碱基数目占该读段序列总碱基数目的比例；所述比对得分是按照每个读段序列与参考基因的匹配程度以及读段序列的长度而确定的数值；

对于已过滤的读段序列，使用每百万读段序列中来自目标基因每千碱基长度的读段序列数目RPKM对所述目标基因表达进行定量，定义为RPKM=（比对到目标基因对应的参考基因的读段序列的数目）*10⁹/（比对到所有参考基因的读段序列的数目*目标基因的长度）。

2.如权利要求1中所述的方法，其特征在于，所述第一预设值为[10%，30%]，所述第二预设值为[15，25]；其中，第一预设值优选为20%，第二预设值优选为20。

3.如权利要求1中所述的方法，其特征在于，所述比对到目标基因对应的参考基因的读段序列的数目是指只能比对到目标基因对应的参考基因上，而且能够比对到所述参考基因的至少一个转录本的读段序列的数目。

4.如权利要求3中所述的方法，其特征在于，所述目标基因的长度是指目标基因的所有转录本中最长的转录本的长度。

5.如权利要求1中所述的方法，其特征在于，所述读段序列与参考基因进行比对的步骤之前还包括：

对读段序列进行修剪，去除接头、低质量的读段序列和3’端位置相同的读段序列，所述低质量的读段序列是指测序准确度低于第四预设值的碱基的数目超过该读段序列总碱基数目的第五预设值的读段序列。

6.如权利要求5中所述的方法，其特征在于，所述第四预设值取值范围为[10%，30%]，所述第五预设值取值范围为[15%，25%]；

优选地，当读段序列获自Ion Torrent测序平台，所述第四预设值为17%，所述第五预设值取值范围为20%。

7.如权利要求5中所述的方法，其特征在于，所述读段序列的长短不一。

8.如权利要求5中所述的方法，其特征在于，所述读段序列长度均一，所述第四预设值为20%，所述第五预设值取值范围为20%。

9.如权利要求1中所述的方法，其特征在于，所述核酸序列信息来源于核糖核酸。

10.一种基因定量表达装置，其特征在于，包括：

数据输入单元，用于输入数据；

数据输出单元，用于输出数据；

存储单元，用于存储数据，其中包括可执行的程序；

处理器，与所述数据输入单元、数据输出单元及存储单元数据连接，用于执行所述可执行的程序，所述程序的执行包括完成如权利要求1-9任意一项所述的方法。