CN104039982B - 一种分析微生物群落组成的方法和装置 - Google Patents

一种分析微生物群落组成的方法和装置 Download PDF

Info

Publication number
CN104039982B
CN104039982B CN201280064063.2A CN201280064063A CN104039982B CN 104039982 B CN104039982 B CN 104039982B CN 201280064063 A CN201280064063 A CN 201280064063A CN 104039982 B CN104039982 B CN 104039982B
Authority
CN
China
Prior art keywords
sequencing
stack
fragments
module
reference set
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201280064063.2A
Other languages
English (en)
Other versions
CN104039982A (zh
Inventor
李胜辉
冯强
覃俊杰
朱剑锋
官远林
王俊
汪建
杨焕明
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
BGI Shenzhen Co Ltd
Original Assignee
BGI Shenzhen Co Ltd
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by BGI Shenzhen Co Ltd filed Critical BGI Shenzhen Co Ltd
Publication of CN104039982A publication Critical patent/CN104039982A/zh
Application granted granted Critical
Publication of CN104039982B publication Critical patent/CN104039982B/zh
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Classifications

    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6869—Methods for sequencing
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/20—Sequence assembly
    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6809—Methods for determination or identification of nucleic acids involving differential detection
    • C—CHEMISTRY; METALLURGY
    • C12—BIOCHEMISTRY; BEER; SPIRITS; WINE; VINEGAR; MICROBIOLOGY; ENZYMOLOGY; MUTATION OR GENETIC ENGINEERING
    • C12Q—MEASURING OR TESTING PROCESSES INVOLVING ENZYMES, NUCLEIC ACIDS OR MICROORGANISMS; COMPOSITIONS OR TEST PAPERS THEREFOR; PROCESSES OF PREPARING SUCH COMPOSITIONS; CONDITION-RESPONSIVE CONTROL IN MICROBIOLOGICAL OR ENZYMOLOGICAL PROCESSES
    • C12Q1/00—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions
    • C12Q1/68—Measuring or testing processes involving enzymes, nucleic acids or microorganisms; Compositions therefor; Processes of preparing such compositions involving nucleic acids
    • C12Q1/6876—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes
    • C12Q1/6888—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for detection or identification of organisms
    • C12Q1/689—Nucleic acid products used in the analysis of nucleic acids, e.g. primers or probes for detection or identification of organisms for bacteria
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B15/00—ICT specially adapted for analysing two-dimensional or three-dimensional molecular structures, e.g. structural or functional relations or structure alignment
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B20/00—ICT specially adapted for functional genomics or proteomics, e.g. genotype-phenotype associations
    • G16B20/20—Allele or variant detection, e.g. single nucleotide polymorphism [SNP] detection
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B30/00—ICT specially adapted for sequence analysis involving nucleotides or amino acids
    • G16B30/10—Sequence alignment; Homology search

Landscapes

  • Life Sciences & Earth Sciences (AREA)
  • Chemical & Material Sciences (AREA)
  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Proteomics, Peptides & Aminoacids (AREA)
  • Analytical Chemistry (AREA)
  • General Health & Medical Sciences (AREA)
  • Biotechnology (AREA)
  • Biophysics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Evolutionary Biology (AREA)
  • Medical Informatics (AREA)
  • Theoretical Computer Science (AREA)
  • Organic Chemistry (AREA)
  • Molecular Biology (AREA)
  • Genetics & Genomics (AREA)
  • Wood Science & Technology (AREA)
  • Zoology (AREA)
  • Immunology (AREA)
  • Microbiology (AREA)
  • Biochemistry (AREA)
  • General Engineering & Computer Science (AREA)
  • Crystallography & Structural Chemistry (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本发明提供了一种用于分析环境样品中的微生物群落组成的方法和装置。该方法包括测序、初级组装、分栈、基于栈的高级组装和鉴定等步骤。

Description

一种分析微生物群落组成的方法和装置
技术领域
本发明涉及宏基因组学(metagenomics)和生物信息学领域。特别地，本发明涉及用于分析环境样品中的微生物群落组成的方法和装置。
背景技术
宏基因组学又称为环境基因组学，元基因组学，生态基因组学，或者群落基因组学，其是一门直接研究各种环境(例如自然环境)中的微生物群落(包含了可培养的和不可培养的细菌、真菌和病毒等的总和)的学科。研究各种环境中的微生物群落和物种多样性具有特别的益处。例如，对人肠道环境中的微生物群落和物种多样性的研究，对于菌群的临床药物开发以及人菌代谢途径的了解是非常有用的。然而，由于传统研究方法的限制，我们对环境(例如肠道环境)中的微生物组成知之甚少。特别地，由于环境中可能包含了不可培养的细菌、真菌或病毒，许多物种无法通过涉及培养的传统研究方法进行鉴定。
如今，全基因组鸟枪法(whole genome shotgun，WGS)已经在宏基因组学研究中逐渐兴起。这种方法一般通过高通量测序得到大量的测序片段(reads)，然后通过组装得到较大的连接片段(contigs)、拼接片段(scaffolds)、或甚至是全基因组。与此同时，新一代高通量测序技术得到了长足的发展，这为利用WGS策略来认识群落结构、研究群落差异和功能提供了良好契机。例如，最近的宏基因组学研究已经在多种环境中，在发现新物种、解析微生物群落多样性与交互关系方面取得了初步的成果：关于海洋环境，参见例如Venter et al.2004；关于石油环境，参见例如Daniel,2005；关于人体环境，参见例如Gillet al.2006。
然而，当利用宏基因组学研究(例如，WGS策略)来分析环境样品中的微生物群落组成时，仍然存在两个巨大的挑战，即，大量的短基因片段(例如，测序片段)的组装以及不同物种的识别。由于宏基因组学研究所收集到的是，一个特定环境中的所有物种的基因信息，因此，如何将这些大量的、混合的短基因片段组装成连接片段或拼接片段，是一个巨大的难题与挑战。同时，在组装得到较长的连接片段或拼接片段后，如何判别这些长片段的物种来源，又是一个巨大的难题与挑战。
目前，已经开发出一些程序，用于组装混合的短基因片段，例如Velvet(Zerbino and Birney2008),EULER-SR(Chaisson and Pevzner2008),Newbler(Mergulies et al.2006)和Soapdenovo(Li et al.2009)。此外，分栈(binning)方法已被广泛用于判别连接片段或拼接片段的归属物种，其包括但不限于，基于相似度(similarity-based)的MEGAN(Husonet al.2007)和CARMA(Tzahoret al.2009)，这类方法是通过与参考基因组进行序列比对来进行片段分类的；基于组成特征(composition-based)的分栈方法，例如基于GC含量、k-mer频率(Schbath et al.1995)或者四核苷酸频率(Teeling et al.2004)等的分栈方法，这类方法在很大程度上受限于片段长度和序列特征的辨别能力；以及，基于片段丰度(abundance-based)的AbundanceBin(Wu and Ye2011)，这类方法根据环境中不同物种的丰度进行片段分类，只适合用于短的测序片段。
然而，宏基因组学的研究目的是重建环境样品中的各种微生物的基因组，以分析环境样品中的微生物群落组成。上述方法把组装和分栈分开，而各自只关注于一个方面。因此，上述方法并不能充分达到宏基因组学的研究目的。另外，即使将上述的组装方法和分栈方法简单组合在一起，由于不同方法所采用的算法、步骤、兼容性并不一定匹配，其最终结果是否能达到宏基因组学的研究目的，以及最终结果的精确度和有效性也都难以预料。
因此，本领域仍然需要一种高效率、高精度的分析环境样品中的微生物群落组成的方法。
发明内容
在本发明中，除非另有说明，否则本文中使用的科学和技术名词具有本领域技术人员所通常理解的含义。并且，本文中所使用的各种实验室操作步骤均为相应领域内广泛使用的常规步骤。同时，为了更好地理解本发明，下面提供相关术语的定义和解释。
如本文中所使用的，术语“环境”是指广义上的各种环境，其包括但不限于，自然环境(例如土壤环境，海洋环境，河流环境)和体内环境(例如口腔环境，肠道环境)。更确切而言，术语“环境”是指，可能存在微生物/微生物群落的任何区域。
如本文中所使用的，术语“环境样品”是指，来自各种环境的可能含有微生物/微生物群落的样品。
如本文中所使用的，术语“微生物”具有本领域技术人员所通常理解的含义，其包括但不限于，细菌、真菌和病毒。
如本文中所使用的，术语“微生物群落”是指，在特定环境中生活在一起的各类微生物的联合。通常，同一微生物群落中的各种微生物不仅相互之间具有直接或间接的交互关系，而且与它们所生存的环境之间也具有交互关系：环境的改变会导致微生物群落的组成(包括，微生物的种类和/或丰度)发生改变；反过来，微生物群落组成的改变也影响着环境。
如本文中所使用的，术语“宏基因组”是指，群落中的各种生物的基因组的总和。特别地，在本发明的方法和装置的背景下，术语“宏基因组”是指，微生物群落中的各种微生物的基因组的总和。相应地，术语“宏基因组测序数据”是指，对整个宏基因组进行测序所获得的数据。由于宏基因组所包含的DNA信息十分巨大，因此，通常使用高通量测序技术(例如，第二代测序技术或第三代测序技术)进行测序。然而，也可以通过其他方法或其他来源，获得所需的宏基因组测序数据。测序数据通常由大量的测序片段(read)构成。
第二代测序技术是本领域技术人员公知的，其包括例如，454测序法(Roche)，Solexa测序法(Illumina)，SOLiD测序法(ABI)和单分子测序法。关于第二代测序技术的详细综述，参见例如，MichaelMetzker(2010),Sequencing technologies-the next generation,NatureGenetics。关于第三代测序技术，可参见例如，Eric E.Schadt等人，A window into third-generation sequencing，Human MolecularGenetics,2010,Vol.19,Review Issue2,R227-240。
表述“测序质量低的序列”的含义是本领域技术人员已知的，其例如可在测序过程中由测序平台和测序软件确定(参见，Quality Scoresfor Next-Generation Sequencing,Technical Note:Sequencing,Illumina)。
如本文中所使用的，表述“去冗余”是指，对于彼此相似度达到95％或以上的序列，只保留一个，例如，将重复的连接片段和拼接片段去除。
如本文中所使用的，表述“参考集”是广义上的组装片段集或基因集，其中，组装片段是指由测序片段组装得到的长片段，例如连接片段(contigs)、拼接片段(scaffolds)；基因集是指在组装片段上预测出来的基因的集合。所述组装片段或基因构成且被称为参考集的“元素”。
如本文中所使用的，术语“分栈(binning)”和“聚类”具有相同的含义，“栈(bin)”和“类”具有相同的含义。它们可互换使用。
如本文中所使用的，术语“多元正态分布模型”和“最大似然函数法”具有本领域技术人员所通常理解的含义。关于这2个术语的详细描述，可参见例如Fraley and Raftery,1998。
如本文中所使用的，术语“基于相似度的聚类方法”是指，通过比较两两序列之间的序列同一性来度量序列之间的相似度(或距离)，并基于这个相似度(或距离)进行聚类；术语“基于组成特征的聚类方法”是指，通过比较两个序列自身组成特征的相似性，如寡核苷酸频率，GC含量等，来度量序列之间的相似度(或距离)，并基于这个相似度(或距离)进行聚类。基于相似度的聚类方法例如但不限于，基于相似度(similarity-based)的MEGAN(Husonet al.2007)和CARMA(Tzahoret al.2009)。基于组成特征的聚类方法例如但不限于，基于GC含量、k-mer频率(Schbath et al.1995)或者四核苷酸频率(Teelinget al.2004)的聚类方法。
本发明所要解决的一个技术问题是，提供一种能有效分析环境样品中的微生物群落组成的方法和装置。基于此，发明人创造性地将组装方法和分栈方法结合在一起，开发了能够高效率且高精度地分析从环境样品获得的宏基因组数据，并进而确定环境样品的微生物群落组成的方法和装置。特别地，本发明的方法也被命名为Soap series ofMetagenome analysis(在下文中简称为SoapMeta)。
因此，在一个方面，本发明提供了一种用于分析环境样品中的微生物群落组成的方法，其包括以下步骤：
1)测序：
对来自环境样品的基因组DNA进行构建文库和测序，从而获得由测序片段池(reads pool)构成的宏基因组测序数据；
2)初级组装：
2a)构建或完善参考集：对测序片段进行组装以得到组装片段，然后去冗余，从而构建非冗余的参考集(即，组装片段集)；任选地，可在所获得的组装片段上预测基因，并将预测出来的基因的集合作为参考集(即，基因集)；或者，如果针对所述环境样品，存在已知的参考集，那么直接将它作为参考集，或者将该已知的参考集与如上所述构建的参考集组合并去冗余，从而获得最终的参考集
2b)构建元素相对丰度谱矩阵：将所述测序片段与参考集进行比对，并计算参考集中的各个元素在样品中的相对丰度；
3)分栈，即，通过下述步骤确定参考集中的每一个元素所归属的栈，得到聚类的栈：
3a)基于丰度的分栈：基于元素在样品中的相对丰度，使用聚类算法，如自底向上的层次聚类方法(HIERARCHICALCLUSTERING SCHEMES，STEPHEN C.JOHNSON，1967)，确定各个元素的初始栈；和
3b)基于模型的分栈：
(i)将每一个初始栈作为一个独立的多元正态分布模型,并基于丰度矩阵，利用最大似然函数法计算所述模型的参数；
(ii)构建一个软矩阵(fuzzy matrix)，用于存储每一个元素归属某一个栈的概率；和
(iii)迭代运算E步和M步，直至似然函数达到最大化：
E步，根据每一个栈的模型参数，分别计算每一个元素属于某一个栈的后验概率，并且修改软矩阵中所述元素属于所述栈的概率；
M步：根据软矩阵，用最大似然函数法计算每一个栈的模型参数；
4)基于栈的高级组装：
4a)通过将测序片段与已分栈的元素进行比对，从宏基因组测序数据中寻找对应到之前确定的各个栈的测序片段；
4b)使用SOAPdenovo或者使用其他针对微生物测序数据的组装软件，将对应到各个栈的测序片段分别进行组装；
4c)使用基于相似度的聚类方法和/或基于组成特征的聚类方法，对每一个栈所包含的元素的分栈进行校正；任选地，还在已获得的栈内部进行再次聚类，然后根据聚类的结果，对已获得的栈进行拆分或保持不变，从而使结果更加准确可信；
4d)重复步骤4a)－4c)，直到各个栈的基因组序列的大小无明显变化为止(总长度增长率小于5％)；
5)鉴定：
利用各个栈的基因组序列，确定各个栈所对应的微生物的类别，从而确定所述环境样品中的微生物群落组成。
关于测序
在一个优选的实施方案中，环境样品来源于自然环境，例如土壤环境，海洋环境和河流环境。在另一个优选的实施方案中，环境样品来源于体内环境，例如口腔环境和肠道环境。
在一个优选的实施方案中，使用第二代测序技术(例如，454测序法，Solexa测序法，SOLiD测序法或单分子测序法)或第三代测序技术对环境样品所包含的微生物群落的宏基因组进行测序，从而提供来自环境样品的宏基因组测序数据。
在一个优选的实施方案中，通过下列步骤来获得宏基因组测序数据：
1a)提供环境样品；
1b)从所述环境样品中提取宏基因组DNA；
1c)利用所述宏基因组DNA构建宏基因组文库；
1d)对所述宏基因组文库进行测序，优选使用Solexa测序法进行测序，从而提供所述环境样品的宏基因组测序数据。
在一个优选的实施方案中，宏基因组测序数据是由测序片段构成的测序片段池(reads pool)。此类测序片段通常通过第二代测序技术(例如Solexa测序法)或第三代测序技术获得。
在一个优选的实施方案中，测序片段是末端配对的测序片段(paired end reads)。
测序片段中可能包含测序过程中所使用的接头(adapter)的序列，测序质量低的序列和/或在分析来自体内环境的样品的情况下，来自宿主基因组的序列。此类序列可能会影响后续的处理和分析，因此，此类序列的去除可能是有利的。
因此，在一个优选的实施方案中，在进行步骤2)之前，对测序数据进行预处理，即，去除接头序列、测序质量低的序列和/或宿主基因组序列。
在一个优选的实施方案中，对来自相同或相似环境的多个样品进行测序，并将所有样品的测序数据组合在一起，构成宏基因组测序数据。
在一个优选的实施方案中，宏基因组的测序深度为至少10×，优选至少20×，优选至少30×，优选至少40×，更优选至少50×。
关于初级组装
在一个优选的实施方案中，使用Soapdenovo将所述测序片段组装成组装片段(例如，连接片段和/或拼接片段)。此类组装方法是本领域技术人员已知的，参见例如，Li et al.2009。
在一个优选的实施方案中，使用多个环境样品来进行本发明的方法，并且针对每个样品分别获得了各自的参考集。在此情况下，将所有样品的参考集组合在一起，并去冗余，从而构建最终的非冗余的参考集。也即，将来自多个样品的参考集组合在一起，并去冗余，从而构建最终的非冗余的参考集。
在一个优选的实施方案中，如果针对所述环境样品，存在已知的参考集，那么可以直接将它作为参考集，也可以将该已知的参考集与步骤2a)中利用测序片段构建的参考集组合并去冗余，从而提供最终的参考集。
例如，在人肠道微生物群落的MWAS研究中，Junjie Qin et al.(2010)A human gut microbial gene catalogue established bymetagenomic sequencing.Nature,464:59-65已构建并公开了3.3M欧洲人肠道微生物群落的非冗余基因集(即，参考集)。因此，在一个优选的实施方案中，所述环境样品是人肠道样品，并且将所述3.3M欧洲人肠道微生物群落的非冗余基因集与步骤2a)所构建的参考集组合并去冗余，从而提供最终的参考基因集。
在一个优选的实施方案中，通过使用SOAP2或MAQ比对软件，将所述测序片段与参考集进行比对。SOAP2和MAQ是本领域技术人员是已知的，参见例如，R Li et al.2009和Li et al.2008。
在一个优选的实施方案中，使用SOAP2将测序片段与参考集进行比对，并按照下列公式计算出参考集中的各元素的相对丰度：
α i = x i / L i Σ j ( x i / L i ) ,
其中
αi：元素i在样品中的相对丰度；
Li：元素i的长度；
xi：元素i在样品中被检测到的次数。
关于分栈
在一个优选的实施方案中，通过下列步骤来确定元素的初始栈：首先，基于元素在样品中的相对丰度，计算两两元素之间的相关性，例如pearson相关系数，spearman相关系数，kendall相关系数，欧几里得距离，曼哈顿距离等；然后，根据两两元素之间的相关性，通过聚类算法，如自底向上层次聚类等，将相关性密切的元素聚到一个类中，从而确定各个元素的初始栈。
在步骤3)的分栈之后，同一个栈里面的各个元素在所有样品中的丰度符合一定的分布模型，如正态分布。因此，聚到同一个栈里面的多个元素具有以下几种可能：(1)这些元素属于同一个物种；(2)这些元素来自共生的物种，因为共生物种的丰度分布相似；(3)这些元素是几个物种共有的，因为几个物种共有的元素的丰度不同于每一个物种各自的丰度。
关于基于栈的高级组装
在一个优选的实施方案中，使用SOAP2来将测序片段与已分栈的元素进行比对。
在一个优选的实施方案中，使用GC-depth spectra classifier和/或tetranucleotide frequencies(TNFs)classifier(Teeling et al.2004)进行校正。
关于鉴定
在一个优选的实施方案中，通过将各个栈的基因组序列与已知的基因组数据库进行比对，从而确定各个栈所对应的微生物的类别。
在一个优选的实施方案中，所述基因组数据库包括但不限于，NCBI/IMG已测序细菌库，NCBI的NR库等。
在一个优选的实施方案中，所述比对是核酸水平和/或蛋白水平的比对。
在另一个方面，本发明提供了一种用于分析环境样品中的微生物群落组成的装置，其包括以下模块：
1)测序模块，其用于对来自环境样品的宏基因组DNA进行测序，提供由测序片段池构成的宏基因组测序数据；
2)初级组装模块，其与测序模块相连，且包括彼此相连的下列模块：
2a)组装构建模块，其用于对测序片段进行组装以得到组装片段，然后去冗余，从而构建非冗余的参考集(即，组装片段集)；任选地，所述组装构建模块还可在所获得的组装片段上预测基因，并将预测出来的基因的集合作为参考集(即，基因集)；和
2b)比对计算模块，其用于将测序片段与参考集进行比对，并计算参考集中的各个元素在样品中的相对丰度；
3)分栈模块，其与初级组装模块相连，用于确定参考集中的每一个元素所归属的栈，得到聚类的栈，且包括彼此相连的下列模块：
3a)丰度分栈模块，其基于丰度确定各个元素的初始栈；和
3b)模型分栈模块，其基于模型确定各个元素所归属的栈；
4)高级组装模块，其与测序模块和分栈模块相连，其用于从宏基因组测序数据中寻找对应到各个栈的测序片段，并将对应到各个栈的测序片段分别进行组装，且对组装结果进行验证及调整；和
5)鉴定模块，其与高级组装模块相连，用于通过各个栈的基因组序列，确定各个栈所对应的微生物的类别，从而确定所述环境样品中的微生物群落组成。
在一个优选的实施方案中，环境样品来源于自然环境，例如土壤环境，海洋环境和河流环境。在另一个优选的实施方案中，环境样品来源于体内环境，例如口腔环境和肠道环境。
在一个优选的实施方案中，所述测序模块使用第二代测序技术(例如，454测序法，Solexa测序法，SOLiD测序法或单分子测序法)或第三代测序技术对环境样品所包含的微生物群落的宏基因组进行测序，从而提供来自环境样品的宏基因组测序数据。
在一个优选的实施方案中，所述装置还包括彼此相连的DNA提取模块和文库构建模块，其中，所述DNA提取模块用于从所述环境样品中提取宏基因组DNA，并且，所述文库构建模块与测序模块相连，且利用所述宏基因组DNA构建基因组文库。
在一个优选的实施方案中，所述测序模块所获得的测序片段是末端配对的测序片段(paired end reads)。
在一个优选的实施方案中，所述装置还包含过滤模块，其与测序模块和初级组装模块相连，用于在进行初级组装之前，去除测序片段中的接头序列、测序质量低的序列和/或宿主基因组序列。
在一个优选的实施方案中，所述测序模块对宏基因组的测序深度为至少10×，优选至少20×，优选至少30×，优选至少40×，更优选至少50×。
在一个优选的实施方案中，所述组装构建模块使用Soapdenovo将测序片段组装成连接片段和/或拼接片段。
在一个优选的实施方案中，所述组装构建模块还包含接收亚模块，其用于接收已知的参考集。在一个优选的实施方案中，所述组装构建模块将所接收的已知参考集作为最终的参考集。在另一个优选的实施方案中，所述组装构建模块将所接收的已知参考集与利用测序片段构建的参考集组合并去冗余，从而提供最终的参考集。
在一个优选的实施方案中，所述组装构建模块能够将来自多个样品的参考集组合在一起，并去冗余，从而构建最终的非冗余的参考集。
在一个优选的实施方案中，所述比对计算模块通过使用SOAP2或MAQ，将测序片段与参考集进行比对。
在一个优选的实施方案中，所述比对计算模块使用SOAP2将测序片段与参考集进行比对，并按照下列公式计算出参考集中各元素的相对丰度：
α i = x i / L i Σ j ( x i / L i ) ,
其中
αi：元素i在样品中的相对丰度；
Li：元素i的长度；
xi：元素i在样品中被检测到的次数。
在一个优选的实施方案中，所述丰度分栈模块基于元素在样品中的相对丰度，计算两两元素之间的相关性，然后通过聚类算法，确定各个元素的初始栈。
在一个优选的实施方案中，所述模型分栈模块通过下列来确定元素所归属的栈：
(i)将每一个初始栈作为一个独立的多元正态分布模型,并基于丰度矩阵，利用最大似然函数法计算所述模型的参数；
(ii)构建一个软矩阵(fuzzy matrix)，用于存储每一个元素归属某一个栈的概率；和
(iii)迭代运算E步和M步，直至似然函数达到最大化：
E步，根据每一个栈的模型参数，分别计算每一个元素属于某一个栈的后验概率，并且修改软矩阵中所述元素属于所述栈的概率；
M步：根据软矩阵，用最大似然函数法计算每一个栈的模型参数。
在一个优选的实施方案中，所述高级组装模块通过下列来实现其功能：
(a)通过将测序片段与已分栈的元素进行比对，从宏基因组测序数据中寻找对应到所述分栈模块所确定的各个栈的测序片段；
(b)使用SOAPdenovo或者使用其他针对微生物数据的组装软件，将对应到各个栈的测序片段分别进行组装；
(c)使用基于相似度的聚类方法和/或基于组成特征的聚类方法，对每一个栈所包含的元素的分栈进行校正；任选地，还在已获得的栈内部进行再次聚类，然后根据聚类的结果，对已获得的栈进行拆分或保持不变，从而使结果更加准确可信；
(d)重复步骤(a)－(c)，直到各个栈的基因组序列的大小没有明显变化为止(总长度增长率小于5％)。
在一个优选的实施方案中，所述高级组装模块使用SOAP2来将测序片段与已分栈的元素进行比对。
在一个优选的实施方案中，所述高级组装模块使用GC-depthspectra classifier和/或tetranucleotide frequencies(TNFs)classifier进行校正。
在一个优选的实施方案中，所述鉴定模块通过将各个栈的基因组序列与已知的基因组数据库进行比对，从而确定各个栈所对应的微生物的类别。
在一个优选的实施方案中，所述基因组数据库包括，但不限于，NCBI/IMG已测序细菌库，NCBI的NR库等。
在一个优选的实施方案中，所述鉴定模块在核酸水平和/或蛋白水平上进行比对。
在另一个方面，还提供了本发明的装置用于分析环境样品中的微生物群落组成的用途。在一个优选的实施方案中，所述环境样品来源于自然环境，例如土壤环境，海洋环境和河流环境。在另一个优选的实施方案中，环境样品来源于体内环境，例如口腔环境和肠道环境。
发明的有益效果
本发明的方法和装置基于高通量测序技术，利用相同或相似环境下多个样品的测序数据进行组装，聚类和再组装，从而得到微生物群落的物种组成信息和物种的基因组信息，有着非常广泛的应用前景。与现有技术中的传统组装方法相比较，本发明的方法和装置有如下优点：
1、系统地将各种测序序列的属性结合起来，用于构建微生物群落的宏基因组的参考集，这特别适合于微生物物种分类，以及从来自同一环境的多个样品的测序数据重建基因组；
2、创造性地将分栈和组装有效地结合在一起，使物种基因组的组装结果更加精确，从而能够实现高效率、高精度地确定微生物群落的组成；
3、首次基于多个样品进行聚类分析，并进行了迭代高级组装。
利用多个样品进行聚类分析具有有两个显著的优点：a)可以覆盖更多的低丰度物种，从而更全面地研究微生物群落；b)由于环境因素，不同的样品可能具有不同的物种组成和丰度，从而可以有利地进行比较研究。相比之下，利用单一样本进行的宏基因组学分析通常只能获得精确的优势物种，而无法全面地分析微生物群落，特别是低丰度物种(参见例如，Hess et al.2011)。
下面将结合附图和实施例对本发明的实施方案进行详细描述，但是本领域技术人员将理解，下列附图和实施例仅用于说明本发明，而不是对本发明的范围的限定。根据附图和优选实施方案的下列详细描述，本发明的各种目的和有利方面对于本领域技术人员来说将变得显然。
附图说明
图1示意性地描述了本发明的SoapMeta方法的流程图，其中，虚线空心框、实线空心框和实心框示意性表示源自三个不同的物种。
图2示意性地描述了本发明的SoapMeta方法的初级组装的流程图。
图3是示意性地描述了本发明的SoapMeta方法的分栈的流程图。
图4是示意性地描述了本发明的SoapMeta方法的高级组装的流程图。
图5是描述了用于实施本发明的SoapMeta方法的装置的结构示意图。
图6-8展示了实施例2中利用第一种策略获得的3个样品(样品A-C)的GC含量-测序深度谱图。图6：样品A；图7：样品B；图8：样品C。结果显示，样品B和样品C中的一些细菌很难区分，因为他们的GC含量和测序深度非常接近。
图9展示了本申请实施例3中通过16S rRNA测序获得的物种分类的信息图。
图10展示了利用16S rRNA测序法获得的Akkermansia属的16SrRNA标签的数量与利用本发明的Soapmeta方法组装出来的相应基因组的测序深度的相关性。
图11展示了利用16S rRNA测序法获得的Lactobacillus属的16SrRNA标签的数量与利用本发明的Soapmeta方法组装出来的相应基因组的测序深度的相关性。
图10-11的结果显示，利用16S rRNA测序法获得的rRNA标签的数量与利用本发明的Soapmeta方法组装出来的相应基因组的测序深度之间具有很强的相关性。这些结果表明，本发明的Soapmeta方法的结果与16S rRNA测序法的结果是基本上一致的，再次证实了本发明的SoapMeta方法的可靠性、准确性和高效性。
具体实施方式
现参照下列意在举例说明本发明(而非限定本发明)的实施例来描述本发明。
除非特别指明，本发明中所使用的分子生物学实验方法，基本上参照J.Sambrook等人，分子克隆：实验室手册，第2版，冷泉港实验室出版社，1989，以及F.M.Ausubel等人，精编分子生物学实验指南，第3版，John Wiley&Sons,Inc.，1995中所述的方法进行；并且各种酶的使用依照产品制造商推荐的条件。那些在实施例中未详细描述的过程和方法是本领域中公知的常规方法。本领域技术人员知晓，实施例以举例方式描述本发明，且不意欲限制本发明所要求保护的范围。
实施例1.模拟环境样品的分析
1、数据模拟
为了模拟环境样品，我们从NCBI基因组数据库(Wheeler et al.2007)中选取了100个不同的物种，这些物种的基因组从变形菌门中随机选择。另外，为了简化模型，不选择同一物种的不同品系。
我们一共模拟了10例样品，每个样品的测序量均为720M。模拟的末端配对的测序片段的长度为90bp，插入片段的大小为500±20bp(均值±标准差)，测序错误率为0.1％。通过Broken-Stick模型(MacArthur1957)的相对物种丰度(relative species abundance，RSA)，来确定每一个样品的物种丰度组成比例。每一个样品所包含的大多数细菌的测序量是比较低的(64％的细菌的RSA<0.01)。将10个样品的测序数据合并后，这些低丰度细菌的测序量达13.6-182.0Mbp，且测序深度为2.7－160.4X。
2、初级组装
我们将所有样品的测序数据(测序片段)合并在一起，并使用组装软件Soapdenovo(Li et al.2009)进行初步的组装(即，不单独对每一个样品的测序数据进行分别的组装)。在组装后，对组装结果进行去冗余，从而得到非冗余的参考集。
特别地，在本实验中，混合样品的初级组装结果(即，参考集)共包含41754条连接片段(contigs)，且连接片段的长度范围为200-2,001,157bp(N50＝93,353bp)(N50是衡量基因组图谱质量的一个判断标准，其是指，当将所有的组装得到的序列按照长度从大到小排列，并从大到小将序列的长度相加，直至相加得到的总长度为所有组装得到的序列的总长度的百分之五十时，那条组装序列的长度，参见例如，Miller et al.2010.Assembly algorithms for next generationsequencing data.Genomics.95(6):315-327)。将这些连接片段与原始细菌基因组进行BLASTN比对。结果显示，组装后的连接片段对原始细菌基因组的平均覆盖度为88.7％，并且每个细菌的覆盖度与测序深度呈现正相关，但是，当测序深度高于20x时，参考集的覆盖度不再发生显著的变化。
使用SOAP2，将测序片段与非冗余参考集进行比对，并通过下式计算出参考集中的各连接片段的相对丰度：
&alpha; i = x i / L i &Sigma; j ( x i / L i ) ,
其中，
αi：连接片段i在样品中的相对丰度；
Li：连接片段i的长度；
xi：连接片段i在样品中被检测到的次数。
3、分栈(bin)
3.1基于丰度的分栈(初始分栈)
首先计算丰度矩阵中各连接片段的两两Kendall's tau秩相关系数；然后根据连接片段两两之间的相关性，采用自底向上层次聚类算法，将相关性比较密切的片段聚到一个类中，从而获得初始的栈。
在本实验中，我们还使用默认的聚类参数，过滤掉了包含小于10个连接片段的初始栈，最终得到343个初始栈。这些栈覆盖了96.8％的连接片段(40,438/41,754)。
对于每一个初始栈，我们还给它定义一个属性，“最优的比对细菌”。也即，如果栈里面大部分的连接片段来自于某一个特定的细菌，那么这个细菌就是这个初始栈的最优的比对细菌。另外，还将栈的精度定为，来自最优的比对细菌的连接片段的总长度占栈里面的连接片段的总长度的百分比。在本实验中，初始栈的精度为50.3％－100.0％(平均值为95.1％)。
3.2基于模型的分栈
我们接着用基于模型的分栈方法来最优化初始分栈的结果。简言之，1)将每一个初始栈作为一个独立的多元正态分布模型,并基于丰度矩阵，利用最大似然函数法计算所述模型的参数；
2)构建一个软矩阵(fuzzy matrix)，用于存储每一个连接片段归属某一个栈的概率；
3)迭代运算E步和M步，直至似然函数达到最大化：
E步，根据每一个栈的模型参数，分别计算每一个连接片段属于某一个栈的后验概率，并且修改软矩阵中所述连接片段属于所述栈的概率；
M步：根据软矩阵，用最大似然函数法计算每一个栈的模型参数。
在该步骤后，所获得的栈减少到135个。与初始分栈相比，这些栈的覆盖度下降到91.9％(38,364/41,754个连接片段)，且精度下降到33.2％－100.0％(平均值92.3％)。在这135个栈中，每一个栈代表一个物种。基于各个栈中的连接片段的序列，我们鉴定到了86个物种(86％)，且每一个物种的基因组覆盖度超过50％。
4、高级组装
高级组装分成以下3步：
1)使用SOAP2，通过序列比对，在模拟的测序数据中寻找对应到之前确定的各个栈的测序片段；
2)使用SOAPdenovo分别将对应到各个栈的测序片段进行深度组装；
3)使用基于相似度和组成特征的聚类方法，对每一个栈所包含的连接片段的分栈进行校正，并且在已有的栈内部进行再次聚类，然后根据再次聚类的结果，对已有的栈进行拆分或保持不变，从而使结果更加准确可信；
4)重复步骤1)-3)，直到各个栈的基因组序列的大小没有明显变化为止(总长度增长率小于5％)。
对之前获得的135个栈进行高级组装之后，得到148个经组装的栈。栈的数目的增加是因为，我们使用了基于组成特征的聚类方法，根据GC含量，测序深度等特征，将一个栈里面的一些可以明显再细分的栈拆开了。
在高级组装后，栈的平均精度达到94.2％(参见，表1)，略微高于前一步的结果。另外，当用原始细菌基因组覆盖组装的栈的基因组时，结果显示，覆盖度为95.5％；反之，当用组装的栈的基因组覆盖原始细菌基因组时，覆盖度为57.4％。
在这148个栈中，基于各个栈的组装的基因组序列，我们鉴定到了100个初始细菌物种中的95个(95％)，且如上所述，每一个物种的基因组覆盖度超过50％。
上述结果表明，本发明的SoapMeta方法的特异度较好，且能够有效地鉴别出模拟样品中所包含的绝大部分物种(95％)。
表1.每一步骤得到的栈的比较
实施例2.简单环境样品(纤维素降解菌群)的分析
本实施例以一个真实的简单环境为例，对本发明的SoapMeta方法进行了进一步的解释说明，并且通过与传统的分析方法相比较，证实了本发明的SoapMeta方法的优势。
在本实施例中，我们收集了三个样品(样品A、B、C)，它们分别来自不同培养条件下的纤维素降解菌群：从同一沼泽的土壤采集3个样品，并且分别用三种包含不同碳源(滤纸、纤维二糖、葡萄糖)的培养基在37℃下培养52小时，然后分别收获菌体，从而获得样品A、B、C。针对每一个样品，我们分别构建了一个测序文库(参数设置：末端配对的测序片段的长度为90bp，插入片段的大小为500±20bp)：首先用HiSeq2000对样品进行测序，从而得到原始测序片段(rawreads)；然后，过滤掉其中的低质量序列和接头序列，从而提供3.88Gbp的用于分析的宏基因组测序数据(3个样品的测序数据的总和)。
在本实施例中，我们应用了两种策略来构建微生物的基因组。第一种策略是，使用传统的分析方法，对每个样品分别进行测序数据的组装，从而构建微生物的基因组(参见，MEGAN(Husonet al.2007))；第二种策略是，使用本发明的SoapMeta方法，将所有样品的测序数据混合在一起，然后进行初级组装，分栈和高级组装，从而构建微生物的基因组。将第一种策略用作对照，以证实本发明的SoapMeta方法在多个样品的混合组装方面的优势。
在第一种策略下，用基于组成特征的聚类方法对来自单个样品的测序片段进行聚类，以判别样品中潜在的微生物。对于所使用的3个样品，我们分别得到了6个类(样品A)，2个类(样品B)，和3个类(样品C)。这3个样品各自的GC图(参见图6-8)显示，样品B和样品C中的一些细菌很难区分，因为他们的GC含量和测序深度非常接近。
在第二种策略下，我们首先在初级组装中得到了连接片段的相对丰度。进一步，通过使用本发明的SoapMeta方法，我们从3个样品的混合测序数据中鉴定到了10个栈，其中有9个栈的组装的基因组序列大于1Mbp，并且这10个栈的基因组序列总长覆盖了所有样品测序数据的89.5％。在这10个栈中，每一个栈对应一个潜在的物种。随后，我们对每个栈的组装的基因组序列进行了TBLASTX比对，以确定各个栈所对应的潜在的物种，结果见表2。
表2的结果显示，在这10个栈中，有6个栈的组装的基因组序列很纯(即，基本上对应至同一个微生物物种的基因组)：短短芽孢杆菌NBRC 100599(Brevibacillus brevis NBRC 100599)、凝结芽孢杆菌2-6(Bacillus coagulans 2-6)、耐盐芽孢杆菌C-125(Bacillus haloduransC-125)、肉毒梭菌A2Kyoto(Clostridium botulinum A2 Kyoto)、热解纤维梭菌ATCC27405(Clostridium thermocellum ATCC 27405)、热解纤维梭菌ATCC27405(Clostridium thermocellum ATCC 27405)，并且其中的热解纤维梭菌(Clostridium thermocellum)是众所周知的纤维素降解菌(Weimer and Zeikus1977；Bayer et al.1983；和Schwarz2001)。此外，其中的短芽孢杆菌(Brevibacillus)和芽孢杆菌(Bacillus)也已知具有纤维降解能力(Liang et al.2009；Li et al.2006；和Rastogi et al.2009)。
从上面的结果可知，本发明的SoapMeta策略不仅在精度和覆盖度上显著优于第一种策略(即，基因组覆盖度更全，分类准确度更高)，而且能够更有效、更精确地鉴定环境样品的微生物组成。
表2、纤维素降解菌群的组装基因组总表
注：图中的*表示，该栈包含有多个物种的序列，并且无法进一步明确区分。例如，B1*表示，栈B1中含有无法进一步区分的多个物种的序列(在使用第二种策略的方法中，栈B1中的这些物种被进一步区分为短短芽孢杆菌NBRC 100599和热解纤维梭菌ATCC 27405)。
实施例3.复杂环境样品(小鼠肠道菌群)的分析
本实施例以一个真实的复杂环境为例，示例性地展示了本发明的SoapMeta方法在小鼠肠道菌群的探测中的应用。
本实验采用了两种常见的小鼠，SV-129和C57Black/6(Fujii etal.1997)。在现实中，小鼠肠道的菌群的相对丰度会随着年龄，性别，饮食等等因素的变化而变化，但是如果小鼠的饮食固定，且环境固定的话，这些菌群的微生物组成一般不会有太大的变动。因此，可以利用本发明的SoapMeta方法来研究特定环境、特定饮食下小鼠的肠道菌群的微生物组成，并构建菌群物种的基因组。
收集了13个粪便样品(其中6个样品来自SV-129小鼠，7个样品来自C57Black/6小鼠)，并构建了测序文库(参数设置：末端配对的测序片段的长度为90bp，插入片段的大小为350±15bp)：首先用HiSeq2000对样品进行测序，从而得到原始测序片段(raw reads)；然后，过滤掉其中的低质量序列、接头序列以及小鼠基因组序列，从而获得3.96±0.55Gbp(每个样品的平均测序数据)的用于分析的宏基因组测序数据。
根据本发明的SoapMeta方法：
首先，对样品的宏基因组测序数据进行了初级组装，得到246.1Mbp的连接片段集(n＝180,056个，N50＝2,613bp)；
然后，进行了分栈，得到325个栈(将序列含量低于100Kbp的栈过滤掉)，这些栈的总序列含量为213.6Mbp(86.8％)，并且其中有56个栈的序列含量大于1Mbp；
最后，对上述序列含量大于1Mbp的56个栈进行了高级组装，最终得到57个基因组(栈)，其总序列含量达141.6Mbp(每个基因组的平均序列含量为2.48Mbp)，并且覆盖了49.5％的测序片段。结果概述于表3中。
使用BLASTN(核酸水平)和TBLASTX(蛋白水平)，将高级组装得到的栈与已知的基因组数据库进行比对。结果显示，有8个栈在核酸水平上与已知的物种十分接近：它们均具有高于90％的精度和高于95％的序列相似度。此外，还有48个栈在蛋白水平上与已知的物种高度同源：它们均具有高于70％的精度和高于50％的序列相似度。另外，还有1个栈比对到未知的物种。
为了验证上述结果，我们通过Solexa测序法对这些样品的16SrRNA(V6高变区)进行了测序，得到高质量的3.63±0.68M(均值±标准差)的16S rRNA标签(tags)(已过滤掉接头序列，低质量序列，重叠序列和引物序列)。利用BLASTN，将这些16S rRNA标签与RefSSU数据库(Huse et al.2010)进行比对。结果示于图9中。结果显示，小鼠肠道菌群中，丰度较高的微生物是：毛螺菌科(Lachnospiraceae)、乳杆菌属(Lactobacillus)、别样棒菌属(Allobaculum)、阿克曼氏菌属(Akkermansia)、Ruminococcaeae、乳头杆菌属(Papillibacter)、拟杆菌属(Bacteroides)和脱硫弧菌科(Desulfovibrionaceae)。这些细菌大部分能够被本发明的SoapMeta方法组装出来的基因组覆盖，这充分表明，本发明的SoapMeta方法能够高效、精确地鉴定环境样品中的微生物组成。
另外，我们还将Akkermansia属和Lactobacillus属的16S rRNA标签的数量与用Soapmeta方法组装出来的基因组的测序深度做比较。结果显示，它们之间具有很强的相关性(参见图10-11)。这再次表明了本发明的SoapMeta方法的准确性和高效性。
尽管本发明的具体实施方式已经得到详细的描述，但本领域技术人员将理解：根据已经公开的所有教导，可以对细节进行各种修改和变动，并且这些改变均在本发明的保护范围之内。本发明的全部范围由所附权利要求及其任何等同物给出。

Claims (61)

1.一种用于分析环境样品中的微生物群落组成的方法，其包括以下步骤：
1)测序：
对来自环境样品的基因组DNA进行构建文库和测序，从而获得由测序片段池构成的宏基因组测序数据；
2)初级组装：
2a)构建或完善参考集：对测序片段进行组装以得到组装片段，然后去冗余，从而构建非冗余的参考集；或者，在所获得的组装片段上预测基因，并将预测出来的基因的集合作为参考集；或者，如果针对所述环境样品，存在已知的参考集，那么直接将它作为参考集，或者将该已知的参考集与如上所述构建的参考集组合并去冗余，从而获得最终的参考集；
2b)构建元素相对丰度谱矩阵：将所述测序片段与参考集进行比对，并计算参考集中的各个元素在样品中的相对丰度；
3)分栈，即，通过下述步骤确定参考集中的每一个元素所归属的栈，得到聚类的栈：
3a)基于丰度的分栈：基于元素在样品中的相对丰度，使用聚类算法，确定各个元素的初始栈；和
3b)基于模型的分栈：
(i)将每一个初始栈作为一个独立的多元正态分布模型,并基于丰度矩阵，利用最大似然函数法计算所述模型的参数；
(ii)构建一个软矩阵，用于存储每一个元素归属某一个栈的概率；和
(iii)迭代运算E步和M步，直至似然函数达到最大化：
E步，根据每一个栈的模型参数，分别计算每一个元素属于某一个栈的后验概率，并且修改软矩阵中所述元素属于所述栈的概率；
M步：根据软矩阵，用最大似然函数法计算每一个栈的模型参数；
4)基于栈的高级组装：
4a)通过将测序片段与已分栈的元素进行比对，从宏基因组测序数据中寻找对应到之前确定的各个栈的测序片段；
4b)使用SOAPdenovo或者使用其他针对微生物测序数据的组装软件，将对应到各个栈的测序片段分别进行组装；
4c)使用基于相似度的聚类方法和/或基于组成特征的聚类方法，对每一个栈所包含的元素的分栈进行校正；任选地，还在已获得的栈内部进行再次聚类，然后根据聚类的结果，对已获得的栈进行拆分或保持不变；
4d)重复步骤4a)－4c)，直到各个栈的基因组序列的大小无明显变化为止，即，总长度增长率小于5％；
5)鉴定：
利用各个栈的基因组序列，确定各个栈所对应的微生物的类别，从而确定所述环境样品中的微生物群落组成；
其中，所述方法用于非诊断目的。
2.权利要求1的方法，其中，所述环境样品来源于自然环境；或者所述环境样品来源于体内环境。
3.权利要求2的方法，其中所述自然环境选自土壤环境，海洋环境和河流环境。
4.权利要求2的方法，其中所述体内环境选自口腔环境和肠道环境。
5.权利要求1的方法，其中在步骤1)中使用第二代测序技术或第三代测序技术对环境样品所包含的微生物群落的宏基因组进行测序，从而提供来自环境样品的宏基因组测序数据。
6.权利要求5的方法，其中所述第二代测序技术选自454测序法，Solexa测序法，SOLiD测序法和单分子测序法。
7.权利要求1的方法，其中，在步骤1)中通过下列步骤来获得宏基因组测序数据：
1a)提供环境样品；
1b)从所述环境样品中提取宏基因组DNA；
1c)利用所述宏基因组DNA构建宏基因组文库；
1d)对所述宏基因组文库进行测序，从而提供所述环境样品的宏基因组测序数据。
8.权利要求7的方法，其中，在步骤1d)中使用Solexa测序法进行测序。
9.权利要求1的方法，其中所述测序片段是末端配对的测序片段。
10.权利要求1的方法，其中在进行步骤2)之前，对测序数据进行预处理，即，去除接头序列、测序质量低的序列和/或宿主基因组序列。
11.权利要求1的方法，其中，对来自相同或相似环境的多个样品进行测序，并将所有样品的测序数据组合在一起，构成宏基因组测序数据。
12.权利要求1的方法，其中所述宏基因组的测序深度为至少10×。
13.权利要求1的方法，其中所述宏基因组的测序深度为至少20×。
14.权利要求1的方法，其中所述宏基因组的测序深度为至少30×。
15.权利要求1的方法，其中所述宏基因组的测序深度为至少40×。
16.权利要求1的方法，其中所述宏基因组的测序深度为至少50×。
17.权利要求1的方法，其中，在步骤2)中，使用Soapdenovo将所述测序片段组装成组装片段。
18.权利要求17的方法，其中，使用Soapdenovo将所述测序片段组装成连接片段和/或拼接片段。
19.权利要求1的方法，其中，在步骤2)中，将来自多个样品的参考集组合在一起，并去冗余，从而构建最终的非冗余的参考集。
20.权利要求1的方法，其中，在步骤2)中，使用SOAP2或MAQ比对软件，将所述测序片段与参考集进行比对。
21.权利要求1的方法，其中，在步骤2)中，使用SOAP2将测序片段与参考集进行比对，并按照下列公式计算出参考集中的各元素的相对丰度：
&alpha; i = x i / L i &Sigma; j ( x i / L i ) ,
其中
αi：元素i在样品中的相对丰度；
Li：元素i的长度；
xi：元素i在样品中被检测到的次数。
22.权利要求1的方法，其中，在步骤3)中，通过下列步骤来确定元素的初始栈：首先，基于元素在样品中的相对丰度，计算两两元素之间的相关性；然后，根据两两元素之间的相关性，通过聚类算法，将相关性密切的元素聚到一个类中，从而确定各个元素的初始栈。
23.权利要求22的方法，其中所述相关性选自：pearson相关系数，spearman相关系数，kendall相关系数，欧几里得距离，和曼哈顿距离。
24.权利要求22的方法，其中所述聚类算法为自底向上层次聚类。
25.权利要求1的方法，其中，在步骤4)中，使用SOAP2来将测序片段与已分栈的元素进行比对。
26.权利要求1的方法，其中，在步骤4)中，使用GC-depthspectra classifier和/或tetranucleotide frequenciesclassifier进行校正。
27.权利要求1的方法，其中，在步骤5)中，通过将各个栈的基因组序列与已知的基因组数据库进行比对，从而确定各个栈所对应的微生物的类别。
28.权利要求27的方法，其中，所述基因组数据库选自下列：NCBI/IMG已测序细菌库，和NCBI的NR库。
29.权利要求27的方法，其中，所述比对是核酸水平和/或蛋白水平的比对。
30.权利要求1的方法，其中，所述步骤3a)中的聚类算法为自底向上的层次聚类方法。
31.一种用于分析环境样品中的微生物群落组成的装置，其包括以下模块：
1)测序模块，其用于对来自环境样品的宏基因组DNA进行测序，提供由测序片段池构成的宏基因组测序数据；
2)初级组装模块，其与测序模块相连，且包括彼此相连的下列模块：
2a)组装构建模块，其用于对测序片段进行组装以得到组装片段，然后去冗余，从而构建非冗余的参考集；任选地，所述组装构建模块还在所获得的组装片段上预测基因，并将预测出来的基因的集合作为参考集；和
2b)比对计算模块，其用于将测序片段与参考集进行比对，并计算参考集中的各个元素在样品中的相对丰度；
3)分栈模块，其与初级组装模块相连，用于确定参考集中的每一个元素所归属的栈，得到聚类的栈，且包括彼此相连的下列模块：
3a)丰度分栈模块，其基于丰度确定各个元素的初始栈；和
3b)模型分栈模块，其基于模型确定各个元素所归属的栈；
4)高级组装模块，其与测序模块和分栈模块相连，其用于从宏基因组测序数据中寻找对应到各个栈的测序片段，并将对应到各个栈的测序片段分别进行组装，且对组装结果进行验证及调整；和
5)鉴定模块，其与高级组装模块相连，用于通过各个栈的基因组序列，确定各个栈所对应的微生物的类别，从而确定所述环境样品中的微生物群落组成。
32.权利要求31的装置，其中，所述环境样品来源于自然环境；或者，所述环境样品来源于体内环境。
33.权利要求32的装置，其中，所述自然环境选自土壤环境，海洋环境和河流环境。
34.权利要求32的装置，其中，所述体内环境选自口腔环境和肠道环境。
35.权利要求31的装置，其中所述测序模块使用第二代测序技术或第三代测序技术对环境样品所包含的微生物群落的宏基因组进行测序，从而提供来自环境样品的宏基因组测序数据。
36.权利要求35的装置，其中所述第二代测序技术选自454测序法，Solexa测序法，SOLiD测序法和单分子测序法。
37.权利要求31的装置，其中所述装置还包括彼此相连的DNA提取模块和文库构建模块，其中，所述DNA提取模块用于从所述环境样品中提取宏基因组DNA，并且，所述文库构建模块与测序模块相连，且利用所述宏基因组DNA构建基因组文库。
38.权利要求31的装置，其中所述测序模块所获得的测序片段是末端配对的测序片段。
39.权利要求31的装置，其中所述装置还包含过滤模块，其与测序模块和初级组装模块相连，用于在进行初级组装之前，去除测序片段中的接头序列、测序质量低的序列和/或宿主基因组序列。
40.权利要求31的装置，其中所述测序模块对宏基因组的测序深度为至少10×。
41.权利要求31的装置，其中所述测序模块对宏基因组的测序深度为至少20×。
42.权利要求31的装置，其中所述测序模块对宏基因组的测序深度为至少30×。
43.权利要求31的装置，其中所述测序模块对宏基因组的测序深度为至少40×。
44.权利要求31的装置，其中所述测序模块对宏基因组的测序深度为至少50×。
45.权利要求31的装置，其中，所述组装构建模块使用Soapdenovo将测序片段组装成连接片段和/或拼接片段。
46.权利要求31的装置，其中所述组装构建模块还包含接收亚模块，其用于接收已知的参考集。
47.权利要求46的装置，其中所述组装构建模块将所接收的已知参考集作为最终的参考集，或者将所接收的已知参考集与利用测序片段构建的参考集组合并去冗余，从而提供最终的参考集。
48.权利要求31的装置，其中所述组装构建模块能够将来自多个样品的参考集组合在一起，并去冗余，从而构建最终的非冗余的参考集。
49.权利要求31的装置，其中所述比对计算模块通过使用SOAP2或MAQ，将测序片段与参考集进行比对。
50.权利要求49的装置，其中所述比对计算模块使用SOAP2将测序片段与参考集进行比对，并按照下列公式计算出参考集中各元素的相对丰度：
&alpha; i = x i / L i &Sigma; j ( x i / L i ) ,
其中
αi：元素i在样品中的相对丰度；
Li：元素i的长度；
xi：元素i在样品中被检测到的次数。
51.权利要求31的装置，其中，所述丰度分栈模块基于元素在样品中的相对丰度，计算两两元素之间的相关性，然后通过聚类算法，确定各个元素的初始栈。
52.权利要求31的装置，其中，所述模型分栈模块通过下列来确定元素所归属的栈：
(i)将每一个初始栈作为一个独立的多元正态分布模型,并基于丰度矩阵，利用最大似然函数法计算所述模型的参数；
(ii)构建一个软矩阵，用于存储每一个元素归属某一个栈的概率；和
(iii)迭代运算E步和M步，直至似然函数达到最大化：
E步，根据每一个栈的模型参数，分别计算每一个元素属于某一个栈的后验概率，并且修改软矩阵中所述元素属于所述栈的概率；
M步：根据软矩阵，用最大似然函数法计算每一个栈的模型参数。
53.权利要求31的装置，其中，所述高级组装模块通过下列来实现其功能：
(a)通过将测序片段与已分栈的元素进行比对，从宏基因组测序数据中寻找对应到分栈模块所确定的各个栈的测序片段；
(b)使用SOAPdenovo或者使用其他针对微生物数据的组装软件，将对应到各个栈的测序片段分别进行组装；
(c)使用基于相似度的聚类方法和/或基于组成特征的聚类方法，对每一个栈所包含的元素的分栈进行校正；任选地，还在已获得的栈内部进行再次聚类，然后根据聚类的结果，对已获得的栈进行拆分或保持不变；
(d)重复步骤(a)－(c)，直到各个栈的基因组序列的大小没有明显变化为止，即，总长度增长率小于5％。
54.权利要求53的装置，其中，所述高级组装模块使用SOAP2来将测序片段与已分栈的元素进行比对。
55.权利要求53的装置，其中，所述高级组装模块使用GC-depthspectra classifier和/或tetranucleotide frequenciesclassifier进行校正。
56.权利要求31的装置，其中，所述鉴定模块通过将各个栈的基因组序列与已知的基因组数据库进行比对，从而确定各个栈所对应的微生物的类别。
57.权利要求56的装置，其中，所述基因组数据库选自NCBI/IMG已测序细菌库和/或NCBI的NR库。
58.权利要求56的装置，其中，所述鉴定模块在核酸水平和/或蛋白水平上进行比对。
59.权利要求31-58任一项的装置用于分析环境样品中的微生物群落组成的用途，其中，所述环境样品来源于自然环境；或者，所述环境样品来源于体内环境。
60.权利要求59的用途，其中，所述自然环境选自土壤环境，海洋环境和河流环境。
61.权利要求59的用途，其中，所述体内环境选自口腔环境和肠道环境。
CN201280064063.2A 2012-08-01 2012-08-01 一种分析微生物群落组成的方法和装置 Active CN104039982B (zh)

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
PCT/CN2012/079492 WO2014019164A1 (zh) 2012-08-01 2012-08-01 一种分析微生物群落组成的方法和装置

Publications (2)

Publication Number Publication Date
CN104039982A CN104039982A (zh) 2014-09-10
CN104039982B true CN104039982B (zh) 2015-09-09

Family

ID=50027091

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201280064063.2A Active CN104039982B (zh) 2012-08-01 2012-08-01 一种分析微生物群落组成的方法和装置

Country Status (4)

Country Link
US (1) US20150242565A1 (zh)
CN (1) CN104039982B (zh)
HK (1) HK1196642A1 (zh)
WO (1) WO2014019164A1 (zh)

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105925664A (zh) * 2016-03-30 2016-09-07 广州精科生物技术有限公司 一种确定核酸序列的方法及系统
CN105950707A (zh) * 2016-03-30 2016-09-21 广州精科生物技术有限公司 一种确定核酸序列的方法及系统
CN108197434A (zh) * 2018-01-16 2018-06-22 深圳市泰康吉音生物科技研发服务有限公司 去除宏基因组测序数据中人源基因序列的方法
TWI629607B (zh) * 2017-08-15 2018-07-11 極諾生技股份有限公司 建立腸道菌數據庫的方法和相關檢測系統

Families Citing this family (28)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2015048595A1 (en) 2013-09-27 2015-04-02 Jay Shendure Methods and systems for large scale scaffolding of genome assemblies
CN105095688A (zh) * 2014-08-28 2015-11-25 吉林大学 检测人体肠道宏基因组的细菌群落及丰度的方法
CN104278091A (zh) * 2014-09-26 2015-01-14 上海交通大学 以废水处理样品微生物元基因组序列拼接细菌基因组的方法
EP3250708B1 (en) * 2015-01-30 2021-08-11 BGI Shenzhen Biomarkers for colorectal cancer related diseases
WO2017156739A1 (zh) * 2016-03-17 2017-09-21 上海锐翌生物科技有限公司 分离的核酸及应用
WO2018080477A1 (en) * 2016-10-26 2018-05-03 The Joan & Irwin Jacobs Technion-Cornell Institute Systems and methods for ultra-fast identification and abundance estimates of microorganisms using a kmer-depth based approach and privacy-preserving protocols
CN106778078B (zh) * 2016-12-20 2019-04-09 福建师范大学 基于kendall相关系数的DNA序列相似性比对方法
US10733214B2 (en) 2017-03-20 2020-08-04 International Business Machines Corporation Analyzing metagenomics data
CN107028606A (zh) * 2017-04-21 2017-08-11 上海耐相智能科技有限公司 医用智能监测环系统
US20200160936A1 (en) * 2017-06-28 2020-05-21 Icahn School Of Medicine At Mount Sinai Methods for high-resolution microbiome analysis
CN107287332A (zh) * 2017-08-03 2017-10-24 华子昂 利用smrt测序技术进行液体酵素菌种鉴定的方法
CN109587001B (zh) * 2018-11-15 2020-11-27 新华三信息安全技术有限公司 一种性能指标异常检测方法及装置
CN111455021B (zh) * 2019-01-18 2024-06-04 广州微远医疗器械有限公司 去除宏基因组中宿主dna的方法及试剂盒
US20220230704A1 (en) * 2019-06-13 2022-07-21 Icahn School Of Medicine At Mount Sinai Dna methylation based high resolution characterization of microbiome using nanopore sequencing
CN110277139B (zh) * 2019-06-18 2023-03-21 江苏省产品质量监督检验研究院 一种基于互联网的微生物限度检查系统及方法
CN110349629B (zh) * 2019-06-20 2021-08-06 湖南赛哲医学检验所有限公司 一种利用宏基因组或宏转录组检测微生物的分析方法
CN111261231A (zh) * 2019-12-03 2020-06-09 康美华大基因技术有限公司 肠道菌群宏基因组数据库构建方法、分析方法及装置
CN111161798B (zh) * 2019-12-31 2024-03-19 余珂 宏基因组的重组装方法、重组装装置及终端设备
CN111477267B (zh) * 2020-03-06 2022-05-03 清华大学 微生物的多关联网络计算方法、装置、设备及存储介质
CN111627500A (zh) * 2020-04-16 2020-09-04 中国科学院生态环境研究中心 一种基于宏基因组技术识别水体中携带毒性因子病原菌的方法
CN114067911B (zh) * 2020-08-07 2024-02-06 西安中科茵康莱医学检验有限公司 获取微生物物种及相关信息的方法和装置
CN112071366B (zh) * 2020-10-13 2024-02-27 南开大学 一种基于二代测序技术的宏基因组数据分析方法
CN112786102B (zh) * 2021-01-25 2022-10-21 北京大学 一种基于宏基因组学分析精准识别水体中未知微生物群落的方法
WO2022222936A1 (en) * 2021-04-20 2022-10-27 Hangzhou Matridx Biotechnology Co., Ltd. Methods, computer-readble media, and systems for filtering noises for dna sequencing data
CN113362890B (zh) * 2021-04-28 2023-09-08 中国科学院生态环境研究中心 一种评价生物滤料降解有机物潜力的方法
CN113284560B (zh) * 2021-04-28 2022-05-17 广州微远基因科技有限公司 病原检测背景微生物判断方法及应用
CN113611359B (zh) * 2021-08-13 2022-08-05 江苏先声医学诊断有限公司 一种提高宏基因组纳米孔测序数据菌种组装效率的方法
CN114999574B (zh) * 2022-08-01 2022-12-27 中山大学 一种肠道菌群大数据的并行识别分析方法及系统

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102517392A (zh) * 2011-12-26 2012-06-27 深圳华大基因研究院 基于宏基因组16s高可变区v3的分类方法和装置

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102517392A (zh) * 2011-12-26 2012-06-27 深圳华大基因研究院 基于宏基因组16s高可变区v3的分类方法和装置

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Environmental microbiology through the lens of high-throughput DNA sequencing:Synopsis of current platforms and bioinformatics approaches;Rmiro Logares等;《JOURNAL OF MICROBIOLOGICAL METHODS》;20120728;第81卷;106-113 *
Taxonomic classification of metagenomic sequences;Wolfgang Gerlach;《德国Bielefeld大学博士学位论文》;20120229;全文 *

Cited By (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN105925664A (zh) * 2016-03-30 2016-09-07 广州精科生物技术有限公司 一种确定核酸序列的方法及系统
CN105950707A (zh) * 2016-03-30 2016-09-21 广州精科生物技术有限公司 一种确定核酸序列的方法及系统
TWI629607B (zh) * 2017-08-15 2018-07-11 極諾生技股份有限公司 建立腸道菌數據庫的方法和相關檢測系統
CN108197434A (zh) * 2018-01-16 2018-06-22 深圳市泰康吉音生物科技研发服务有限公司 去除宏基因组测序数据中人源基因序列的方法
CN108197434B (zh) * 2018-01-16 2020-04-10 深圳市泰康吉音生物科技研发服务有限公司 去除宏基因组测序数据中人源基因序列的方法

Also Published As

Publication number Publication date
HK1196642A1 (zh) 2014-12-19
WO2014019164A1 (zh) 2014-02-06
US20150242565A1 (en) 2015-08-27
CN104039982A (zh) 2014-09-10

Similar Documents

Publication Publication Date Title
CN104039982B (zh) 一种分析微生物群落组成的方法和装置
Wu et al. A novel abundance-based algorithm for binning metagenomic sequences using l-tuples
Gruber-Vodicka et al. phyloFlash: rapid small-subunit rRNA profiling and targeted assembly from metagenomes
Bickhart et al. Generating lineage-resolved, complete metagenome-assembled genomes from complex microbial communities
Eren et al. Oligotyping: differentiating between closely related microbial taxa using 16S rRNA gene data
Kermarrec et al. Next‐generation sequencing to inventory taxonomic diversity in eukaryotic communities: a test for freshwater diatoms
Antipov et al. plasmidSPAdes: assembling plasmids from whole genome sequencing data
Brbić et al. The landscape of microbial phenotypic traits and associated genes
Rigaill et al. Synthetic data sets for the identification of key ingredients for RNA-seq differential analysis
Jin et al. Hybrid, ultra-deep metagenomic sequencing enables genomic and functional characterization of low-abundance species in the human gut microbiome
CN113744807B (zh) 一种基于宏基因组学的病原微生物检测方法及装置
US20160364523A1 (en) Systems and methods for identifying microorganisms
Meiser et al. Sequencing genomes from mixed DNA samples-evaluating the metagenome skimming approach in lichenized fungi
Brealey et al. Dental calculus as a tool to study the evolution of the mammalian oral microbiome
Jeraldo et al. Capturing one of the human gut microbiome’s most wanted: reconstructing the genome of a novel butyrate-producing, clostridial scavenger from metagenomic sequence data
Scott et al. Optimization and performance testing of a sequence processing pipeline applied to detection of nonindigenous species
US20190119717A1 (en) Ribosomal rna sequence extraction method and microorganism identification method using extracted ribosomal rna sequence
US20140288844A1 (en) Characterization of biological material in a sample or isolate using unassembled sequence information, probabilistic methods and trait-specific database catalogs
Zhang et al. A comprehensive investigation of metagenome assembly by linked-read sequencing
CN113260710A (zh) 用于通过多个定制掺合混合物验证微生物组序列处理和差异丰度分析的组合物、系统、设备和方法
Wu et al. Constructing metagenome-assembled genomes for almost all components in a real bacterial consortium for binning benchmarking
Yuan et al. RNA-CODE: a noncoding RNA classification tool for short reads in NGS data lacking reference genomes
Wu et al. A novel abundance-based algorithm for binning metagenomic sequences using l-tuples
Chandrasiri et al. CH-Bin: A convex hull based approach for binning metagenomic contigs
Taghavi et al. Distilled single-cell genome sequencing and de novo assembly for sparse microbial communities

Legal Events

Date Code Title Description
C06 Publication
PB01 Publication
C10 Entry into substantive examination
SE01 Entry into force of request for substantive examination
REG Reference to a national code

Ref country code: HK

Ref legal event code: DE

Ref document number: 1196642

Country of ref document: HK

C14 Grant of patent or utility model
GR01 Patent grant
REG Reference to a national code

Ref country code: HK

Ref legal event code: GR

Ref document number: 1196642

Country of ref document: HK

CP01 Change in the name or title of a patent holder
CP01 Change in the name or title of a patent holder

Address after: 518083 comprehensive building, Beishan Industrial Zone, Yantian District, Guangdong, Shenzhen

Patentee after: BGI SHENZHEN

Patentee after: Shenzhen Huada Gene Technology Co., Ltd.

Address before: 518083 comprehensive building, Beishan Industrial Zone, Yantian District, Guangdong, Shenzhen

Patentee before: BGI SHENZHEN

Patentee before: Shenzhen Huada Gene Technology Co., Ltd.