基于宏基因组 16S 高可变区 V3 的分类方法和装置 【技术领域】
本发明涉及生物信息学分析技术领域， 尤其涉及一种基于宏基因组 16S 高可变区 V3 的分类方法和装置。背景技术
为了研究生物环境中微生物群体的种类， 一般传统的方法包括 ： 直接对微生物 进行培养， 变性梯度凝胶电泳 (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(RiboNucleicAcid)) 的序列高度保守， 可精确指示细菌之间的亲缘关系 ； 16S rRNA 的大小为 1500bp( 碱基对， Base Pair) 左右， 所含信息能反映生物界进化关系， 易操作， 适用于各级分类单元 ； 所以在 宏基因组的研究中， 16S 区测序是最常用的聚类和分类方法。 传统的宏基因组的测序是通过 Sanger 技术测序 16S rRNA gene(16S rDNA) 得到至少 500bp 的读长， 这个读长的长度足够 长， 能够装配出近乎完整的 16S rDNA 序列， 帮助我们去精准地研究每一条序列的物种来源， 但它容易产生嵌合体， 而且测序成本比较高， 费时又费力。
随着新开发出的测序技术以及测序成本的逐步降低， 宏基因组的研究变得越来越 实用， 所涉及的技术包括 Pyrosequencing、 Solexa 等。对于这些革命性的技术的一个主要 挑战就是读长太短， 无法对每个个体的 16S rDNA 进行测序， 因而它的测序信息不足以让我 们去精准地对微生物进行分类。为了解决读长的问题， 有研究 (Bacterial flora-typing with targeted， chip-based Pyrosequencing， BMC Microbiology 2007， 7： 108doi ： 10.1186/1471-2180-7-108，公 开 于 2007 年 11 月 30 日 ) 通 过 Genome Sequencer 20 system(454 Life Sciences) 测序 16S rDNA 可变区来对微生物进行分类， 通过设计特定的 通用引物对 16S 可变区进行特定的 PCR( 聚合酶链式反应， Polymerase Chain Reaction)， 然后用 454 序仪测序， 建立在这种方法上的系统树显示了很好的生物多样性， 但它的测序 成本高， 虽然是传统毛细管测序法费用的 1/10， 但却是其他新一代测序仪测序费用的 10 倍 左右。
综上所述， 提供一种更加准确地对微生物进行聚类分析的方法且方便快捷、 成本 低廉成为本领域亟待解决的技术问题。 发明内容
本发明要解决的一个技术问题是提供一种基于宏基因组 16S 高可变区 V3 的分类方法和装置， 通过对 16S 的高可变区 V3 区进行 solexa 测序， 并通过对这些 16S 可变区的短 序列进行系统分类， 可以在成本低廉的基础上准确反映物种的丰度信息。
本发明的第一方面提供了一种基于宏基因组 16S 高可变区 V3 的分类方法， 该方法 包括 ： 提取微生物样品中的脱氧核糖核酸 (DNA) ； 对提取 DNA 的宏基因组 16S 核糖体脱氧核 糖核酸 (rDNA) 的高可变区 (V3) 进行扩增， 得到作为扩增产物的 DNA 片段 ； 对 DNA 片段进 行 PCR-FreeSolexa 建库， 建库过程中在 DNA 片段上加上标签序列以对每个样品进行标记 ； 将各个样品的带有标签序列的 DNA 片段进行混合， 使用 Solexa 测序工具对混合后的 DNA 片 段进行测序， 得到按照标签区分的测序序列 reads ； 利用 reads 的重叠关系组装得到高可变 区 V3 的全长序列 unique reads ； 对 unique reads 进行分类分析， 以实现对微生物群体的 分类。
优选地， 该方法还包括 ： 在步骤 “提取微生物样品中的脱氧核糖核酸 DNA” 之前， 执 行微生物群体的取样。
优选地， 所述对 unique reads 进行分类分析包括 ： 计算 unique reads 之间的序列 差异度 ； 根据序列差异度执行操作分类学单元 OTU 的分类， 将 unique reads 分配到 OTU 中 ； 将每一个 OTU 分类中的 unique reads 比对到 16S rDNA 的 v3 数据库中， 将比对结果根据众 数原则对 OTU 进行物种注释。 优选地， 根据序列差异度执行操作分类学单元 (OTU) 的分类是指根据本领域公知 的 OTU 分类中 “种” 水平之间的差异度将 unique reads 分配到相应的 OTU 中。在本发明的 一个实施方案中， 将序列差异度在 3％以内的 unique reads 分配到一个 OTU 中。
优选地， 将比对结果根据众数原则对 OTU 进行物种注释是指如果一个 OTU 中 66％ 以上的比对结果均为同一个物种， 则将该 OTU 注释为该物种 ； 如果未达到该比例， 则将物种 分类信息上移一个水平 ( 例如从 “种” 上移到 “属” ， 或从 “属” 继续上移到 “科” ) 再进行统 计， 直到达到 66％的比例标准为止。
优选地， 该方法还包括 ： 在步骤 “对 unique reads 进行分类分析” 之后， 基于分类 分析结果， 进行种群多样性分析和 / 或统计得到微生物群体的相对丰度值。
优选地， 步骤 “对宏基因组 16S rDNA 的高可变区 V3 进行扩增” 是指利用本领域公 知的方法扩增 DNA 序列， 在本发明的一个实施方案中， 采用聚合酶链式反应 (PCR) 扩增 16S rDNA 的高可变区 V3， 所述 PCR 反应的引物为引物 338F ： ACTCCTACGGGAGGCAGCAG 和 533R ： TTACCGCGGCTGCTGGCAC。
优选地， 步骤 “对 DNA 片段进行 PCR-Free Solexa 建库， 建库过程中在 DNA 片段上 加上标签序列， 对每个样品进行标记” 进一步包括 ： 将所述 DNA 片段进行纯化， 对纯化后的 DNA 片段进行浓度定量， 定量后不同样品取等浓度的量分别进行末端修复， 在 3’ 端加上碱 基 A， 然后加上标签序列， 再进一步加上 PCR-Free 的接头， 最后对样品进行纯化。
优选地， 在得到按照标签区分的原始的测序序列 reads 后， 还包括对所述测序序 列进行筛选的步骤， 以过滤掉低质量的测序序列 ； 所述低质量的测序序列选自以下序列中 的任意一种或数种 ： 接头污染序列， 含有多个 poly(A|T|C|G) 的序列、 以及含有连续 2 个以 上的 N 的序列 ；
优选地， 步骤 “利用 reads 的重叠关系组装得到高可变区 V3 的全长序列 unique reads 是指按照本领域公知的条件进行序列的拼接， 例如运用拼接软件， 根据序列两端的重
叠关系对 reads 进行拼接， 将其组装成 V3 的全长序列 unique reads。 在本发明的一个实施 方案中， 拼接的条件是最小匹配长度为 5bp， 重叠区域不允许错配， 重叠区域 N 所占最大百 分比是 0.4％ ； 为了更多的利用序列， 不满足以上结果的序列将各切除 5bp 继续组装， 如此 重复多次， 最终产生的就是 V3 的序列， 如果最终的拼接结果小于 50bp 也不用于后续分析。
本发明的第二方面提供了一种基于宏基因组 16S 高可变区 V3 的分类装置， 所述 装置包括 ： 脱氧核糖核酸 DNA 提取设备， 用于提取微生物样品中的脱氧核糖核酸 DNA ； 扩增 设备， 用于对宏基因组 16S rDNA 的高可变区 V3 进行扩增， 得到作为扩增产物的 DNA 片段 ； Solexa 建库设备， 用于对 DNA 片段进行 PCR-Free Solexa 建库， 建库过程中在 DNA 片段上 加上标签序列， 对每个样品进行标记 ； Solexa 测序设备， 将各个样品的带有标签序列的 DNA 片段进行混合， 使用 Solexa 测序工具对混合后的 DNA 片段进行测序， 得到按照标签区分的 原始的测序序列 reads ； 全长序列组装设备， 用于利用 reads 的重叠关系组装得到高可变区 V3 的全长序列 unique reads ； 分类设备， 用于对 unique reads 进行分类分析， 以实现对微 生物群体的分类。
优选地， 该装置还包括取样设备， 用于执行微生物群体的取样。
优选地， 分类设备包括 ： 序列差异度计算单元， 用于计算 unique reads 之间的序 列差异度 ； OTU 分类单元， 用于根据序列差异度执行操作分类学单元 OTU 的分类， 将 unique reads 分配到 OTU 中 ； 物种注释单元， 用于将每一个 OTU 分类中的 unique reads 比对到 16S rDNA 的 v3 数据库中， 将比对结果根据众数原则对 OTU 进行物种注释。
优选地， 根据序列差异度执行操作分类学单元 (OTU) 的分类是指根据本领域公知 的 OTU 分类中 “种” 水平之间的差异度将 unique reads 分配到相应的 OTU 中。在本发明的 一个实施方案中， 将序列差异度在 3％以内的 unique reads 分配到一个 OTU 中。
优选地， 将比对结果根据众数原则对 OTU 进行物种注释是指如果一个 OTU 中 66％ 以上的比对结果均为同一个物种， 则将该 OTU 注释为该物种 ； 如果未达到该比例， 则将物种 分类信息上移一个水平 ( 例如从 “种” 上移到 “属” ， 或从 “属” 继续上移到 “科” ) 再进行统 计， 直到达到 66％的比例标准为止。
优选地， 还可以包括数据分析设备， 用于在对 unique reads 进行分类分析之后， 对 所得到的数据结果进行进一步分析 ； 所述数据分析设备包括种群多样性分析单元， 用于分 析种群多样性 ； 和 / 或相对丰度统计单元， 用于统计得到微生物群体的相对丰度值。
优 选 地， 对 宏 基 因 组 16S 核 糖 体 脱 氧 核 糖 核 酸 rDNA 的 高 可 变 区 V3 进 行 扩 增是指利用本领域公知的方法扩增 DNA 序列， 在本发明的一个实施方案中， 采用聚合 酶 链 式 反 应 (PCR) 扩 增 16S rDNA 的 高 可 变 区 V3， 所 述 PCR 反 应 的 引 物 为 引 物 338F ： ACTCCTACGGGAGGCAGCAG 和 533R ： TTACCGCGGCTGCTGGCAC。
优选地， Solexa 建库设备包括 ： DNA 片段纯化单元， 用于将所获得的 DNA 片段进行 纯化 ； 定量单元， 用于对纯化后的 DNA 片段进行浓度定量 ； DNA 片段修饰单元， 用于在定量后 不同样品取等浓度的量分别进行末端修复， 在 3’ 端加上碱基 A， 然后加上标签序列， 再进一 步加上 PCR-Free 的接头 ； 修饰后产物纯化单元， 用于对修饰后的 DNA 片段进行纯化。
优选地， 还包括测序序列筛选设备， 用于在得到按照标签区分的原始的测序序列 reads 后， 对所述测序序列进行筛选， 以过滤掉低质量的测序序列 ； 所述低质量的测序序列 选自以下序列中的任意一种或数种 ： 接头污染序列， 含有多个 poly(A|T|C|G) 的序列、 以及含有连续 2 个以上的 N 的序列。
优选地， 所述全长序列组装设备按照本领域公知的条件进行序列的拼接， 例如运 用拼接软件， 根据序列两端的重叠关系对 reads 数据进行拼接， 将其组装成 V3 的全长序列 unique reads。在本发明的一个实施方案中， 拼接的条件是最小匹配长度为 5bp， 重叠区域 不允许错配， 重叠区域 N 所占最大百分比是 0.4％ ； 为了更多的利用序列， 不满足以上结果 的序列将各切除 5bp 继续组装， 如此重复多次， 最终产生的就是 V3 的序列， 如果最终的拼接 结果小于 50bp 也不用于后续分析。
本发明提供的基于宏基因组 16S 高可变区 V3 的分类方法， 采用结合了加标签技术 的 Solexa 技术， 对特定环境下的微生物群体进行了高通量测序， 既减少了人力劳动也节省 了经济花费， 使得在研究微生物群落结构与健康、 环境因子等的关系上变得容易可行。 附图说明 图 1 示出本发明实施例提供的一种基于宏基因组 16S 高可变区 V3 进行分类的方 法的流程图。
图 2 示出对 unique reads 进行分类分析的方法的流程图。
图 3 示出肠道样本扩增结果的稀释曲线分析图。 图 4 示出本发明的基于宏基因组 16S 高可变区 V3 的分类装置的一个实施例的结构图。
图 5 示出本发明的基于宏基因组 16S 高可变区 V3 的分类装置的另一个实施例的 结构图。具体实施方式
现在将参照附图来详细描述本发明的各种示例性实施例。应注意到 ： 除非另外具 体说明， 否则在这些实施例中阐述的部件和步骤的相对布置、 数字表达式和数值不限制本 发明的范围。
同时， 应当明白， 为了便于描述， 附图中所示出的各个部分的尺寸并不是按照实际 的比例关系绘制的。
以下对至少一个示例性实施例的描述实际上仅仅是说明性的， 决不作为对本发明 及其应用或使用的任何限制。
对于相关领域普通技术人员已知的技术、 方法和设备可能不作详细讨论， 但在适 当情况下， 技术、 方法和设备应当被视为授权说明书的一部分。
在这里示出和讨论的所有示例中， 任何具体值应被解释为仅仅是示例性的， 而不 是作为限制。因此， 示例性实施例的其它示例可以具有不同的值。
应注意到 ： 相似的标号和字母在下面的附图中表示类似项， 因此， 一旦某一项在一 个附图中被定义， 则在随后的附图中不需要对其进行进一步讨论。
下面介绍几个本发明技术方案涉及的概念。
宏基因组是指 (Metagenome)( 也称微生物环境基因组 Microbial Environmental Genome， 或元基因组 ) 特定环境中全部微小生物遗传物质的总和。它包含了可培养的和未 可培养的微生物的基因， 目前主要指环境样品中的细菌和真菌的基因组总和。PCR-Free Solexa 建库是指对 PCR 产物进行纯化， 然后进行浓度定量。 不同样品取 等浓度的量分别进行末端修复 ( 即通过酶反应使所有 DNA 双链的粘性末端成为平末端 )， 然后加上碱基 “A” ， 再加上 PCR-Free 的接头， 加完接头后， 对样品进行纯化。所述 PCR-Free 的接头是指测序用的引物序列。
标签序列 (barcode) 是指加在引物 5’ 端前面的一段碱基序列， 用于区分不同的样 本。标签序列可以是由 6 个碱基组成条形码序列， 标签序列的设计要符合一定规则， 比如 碱基含量和不同碱基数目等， 目的是防止因为个别测序错误等原因导致标签相互之间的混 淆， 例如可以参考美国专利申请 US20100267043A1 中公开的方法和原则。
Read(s) 是指 Solexa 测序后产生的测序片段的序列， 也称为测序序列。
Unique read(s) 是指 Read(s) 经过拼接后产生的 v3 的全长序列。 Unique read(s) 的数目是指将序列相同的全长序列合并后得到的全长序列数。
序列差异度是指两条序列比对时不同碱基的个数所占的百分比。
图 1 示出本发明实施例提供的一种基于宏基因组 16S 高可变区 V3 进行分类的方 法的流程图。
如图 1 所示， 基于宏基因组 16S 高可变区 V3 进行分类的方法流程包括 ： 步骤 102， 提取微生物样品中的脱氧核糖核酸 DNA。例如， 采用 Ultraclean Soil DNA kit 试剂盒 (MoBio， USA) 从样品沉积物中提取微生物的 DNA。
步骤 104， 对宏基因组 16S 核糖体脱氧核糖核酸 rDNA 的高可变区 V3( 该区的两端 各有 20 碱基对 bp 左右的保守区， 中间的可变区为 130-160bp 左右 ) 进行扩增， 得到作为扩 增产物的 DNA 片段。例如利用聚合酶链式反应 PCR， 使用引物 338F ： ACTCCTACGGGAGGCAGCAG 和 533R ： TTACCGCGGCTGCTGGCAC 去扩增微生物群体中细菌的 16S 高可变区 V3 区片段。
步 骤 106， 对获 得的 DNA 片段 进行 PCR-Free Solexa 建 库 法建库， 建 库过程 中 在 DNA 片段上加上标签序列， 对每个样品进行标记。例如， 把 PCR 产物用 QIAquick PCR purification Kit(Qiagen) 进行纯化， 用分光光度计对所述 16S 的高可变区 V3 的 PCR 产 物进行浓度定量。不同样品取等浓度的量分别进行末端修复 ( 即通过酶反应使所有 DNA 双链的粘性末端成为平末端 )， 加 “A” ， 加上标签序列， 再加上 PCR-Free 的接头 (Pair-end library preparation kit， Illumina) ； 加完接头后， 对样品进行纯化 ( 用 1.8x ampureXP 磁珠进行纯化 )。有些试剂盒中 PCR-Free 的接头是带有标签序列的接头， 这样可以将加标 签序列和加作为引物的接头一步完成。
步骤 108， 将各个样品的带有标签序列的 DNA 片段进行混合， 使用 Solexa 测序工 具对混合后的 DNA 片段进行测序， 得到按照标签区分的原始的测序序列 reads。例如， 直接 用 Illumina GA II(150bp pair-end 策略 ) 进行测序。Solexa 测序仪 (Illumina genome analyzer) 是新一代的高通量测序仪， 测序价格低廉， 数据读取量大， 相同的测序量的情况 下， Solexa 测序费用是 454 测序费用的十分之一， 而且错误率低 ( 如单碱基测序错误率 -5 ＜ 10 )， 测序无偏性， 对于宏基因组， 可以真实反映物种的丰度信息。而且得到的测序结果 是已根据标签序列进行区分的测序序列 reads。
步骤 110， 利用 reads 的重叠关系组装得到高可变区 V3 的全长序列 unique reads。 例如， 运用拼接软件对 reads 数据进行拼接， 得到拼接的结果。通过序列两端的重叠关系将 两端测序得到的序列组装成 V3 的全长序列 unique reads。拼接的条件是最小匹配长度为
5bp， 重叠区域不允许错配， N 所占最大百分比是 0.4％。为了更多的利用序列， 不满足以上 结果的序列将各切除 5bp 继续组装， 如此重复多次。最终产生的就是 V3 的序列。如果最终 的拼接结果小于 50bp 也不用于后续分析。所述拼接软件可以为 merger、 CABOG、 ARACHNE、 RePS、 phrap 及 newbler 等软件， 在本发明的一个实施例中， 应用了 merger 拼接软件。根据 标签序列即可以把 unique reads 分配到对应的样品上。
步 骤 112， 对 unique reads 进 行 分 类 分 析， 以实现对微生物群体进行高通量 的 分 类。 本 发 明 采 用 结 合 了 标 签 技 术 的 Solexa 技 术， 分 辨 率 大 大 提 高， 单 个 Run 上 Solexa(Illumina) 能产生比 454 多 100 倍的 reads， 因此， 仅仅通过测序这么短的长度就能 得到很好的分类效果， 另外由于结合了标签技术， 能够在单个 Lane(Illumina 高通量测序 仪一张芯片有 8 个通道， 每个通道被称为 “lane” ) 上点更多的样， 大大节约了每个样品的测 序成本。
在本发明的一个实施例中， 对 unique reads 进行分类分析的基本思想为根据 unique reads 的序列差异度将其分配到各个 OTU 中， 再将每一个 OTU 中的序列比对到 16S rDNA 的 v3 数据库中， 得到每一个 OTU 的物种分类。具体流程如图 2 所示 ：
步骤 202， 计算 unique reads 之间的序列差异度。 步骤 204， 根据序列差异度执行操作分类学单元 OTU 的分类， 将 unique reads 分配 到 OTU 中。分配到 OTU 的软件可以例如为 Mothur、 RDP classifier、 qiime 等软件， 在本发 明中采用 Mothur 软件， 其版本为 v.1.6.0， 下载网址为 http://www.mothur.org/wiki/Main Page。
步骤 206， 将每一个 OTU 分类中的 unique reads 比对到 16S rDNA 数据库中。
步骤 208， 将比对结果根据众数原则对 OTU 进行物种注释。
在本发明的一个具体实施例中， 将序列差异度在 3％以内的 unique reads 分配到 一个 OTU 中。
在本发明的一个具体实施例中， 将比对结果根据众数原则对 OTU 进行物种注释是 指如果一个 OTU 中 66％以上的比对结果均为同一个物种， 则将该 OTU 注释为该物种 ； 如果 未达到该比例， 则将物种分类信息上移一个水平 ( 例如从 “种” 上移到 “属” ， 或从 “属” 继续 上移到 “科” ) 再进行统计， 直到达到 66％的比例标准为止。
本发明首先根据序列的差异度进行聚类分析， 将拼接后得到的全长序列按照种的 水平分配到不同的 OTU 中， 再将 OTU 中的序列比对到 16SrDNA v3 数据库中， 在最好匹配的 基础上进行物种的分类， 实现了对复杂微生物样品的准确注释 ； 而且采用大量平行测序能 够发现更多的稀有微生物种类。
另外， 在对 unique reads 进行分类的基础上， 还可以基于分类分析结果， 进行群多 样性分析和 / 或统计得到微生物群体的相对丰度值。
在本发明的一个实施方案中， 利用 Mothur Canoco 软件进行 Chao1 分析， 计算样品 OTU 的 alpha 多样性， 可以得到该环境样品的物种丰富度 ； 同时， 对比对后得到的物种注释 结果进行统计， 可以得到各种微生物的相对丰度。其与 16S 测全长序列在微生物分类和测 量群体的相对丰度上具有等同的技术效果。
接下来详细描述本发明提供的基于宏基因组 16S 高可变区 V3 进行分类的方法的 一个具体实施方式。
步骤 1、 进行微生物群体的取样。 具体来说， 提取深圳北大医院肠道样品共 20 个样本。( 命名以标签序列的数字代码命名 ) 步骤 2、 提取微生物样本的基因组 DNA。
具体来说， 所有样本的 DNA 都采用 Ultraclean Soil DNA kit(MoBio， USA) 从肠道 样品中提取出来。
步骤 3、 使用特定的引物进行 PCR 扩增。
具体来说， 使用引物 338F ： ACTCCTACGGGAGGCAGCAG 和 533R ： TTACCGCGGCTGCTGGCAC 去扩增微生物群体中细菌的 16S 高可变区 V3 区片段。
步骤 4、 对 PCR 产物进行 PCR-Free Solexa 建库法建库。具体来说， 把 PCR 产物用 QIAquick PCR purification Kit(Qiagen) 进行纯化， 用分光光度计对所述 16S 的高可变区 V3 的 PCR 产物进行浓度定量。 20 个样本取等浓度的量分别进行末端修复 ( 即通过酶反应使 所有 DNA 双链的粘性末端成为平末端 )， 接着 3′端加 “A” ， 然后 3′端加上标签序列， 再在 3′端加上 PCR-Free 的接头 (Pair-end library preparation kit， Illumina) ； 加完接头 后， 对样品进行纯化 ( 用 1.8x ampureXP 磁珠进行纯化 )。用 12μL 的 EB( 洗脱液， Elution buffer) 进行溶解。
步骤 5、 Solexa 测序。 具体来说， 可以按照厂家 Illumina 的说明书直接用 Illumina GA II 进行测序 (150bp pair-end 策略， 即读长为 150 个碱基的双末端测序 )， 得到原始的 测序序列 reads。而且这些 reads 已根据标签序列进行了区分。如表 1 所示。
表 1 Solexa 原始测序序列 reads 和 unique reads 数
步骤 6、 得到原始的测序数据后， 过滤掉那些低质量的数据。具体来说， 低质量的 数据是指接头污染序列， 含有多个 poly(A|T|C|G) 的序列、 以及含有连续 2 个以上的 N 的序 列； 如表 2 所示。接头污染序列的判断标准为 ： 若 reads 与接头序列可以连续比对上 15bp 的长度， 则认为该 reads 有接头污染。含有多个 poly(A|T|C|G) 的序列， 即低复杂度序列的 判断标准为 ： 若 reads 中 poly(A|T|C|G) 序列的长度≥ 10bp， 则为低复杂度序列的 reads。
表 2 宏基因组数据
11CN 102517392 A Item 项目 原始测序数据 过滤掉低质量数据后的序列 序列组装条数 全长序列 (Unique read)
说明书产出数据 (Mb) 3000 2819 2184 2269/17 页Reads 数量 20,000,000 18,799,024 12,864,308 1,673,626其中产出数据＝ reads 数 ×150bp。序列组装条数是满足序列拼接条件的序列统 计得到的个数。
步骤 7、 利用重叠关系来组装 V3 的全长序列。
具体来说， 通过所述 Pair-end 的 reads 重叠区来组装高可变区 V3 的序列。运 用 merger 拼接软件 (http://emboss.sourceforge.net/apps/release/6.2/emboss/apps/ merger.html) 对 reads 数据进行拼接， 得到拼接的结果， 即通过重叠关系将两条两端测序 得到的序列组装成一条序列。拼接的条件是最小匹配长度为 5bp， 重叠区域不允许错配， 重 叠区域 N 所占最大百分比是 0.4％。为了更多的利用序列， 不满足以上结果的序列将各切 除 5bp 继续组装， 如此重复多次。最终产生的就是 V3 的全长序列。如果最终的拼接结果小 于 50bp 也不用于后续分析。利用 Mothur 所带 unique 程序去冗余， 统计就能得到 unique reads 的个数。所述去冗余是指将序列相同的 unique reads 合并为一条。通过条形码标签 序列把 unique reads 都对应到相应的样品上。
步骤 8、 执行 OTU(operational taxonomic unit， 操作分类学单元 ) 的分类。使用 SLP 软件 (Huse SM， D Mark Welch et al.(2010).Ironing out the wrinkles in the rare biosphere through improved OTU clustering.Environ Microbiol 12 ： 1889-1898.) 对 样品的 unique reads 序列进行差异度计算， 根据 unique reads 之间的差异， 利用 Mothur 软 件 ( 软 件 mothur(v.1.6.0) 的 下 载 途 径 为 http://www.mothur.org/wiki/Main_Page) 对 unique reads 进行 OTU 分类。差异度在 3％以内 ( 相当于分类水平中 “种” 的水平 ) 的 unique reads 被聚为一个 OTU。
步骤 9、 对各 OTU 进行物种注释。使用 BLASTN 将 OTU 中的 Tags 序列比对到 16S rDNA 数据库 refhvr_V3(http://vamps.mbl.edu/resources/databases.php) 中。参数设 定为 ： -F F-e 1e-5。将得到的 Blast 结果根据众数原则， 对 OUT 进行注释， 即如果一个 OTU 中 66 ％ ( 三分之二 ) 的比对结果均为同一个物种， 则将该 OTU 注释为该物种 (Huse SM， Dethlefsen L， Huber JA， Mark Welch D， Relman DA， Sogin ML.Exploring microbial diversity and taxonomy using SSU rRNA hypervariable tag sequencing.PLoS Genet.2008Nov ； 4(11) ： e1000255.Epub 2008 Nov 21)。如果未达到该比例， 则将物种分类 信息上移一个水平 ( 从 “种” 上移到 “属” ， 或从 “属” 继续上移到 “科” ) 再进行统计， 直到达 到 66％的比例标准为止。
步骤 10、 数据分析。
具体来说， 基于步骤 8 的 OTU 分类的结果， 采用 Mothur. 软件进行 Chao1， 即群落内 的多样性分析， 计算样品 OTU 的 alpha 多样性。 Alpha 多样性是一个环境中物种的多样性分析的结果。Shannon， npshannon 模型计算的结果越大表示该环境的物种越丰富。Simpson 模型的结果越接近于 0 表示物种越丰富， 结果如表 3 所示。同时对步骤 9 中得到的 OTU 注 释结果进行统计， 得到表 4 所示的微生物群体的相对丰度值。并根据各样品的 OTU 数量及 其序列总数量绘制其稀释 (Rarefaction) 曲线， 见图 3。
表 3 微生物 Alpha 多样性指数
样品名称 肠道样品 57 肠道样品 63 肠道样品 93 肠道样品 101 肠道样品 115 肠道样品 123 肠道样品 100 肠道样品 103 肠道样品 124 肠道样品 137 肠道样品 153 肠道样品 35 肠道样品 74 肠道样品 158 肠道样品 89 肠道样品 92 肠道样品 97 肠道样品 1Shannon 5.62 5.70 5.86 5.93 5.93 5.97 5.64 5.76 5.95 5.85 5.88 5.30 5.40 5.47 5.63 5.38 5.57 5.63Npshannon 5.69 5.78 5.95 6.01 6.02 6.06 5.72 5.84 6.04 5.93 5.96 5.36 5.46 5.53 5.70 5.44 5.63 5.71Simpson 0.031 0.030 0.028 0.026 0.027 0.027 0.030 0.029 0.026 0.028 0.027 0.034 0.030 0.030 0.029 0.033 0.030 0.02913CN 102517392 A 肠道样品 2 肠道样品 3
物种 Acidaminococcus Akkermansia Alistipes 样品 57 7436 1627 10404 样品 63 7762 1354 10230说5.37 5.48明书0.035 0.03111/17 页5.44 5.54表 4 肠道样品扩增物种 reads 丰度值样品 93 6292 1192 9615 样品 101 7368 1648 10387 样品 115 6739 1350 8985 样品 123 6121 1302 8827 样品 100 7354 1339 10449
14Bacteroides 2154 1001 821 96475 811 81 1588 7347 2 19679 5159 2453 1409 样品 103 样品 124 1227 2653 2885 1135 样品 137 4716 4280 19787 18416 20092 4585 2524 1181 样品 153 4 2 1 6766 6307 6912 1690 1410 1464 1499 6045 1 18918 4234 2561 1013 样品 1 75 60 54 47 844 753 774 729 767 39 1433 5725 4 18414 3793 2624 1038 样品 2 111236 120329 105584 111917 119035 1076 993 710 711 614 924 761 817 793 718 954 1185 114044 914 56 1952 2277 2123 1877 2041 2292340196324361313269325480332439324911328043BifidobacteriumCN 102517392 ABilophilaEscherichiaFaecalibacteriumOdoribacterOxalobacter说Papillibacter1481 7237 2 20440 4691 2848 1233 样品 3明15Parabacteroides书PeptostreptococcusRoseburiaRuminococcusSubdoligranulumSutterella物种12/17 页
6289 1207 9807 316047 2462 888 733 128984 830 65 1636 6198 3 20193 20217 5 1 18114 6768 5507 1716 1495 50 55 46 1315 5840 6 19585 870 688 619 108252 115650 127634 958 633 644 1529 96158 2180 165 1129 3273 1 31410 785 740 744 334 2086 1835 2141 1831 321657 342553 335685 269333 279904 2338 201 1529 115933 1908 232 1097 2477 3 35069 10159 10003 7448 10773 10286 1682 1377 1430 4262 4564 4256 13287 251964 2878 185 7098 6819 6211 21338 21984 26085 1819 123340 1956 115 1869 1998 4 38978AcidaminococcusCN 102517392 AAkkermansiaAlistipesBacteroidesBifidobacteriumBilophila说Escherichia明16Faecalibacterium书OdoribacterOxalobacterPapillibacterParabacteroidesPeptostreptococcusRoseburia13/17 页Ruminococcus 3056 1178 样品 35 15255 3061 10901 259635 2036 218 1310 104589 1753 47 1365 1288 146 124 1237 1993 1808 117416 105177 1153 1279 1515 98727 2062 197 1267 348 369 419 2413 2180 2023 232951 229491 262062 265752 2344 243 1636 121837 1872 225 1156 8762 8864 10872 9834 4355 3375 4595 4777 17802 16315 22112 20783 26141 3946 13672 样品 74 样品 158 样品 89 样品 92 样品 97 1102 1063 975 1544 1806 1690 2642 2532 2802 3445 3844 435844314839421940841245603860SubdoligranulumCN 102517392 ASutterella物种AcidaminococcusAkkermansiaAlistipes说Bacteroides257709 2750 206 1696 116235 71 1906 122明17Bifidobacterium书BilophilaEscherichiaFaecalibacteriumOdoribacterOxalobacterPapillibacter14/17 页
3150 2 32633 845 3712 1114 1454 1287 1659 4562 3874 3573 1127 1284 1288 36911 32773 32183 37113 681 4304 1786 4 2 3 4 2627 2882 3446 2474 1764 2052 37129 778 4092 1699CN 102517392 AParabacteroidesPeptostreptococcus说Roseburia明从表 3、 表 4 和图 3 中可看出， 拟杆菌属 (Bacteroides) 是占绝对优势的物种， 占总 tags 数的 60％左右， 而对突柄杆菌属 (Akkermansia) 扩增的数量都很少。对于柔嫩梭菌属18Ruminococcus书SubdoligranulumSutterella15/17 页CN 102517392 A说明书16/17 页(Faecalibacterium) 的扩增效果， 肠道样品几乎相同， 该菌是肠道中的有益菌群之一， 据报 道其分泌产物能够帮助人体抵御炎症的侵袭。汇总各样品的物种分类信息， 能在属的水平 上进行区分的共有 158 个属。将不同样品间能分到属的物种做比较分析， 找出具有显著差 异的物种。绝大多数的物种在不同样品间都没有显著性差异。扩增结果的差异主要表现在 优势类群上， 特别是拟杆菌的五个属。 其中有十一个样品在三十个属上有显著差异， 但在属 于拟杆菌门的 Odoribacter、 属于厚壁菌门的 Anaerovorax、 Roseburia、 Subdoligranulum 等十个样品间没有差异。
图 4 示出本发明的基于宏基因组 16S 高可变区 V3 的分类装置的一个实施例的结 构图。如图 14 所示， 该装置包括 ： 脱氧核糖核酸 DNA 提取设备 402， 用于提取微生物样品中 的脱氧核糖核酸 DNA ； 扩增设备 404， 用于对宏基因组 16S rDNA 的高可变区 V3 进行扩增， 得 到作为扩增产物的 DNA 片段 ； Solexa 建库设备 406， 用于对 DNA 片段进行 PCR-FreeSolexa 建 库， 建库过程中在 DNA 片段上加上标签序列以对每个样品进行标记 ； Solexa 测序设备 408， 将各个样品的带有标签序列的 DNA 片段进行混合， 使用 Solexa 测序工具对混合后的 DNA 片 段进行测序， 得到按照标签区分的原始的测序序列 reads ； 全长序列组装设备 410， 用于利 用 reads 的重叠关系组装得到高可变区 V3 的全长序列 unique reads ； 分类设备 412， 用于 对 unique reads 进行分类分析， 以实现对微生物群体的分类。
图 5 示出了本发明的基于宏基因组 16S 高可变区 V3 的分类装置的另一个实施 例的结构图。与图 4 相比， 该装置还包括测序序列筛选设备 509， 用于在得到按照标签 区分的原始的测序序列 reads 后， 对所述测序序列进行筛选， 以过滤掉低质量的测序序 列； 所述低质量的测序序列选自以下序列中的任意一种或数种 ： 接头污染序列， 含有多个 poly(A|T|C|G) 的序列、 以及含有连续 2 个以上的 N 的序列 ； 还包括数据分析设备 514， 用于 在对 unique reads 进行分类分析之后， 对所得到的数据结果进行种群多样性分析和 / 或统 计得到微生物群体的相对丰度值。
在本发明的一个实施例中， 分类设备 412 包括 ： 序列差异度计算单元 4121， 用于计 算 unique reads 之间的序列差异度 ； OTU 分类单元 4122， 用于根据序列差异度执行操作分 类学单元 OTU 的分类， 将 unique reads 分配到 OTU 中 ； 物种注释单元 4123， 用于将每一个 OTU 分类中的 unique reads 比对到 16S rDNA 的 v3 数据库中， 将比对结果根据众数原则对 OTU 进行物种注释。
对于图 4、 5 中各个装置或单元的功能， 可以参考上文中关于本发明方法的实施例 中对应部分的说明， 为简洁起见， 在此不再详述。
本领域的技术人员应当理解， 对于图 4、 5 中的各个装置， 可以通过单独的计算处 理设备实现， 或者将其集成为一个独立的设备实现。在图 4、 5 中用框示出以说明它们的功 能。这些功能块可以用硬件、 软件、 固件、 中间件、 微代码、 硬件描述语音或者它们的任意组 合来实现。举例来说， 一个或者两个功能块都可以利用运行在微处理器、 数字信号处理器 (DSP) 或任何其他适当计算设备上的代码实现。代码可以表示过程、 功能、 子程序、 程序、 例 行程序、 子例行程序、 模块或者指令、 数据结构或程序语句的任意组合。代码可以位于计算 机可读介质中。计算机可读介质可以包括一个或者多个存储设备， 例如， 包括 RAM 存储器、 闪存存储器、 ROM 存储器、 EPROM 存储器、 EEPROM 存储器、 寄存器、 硬盘、 移动硬盘、 CD-ROM 或 本领域公知的其他任何形式的存储介质。计算机可读介质还可以包括编码数据信号的载波。 本发明提供的基于宏基因组 16S 高可变区 V3 进行分类的方法和装置， 采用结合了 加标签技术的 Solexa 技术， 对特定环境下的微生物群体进行了高通量测序， 在单个的 Lane 中我们测了来自 20 个样品的大约 2 千万条 16S rRNA V3 标签序列。由此可见 Solexa 测序 16S rRNA v3 可变区对微生物群体进行分类的方法是经济节约型的， 既减少了人力劳动也 节省了经济花费， 使得在研究微生物群落结构与健康， 环境因子等等的关系上变得容易可 行。此外， 不管是总共的 reads 还是 0 个错配的 reads 数量都比以前报道的测序 16S 标签 序列的数量高。
参考前述本发明示例性的描述， 本领域技术人员可以清楚的知晓本发明具有以下 优点 ：
本发明提供的基于宏基因组 16S 高可变区 V3 进行分类的方法和装置， 仅仅用高变 区 v3 测序来对样品中的微生物进行分类， 这种方法在分类和测量微生物群体的相对丰度 上显示了很好的效果， 甚至在可变区 V3 区序列与它们最近的参考序列有一定差异的情况 下同样能够达到很好的效果。结果显示通过测 V3 可变区进行微生物种类的分析， 不仅可以 测到主要的那些微生物， 还可以测到更多的稀有微生物。并发现微生物的多样性不仅仅局 限于以前按照表型来分类的伯吉分类法， 而且微生物群体也远比想象中的复杂。 此外， 在开 发微生物群体的多样性和相对丰度上， 大量平行 Solexa 测序 V3 可变区序列有超越其他技 术很多的优势， 比如微生物多样性的相对水平， 序列的长度， 同聚物的密度， 能够识别到种 水平的能力， 或适应不同扩增引物的优点。
另外， 由于 Solexa 测序的序列比较短， 同样一个 run， 它能提供更多样品的 reads， 识别更多的微生物， 比传统的全长 SSU rRNA 测序每个 read 花费更少。随着技术的进步， 产 生更多的 Reads 数据和更长的序列的 Solexa 测序将会给可变区测序分类微生物提供更广 阔的机会， 比如长测序， 可变区域的应用， 各种各样可变区的结合， 或者更深的测序深度。 可 变区标签测序的最大优势在于它应用了大量平行 Solexa 测序的优势， 比原先所达到的大 好几个数量级的测序深度和宽度， 促进了微生物群体和稀有生物圈广阔多样性的开发。
至此， 已经详细描述了根据本发明的基于宏基因组 16S 高可变区 V3 进行分类的方 法和装置。为了避免遮蔽本发明的构思， 没有描述本领域所公知的一些细节。本领域技术 人员根据上面的描述， 完全可以明白如何实施这里公开的技术方案。
虽然已经通过示例对本发明的一些特定实施例进行了详细说明， 但是本领域的技 术人员应该理解， 以上示例仅是为了进行说明， 而不是为了限制本发明的范围。 本领域的技 术人员应该理解， 可在不脱离本发明的范围和精神的情况下， 对以上实施例进行修改。 本发 明的范围由所附权利要求来限定。