环境微生物群宏基因组测序中序列覆盖度的置信测定方法

1.本发明涉及生物工程技术领域，尤其是一种环境微生物的测定方法。


背景技术：

2.宏基因组(metagenome)的定义为“the genomes of the total microbiota found in nature”,即环境中所有微生物基因组的总和。宏基因组不包含对某个特定微生物种群的靶向，而是针对所有微生物基因组的总和。参见附图1，现有技术中，经典的宏基因组测序分为四步：样品制备、文库构建、上机测序和数据分析。样品制备一般由两步组成，分别为样品收集和dna提取，这一步骤需要极力避免污染，尽量确保使用的所有试剂的“无菌”状态。由于二代测序的高敏感性，在dna文库中即便极低的dna含量也会被扩增并测序，一旦引入污染菌群，就会对样品中真实的信号产生覆盖。文库构建&上机测序可以根据测序平台的选择采用不同的建库方案及上机测序流程。对于数据分析，宏基因组流行的分析方法主要分为两种，即测序序列分类(read classification)和宏基因组组装(metagenomic assembly)。前者为宏基因组测序结果与数据库中已知微生物基因组的比对，按照比对结果对reads进行分类，并根据各微生物reads的相对丰度分析样本中各微生物的相对种群丰度；后者为根据宏基因组测序结果，对微生物基因组数据组装为完整的基因组序列。
3.可见，基因组数据分析通常会报一个覆盖度(coverage)的参数，即基因组被序列覆盖到的区域的比例，对于动植物和单菌来说，覆盖度可以反映出基因组的完整度，覆盖度较低则说明基因组测序深度不够，则需要对基因组进行加测，并观测覆盖度是否得到有效改善，若没有改善则需要用大片段文库进行sanger测序或者nanopore测序补齐其中的缺失区域(gap)。
4.然而，宏基因组则比较特殊，目的是对环境微生物群组(包括人体基因组中的微生物)进行定性，而不是进行定量。对这些标本提取好dna或者rna后，进行高通量基因组测序，然后进行后续分析。如图1所示，在实验环节会引入大量的人工序列(architects)，人工序列主要来自于试剂及文库构建中的空载体序列，这部分序列可以通过通过背景数据库及标准的流程删去，人体细胞测序后得到的序列也可以通过和人类基因组参考序列比对删去。而剩余的序列和微生物基因组参考基因组比对后，如果未有收录，则无法判定其来源，可能是未知的物种，也可能是实验环节带入的污染序列。如参考基因组中有收录这些序列信息，则可以判定为已知物种。本发明主要指向是已知环境微生物序列。


技术实现要素：

5.本发明要解决的技术问题是提供一种环境微生物群宏基因组测序中序列覆盖度的置信测定方法。
6.为解决上述技术问题，本发明所采取的技术方案如下。
7.环境微生物群宏基因组测序中序列覆盖度的置信测定方法，其特征在于：
8.该方法首先对参考基因组进行随机化区域划分，然后将测序序列和基因组进行比
对，对覆盖到的随机区域进行占比统计得到随机覆盖度，最后基于随机覆盖度数据进行概率置信测定。
9.作为本发明的一种优选技术方案，该方法包括如下实施步骤：
10.a、参考基因组的构建和随机化区域划分
11.将参考基因组分类后进行随机化分，然后按照统一的文本格式进行存储，以便于后续数据处理过程中随时调取；
12.b、随机覆盖度的测定
13.将测定序列和基因组进行比对，统计覆盖到的随机区域的占比，并按照随机批次的划分进行若干次重复计算，最后取平均值，得到随机覆盖度；
14.c、随机置信度的测定
15.确认环境样品中微生物的测序序列所在基因组坐标位置的随机分布性，基于此，对于随机性区域进一步划分为去除重复的区域和所有区域，并对这两种区域分别进行计算测定，最后设定量化可信度指标pr，pr的范围为0-1，越接近1，说明环境样品中存在目标微生物的可信度越高；pr的表达式为：
16.pr＝(cov～{mcov}i)*sigmoid(cov)；
17.其中，cov为当次抽样随机覆盖度，{mcov}i为在序列数为i的情况下通过蒙特卡洛模拟得到的模拟值集合，(cov～{mcov}i)为密度函数，即此次测序的随机覆盖度，落入{mcov}i中的可能性，为0-1之间的值；sigmoid(cov)为sigmoid惩罚函数。
18.作为本发明的一种优选技术方案，该方法包括如下实施步骤：
19.a、参考基因组的构建和随机化区域划分，包括如下分步骤：
20.a-1、参考基因组包括完整基因组以及碎片基因组，将碎片基因组和完整基因组比对后进行排序；
21.a-2、通过环化处理直接将其首尾相联；
22.a-3、将环化后的基因组随机化分，并照统一的文本格式进行存储，以便于后续数据处理过程中随时调取；
23.b、随机覆盖度的测定
24.将测定序列和基因组进行比对，统计覆盖到的随机区域的占比，并按照随机批次的划分进行若干次重复计算，最后取平均值，得到随机覆盖度；计算公式为：
25.c＝{n}/nm；
26.其中，n为覆盖到的随机区域的个数，n为随机批次的个数，m为随机区域的个数；
27.c、随机置信度的测定
28.测序过程中将遗传物质随机打断成大小一致的片段并对基因组片段进行随机抽样，由此确认环境样品中微生物的测序序列所在基因组坐标位置的随机分布性，基于此，对于随机性区域进一步划分为去除重复的区域和所有区域，并对这两种区域分别进行计算测定，并基于随机覆盖度集合不符合正态分布，采用非参数检验的落入区间概率,同时对随机覆盖度较低的随机覆盖度引入sigmoid惩罚打分，最后得到量化可信度指标pr，pr的范围为0-1，越接近1，说明环境样品中存在目标微生物的可信度越高；pr的表达式为：
29.pr＝(cov～{mcov}i)*sigmoid(cov)；
30.其中，cov为当次抽样随机覆盖度，{mcov}i为在序列数为i的情况下通过蒙特卡洛
模拟得到的模拟值集；(cov～{mcov}i)为密度函数，即此次测序的随机覆盖度，落入{mcov}i中正常区间可能性，为0-1之间的值，sigmoid(cov)为sigmoid函数。
31.我们假设某次测序发现，环境中物种的序列数为28，对于这个公式的详细描述如下所示：
32.{mcov}i即为{mcov}
28
,为一个500个模拟随机覆盖度的集合。如下表所示，平均值为0.76，最大值为0.95，最小值为0.15，标准偏差为0.07，出现频率最高的区间为0.75-0.80。异常区间的检测使用turkey’s test,即为q3-k(q3-q1),q1为下四分位数，q3为上四分位数，k取为3(极度异常)。如下表所示，其中0.15为异常值。
[0033][0034]
(cov～{mcov}i)为，此次的随机覆盖度，落入{mcov}
28
这个数据集合中的可能性，返回值为0到1之间的小数。0.15-0.55为异常区间。计算方式如下表所示：
[0035][0036]
sigmoid(cov)为固定的函数值，如下表所示。
[0037][0038]
作为本发明的一种优选技术方案，步骤a-3中，所述文本格式统一按照下表进行设定：
[0039][0040]
作为本发明的一种优选技术方案，文本内容提前存储在一个文本文件中，所述文本文件命名为：randomregion.txt。
[0041]
作为本发明的一种优选技术方案，开始坐标和结束坐标都以线性基因组的原始标尺作为为基线标尺。
[0042]
采用上述技术方案所产生的有益效果在于：
[0043]
本发明采用蒙特卡洛模拟研发了一种全新的置信参考，能够用于展示宏基因组测
序报告中各类环境微生物基因组的随机覆盖度，随机覆盖度越高，则说明结果可靠度越高。
[0044]
本发明提出了随机覆盖度和随机置信度的概念，能够用于矫正环境微生物群组宏基因组测序中的假阳性问题，为环境微生物检测的可靠度提供了高置信参考。
附图说明
[0045]
图1为现有技术中常规性的宏基因组测序步骤。
[0046]
图2为本发明参考基因组序列随机区域划分示意图。
[0047]
图3：随机覆盖度测算过程示意图。
[0048]
图4为采用本发明的方法进行模拟测算的示意性结果。
[0049]
图5为采用本发明的方法进行模拟测算的示意性结果。
[0050]
图6为本发明随机区域分类方法示意图。
[0051]
图7为随机区域的分布示意图,包括正态分布和sigmoid函数。
[0052]
图8为随机度和序列数的分布示意图。
具体实施方式
[0053]
以下实施例详细说明了本发明。本发明所使用的各种原料及各项设备均为常规市售产品，均能够通过市场购买直接获得。
[0054]
在以下实施例的描述中，为了说明而不是为了限定，提出了诸如特定系统结构、技术之类的具体细节，以便透彻理解本技术实施例。然而，本领域的技术人员应当清楚，在没有这些具体细节的其它实施例中也可以实现本技术。
[0055]
应当理解，当在本技术说明书和所附权利要求书中使用时，术语“包括”指示所描述特征、整体、步骤、操作、元素的存在，但并不排除一个或多个其它特征、整体、步骤、操作、元素和/或其集合的存在或添加。
[0056]
还应当理解，在本技术说明书和所附权利要求书中使用的术语“和/或”是指相关联列出的项中的一个或多个的任何组合以及所有可能组合，并且包括这些组合。
[0057]
如在本技术说明书和所附权利要求书中所使用的那样，术语“如果”可以依据上下文被解释为“当...时”或“一旦”或“响应于确定”或“响应于检测到”。类似地，短语“如果确定”或“如果检测到[所描述条件或事件]”可以依据上下文被解释为意指“一旦确定”或“响应于确定”或“一旦检测到[所描述条件或事件]”或“响应于检测到[所描述条件或事件]”。
[0058]
另外，在本技术说明书和所附权利要求书的描述中，术语“第一”、“第二”、“第三”等仅用于区分描述，而不能理解为指示或暗示相对重要性。
[0059]
在本技术说明书中描述的参考“一个实施例”或“一些实施例”等意味着在本技术的一个或多个实施例中包括结合该实施例描述的特定特征、结构或特点。由此，在本说明书中的不同之处出现的语句“在一个实施例中”、“在一些实施例中”、“在其他一些实施例中”、“在另外一些实施例中”等不是必然都参考相同的实施例，而是意味着“一个或多个但不是所有的实施例”，除非是以其他方式另外特别强调。术语“包括”、“包含”、“具有”及它们的变形都意味着“包括但不限于”，除非是以其他方式另外特别强调。
[0060]
实施例1、参考基因组的构建和随机化区域划分
[0061]
参考基因组一般分为两种，包括完整基因组以及碎片基因组。1)首先要将碎片基
因组和参考基因组比对后进行排序。2)直接将首尾相联。3)将基因组随机化分，并按照如下表1的格式，提前存储在一个randomregion.txt的文件中，方便随时调取。
[0062]
表1随机区域的格式
[0063][0064]
因为每个基因组会进行多次随机化区域处理，因此不同随机批次区域会有交叉区域。尽管随机区域的名称相同，但是实际上是不同的区域，另外，开始坐标和结束坐标都以线性基因组的原始作为为基线标尺，参见图2。
[0065]
实施例2、随机覆盖度的计算
[0066]
将序列和基因组进行比对后，统计覆盖到的随机区域的占比，并按照随机批次的划分计算多次，最后取平均值为随机覆盖度，计算过程如图3所示。公式为：
[0067]
c＝{n}/nm
[0068]
其中，n为覆盖到的随机区域的个数，n为随机批次的个数，m为随机区域的个数。
[0069]
我们采用了大肠杆菌escherichia coli o157:h7 str.sakai,新冠病毒sars-cov-2,白假丝酵母candida albicans sc5314，猪肉绦虫taenia solium mex_genome_complete.1-6-13，分别进行模拟数据测试，如表二所示。将基因组分为30个线性区域，并按照2-100的序列梯度，进行随机抽样，每个梯度重复500次。模拟结果如表2所示。
[0070]
表2模拟结果ⅰ[0071][0072]
表3模拟结果ⅱ[0073][0074]
从表3可以看出，序列数达到50时，就可以覆盖92％以上的区域，随机覆盖度为92
以上。序列数到达100，就可以覆盖99％以上的区域，随机覆盖的为99以上，且和基因组的大小无相关性，如图4所示。
[0075]
在本次宏基因组样品中，人体的基因组占据绝大多数，因此，我们的假设是，微生物体的大部分基因组应该存在总体数据集中，而测序的随机抽取值有限，只能随机挑选部分微生物体序列进行测序，如果随机抽取值加大，基因组的真实覆盖度也应该会上升。我们使用了真实数据集来验证这一猜想(真实实验的数据集说明如上表3所示)。
[0076]
同时，我们从测序序列中，按照10％-100％的比例进行梯度抽取，观测真实覆盖度的增长区间，结果如图5所示。
[0077]
综上可见，本次真实实验的模拟结果能够验证我们的猜想，即随机抽取值加大，基因组的真实覆盖度也应该会上升，说明微生物体全基因组应该是在宏基因组中存在的。标本c和标本d中，因为序列较多，因此尽管抽取10％的序列，随机覆盖度仍然是1，即可以完全覆盖整个基因组。
[0078]
实施例3、随机置信度的计算
[0079]
基因组测序是将遗传物质随机打断成大小一致的片段，对基因组片段进行随机抽样，因此我们的假设为：如果环境样品中存在这个微生物体，那么测序序列所在的基因组坐标位置，应是随机分布的。而随机性区域又可以分为两部分，即去除重复的区域，以及所有区域，这里我们分别进行计算。如图6所示。
[0080]
我们在这里定义了pr为，在抽取序列数固定为n的情况下本次抽样得到的随机覆盖度落入模拟随机覆盖度区间的概率；模拟随机覆盖度集合既为通过蒙特卡洛模拟得到的模拟值。pr的范围为0-1，越接近1，说明可信度越高。
[0081]
p＝(cov～{mcov}i)*sigmoid(cov)
[0082]
其中，cov为当次抽样随机覆盖度，{mcov}i为在序列数为i的情况下通过蒙特卡洛模拟得到的模拟值集。(cov～{mcov}i)为密度函数，即此次测序的随机覆盖度，落入{mcov}i中的可能性，为0-1之间的值，sigmoid(cov)为sigmoid函数。
[0083]
我们假设某次测序发现，环境中物种的序列数为28，对于这个公式的详细描述如下所示：
[0084]
{mcov}i即为{mcov}
28
,为一个500个模拟随机覆盖度的集合。如下表所示，平均值为0.76，最大值为0.95，最小值为0.15，标准偏差为0.07，出现频率最高的区间为0.75-0.80。异常区间的检测使用turkey’s test,即为q3-k(q3-q1),q1为下四分位数，q3为上四分位数，k取为3(极度异常)。如下表所示，其中0.15为异常值。
[0085][0086]
(cov～{mcov}i)为，此次的随机覆盖度，落入{mcov}
28
这个数据集合中的可能性，返回值为0到1之间的小数。0.15-0.55为异常区间。计算方式如下表所示：
[0087][0088]
sigmoid(cov)为固定的函数值，如下表所示。
[0089][0090]
我们使用表3中的标本b，首先构建梯度抽样(2到84)下的随机覆盖度集合，如图7所示。可以看出随机覆盖度集合，都不符合正态分布。因此我们使用了非参数检验的落入区间概率,同时我们需要对随机覆盖度较低的随机覆盖度，引入sigmoid惩罚打分。
[0091]
我们仍然使用了表3中的标本b，按照公式计算，分布如图8所示，可以看出，随着抽样序列数的增多，随机置信度和随机覆盖度都协同增长，符合预期。
[0092]
综上实施例可见，本发明提出了随机覆盖度和随机置信度的概念，为首个可以直观量化环境微生物群体宏基因组测序中微生物个体可信度的算法，基于参考基因组的管理和更新以及随机值的构建，可以有效矫正环境微生物群体宏基因组测序中的假阳性问题。
[0093]
在上述实施例中，对各个实施例的描述都各有侧重，某个实施例中没有详述或记载的部分，可以参见其它实施例的相关描述。
[0094]
以上所述实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分技术特征进行等同替换；而这些修改
或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的精神和范围，均应包含在本发明的保护范围之内。