专利名称：对宏基因组16s高可变区v6进行测序聚类分析的方法
技术领域：
本发明涉及微生物基因测序分析技术领域，尤其涉及一种对宏基因组16S高可变区V6进行测序聚类分析的方法。
背景技术：
为了研究生物环境中微生物群体的种类，一般传统的方法包括直接对微生物进行培养，变性梯度凝胶电泳(DGGE，Denaturing Gradient Gel Electrophoresis)，末端限制性内切酶片段长度多态性(T-RFLP，^Terminal Restriction Fragment Length Polymorphism),焚光原位杂交(FISH, Fluorescence In Situ Hybridization),对可能的微生物种类进行PCR(聚合酶链式反应，Polymerase Chain Reaction)；但这些方式都只能揭露环境中很小一部分微生物种类。如果能进行宏基因组的分析，通过直接对环境中的微生物群体进行基因组研究，得到一个比较全面的微生物种类目录，有助于对微生物群体的后续研究和应用。
由于原核生物中16S rRNA (核蛋白核糖核酸，ribosomal RNA (RiboNucleic Acid))的序列高度保守，可精确指示细菌之间的亲缘关系；16S rRNA的大小为1500bp (碱基对，Base Pair)左右，所含信息能反映生物界进化关系，易操作，适用于各级分类单元；所以在宏基因组的研究中，16S区测序是最常用的聚类和分类方法。传统的宏基因组的测序是通过Sanger技术测序16S rRNA gene(16S rDNA)得到至少500bp的读长，这个读长的长度足够长，能够装配出近乎完整的16S rDNA序列，帮助我们去精准地研究每一条序列的物种来源，但它容易产生嵌合体，而且测序成本比较高，费时又费力。
随着新开发出的测序技术以及测序成本的逐步降低，宏基因组的研究变得越来越实用，所涉及的技术包括Pyrosequencing、Solexa等。对于这些革命性的技术的一个主要挑战就是读长太短，无法对每个个体的16S rDNA进行测序，因而它的测序信息不足以让我们去精准地对微生物进行分类。
综上所述，提供一种更加准确地对微生物进行聚类分析的方法且方便快捷、成本低廉成为本领域亟待解决的技术问题。发明内容
本发明要解决的一个技术问题是提供一种对宏基因组16S高可变区V6进行测序聚类分析的方法，通过对16S的高可变区V6区进行solexa测序，并通过对这些16S可变区的短序列进行系统分类，可以在成本低廉的基础上准确反映物种的丰度信息。
本发明的一个方面提供了一种对宏基因组16S高可变区V6进行测序聚类分析的方法，该方法包括提取微生物的脱氧核糖核酸DNA ；通过引物对宏基因组16S核糖体脱氧核糖核酸rDNA的高可变区V6进行聚合酶链式反应PCR，并为每个样品加上标签序列；把不同样品的PCR产物进行混合；对混合后的PCR产物进行Solexa建库法建库；使用Solexa测序工具对高可变区V6的文库进行双末端pair-end测序，得到原始的测序数据；对测序数据进行筛选，以过滤掉低质量的数据；利用重叠群的关系对高可变区V6的全长序列进行组装；通过标签序列把reads分配到对应的样品上；通过对reads进行分类分析，以实现使用高可变区的测序对微生物群体进行高通量的分类。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，该方法还包括在步骤“提取微生物的脱氧核糖核酸DNA”之前，执行微生物群体的取样。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，该方法还包括在步骤“通过对reads进行分类分析”之后，对不同差异度的序列进行操作分类学单元OTU的分类；根据标签序列和reads，进行种群多样性估计Chaol算法和血管紧张素转化酶ACE的多样性分析。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，在进行种群多样性估计Chaol算法和血管紧张素转化酶ACE的多样性分析之后，输出微生物群体的多样性分析图和相对丰度图。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，步骤“通过引物对宏基因组16S核糖体脱氧核糖核酸rDNA的高可变区V6 进行聚合酶链式反应PCR，并为每个样品加上标签序列”进一步包括使用引物967f CNACGCGAAGAACCTTANCGeq ID NO 1) ^P 1406R :GACAGCCATGCANCACCT(Seq ID NO 2)去复制微生物群体中细菌的16S高可变区V6区片段；对每个微生物样品加标签序列，标签序列被加到引物967f的5，端的前面，以及在标签序列和引物967f之间加上碱基GT。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，该方法还包括对于古生菌的高可变区V6的聚合酶链式反应PCR，使用引物958AR: AATTGGANTCAACGCCGG(Seq ID NO 3)和 1048AR :CGRCGGCCATGCACCWC(kq ID NO :4)。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，步骤“把不同样品的PCR产物进行混合”进一步包括对16S的高可变区V6的PCR产物进行浓度定量；以及按照等摩尔的量混合在一起。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，步骤“对混合后的PCR产物进行Solexa建库法建库”进一步包括把混合产物进行纯化，末端修复，在3’端加上碱基A，加上双末端I^ir-end测序接头；加完接头后，对样品进行纯化；对纯化后的样品进行溶解，并作为模板进行聚合酶链式反应PCR扩增；以及对聚合酶链式反应PCR产物进行胶纯化。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，低质量的数据与最邻近的引物不匹配的序列、小于50碱基对的序列，或者具有至少一个引起歧义碱基的序列。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，步骤“利用重叠群的关系对高可变区V6的全长序列进行组装”进一步包括采用高可变区V6的PCR产物5’端的前75、70、65、60和55碱基对来进行重叠从而组装；其中，组装的标准是一对序列具有大于5个碱基对的重叠长度和在重叠区域小于10%的不匹配度。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，步骤“通过对reads进行分类分析”进一步包括将分配到对应样品上的reads比对5到现有16S v6数据库中，来达到使用高可变区的标签测序对微生物群体进行高通量的分类分析，进而研究微生物群体的结构。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法，采用结合了加标签技术的Solexa技术，对特定环境下的微生物群体进行了高通量测序，既减少了人力劳动也节省了经济花费，使得在研究微生物群落结构与健康、环境因子等等的关系上变得容易可行。
图1示出本发明实施例提供的一种对宏基因组16S高可变区V6进行测序聚类分析的方法的流程图2示出本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的另一个实施例的流程图3示出不同环境下微生物群体在0. 03和0. 3的差异度情况下的OTOs的数量。
具体实施方式
下面参照附图对本发明进行更全面的描述，其中说明本发明的示例性实施例。
图1示出本发明实施例提供的一种对宏基因组16S高可变区V6进行测序聚类分析的方法的流程图。
如图1所示，对宏基因组16S高可变区V6进行测序聚类分析的方法流程100包括
步骤102，提取微生物的脱氧核糖核酸DNA。例如，采用Ultraclean Soil DNA kit 试剂盒(MoBio，USA)从样品沉积物中提取微生物的DNA。
步骤104，通过引物对宏基因组16S核糖体脱氧核糖核酸rDNA的高可变区V6 (该区的两端各有20碱基对bp左右的保守区，中间的可变区为60-90bp左右)进行聚合酶链式反应PCR，并为每个样品加上标签序列。例如，使用引物967f :CNACGCGAAGAACCTTANC(Seq ID NO 1)和 1406R :GACAGCCATGCANCACCT (Seq ID NO :2)去复制微生物群体中细菌的 16S 高可变区V6区片段；并对每个微生物样品加标签序列，标签序列被加到引物967f的5’端的前面，以及在标签序列和引物967f之间加上碱基GT (即碱基G和T)。其中，标签序列可以是由8个碱基组成条形码序列，标签序列的设计要符合一定规则，比如碱基含量和不同碱基数目等，目的是防止因为个别测序错误等原因导致标签相互之间的混淆，例如可以参考美国专利申请US20100267043A1中公开的方法和原则。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个实施例中，对于古生菌的高可变区V6的聚合酶链式反应PCR，使用引物958AR: AATTGGANTCAACGCCGG(Seq ID NO 3)和 1048AR =CGRCGGCCATGCACCffC(Seq ID NO 4)。
步骤106，把不同样品的PCR产物进行混合。例如，采用分光光度计(如Nanodrop) 对16S的高可变区V6的PCR产物进行浓度定量，之后再按照等摩尔的量混合在一起。
步骤108，对混合后的PCR产物进行Solexa建库法建库。例如，把混合产物用QIAquick PCR purification Kit(Qiagen)进行纯化，末端修复(即通过酶反应使所有DNA双链的粘性末端成为平末端)，加“A”，加上I^air-end的接头(Pair-end librarypreparation kit, Illumina)；加完接头后，对样品进行纯化；对纯化后的样品进行溶解，并作为模板进行聚合酶链式反应PCR扩增(12CyCleS)；以及用OlIAquick gel extraction kit,Qiagen)对聚合酶链式反应PCR产物进行胶纯化(即点样电泳，在DNA位置切胶，用试剂盒纯化)。
步骤 110，使用 Solexa 测序工具(如 Illumina GA, illumina GA2, illumina Hiseq2000, illumina HiseqlOOO等)对高可变区V6的文库进行pair-end测序，得到原始的测序数据。例如，直接用Illumina GA II(75bp pair-end策略)进行测序。Solexa测序仪(Illumina genome analyzer)是新一代的高通量测序仪，测序价格低廉，数据读取量大， 相同的测序量的情况下，Solexa测序费用是妨4测序费用的十分之一，而且错误率低(如单碱基测序错误率< 10_5)，测序无偏性，对于宏基因组，可以真实反映物种的丰度信息。
步骤112，对测序数据进行筛选，以过滤掉低质量的数据。例如，低质量的数据选自以下序列中的任意一种与最邻近的引物不匹配的序列、小于50碱基对bp的序列，或者具有至少一个引起歧义碱基的序列。
步骤114，利用重叠群的关系对高可变区V6的全长序列进行组装。例如，采用高可变区V6的PCR产物5’端的前75、70、65、60或55碱基对bp来进行重叠从而组装；其中，组装的标准可以是一对序列具有大于^P的重叠长度，并且在重叠区域小于10%的不匹配度 (即高于90%的匹配度)。
步骤116，通过标签序列把reads分配到对应的样品上。
步骤118，通过对reads进行分类分析，以实现使用高可变区的测序对微生物群体进行高通量的分类。例如，通过GAST软件将分配到对应样品上的reads比对到数据库16S v6数据库refhvr_V6中，来达到使用高可变区的标签测序对微生物群体进行高通量的分类分析，进而研究微生物群体的结构。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法，通过把16S 高可变区的序列与rRNA数据库进行比较，在最好匹配的基础上进行分类。该分类方法在微生物群体的构成和多样化上都可以提供信息，其与16S测全长在微生物分类和测量群体的相对丰度上具有等同的技术效果；此外，本发明采用大量平行测序能够发现更多的稀有微生物种类。
进一步的，由于Solexa的读长为75bp左右，它的通量大，产出数据多，这种方法在探究微生物群落(包括稀薄生物圈)结构的变化上具有很好的成本效益。
图2示出本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的另一个实施例的流程图。
如图2所示，对宏基因组16S高可变区V6进行测序聚类分析的方法流程200包括步骤201,202-218,219和220，其中步骤202_218、204、206和208可以分别执行与图1 所示的步骤102-118相同或相似的技术内容，为简洁起见，这里不再赘述其技术内容。
如图2所示，在步骤202 “提取微生物的脱氧核糖核酸DNA”之前，执行步骤201， 微生物群体的取样。例如，从湖泊等水域中提取沉淀物作为取样样本。
在步骤218 “通过对reads进行分类分析”之后，执行步骤219，对不同差异度的序列进行操作分类学单元(OTU)的分类。例如，利用V. 1. 6. 0版本的Mothur软件(下载网址为http://www. mothur. org/wiki/Main_Page)，采用GAST-0TU策略对不同差异度的序列进行OTU的分类。
步骤220，根据标签序列和reads，利用Mothur. Canoco软件进行Chaol和血管紧(ACE, Angiotensin Converting Enzyme)白勺个生 ^。
本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法，尽管Solexa 测得的16S高可变区V6的reads很短并不能包含足够的进化信息去推论出系统分类，但本发明利用搜索软件如GAST，Mothur软件等，通过比对每个样品的reads到数据库16S v6 region database refhvr_V6中，来达到使用高可变区的标签测序对微生物群体进行高通量的分类分析。总之，使用Solexa测序技术对微生物样本进行测序能够在通量、成本和有效的分类效果上达到很好的平衡，此外，本发明采用的测序技术结合了标签序列，在分辨率上大大提高，单个Run Ifelexa(Illumina)能产生比妨4多100倍的reads。因此，仅仅通过测序16S rRNA V6区这么短的长度就能得到很好的分类效果，另外由于结合了标签技术， 测的长度相对比较短，因而能够在单个Lanedllumina高通量测序仪一张芯片有8个通道， 每个通道被称为“lane”)上点更多的样，大大节约了每个样品的测序成本。
接下来详细描述本发明提供的对宏基因组16S高可变区V6进行测序聚类分析的方法的一个具体实施方式
。
步骤1、进行微生物群体的取样。
具体来说，分别提取深圳-北山水库沉积物、深圳-仙湖植物园沉积物、深圳-红树林沉积物、深圳-大梅沙沉积物、深圳-龙岗河沉积物、深圳-污水处理厂沉积物、深圳-东湖公园沉积物，总共65份样品。
步骤2、提取微生物样本的DNA。
具体来说，所有的沉积物的DNA都采用Ultraclean Soil DNA kit (MoBio，USA)从新鲜或深藏冷冻的沉积物样品中提取出来。
步骤3、使用特定的引物进行PCR扩增，同时对每一个样品加上一个序列标签。
具体来说，使用引物967f :CNACGCGAAGAACCTTANC(Seq ID NO :1)禾Π 1406R GACAGCCATGCANCACCT(Seq ID NO :2)去复制微生物群体中细菌的16S V6区片段。由于随后需要对所有微生物进行混合测序，可以对每个样品加一个标签序列，这个序列可以是由8 个碱基组成的经过修正了错误的条形码序列，这个序列被加到引物967f的5’端的前面，在标签序列(条形码序列)和引物967f之间加了一个连接物“GT”。
此外，对于古生菌的V6区域的聚合酶链式反应PCR产物，可以使用引物958AR AATTGGANTCAACGCCGG (Seq ID NO 3)和 1048AR :CGRCGGCCATGCACCWC (kq ID NO :4)，随后采用相同的方式对微生物样品加上条形码序列和“GT”连接物。
步骤4、样品的PCR产物混合并对混合的PCR产物使用优化的So 1 exa建库。
具体来说，对于加上条形码标签的16S V6区的PCR产物，采用分光光度计 Nanodrop进行浓度的定量，然后按照等摩尔的量混合在一起。该具体实施方式
中是将52个细菌V6的PCR产物和13个古生菌的V6的PCR产物共65个样品混合在一起。
把这些混合产物用试剂盒QIAquick PCR purification Kit(Qiagen)进行纯化， 末端修复，在3，端加上碱基A，加上双末端I^air-end测序接头(用试剂盒I^air-end 1 ibrary preparation kit, Illumina)。加完接头后，对样品进行纯化，用30 μ L的EB (洗脱液， Elution buffer)进行溶解。然后取1 μ L的溶液作为模板进行PCR扩增(12CyCleS)。使用试剂盒(QIAquick gel extraction kit, Qiagen)对 PCR 产物进行胶纯化。
步骤5、Solexa测序。具体来说，可以按照厂家11 Iumina的说明书直接用11 Iumina GA II进行测序(7^p pair-end策略，即读长为75个碱基的双末端测序)，如表1所示。
表1样品名-标签序列-Solexa reads
样品名称标签序列备注总 reads 射· 数特异 reads数深圳-北山水库4样品以细菌v6区引物扩增北山水库一 1AACGGCAA76,01930,309北山水库_2AAGGAACC71,05124,714北山水库_3AATTGCGC60,44123,558北山水库4ACAGACTC83,15632,553北山水库_5ACTCAGAC74,63926,703北山水库_6CACTACTC84,39031,436北山水库_7CACTAGAC77,87428,372北山水库_8CAGTGTCA74,47530,909北山水库_9GGAAGCAT38,81715,630北山水库jOGTAGCATC82,99430,159北山水库—11GTCTTGAG66,02925,17权利要求
1 一种对宏基因组16S高可变区V6进行测序聚类分析的方法，其特征在于，所述方法包括提取微生物的脱氧核糖核酸DNA ；通过引物对宏基因组16S核糖体脱氧核糖核酸rDNA的高可变区V6进行聚合酶链式反应PCR，并为每个样品加上标签序列；把不同样品的PCR产物进行混合；对混合后的PCR产物进行Solexa建库法建库；使用Solexa测序工具对所述高可变区V6的文库进行双末端pair-end测序，得到原始的测序数据；对所述测序数据进行筛选，以过滤掉低质量的数据；利用重叠群的关系对所述高可变区V6的全长序列进行组装；通过标签序列把reads分配到对应的样品上；通过对所述reads进行分类分析，以实现使用所述高可变区的测序对微生物群体进行高通量的分类。
2.根据权利要求1所述的方法，其特征在于，所述方法还包括在步骤“提取微生物的脱氧核糖核酸DNA”之前，执行微生物群体的取样。
3.根据权利要求1所述的方法，其特征在于，所述方法还包括在步骤“通过对所述 reads进行分类分析”之后，对不同差异度的序列进行操作分类学单元OTU的分类；根据所述标签序列和reads，进行种群多样性估计Chaol算法和血管紧张素转化酶ACE 的多样性分析。
4.根据权利要求3所述的方法，其特征在于，在进行种群多样性估计Chaol算法和血管紧张素转化酶ACE的多样性分析之后，输出微生物群体的多样性分析图和相对丰度图。
5.根据权利要求1所述的方法，其特征在于，步骤“通过引物对宏基因组16S核糖体脱氧核糖核酸rDNA的高可变区V6进行聚合酶链式反应PCR，并为每个样品加上标签序列”进一步包括使用引物 967f CNACGCGAAGAACCTTANC 和 1406R :GACAGCCATGCANCACCT 去复制微生物群体中细菌的16S高可变区V6区片段；对每个微生物样品加标签序列，所述标签序列被加到所述引物967f的5’端的前面，以及在所述标签序列和所述引物967f之间加上碱基GT。
6.根据权利要求5所述的方法，其特征在于，所述方法还包括对于古生菌的高可变区V6的聚合酶链式反应PCR，使用引物958AR AATTGGANTCAACGCCGG 和 1048AR =CGRCGGCCATGCACCffC0
7.根据权利要求1所述的方法，其特征在于，步骤“把不同样品的PCR产物进行混合” 进一步包括对所述16S的高可变区V6的PCR产物进行浓度定量；以及按照等摩尔的量混合在一起。
8.根据权利要求1所述的方法，其特征在于，步骤“对混合后的PCR产物进行Solexa 建库法建库”进一步包括把混合产物进行纯化，末端修复，在3’端加上碱基A，加上双末端I^ir-end测序接头；加完接头后，对样品进行纯化；对纯化后的样品进行溶解，并作为模板进行聚合酶链式反应PCR扩增；以及对所述聚合酶链式反应PCR产物进行胶纯化。
9.根据权利要求1所述的方法，其特征在于，所述低质量的数据包括与最邻近的引物不匹配的序列、小于50碱基对的序列，或者具有至少一个不同碱基的序列。
10.根据权利要求1所述的方法，其特征在于，步骤“利用重叠群的关系对所述高可变区V6的全长序列进行组装”进一步包括采用所述高可变区V6的PCR产物5’端的前75、70、65、60和55碱基对来进行重叠从而组装；其中，组装的标准是一对序列具有大于5碱基对的重叠长度和在重叠区域小于10% 的不匹配度。
11.根据权利要求1所述的方法，其特征在于，步骤“通过对所述reads进行分类分析” 进一步包括将分配到对应样品上的所述reads比对到现有16s v6数据库中，来达到使用高可变区的标签测序对微生物群体进行高通量的分类分析，进而研究微生物群体的结构。
全文摘要
本发明公开一种对宏基因组16S高可变区V6进行测序聚类分析的方法，该方法包括提取微生物DNA；通过引物对宏基因组16S rDNA的高可变区V6进行PCR，并为每个样品加上标签序列；把不同样品的PCR产物进行混合；对混合后的PCR产物进行Solexa建库法建库；使用Solexa测序工具对高可变区V6的文库进行双末端pair-end测序，得到原始的测序数据；对测序数据进行筛选，以过滤掉低质量的数据；利用重叠群的关系对高可变区V6的全长序列进行组装；通过标签序列把reads分配到对应的样品上；通过对reads进行分类分析，以实现使用高可变区的测序对微生物群体进行高通量和精准的分类。
文档编号C12Q1/68GK102477460SQ201010557119
公开日2012年5月30日 申请日期2010年11月24日 优先权日2010年11月24日
发明者刘晓, 周宏伟, 栗东芳 申请人:深圳华大基因研究院, 深圳华大基因科技有限公司