CN103186716A

CN103186716A - 基于元基因组学的未知病原快速鉴定系统及分析方法

Info

Publication number: CN103186716A
Application number: CN2011104526667A
Authority: CN
Inventors: 韦朝春; 蔡锴晔; 宣黎明; 贾鹏; 刘雷; 李亦学
Original assignee: SHANGHAI CENTER FOR BIOINFORMATION TECHNOLOGY
Current assignee: SHANGHAI CENTER FOR BIOINFORMATION TECHNOLOGY
Priority date: 2011-12-29
Filing date: 2011-12-29
Publication date: 2013-07-03
Anticipated expiration: 2031-12-29
Also published as: CN103186716B

Abstract

本发明提供一种基于元基因组学的未知病原的快速鉴定系统及分析方法，运用元基因组学的方法进行未知病原鉴定及分析：包括：测序质量控制系统：对测序数据进行质量分析，把测序质量不合要求的部分去除；Meta-All分析系统：本系统的主体部分,从测序数据分析元基因组中的物种组成结构；NeSSM第二代测序模拟系统：独立部分可以产生模拟的测序数据；病原数据库系统自动更新Meta-All分析系统需要用到的数据库资源；病原组成比较分析单元：根据物种组成表比较不同来源的元基因组病原组成结构的异同，找出与疾病相关性高的病原。本发明的优点在运行元基因组学的方法可以迅速缩小未知病原鉴定的范围，从而加速鉴定的速度，同时对完全未知的病原可以从相似病原进行分析。

Description

基于元基因组学的未知病原快速鉴定系统及分析方法

技术领域

本发明涉及未知病原的分析，具体涉及一种基于元基因组学的未知病原分析系统及分析方法。

背景技术

由于气候、生态和人类行为等各方面因素的变化，出现了许多未知的病原体引发的传染病，例如尼帕病毒、SARS等。甚至在常见的传染病爆发时，仍有相当一部分病人无法确定病原。因此，快速检测未知病原的平台显得尤为重要，它能在传染病爆发的早期帮助确定或者是缩小可疑病原微生物的范围，以便更迅速做出反应，部署有针对性的防控措施，减少疾病流行对社会和经济造成的损失。

相对于目前已有的病原检测方法，用元基因组学的方法来推断微生物的组成，显著的特点是不依赖于培养技术和对微生物的先验知识。在自然界中，不能培养的微生物占绝大多数，因此元基因组学方法有希望揭示真实的微生物群落组成和各微生物之间相互作用。随着核酸测序技术和生物信息学的发展，用基于元基因组学的方法来快速鉴定样本的物种组成的技术逐渐成熟，使得基于元基因组学的未知病原寻找系统成为可能。

关于Meta-All分析系统，主要包括以下功能单元，测序数据的两种处理方式之一Meta-All，该方式是基于测序序列与已知基因组序列比对得到结果；测序数据的另一种处理方式Meta-BinG，该方式是分析测序数据的K-mer组成得到结果，基于比对的方式针对常见的两种不同测序技术应用不同的序列比对程序，solexa序列用Bowtie做比对，基于比对的方式针对常见的两种不同测序技术应用不同的序列比对程序，454序列用Blat做比对，比对结果修正单元，若为16s rRNA则基于拷贝数进行修正，若为全基因组WGS则基于基因组大小进行修正，结果统计输出单元，在门、纲、目、科、属、种这6个层次输出相应的物种组成表和饼图

发明内容

本发明的第一目的是提供一种基于元基因组的未知病原快速鉴定系统。

本发明的第二目的是提供一种基于元基因组的未知病原快速鉴定系统的分析方法。

一种基于元基因组的未知病原快速鉴定系统，其包括：

质量控制系统：对测序数据进行质量分析，把测序质量不合要求的部分去除；

Meta-All分析系统：本系统的主体部分；

NeSSM第二代测序模拟系统：独立部分，可以产生模拟的测序数据；

数据库更新组件：自动更新Meta-All分析系统需要用到的数据库资源；

病原组成比较分析单元：根据物种组成表，比较不同来源的元基因组病原组成的异同，找出与疾病相关性高的病原。

进一步地，所述数据库更新组件包括以下单元：

下载单元：查询公共数据库里新增和修改过的数据，并下载到本地；

数据处理单元：把已下载的数据根据类型整理并转化为本地数据库需要的格式；

数据记录单元：把格式化后的数据写入对应的本地数据库，把更新时间记入日志；

更新控制单元：检查日志文件，判断是否需要下载数据更新本地数据库。

一种使用如所述的基于元基因组的未知病原快速鉴定系统的未知病原分析方法，所述方法包括如下步骤：

首先，利用数据库更新组件进行以下数据更新步骤

通过下载单元，查询公共数据库里新增和修改过的数据，并下载到本地；

通过数据处理单元，把已下载的数据根据类型整理并转化为本地数据库需要的格式；

通过数据记录单元，把格式化后的数据写入对应的本地数据库，把更新时间记入日志；

通过更新控制单元，检查日志文件，判断是否需要下载数据更新本地数据库。

然后，进行以下分析步骤

步骤I：质量控制；

用测序质量值Q表示某个碱基测序质量，若该碱基的错误概率p，p为(0，1)之间的一个实数，则Q为p的一个整数映射，使用Q值作为质量控制的依据，设定碱基水平和序列水平的质量控制的阈值Q均为20，所对应的碱基错误概率p大约为0.01；

步骤II：物种分析；

经过质量控制的序列即可进行下一步的物种分析，这里的物种分析基于序列相似性的分析；其中包括：

对比：采用了Bowtie作为Meta-All的比对软件，Bowtie采用了Burrows-Wheeler转换方法来处理基因组序列，可以快速地将短序列比对到基因组上；

物种注释：将短序列和基因组序列用bowtie比对，比对结果中包含着基因组序列的核酸GI号，GI号可以唯一地映射到一个NCBI的物种号，在NCBI的Taxonomy数据库中可以得到该物种号在门纲目科属种等不同层次上的物种名称，GI号和物种号的映射关系以及物种号的注释信息都可以在NCBI的FTP下载；

步骤III：元基因组模拟测序系统；

设计一个元基因组学模拟测序系统来产生模拟的元基因组学测序数据，模拟测序系统的基础是目前已经被完整测序的细菌基因组，给定细菌的组成(细菌的名称和细菌在种群中所占的百分比)、测序方式(全基因组/16S)、测序技术(454/Solexa)和测序数量，模拟测序系统会产生对应的模拟测序数据集；

步骤IV：基于16S核糖体RNA基因序列的物种分析的校正；

对于元基因组学16S基因序列的物种分析结果需要考虑不同基因组的16S拷贝数差异：

P_{i}^{'} = \frac{P_{i}}{C_{i}}

这里，Pi表示第i个物种在种群中所占的百分比，Ci表示第i个物种的16S基因拷贝数，Pi’表示校正后第i个物种在种群中所占的百分比，校正后，PI’之和可能并不为1，需要重新进行标准化；

步骤V：全基因组随机测序序列物种分析的校正；

P_{i}^{'} = \frac{P_{i}}{U_{i} * G_{i}}

这里，Pi表示第i个物种在种群中所占的百分比，Ui表示第i个细菌基因组的物种特异片段占整个细菌基因组的百分比，Gi表示第i个细菌基因组的大小，Pi’表示校正后第i个物种在种群中所占的百分比，校正后，PI’之和可能并不为1，需要重新进行标准化；

步骤VI：全基因组随机测序序列的Binning；

利用GPU加速元基因组学分析的系统MetaBinG，MetaBinG使用了k阶马尔科夫模型分析元基因组学数据，k阶马尔科夫的状态是指k个寡核苷酸，一个状态之后连接着四个状态，前一个状态的后k-1个字符和后一个状态的前k-1个字符是相同的，因此每一个字符都被视为是由前k个字符决定的，具体地说，在第i个基因组中，从状态m到状态n的概率可以由如下的公式计算，该概率也可以理解为观测到状态m的k个字符后下一个字符出现状态n的最后一个字符的概率，

{kMM}_{i, mn} = P_{i} (O_{m} | O_{n}) = \frac{F_{i} (O_{m} | O_{n})}{F_{i} (O_{m})}

这里Om和On是长度为k的寡核苷酸，Pi(Om|On)在第i个基因组中是从状态Om转移到状态On的概率，Fi(Om|On)是在第i个基因组中观测到从状态Om转移到状态On的频数，Fi(Om)是在第i个基因组中观测到状态Om的频数；

步骤VII：多样本物种比较；

找到在不同表型间分布有着显著差异的物种，实现未知病原的鉴定。

在基于元基因组学的未知病原分析方法的步骤VII多样本物种比较中，具体包括：

种群自身复杂程度比较；

种群自身的复杂程度属于α-多样性的范畴，常用的α-多样性的研究方法是计算种群的Shannon系数：

H = - Σ_{i = 1}^{S} (p_{i} \ln p_{i})

其中，Pi是指第i个物种在种群中所占的百分比，种群中，物种数量S越大，物种分布越均匀，Shannon系数也越大，特别的，当物种数量S一定时，完全均匀分布的种群(即每个物种占种群的百分比都是1/S)的Shannon系数有最大值lnS，

通过Shannon系数的计算，可以比较不同样本的微生物种群在复杂程度上的差异；

种群间距离的计算；

种群间的距离研究属于β-多样性的范畴，这里我们用Bary-Curtis距离度量种群间的距离，计算i种群和j种群的Bary-Curtis距离：

{BC}_{ij} = \frac{Σ_{k = 1}^{n} | x_{ik} - x_{jk} |}{Σ_{k = 1}^{n} (x_{ik} + x_{jk})}

这里取i种群和j种群物种的并集作为全集A，Xik，Xjk分别表示i种群和j种群中k物种所占的百分比，k∈A，

Bary-Curtis距离作为一种种群间距离的度量，可以定量地度量样本和样本之间的相似程度，用于样本的聚类；

寻找分布显著差异的物种；

样本中每个物种的多少是以匹配到该物种的短序列条数来度量的，由于每个样本的测序深度不同，为了比较不同样本的物种差异，需要将各个样本的种群结构从序列数量形式转换为百分比形式，即用每个物种对应的短序列数量除以总体序列数量，随后可以对每一个物种进行假设检验，寻找在不同组别里物种丰度存在显著差异的物种，致病的原因可能是由于某个物种的异常增多造成，也可能是由于某物种的异常减少引起，可以用oneway-ANOVA来检验。

本发明的技术原理如下：要利用元基因组学技术来研究未知病原，首先需要判定元基因组测序序列的物种来源。这里我们利用相似性比对的方法对元基因组测序序列进行物种分析。目前，NCBI的Genome数据库已经收录了超过一千个已经被完整测序的细菌基因组(见下文数据集1)。我们将待分析的元基因组序列比对到细菌全基因组数据库，如果某条元基因组序列在预先设定的匹配标准下能够唯一地匹配到某个物种，我们认为该序列来自于这个物种。如果某条元基因组序列能够匹配到多个物种，我们将该序列归类到这多个物种的在物种树上的最低父节点。这里采用的物种分类体系是NCBI的物种分类体系。

根据测序类型的不同，元基因组序列可分为全基因组随机测序和16S核糖体RNA特征基因测序。我们也提出了两个方法Meta-All-MSS(Metagenome analysis system based on Allavailable genomes using Metagenome Shotgun Sequencing)和Meta-All-16S(Metagenomeanalysis system based on All available genomes using 16S rRN Agene sequencing)，分别对应上述两种元基因组学测序方式。我们专门设计了3个数据集来测试Meta-All是否能准确还原出样本的细菌种群结构。从结果来看，Meta-All-MSS和Meta-All-16S都能显著的改善样本细菌种群还原的准确性。

但另一方面，通常元基因组中包含着大量的未知物种，用比对的办法很难将所有序列很好地归类。而基于k-mer的方式通过分析元基因组，可能并不能将元基因组序列准确地每一条序列都分类到种和属的水平，但是却可以在门和纲的水平上大致地推测出每一条序列的物种来源。

目前最好的基于组成的元基因组学物种分析系统Phymm使用的是插值马尔科夫模型，在模拟的测试中Phymm被证明有着最高的敏感性和特异性。但是Phymm的计算量非常大，处理新一代测序技术产生的海量数据需要相当长的时间。

我们提出了一种基于元基因组学的未知病原分析系统，其中包括数据库更新组件，该数据库更新组件可以自动更新Meta-All分析系统需要用到的数据库资源，本发明之所以可以用作未知病原鉴定，是因为我们整理、格式化并维持了数据库的更新，使之后的比对分析有了基础。

附图说明

图1为本发明的整体架构图。

图2为Meta-All分析系统的架构图。

图3为本发明的数据库更新组件的架构图。

具体实施方式

以下通过特定的具体实例说明本发明的实施方式，本领域技术人员可由本说明书所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。

数据集：

数据集1.现有细菌基因组数据库

在2011年10月31日，我们从NCBI的FTP站点下载的1483个细菌基因组，2422个病毒基因组和28个真菌基因组。

数据集2.细菌的16S核糖体RNA的基因序列

下载GenBank格式的细菌基因组，用正则表达式查找其中的16S核糖体RNA基因注释，根据注释信息从数据集1中将16S序列提取出来。

数据集3.口腔元基因组学16S核糖体基因测序序列

该数据集对应的样本来自60个杭州幼儿园的3岁到6岁的儿童，其中34个男孩，26个女孩。这些样本分为4个组：MN(男孩，无龋齿，n＝17)，MC(男孩，有龋齿，n＝17)，FN(女孩，无龋齿，n＝11)，FC(女孩，有龋齿，n＝15)。每个儿童均采集了唾液样本和牙菌斑样本。该数据集共包含186,787条16S核糖体RNA基因的V3区序列。

首先，参见图3，数据库更新组件自动更新Meta-All分析系统需要用到的数据库资源，本发明之所以可以用作未知病原鉴定，是因为我们整理、格式化并维持了数据库的更新，使之后的比对分析有了基础；其中，数据库更新组件包括以下单元：

然后，根据以下步骤进行分析：

I.质量控制

通常用测序质量值Q表示某个碱基测序质量。若该碱基的错误概率p，p为(0，1)之间的一个实数，则Q为p的一个整数映射。对于Sanger测序，这个映射关系是：

Q_ssnger＝-10log₁₀p 公式(1)

而对于Solexa/Illumina测序这个映射关系是：

Q_{solexa} = - {10 \log}_{10} \frac{p}{1 - p}

公式(2)

尽管映射关系不同，但是当p＜0.05的时候，两个公式计算出来的Q几乎是相同的。为了方便起见，后面都使用Q值作为质量控制的依据。这里设定碱基水平和序列水平的质量控制的阈值Q均为20，所对应的碱基错误概率p大约为0.01。

由于目前有多个主流的测序平台，各个测序平台产生的序列并不相同，对不同平台的测序数据进行的质量控制流程和方法并非完全相同。

1.Solexa/Illumina序列

常见的存储生物序列的格式是FASTA格式，而Solexa/Illumina序列采用FASTQ格式存储序列信息。FASATA序列仅包含序列信息，FASTQ序列则还包含各位点的质量信息。下面文本框展示了一条典型的FASTQ序列。

FASTQ格式中第一行和第三行是序列的标识符信息，第二行的GGGTGATGGCCGCTGCCGATGGCGTCAAATCCCACC表示该序列的碱基信息，而第四行的IIIIIIIIIIIIIIIIIIIIIIIIIIIIII9IG9IC表示这些碱基的测序质量。质量值用字符表示，某个字符ch和质量Q的换算公式如下：

Q = 10 * \log_{10} (1 + 10^{\frac{ord (ch) - 64}{10}})

公式(3)

其中，ord(ch)表示字符ch的ASCII值。按照上述公式可以将FASTQ格式中的碱基质量由字符形式转换为数字形式，进行下一步的质量控制。

2.454序列

454序列存储为两个文件：一个FASTA文件和一个QUAL文件。FASTA文件存储碱基序列，QUAL文件存储碱基质量，如下所示：

Example.fa

Example.qual

由于454序列的碱基质量已经是数值形式，可以直接进行质量控制。

这里我们将质量分为碱基水平的质量控制和序列水平的质量控制，前者是指去除原始序列头部和尾部质量低于设定阈值的碱基，后者是指去除平均碱基质量低于设定阈值的序列。原始的测序序列需要先通过碱基水平的质量控制，再通过序列水平的质量控制才能进入后续的生物信息学分析，使得后续的分析更加可信。

II.物种组成分析

经过质量控制的序列即可进行下一步的物种组成分析，这里的物种组成分析基于序列相似性的分析。

1.比对

新一代测序技术可以在一次实验中产生海量的测序数据。由于传统的比对程序如BLAST和BLAT出现于Sanger测序技术占统治地位的时代，用BLAST或者BLAT来比对高通量的测序数据需要消耗很长的时间。随着新一代测序技术的发展而产生的短序列比对程序如Bowtie，Maq，和Bwa等能快速的将短序列比对到基因组上。在这些软件中，Bowtie是最快的，出于速度上的考虑，我们采用了Bowtie作为Meta-All的比对软件。

Bowtie采用了Burrows-Wheeler转换方法来处理基因组序列，可以快速地将短序列比对到基因组上。具体而言，Bowtie每小时可以将25,000,000条短序列比对到人类基因组上。Bowtie也可以同时使用多个CPU来加速比对过程。

用Bowtie比对需要先用Bowtie工具包中的build-index工具为基因组文件建立索引。由于测序平台和测序方式的不同，我们需要处理不同存储形式的序列，因此用Bowtie将短序列比对到基因组需要采用不同的参数，具体参数设置如下：

a.单端Solexa/Illumina序列，FASTQ格式

--best--strata--tryhard--all

b.单端Solexa/Illumina序列，FASTA格式

--best--strata--tryhard--all-f

c.双端Solexa/Illumina序列，FASTQ格式

--tryhard--all-1-2

d.双端Solexa/Illumina序列，FASTA格式

--tryhard--all-f-1-2

e.单端454序列

--best--strata--tryhard--all-f

参数说明：

a.--best--strata：如果短序列存在许多匹配结果，则只输出最好的一个结果；

b.--tryhard：尽可能找到最准确的匹配，这个选项会降级比对的速度；

c.-f：输入文件是一个fasta文件

d.参数说明：

Bowtie的比对结果如下所示：

Bowtie的比对结果每一行表示一个匹配结果；比对结果有多列，每列用制表符隔开，从左到右各列分别表示：

a.短序列名称

b.匹配到参考序列上的是正链还是负链

c.参考序列的名称

d.匹配相对于参考序列的起始位置

e.短序列

f.短序列的质量

g.-M选项阈值(若短序列和基因组满足条件的匹配数量超过阈值，则随机报告其中一个匹配)；

h.该匹配结果的错配碱基

2.物种注释

每一条读序列可以被视为从元基因组中进行的一次随机采样，如果能判断出每一条序列来自于哪种细菌，就可以用匹配到各个物种的Read数目大致推测出元基因组学的物种组成结构。

将短序列和基因组序列用bowtie比对，比对结果中包含着基因组序列的核酸GI号。GI号可以唯一地映射到一个NCBI的物种号，在NCBI的Taxonomy数据库中可以得到该物种号在门纲目科属种等不同层次上的物种名称。GI号和物种号的映射关系以及物种号的注释信息都可以在NCBI的FTP下载(ftp://ftp.ncbi.nih.gov/pub/taxonomy/)。

在某个层次上，如果某条元基因组学序列在预先设定的匹配标准下能够唯一地匹配到某个物种，我们认为该序列来自于这个物种，而如果某条元基因组学序列能够匹配到多个物种，我们将该序列归类到这多个物种的在物种树上的最低父节点。

III.元基因组模拟测序系统

由于现在缺少元基因组学分析的金标准数据，因此我们设计了一个元基因组学模拟测序系统来产生模拟的元基因组学测序数据。模拟测序系统的基础是目前已经被完整测序的一千多细菌基因组。给定细菌的组成(细菌的名称和细菌在种群中所占的百分比)、测序方式(全基因组/16S)、测序技术(454/Solexa)和测序数量，模拟测序系统会产生对应的模拟测序数据集，该元基因组模拟测序系统NeSSM可独立运行，用于本方法的评估和作为补充单元。

IV.16S核糖体RNA基因序列的物种分析的校正

16S核糖体RNA基因是一个在细菌和古菌中都相当保守的基因，并且越相近的物种的16S基因序列越相似，因此，16S核糖体基因序列被作为理想的特征基因运用到系统发育的研究中。

在元基因组学技术发展的早期，由于技术条件的限制，相当多的元基因组学项目采用的方法是用PCR引物扩增16S基因片段，对16S基因片段的分析来推断整个元基因组学的组成。这种方法最初取得了很大的成功，对16S基因的研究揭示了自然界中存在的细菌数量远远大于我们之前所了解的细菌数量。

由于新一代测序技术所提供的序列读长较短，其中最长的454技术目前也仅能提供400bp左右的读长，而16S核糖体RNA大约有1500bp。16S核糖体基因中含有八个保守区，九个可变区。因此目前基于16S的研究思路是使用PCR引物扩增16S的一个或几个可变区。目前研究已经证明，恰当的16S可变区序列可以和全长的16S序列有着相近的分类准确性。在新一代测序技术的帮助下，16S序列的数据正在以前所未有的速度增长。

但是由于16S基因自身的一些性质，其用于元基因组学的研究也受到了一些人的批评。

1.16S基因序列用于物种分类的假设是越相近的物种16S基因序列也越相似，但是在一些情况下并不是这样的，甚至于在一个基因组上的多个16S基因也并不完全相似，因此Pontes等认为不能将16S分配到具体的物种；

2.16S序列在基因组上存在着多个拷贝，并且不同基因组上拷贝数目并不一样；

3.PCR可能会为后续的分析带来误差，由于PCR引物设计的问题，可能某些物种并不能够被PCR引物扩增出来，也有可能另一些物种和PCR引物的亲和力较好而获得了更多的扩增产物。综上所述，通过对16S基因序列的分析推断出的细菌群落组成可能同实际情况存在着差异。

尽管用16S基因做元基因组学研究存在着上述不足，但是16S基因仍然是目前元基因组学研究的主流方法。鉴于此，这里对基于16S核糖体基因的元基因组学研究做了一些改进。

截止2010年12月，NCBI的Genome数据库中共包含1212个已经完整测序的细菌基因组。将这些基因组以Genbank的格式下载(ftp://ftp.ncbi.nih.gov/genomes/Bacteria/all.gbk.tar.gz)，用正则表达式统计每个基因组上的16S注释信息。这1212个基因组上，16S核糖体基因的拷贝数最少1个，最多16个，中值为4个。由于各个基因组上的16S拷贝数差异是如此的大，因此，我们认为不考虑16S拷贝数的元基因组学分析不能体现元基因组学的实际情况。

按照NCBI的物种分类体系，对每个物种的16S rRNA基因拷贝数进行统计，在门的水平上的统计结果见下图，其中Firmicutes，Fusobacteria和Proteobacteria门的平均拷贝数均超过了5，而各门的平均拷贝数仅为2.43个，这个结果说明如果不考虑在不同基因组上的16S rRNA基因的拷贝数差异的话，分析得到的结果是不准确的。

现有的方法通常以16S基因的组成结构来代替物种的组成结构，但是基于上述的结果，这样得到的物种组成是粗糙的。我们认为对于元基因组学16S基因序列的物种分析结果需要考虑不同基因组的16S拷贝数差异：

P_{i}^{'} = \frac{P_{i}}{C_{i}}

公式(4)

这里，Pi表示第i个物种在种群中所占的百分比，Ci表示第i个物种的16S基因拷贝数，Pi’表示校正后第i个物种在种群中所占的百分比。校正后，PI’之和可能并不为1，需要重新进行标准化。

V.全基因组随机测序序列物种分析的校正

由于新一代测序技术提供的读长较短，一条短序列仅能覆盖基因组上很短的一个区间。因此，相当比例的短序列可能匹配到多个细菌基因组上，由于并不能确认这些短序列的确切来源，进行物种统计的时候并没有考虑这一部分的短序列，而仅仅统计那些可以唯一匹配到一个细菌基因组的短序列。换一种说法，某个细菌的在细菌群落的多少，按照现有的分析方法，取决于细菌基因组上物种特异片段的大小，由于每一个细菌基因组上的物种特异片段大小并不相同，这会给分析结果带来一定的误差，导致从元基因组学测序序列中还原出来的细菌种群结构和实际的细菌种群结构并不完全相同。

基于上述的理由，全基因组随机测序序列的物种分析仍然需要进行一定的校正。校正主要考虑两个因素：

1.每一个细菌基因组的物种特异片段占整个细菌基因组的百分比，该百分比越大，能够鉴定到该细菌的短序列就越多；

2.细菌基因组的大小，细菌基因组越大，能鉴定到该细菌的短序列也越多。因此对全基因组随机测序序列的物种分析需要进行的校正如下：

P_{i}^{'} = \frac{P_{i}}{U_{i} * G_{i}}

公式(5)

这里，Pi表示第i个物种在种群中所占的百分比，Ui表示第i个细菌基因组的物种特异片段占整个细菌基因组的百分比，Gi表示第i个细菌基因组的大小，Pi’表示校正后第i个物种在种群中所占的百分比。校正后，PI’之和可能并不为1，需要重新进行标准化。

VI.全基因组随机测序序列的Binning

进行元基因组学分析，最基础也是最重要的一个步骤是准确地将元基因组学序列分配到其来源的物种。这个过程叫做元基因组学物种分析。按照分析的手段不同，可以将元基因组学物种分析分为两类：基于比对的方法和基于组成的方法。当元基因组学序列所来自的物种已经被完整测序，基于比对的方法可以得到更加准确的结果。但是，通常元基因组学中包含着大量的未知物种，用比对的办法很难将所有元基因组学序列很好地归类。

而基于组成的元基因组学分析方法，可能并不能将元基因组学序列准确地将每一条序列都精确地分类到种和属的水平，但是却可以在门和纲的水平上大致地推测出每一条物种的来源。

目前最好的元基因组学物种分析系统Phymm使用的是插值马尔科夫模型，在模拟的测试中Phymm被证明有着最高的敏感性和特异性。但是Phymm的计算量非常大，处理新一代测序技术产生的海量数据需要相当长的时间。

我们之前开发的MetaBinG系统是一种利用GPU加速元基因组学分析的系统。结果显示MetaBinG在准确性上和Phymm是可比的，但是MetaBinG比Phymm快至少两个数量级。我们将MetaBinG作为本方法的一条可选分析途径。

VII.多样本物种比较

上述元基因组学物种分析的工作都是为了后续的多样本物种比较打下基础。只有找到在不同表型间分布有着显著差异的物种，才能实现未知病原的鉴定。

1.种群自身复杂程度比较

H = - Σ_{i = 1}^{S} (p_{i} \ln p_{i})

公式(6)

其中，Pi是指第i个物种在种群中所占的百分比。种群中，物种数量S越大，物种分布越均匀，Shannon系数也越大。特别的，当物种数量S一定时，完全均匀分布的种群(即每个物种占种群的百分比都是1/S)的Shannon系数有最大值lnS。

通过Shannon系数的计算，可以比较不同样本的微生物种群在复杂程度上的差异。

2.种群间距离的计算

种群间的距离研究属于β-多样性的范畴，这里我们用Bary-Curtis距离度量种群间的距离。计算i种群和j种群的Bary-Curtis距离：

{BC}_{ij} = \frac{Σ_{k = 1}^{n} | x_{ik} - x_{jk} |}{Σ_{k = 1}^{n} (x_{ik} + x_{jk})}

公式(7)

这里取i种群和j种群物种的并集作为全集A，Xik，Xjk分别表示i种群和j种群中k物种所占的百分比，k∈A。

Bary-Curtis距离作为一种种群间距离的度量，可以定量地度量样本和样本之间的相似程度，用于样本的聚类等。

3.寻找分布显著差异的物种

样本中每个物种的多少是以匹配到该物种的短序列条数来度量的，由于每个样本的测序深度不同，为了比较不同样本的物种差异，需要将各个样本的种群结构从序列数量形式转换为转换为百分比形式(用每个物种对应的短序列数量除以总体序列数量)。

随后可以对每一个物种进行假设检验，寻找在不同组别里物种丰度存在显著差异的物种。致病的原因可能是由于某个物种的异常增多造成，也可能是由于某物种的异常减少引起。可以用oneway-ANOVA来检验。

实施例1.16S rRNA拷贝数影响的校正

1.所用数据集

我们这里使用了Zhang(Zhang，C.et al.Interactions between gut microbiota，host geneticsand diet relevant to development of metabolic syndromes in mice.Isme J 4，232-241，doi：ismej2009112)等人的小鼠16S核糖体基因序列，包括29343条16S核糖体RNA基因序列，将这些序列用Meta-All进行物种分析，可以得到一个小鼠的肠道细菌群落物种组成结构。

2.还原菌种种群结构

a.传统方式

所用引物对如下表所示：

口腔元基因组学样本分析的实施例

口腔细菌和口腔疾病有着密切的关系，龋齿是儿童口腔疾病中最常见的一种。科学研究认为龋齿的发展和病原细菌有着密切的关系。但是，儿童龋齿环境的细菌种群的整体结构并没有被研究过。这里，我们用高通量的编码焦磷酸测序和PCR-DGGE技术来检测儿童的口腔细菌多样性从而得到了数据集5。数据集5共包含120个样本，表型包含三个特征：性别，有龋/无龋，牙菌斑/唾液。3个因素可以把样本分为8个组。

通过Meta-All的分析，我们发现，口腔的细菌多样性比先前报道的更多。

在门水平上对8个组的物种统计得到口腔细菌主要由8个门的细菌组成Bacteroidetes，Firmicutes，Proteobacteria，Actinobacteria，Fusobacteria，Spirochaetes，TM7和SR1。其中，Bacteroidetes和Firmicutes处于优势地位。

计算各个样本的Shannon系数和Simpson系数，并用ANOVA检验八个组的差异，结果显示各个组的多样性水平是相似的(p＞0.05)。

聚类分析显示3个分组因素中对细菌群落组成影响最大的因素是牙菌斑/唾液。

要找到和龋齿相关的病原细菌，需要在更精确的物种层次进行分析。这里我们选择了属这个层次。在这个数据集中，Meta-All共鉴定到超过70个属的细菌。

其中主要的细菌及对应的物种数量见下表：

由于性别因素对于口腔细菌的组成并没有显著的影响，因此这里只研究有龋/无龋和唾液/牙菌斑因素，这两个因素将样本分为有龋-唾液，有龋-牙菌斑，无龋-唾液，和无龋-牙菌斑四个组。ANOVA分析发现其中六个物种在有龋-牙菌斑和无龋-牙菌斑两组中存在着显著差异(p＜0.05)，这六个物种分别是：Streptococcus，Veillonella，Actinomyces，Granulicatella，Leptotrichia，和Thiomonas，说明龋齿的形成和发展与牙菌斑中的这六种细菌密切相关。六种细菌除Thiomonas外，都能找到文献证据的支持(即该细菌和儿童龋齿是相关的)。这个例子充分证明了Meta-All系统的有效性，也说明了用元基因组学来寻找病原的是一条行之有效的方法。

上述实施例仅例示性说明本发明的原理及其功效，而非用于限制本发明。任何熟悉此技术的人士皆可在不违背本发明的精神及范畴下，对上述实施例进行修饰或改变。因此，举凡所属技术领域中具有通常知识者在未脱离本发明所揭示的精神与技术思想下所完成的一切等效修饰或改变，仍应由本发明的权利要求所涵盖。

Claims

1.一种基于元基因组的未知病原快速鉴定系统，其特征在于，包括：

Meta-All分析系统：本系统的主体部分，从测序数据分析元基因组的物种组成结构；

病原数据库系统：病原分析相关的Meta-All分析系统需要用到的数据库资源，包括数据库自动更新组件；

2.根据权利要求1所述的基于元基因组的未知病原快速鉴定系统，其特征在于，所述病原数据库系统包括以下单元：

3.一种使用如权利要求1所述的基于元基因组的未知病原快速鉴定系统的未知病原分析方法，其特征在于：所述方法包括如下步骤：

首先，利用数据库更新组件进行以下数据更新步骤；

通过更新控制单元，检查日志文件，判断是否需要下载数据更新本地数据库；

然后，进行以下分析步骤

步骤I：质量控制；

用测序质量值Q表示某个碱基测序质量，若该碱基的错误概率p，p为(0，1)之间的一个实数，则Q为p的一个整数映射，使用Q值作为质量控制的依据；设定碱基水平和序列水平的质量控制的阈值Q均为20，所对应的碱基错误概率p大约为0.01；去掉所有质量值低于20的碱基；

步骤II：物种分析；

经过质量控制的序列即可进行下一步的物种组成分析，这里的物种组成分析是基于序列相似性的分析；其中包括：

序列比对：采用了Bowtie作为Meta-All的比对软件，Bowtie采用了Burrows-Wheeler转换方法来处理基因组序列，可以快速地将短序列比对到基因组上；

步骤IV：16S核糖体RNA基因序列的物种分析的校正；

P_{i}^{'} = \frac{P_{i}}{C_{i}}

步骤V：全基因组随机测序序列物种分析的校正；

P_{i}^{'} = \frac{P_{i}}{U_{i} * G_{i}}

步骤VI：全基因组随机测序序列的Binning；

{kMM}_{i, mn} = P_{i} (O_{m} | O_{n}) = \frac{F_{i} (O_{m} | O_{n})}{F_{i} (O_{m})}

步骤VII：多样本物种组成结构比较；

找到在不同样本间分布有着显著差异的物种，实现未知病原的鉴定。

4.根据权利要求3的未知病原分析方法，其特征在于，在步骤II和步骤VI之间，还可以包括如下步骤：

步骤III：元基因组模拟测序系统；

设计一个元基因组学模拟测序系统来产生模拟的元基因组学测序数据，模拟测序系统的基础是目前已经被完整测序的细菌基因组，给定元基因组的组成结构(包括细菌的名称和细菌在种群中所占的百分比)、测序方式(全基因组/16S rRNA)、测序技术平台(454/Solexa)和测序数量，模拟测序系统会产生对应的模拟测序数据集；

5.根据权利要求3的未知病原分析方法，其特征在于，在步骤VII多样本物种比较中，具体包括：

种群自身复杂程度比较；

H = - Σ_{i = 1}^{S} (p_{i} \ln p_{i})

种群间距离的计算；

{BC}_{ij} = \frac{Σ_{k = 1}^{n} | x_{ik} - x_{jk} |}{Σ_{k = 1}^{n} (x_{ik} + x_{jk})}

寻找分布显著差异的物种；

样本中每个物种的多少是以匹配到该物种的短序列条数来度量的，由于每个样本的测序深度不同，为了比较不同样本的物种差异，需要将各个样本的种群结构从序列数量形式转换为百分比形式，即用每个物种对应的短序列数量除以总体序列数量，随后可以对每一个物种进行假设检验，寻找在不同组别里物种丰度存在显著差异的物种，致病的原因可能是由于某个物种的异常增多造成，也可能是由于某物种的异常减少引起，可以用oneway-ANOVA来检验不同样本中的物种分布差异。