WO2013041021A1

WO2013041021A1 - 一种分析基因表达定量的方法

Info

Publication number: WO2013041021A1
Application number: PCT/CN2012/081612
Authority: WO
Inventors: 章文蔚; 张艳艳; 龚梅花; 彭智宇; 韩祖晶; 高欢; 李英睿; 王俊; 汪建; 杨焕明
Original assignee: 深圳华大基因科技有限公司; 深圳华大基因研究院
Priority date: 2011-09-22
Filing date: 2012-09-19
Publication date: 2013-03-28
Also published as: CN103014137A; CN103014137B

Abstract

本发明公开了一种分析基因表达定量的方法，所述方法包括：从总 RNA中纯化mRNA，制备片段化的mRNA；将所述片段化的mRNA逆转录，制备得到 cDNA，将所述cDNA纯化后，制备为平末端 DNA，纯化所述平末端 DNA；将所述平末端 DNA片段制备得到末端加"A"碱基的 DNA片段；在所述末端加"A"碱基的 DNA片段两端加接头序列，得到两端加接头序列的 DNA片段，并进行纯化，对所述两端加接头序列的 DNA片段进行 PCR反应，纯化 PCR反应产物；对所述 PCR反应产物测序；对所述测序得到的数据进行过滤，去除不合格序列，得到干净序列，利用短序列映射程序将所述干净序列与参考序列比对，对所述比对结果进行分析。

Description

一种分析基因表达定量的方法

技术领域

本发明涉及核酸测序技术领域，特别是 RNA-seq技术领域以及测序后信息分析的方法。背景技术

目前，基因表达定量研究领域主要有两种技术：传统的芯片技术和测序技术。其中，芯片技术通量高，自动化，成本低，但是芯片技术依赖于已知基因，信号噪音高，重复性差，检测阈值窄；测序技术又分为 SAGE ( Serial Analysis of Gene Expression) ，数字基因表达谱 ( Digital Gene Expression, DGE) 和数字基因表达谱升级版 RNA-Seq(Quantification)技术，其中， SAGE技术测序准确，但操作繁琐，测序成本高。基于第二代高通量测序平台的 DGE和 RNA-Seq 技术克服了芯片技术和 SAGE技术的缺点，它们通量高，自动化，测序成本低，噪音小，不依赖于已知基因，检测阈值宽。

但是 DGE 由于实验本身的局限性，导致了该项技术不能够检测到不含 CATG位点的基因，并且 DGE技术在研究基因表达定量时对参考基因的依赖性很强，对于一些非模式生物的定量分析也存在一定的局限性。

以 illumina测序平台为代表的第二代高通量测序技术不仅节省了大量的人力和物力，而且还具有测序通量高、准确度高和成本低的众多优点。目前该平台已经广泛应用于：全基因组测序，新物种测序，目标基因组测序，转录组和表观遗传分析等领域。

随着第二代高通量 illumina测序平台的广泛应用，多物种基因组测序和全基因组研究的大规模开展，降低测序成本，减少测序流程，提高劳动效率成为测序技术的一个重要研究方向。而基于 illumina测序平台 RNA-seq的基因表达分析存在步骤多，成本高，操作过程繁琐，不适合用于自动化工作站等缺点。发明内容

本发明的一个方面，提供了一种分析基因表达定量的方法，包括下述步骤： ( 1 ) 从总 RNA中纯化 mRNA，制备片段化 mRNA;

(2) 将所述片段化 mRNA逆转录制备得到 cDNA，将所述 cDNA纯化后制备为平末端 DNA，纯化平末端 DNA;

(3) 将所述平末端 DNA的末端加 "A" 碱基，得到末端加 "A" 碱基的

DNA;

(4) 在末端加 "A"碱基的 DNA两端加接头序列，纯化两端加接头序列的 DNA进行 PCR反应，纯化 PCR反应产物；

(5) 对所述 PCR反应产物测序；

(6) 将测序得到的数据过滤不合格序列得到干净序列，将所述干净序列与参考序列比对，对比对结果进行分析。

在本发明的一个实施方案中，所述的总 RNA的选取量为 0.1 μ g〜2w g。在本发明的一个实施方案中，使用 Invitrogen公司生产的 Oligo (dT) 25(产品号 610.06)磁珠从总 RNA中纯化 mRNA。

在本发明的一个实施方案中，使用 Beckman公司生产的 Ampure XP磁珠 (产品号 A63882)纯化所述 cDNA、两端加接头序列的 DNA、 PCR反应产物。

在本发明的一个实施方案中，在步骤（1) 中，使用试剂 I制备片段化的 mRNA, 所述试剂 I含有： 10-400mM可溶性盐， 200mM-300mM缓冲盐， pH 8.0-8.5，溶剂为水。优选地，试剂 I 中缓冲盐选自： Tris-HCl、磷酸盐。优选地，试剂 I中可溶性盐选自氯化钠，氯化镁。优选地， mRNA与试剂 I混合温度为 65°C〜94°C。

在本发明的一个实施方案中，在步骤（2) 中，使用试剂 II对 cDNA进行末端修复，得到平末端 DNA，所述试剂 II含有： 1.2 LT4DNA 聚合酶 (3U/ L)， 1.2 L T4多聚核苷酸激酶（10U/ L)， 0.2μ1 Klenow DNA 聚合酶（5U/ L)， 0AμL 25mM dNTP; T4多聚核苷酸激酶缓冲液含有 700 mM Tris-HCl, 100 mM氯化镁， 50mMDTT。

在本发明的一个实施方案中，在步骤（3) 中，使用试剂 III对所述平末端 DNA的末端加 " A"碱基，所述试剂 III含有： 100mM-500mM 可溶性盐， 100 mM 缓冲盐， 10mM-50mM 二硫苏糖醇， 5mM dATP， 0.2 μ L Klenow(3 ' -5' exo)酶 (5U/ L)， pH7.6-7.9, 溶剂是水。优选地，试剂 III中缓冲盐选自 Tris-HCl、磷酸盐。优选地，试剂 III中可溶性盐为氯化钠。优选地，样品与试剂 III混合温度为 16°C-37°C o

在本发明的一个实施方案中，在步骤（4 ) 中，使用试剂 IV在末端加 "A " 碱基的 DNA两端加接头序列，所述试剂 IV含有： 100 mM 缓冲盐溶液， 10mM〜 50mM 二硫苏糖醇， 5〜10mM ATP， 1.2 L T4 DNA 连接酶， pH值为 7.6〜7.9，溶剂是水。优选地，缓冲盐溶液为 Tris-HCl、磷酸盐缓冲溶液。

在本发明的一个实施方案中，在步骤（5 ) 对所述 PCR产物测序前，还包括步骤：采用 Agilent Bioanalyzer 2100和 Q-PCR检测 DNA浓度及 DNA片段大小。

在本发明的一个实施方案中，所述的测序为高通量测序技术。优选地，为 illumina solexa测序技术。

在本发明的一个实施方案中，在步骤（6 ) 中，所述的不合格序列包括：测序质量低于预定阈值的碱基个数超过整条序列碱基个数的 50%的序列，序列中测序结果不确定的碱基个数超过整条序列碱基个数的 10%的序列，除样本接头序列外引入的外源序列。

在本发明的一个实施方案中，在步骤（ 6 ) 中，所述比对使用 SOAPaligner/soap2。

在本发明的一个实施方案中，在步骤（6 ) 中，所述对比对结果进行分析包括：高通量测序的质量评估，基因表达量的统计，差异表达基因筛选，实验重复性分析，差异基因表达模式聚类分析， Gene Ontology ( GO ) 功能显著性富集分析，通路（Pathway) 显著性富集分析，蛋白相互作用网络分析。

双末端测序

对基因片段 (包括 DNA片段和 cDNA片段)进行测序，其测序对象都是一段物理连续的碱基序列片段，该片段称为插入片段，其长度称为插入片段长度 (insertsize)。

如本文所用，术语"双末端测序"是对该片段的两侧碱基序列从边缘向内部的测序，测得的序列称为读序 (read) , 长度称为读长 (read-length)。两侧测得的读序是来自于同一个插入片段，并且其末端距离为 insertsize, 故两侧读序的配对关系确定。这两个读序被称为配对读序 (Pair-end reads)。高通量测序

基因组的高通量测序使得人类能够尽早地发现与疾病相关基因的异常变化，有助于对个体疾病的诊断和治疗进行深入的研究。本领域技术人员通常可以采用三种第二代测序平台进行高通量测序，如： 454FLX(Roche公司）、 Solexa Genome Analyzer(Illumina测序平台)禾卩 Applied Biosystems 公司的 SOLID等。这些平台共同的特点是极高的测序通量，相对于传统测序的 96道毛细管测序，高通量测序一次实验可以读取 40万到 400万条序列，根据平台的不同，读取长度从 25bp到 450bp不等，因此不同的测序平台在一次实验中，可以读取 1G到 14G 不等的碱基数。

本发明优选 Illumina测序平台，它包括 DNA簇形成和上机测序两个步骤： PCR扩增产物的混合物与固相载体上固定的测序探针进行杂交，并进行固相桥式 PCR扩增，形成测序簇；对所述测序簇用"边合成 -边测序法"进行测序，从而得到样本中核酸分子的序列。

DNA簇的形成是使用表面连有一层单链引物 (primer)的测序芯片 (flow cell), 单链状态的 DNA片段通过接头序列与芯片表面的引物通过碱基互补配对的原理被固定在芯片的表面，通过扩增反应，固定的单链 DNA变为双链 DNA，双链再次变性成为单链，其一端锚定在测序芯片上，另一端随机和附近的另一个引物互补从而被锚定，形成"桥"；在测序芯片上同时有上千万个 DNA单分子发生以上的反应；形成的单链桥，以周围的引物为扩增引物，在扩增芯片的表面再次扩增，形成双链，双链经变性成单链，再次成为桥，称为下一轮扩增的模板继续扩增；反复进行了 30轮扩增后，每个单分子得到 1000倍扩增，称为单克隆的 DNA簇。

DNA簇在 Solexa测序仪上进行边合成边测序，测序反应中，四种碱基分别标记不同的荧光，每个碱基末端被保护碱基封闭，单次反应只能加入一个碱基，经过扫描，读取该次反应的颜色后，该保护集团被除去，下一个反应可以继续进行，如此反复，即得到碱基的精确序列。在 Solexa多重测序 (Multiplexed Sequencing)过程中会使用 Index(标签)来区分样品，并在常规测序完成后，针对 Index部分额外进行测序，通过 Index的识别，可以在 1条测序甬道中区分多达 12 种不同的样品。本发明的主要优点如下：

NA-Seq应用于基因表达定量研究克服了 DGE技术对 CATG位点和参考基因完整性依赖性很强的缺点同时也克服了芯片技术检测阈值窄，噪音污染大等缺点。从而真正地达到定量准、可重复性高、费用低、检测阈值宽、信号噪音小等优点。附图说明

图 1示出本发明文库构建流程图；

图 2示出本发明的信息分析流程图；

图 3示出图 2所示应用例中样本一 Reads在参考基因组 chrl O上的分布图; 图 4示出图 2所示应用例中两次平行实验的结果相关性分析结果图；图 5示出图 2所示应用例中样本一测序 reads在基因上的分布图。具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图及实施例，对本发明进行进一步详细说明。实施例中未注明具体条件者，按照常规条件或制造商建议的条件进行。所用试剂或仪器未注明生产厂商者，均为可以通过市购获得的常规产品。实施例两个人体组织样本的 RNA-seq分析

组织样本由北京大学提供。

文库构建过程如图 1所示。具体如下：取 0.1 w g〜2 w g的总 RNA样品，用脱氧核糖核酸酶 K Dnasel )进行消化，乙醇沉淀纯化消化后产物，使用 Oligo ( dT ) 25磁珠将所得总 RNA中的 mRNA调取出来并纯化，将所得 mRNA与试剂 I混匀反应，得到片段化的 mRNA, 所得 mRNA与试剂 I混匀反应得到的片段化的 mRNA, 经反转录合成 cDNA , 使用 Ampure XP磁珠纯化产物，所得 cDNA与试剂 II混匀反应，形成平末端的 DNA片段，使用 Ampure XP磁珠纯化产物，所得平末端 DNA片段与试剂 III混匀反应，得到 3 ' 端加上一个" A" 碱基的 DNA片段，与试剂 IV混匀反应，得到两端加接头的 DNA片段，使用 Ampure XP磁珠纯化产物，采用聚合酶链式反应（PCR)扩增所得 DNA片段， Ampure XP磁珠纯化 PCR产物，上机测序。测序使用 Illumina Hiseq2000。

试剂 I为： 10-400mM氯化镁， 200mM-300mM Tris-HCl， pH 8.0-8.5 , 溶剂为水。

试剂 II为： 1.2uLT4 DNA 聚合酶（3U/ L)， 1.2uLT4 多聚核苷酸激酶 (10U/ L),0.2ulKlenow DNA 聚合酶（5U/ L)，0.4uL 25mM dNTP; T4多聚核苷酸激酶缓冲液含有 700 mM Tris-HCl， 100 mM 氯化镁， 50 mM DTT。

试剂 III为： 100 mM -500mM 氯化钠， 100 mM Tris-HCl， 10mM -50mM 二硫苏糖醇， 5mM dATP， 0.2 μ L Klenow(3 ' -5 ' exo)酶（5U/ L)， pH7.6-7.9, 溶剂是水。

试剂 IV为： 100 mM Tris-HCl， 10mM〜50mM 二硫苏糖醇， 5〜 lOmM ATP， 1.2 L T4 DNA 连接酶， pH值为 7.6〜7.9，溶剂是水。

图 2 示出了本发明实施例提供的数字基因表达谱升级版 NA- Seq(Quantification)生物信息学分析方法的实现流程，详述如下：

在步骤 S1中，接收高通量测序技术得到的测序片段。在本发明实施例中，采用 Illumina Hiseq2000测序。接收到原始测序序列后，对原始测序序列进行过滤，去除不合格的序列。不合格序列包括：测序质量值低于 5的碱基个数超过整条序列碱基个数的 50%则认为是不合格序列；序列中测序结果中测序结果不确定的碱基个数超过整条序列碱基个数的 10%则认为是不合格序列；与测序接头序列库进行比对，若序列中存在测序接头序列则认为是不合格序列。

在步骤 S2中，将每个序列中样本接头序列与样本接头序列库进行比对，实现分样本操作，同时将样本接头序列从序列片段中去除。将接头序列（本实施例为 8bp ) 中有测序质量低于 5的碱基个数大于 3个的序列去除。

在步骤 S3中，本发明实施例采用 SOAPaligner/soap2，将高通量测序技术得到的测序片段比对到参考基因组序列上。

在步骤 S4中，本发明实施例主要是以图形的方式概括地给出 Reads在基因组各个位置大致的分布情况，以及该位置基因的分布情况。如图 3画出 Reads 在最长的 1条染色体（或 Scaffold) 上的分布图，样本一 Reads在参考基因组 chrlO上的分布。其中 Gene指每个窗口中 gene的个数， Coverage指每个窗口下被 reads覆盖的区域与窗口长度之比， Reads指每个窗口的平均测序深度，数值取了 log2。

在步骤 S5中，是用来衡量样品的测序量多少的标准，随着测序量（reads 数量）的增多，检测到的基因数也随之上升，当测序量达到某个值时，其检测到的基因数增长速度趋于平缓，说明检测到的基因数趋于饱和。

在步骤 S6 中，本发明用 RPKM法计算基因的表达量，其计算公式为：

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

然后，本发明根据国际标准化的基因功能分类体系 Gene Ontology全面描述基因的属性，其中包括基因的分子功能（molecular function) 、所处的细胞位置 ( cellular component) 、参与的生物过禾呈 ( biological process ) 。

在步骤 S7中，本发明通过比较不同样本间的数据从而筛选出差异表达的基因，后续分析中的差异基因表达模式聚类分析， Gene Ontology功能显著性富集分析， Pathway显著性富集分析，蛋白互作网络分析均是基于差异表达基因。

参照 Audic S.等人发表在 Genome Research上的基于测序的差异基因检测方法 (Audic S. and Claverie J. The Significance of Digital Gene Expression Profiles. Genome Research, 1997 7: 986-995.)，筛选两样本间的差异表达基因。

表达模式相似的基因通常具有相似的功能。我们利用 cluster 软件，以欧氏距离为距离距阵计算公式，对差异表达基因和实验条件同时进行等级聚类分析。

GO功能显著性富集分析提供与参考基因比较后，在差异表达基因中显著富集的 GO功能条目，并筛选出差异表达基因与哪些生物学功能显著相关。该分析首先把所有差异表达基因向 Gene Ontology 数据库 ( http://www.geneontology.org/)的各个 term映射，计算每个 term的基因数目，然后应用超几何检验，找出与整个基因组背景相比，在差异表达基因中显著富集的 GO条目。

在生物体内，不同基因相互协调行使其生物学功能，基于 pathway的分析有助于更进一步了解基因的生物学功能。 KEGG是有关 pathway的主要公共数据库， pathway显著性富集分析以 KEGG pathway为单位，应用超几何检验，找出与整个基因组相比较后差异表达基因中显著性富集的 pathway。

蛋白互作网络分析整合了 BIND、 BioGrid, HPRD等相互作用网络数据库的信息，结果文件中的网络由差异表达基因以及跟差异表达基因有直接相互作用的基因组成。

在步骤 S8中，本发明对两次平行实验的结果相关性分析可获得对实验结果可靠性和操作稳定性的评估。如图 4所示，若同一样本两次平行实验之间的相关性越接近 1，说明可重复性越高。

在步骤 S9中，本发明以 reads在参考基因上的分布情况来评价 mRNA打断的随机程度。由于不同参考基因有不同的长度，我们把 reads在参考基因上的位置标准化到相对位置（reads 在基因上的位置与基因长度的比值），然后统计基因的不同位置比对上的 reads数。如果打断随机性好， reads在基因各部位应分布得比较均匀。图 5给出的是样本一测序 reads在基因上的分布。本发明的描述是为了示例和描述起见而给出的，而并不是无遗漏的或者将本发明限于所公开的形式。很多修改和变化对于本领域的普通技术人员而言是显然的。选择和描述实施例是为了更好说明本发明的原理和实际应用，并且使本领域的普通技术人员能够理解本发明从而设计适于特定用途的带有各种修改的各种实施例。

Claims

权利要求

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

( 1 ) 从总 RNA中纯化 mRNA，制备片段化 mRNA;

( 2) 将所述片段化 mRNA逆转录制备得到 cDNA，将所述 cDNA纯化后制备为平末端 DNA，纯化所述平末端 DNA;

( 3 ) 将所述平末端 DNA片段制备得到末端加 "A" 碱基的 DNA片段；

( 4) 在所述末端加 "A" 碱基的 DNA片段两端加接头序列，得到两端加接头序列的 DNA片段并进行纯化，对所述两端加接头序列的 DNA片段进行 PCR 反应，纯化 PCR反应产物；

( 5 ) 对所述 PCR反应产物测序；

( 6) 将所述测序得到的数据过滤不合格序列得到干净序列，利用短序列映射程序将所述干净序列与参考序列比对，对所述比对结果进行分析。

2、根据权利要求 1所述的方法，其特征在于，在步骤（1 ) 中，所述的总 RNA的选取量为 0.1 μ g〜2 w g。

3、根据权利要求 1所述的方法，其特征在于，在步骤（1 ) 中，使用 Oligo ( dT) 25磁珠从总 RNA中纯化 mRNA。

4、根据权利要求 1所述的方法，其特征在于，使用 Ampure XP磁珠纯化所述 cDNA、两端加接头序列的 DNA片段、 PCR反应产物。

5、根据权利要求 1所述的方法，其特征在于，在步骤（5 ) 的对所述 PCR 反应产物测序前，还包括步骤：采用 Agilent Bioanalyzer 210(^nQ-PClM^IlDNA 浓度及 DNA片段大小。

6、根据权利要求 1所述的方法，其特征在于，在步骤（5 ) 中，所述测序使用高通量测序技术。

7、根据权利要求 1所述的方法，其特征在于，在步骤（6) 中，所述的不合格序列包括：测序质量低于预定阈值的碱基个数超过整条序列碱基个数的 50%的序列，序列中测序结果不确定的碱基个数超过整条序列碱基个数的 10% 的序列，除样本接头序列外引入的外源序列。

8、根据权利要求 1所述的方法，其特征在于，在步骤（6) 中，所述的短序列映射程序选用 SOAPaligner/soap2。

9、根据权利要求 1所述的方法，其特征在于，在步骤（6 ) 中，对比结果所进行的生物信息分析包括：高通量测序的质量评估、基因表达量的统计、差异表达基因筛选、实验重复性分析、差异基因表达模式聚类分析、 Gene Ontology 功能显著性富集分析、通路显著性富集分析、蛋白相互作用网络分析。

10、根据权利要求 1所述的方法，其特征在于，步骤（3 ) 和步骤（4 ) 中，添加的 " A " 碱基的数量为一个；以及，

在步骤（6 ) 中，通过分析，获得基因表达定量的数据。