阅读说明：本技术 一种通过微生物宏基因组对未知样本进行地理定位的方法 (Method for carrying out geographical positioning on unknown sample through microorganism metagenome ) 是由 许灿强 黄丽红 杨文娴 俞容山 于 2019-10-09 设计创作，主要内容包括：本发明涉及一种通过微生物宏基因组对未知样本进行地理定位的方法，其基于微生物宏基因组对未知样本的来源城市进行预测定位，且在对训练样本和未知样本进行数据处理的过程中，对样本中各菌株的丰度进行分级，通过多个门限值将双精度的菌株丰度转换成离散的多元值。对相较于现有的定位方法，本发明预测准确性高。(The invention relates to a method for carrying out geographical positioning on an unknown sample through a microorganism macro genome, which carries out prediction positioning on a source city of the unknown sample based on the microorganism macro genome, grades the abundance of each strain in the sample in the process of carrying out data processing on a training sample and the unknown sample, and converts the abundance of double-precision strains into discrete multivariate values through a plurality of threshold values. Compared with the existing positioning method, the method has high prediction accuracy.)

一种通过微生物宏基因组对未知样本进行地理定位的方法

技术领域

本发明涉及微生物技术领域，具体涉及一种通过微生物宏基因组对未知样本进行地理定位的方法。

背景技术

微生物是地球上数量最多，种类最丰富，分布范围最广泛的生物类群。基于高通量测序的宏基因组学技术不需要进行微生物培养，可以直接对取自环境里的微生物样品进行分析研究。整个流程是先将样品里的DNA提取出来，进行测序，然后通过算法和计算机软件对测序结果进行分析。目前宏基因组学的发展，可以从多个不同的环境样本中快速、准确的获得环境微生物的基因组序列。通过宏基因组测序数据可以进行菌群检测和菌群丰度的定量，分析样本的物种组成和功能组成。宏基因组学技术给病源识别鉴定、溯源分析等带来了新的方法和思路，在食品安全、传染病防控等方面有巨大潜力和发展空间。

未知样本的地理定位是指对于一个未知的微生物样本，通过对该样本的测序数据的分析，定位其地理来源。现有的方法大多是基于16s RNA测序进行微生物来源城市的预测，在小样本量预测的准确性并不理想。

发明内容

针对上述问题，本发明的目的在于提供一种通过微生物宏基因组对未知样本进行地理定位的方法，其准确性高。

为实现上述目的，本发明采用的技术方案是：

一种通过微生物宏基因组对未知样本进行地理定位的方法，其包括以下步骤：

步骤1、训练预测模型

将已知来源的微生物样本作为训练样本，将输入到预测模型中进行训练，而在训练之前需要对训练样本进行数据预处理和特征选择；

步骤1.1、数据预处理

对训练样本的宏基因测序数据进行预处理，包括质控、丰度定量、丰度分级；

对训练样本的各菌株的丰度进行分级，分级方法通过多个门限值将双精度的菌株丰度转换成离散的多元值；

步骤1.2、特征选择

从样本的所有菌株的集合中，挑选出有区分鉴定能力的菌株，作为特征菌株；

步骤1.3、预测模型训练

在这一步骤中，采用训练样本中经过特征选择后的所有特征菌株丰度值分级多元值，和这些训练样本的来源城市作为输入，利用机器学习的办法进行训练，从而得到预测模型；

步骤2、未知样本的地理定位

步骤2.1、对未知样本进行数据预处理

对未知样本的的宏基因测序数据进行预处理，包括质控、丰度定量、丰度分级；

步骤2.2、特征选择

从未知样本的所有菌株的集合中，挑选出一部分有区分鉴定能力的菌株，作为特征菌株；

步骤2.3、未知样本在训练集城市上的概率预测

将未知样本的所有特征菌株的丰度分级多元值做为输入，输入至预测模型中，得到该样本来自n个城市的概率y_i,(i＝1,2,...,n)；

步骤2.4、地理定位

如果该未知样本来自这n个城市，那么，将预测概率最高的城市作为该未知样本的来源城市。

所述步骤2.4还包括：

若未知样本不是来自训练集的n个城市，假设有n个训练集内城市在指定坐标系下的坐标为(x_i,y_i),(i＝1,2,...,n),然后，未知样本在这些城市的概率分别为z_i,(i＝1,2,...,n)，采用插值法对指定坐标系上的所有城市进行概率计算，概率最高的城市即为未知样本的来源城市。

所述指定坐标系为地理坐标系，城市在地理坐标系下的地理坐标为城市所在的经纬度坐标。

所述指定坐标系为生物坐标系，城市在生物坐标系下的生物坐标由其地理坐标经仿射变换得到，具体如下：

将训练样本的所有特征菌株的丰度分级多元值做为输入，通过流形学习的方法TSNE进行降维，从而得到训练集中每个样本的二维坐标；对于训练集内城市，通过训练集中来自该城市的所有样本的二维坐标计算中心点坐标，将其作为该城市的生物坐标；城市所在的经纬度坐标作为其地理坐标，通过仿射变换将训练集内城市的地理坐标转换为其对应的生物坐标；将训练集外城市的地理坐标也通过此仿射变换，转化为生物坐标；

在生物坐标系上，当通过插值法概率最大点的生物坐标时，将该生物坐标通过你仿射变换得到其地理坐标，该地理坐标对应的城市即为未知样本的来源城市。

所述步骤1.2和步骤2.2中，采用结合了递归特征消除和随机森林两种算法的集成学习方法来进行特征选择。

在数据预处理过程中，对各菌株的丰度进行分级具体为：将双精度丰度值转换成-1，0，1的三元值；对一个样本中所含的每一个菌株，其丰度值低于25％的转换为-1，丰度值在25％与75％之间的转换为0，丰度值大于75％的转换为1。

采用上述方案后，本发明基于微生物宏基因组对未知样本的来源城市进行预测定位，且在对训练样本和未知样本进行数据处理的过程中，对样本中各菌株的丰度进行分级，通过多个门限值将双精度的菌株丰度转换成离散的多元值。该分级方法是一种将连续值转换为离散值的量化方法，提取不同菌株的丰度值之间的显著差异，而忽略微小差异。通过该分级方法去噪，从而增加算法的稳定性和鲁棒性。对相较于现有的定位方法，本发明预测准确性高。

此外，本发明通过设定一指定坐标系，并将训练集内的城市和训练集外的城市均用该指定坐标系下的坐标来表示，然后采用插值法对指定坐标系上的所有城市进行概率计算，概率最高的城市即为未知样本的来源城市，而该城市可能不存在于训练样本的来源城市集合中。也就是说，本发明不仅可以预测出属于训练样本来源城市的未知样本，也可以预测出属于训练样本来源城市外的其他城市的未知样本，进一步提高了未知样本地理未知预测的准确度。

附图说明

图1为本发明数据处理流程图；

图2为本发明预测模型训练流程图；

图3为本发明未知样本定位流程图；

图4为本发明地理坐标与生物坐标的仿射变换示意图。

具体实施方式

如图1至图3所示，本发明揭示了一种通过微生物宏基因组对未知样本进行地理定位的方法，其具体包括以下步骤：

步骤1、训练预测模型

将已知来源的微生物菌株作为训练样本，将其输入到预测模型中进行训练。而在训练之前需要对训练样本进行数据预处理和特征选择。

步骤1.1、数据预处理

训练样本的宏基因测序数据在FASTAQ文件中以短读长序列的方式存储，每个短读长序列的内容表示为4行文本：1)短读长序列的头信息，2)短读长序列本身，或称碱基序列，3)预留给其他附加信息，4)对应于碱基序列的质量值序列。

对训练样本的宏基因测序数据进行预处理，包括质控、丰度定量、丰度分级。

首先，首先对宏基因测序数据中的短读长序列进行质控。

在测序建库的过程中，可能有各种物理化学原因或污染，测序技术和测序仪本身的缺陷等，都会造成测序结果里有碱基质量过低，或者含有其他来源的污染序列。为了后续生信分析的可靠性，通过质控把这些序列部分或整条清除，滤除掉不符合质量标准的测序数据。

随后，通过将质控处理后的宏基因测序数据和参考基因组序列集合，检测训练样本中存在的微生物，并对各个菌株的丰度进行定量。

最后，对训练样本的各个菌株的丰度进行分级。

对各个样本中的微生物菌株进行丰度定量后，得到该样本中每个菌株的比例，如果直接采用这些值进行后续的机器学习，可能会导致模型的过拟合，尤其是在样本数量不足的情况下。因此，需要对各菌株的丰度进行分级，通过多个门限值将双精度的菌株丰度转换成离散的多元值。在本实施例中将双精度丰度值转换成-1，0，1的三元值。对一个样本中所含的每一个菌株，其丰度值低于25％的转换为-1，丰度值在25％与75％之间的转换为0，丰度值大于75％的转换为1。至此，每个样本可以用一个取值在{-1,0,1}范围内的向量表示，向量的每一元表示该样本中特定菌株的特征，菌株的集合为样本的所有菌株的并集。

步骤1.2、特征选择

由于训练集的样本上包含了大量的菌株，因此，如果直接把所有的菌株的丰度作为预测模型的输入特征，其特征向量的维度会非常大，不利用后续分析。同时，许多菌株对样本的地理定位也没有帮助。因此，我们首先需要通过特征选择的办法，选定一些有区分鉴定能力的菌株作为特征。

有许多机器学习的算法可以进行特征选择，在本实例中，采用集成学习(ensemblelearning)的方法，结合了递归特征消除和随机森林[4]两种算法，采取这两种算法分别选取的特征的并集作为预测模型的输入。

在递归特征消除的算法中，对于初始的特征集合，通过逻辑回归模型的特征权重去评价每个特征，将特征权重最低的特征从特征集合中剔除，再将新的特征集合做为输入，重新评估，直到选出权重最高的指定数目的特征。在我们的应用实例中，通过递归特征消除选取了50个特征菌株。

在随机森林进行特征选择的算法中，随机森林中每一棵决策树中的每一个节点都是关于某个特征的条件，通过这些决策树可以将数据集按照不同的标签分类。对每棵决策树的每一个节点，可以计算其基尼不纯度(Gini impurity)。节点的基尼不纯度是指，根据节点中样本的分布对样本分类时，从节点中随机选择的样本被分错的概率。在训练决策树森林的过程中，可以计算出每个特征平均减少了多少不纯度，并把它平均减少的不纯度作为特征选择的权重。在我们的应用实例中，通过随机森林选出了243个有区分度的特征菌株。

将递归消除法和随机森林得到的两个特征集取并集，得到的集合做为我们的下一步预测模型的输入。

步骤1.3、预测模型训练

在这一步骤中，采用训练样本中经过特征选择后的所有特征菌株丰度值分级三元值，和这些训练样本的来源城市作为输入，利用机器学习的办法进行训练，从而得到预测模型。

在实例中，采用逻辑回归结合OVR(one vs rest)构造多分类器，预测每个样本来自训练集上不同城市的概率。假设我们的训练样本来自n个城市，在训练过程中，我们将其中某个城市的样本作为一类，将其余的所有城市的样本作为另一类，这样就形成了n个二分类问题。然后，使用逻辑回归算法对这n个二分类问题进行建模，并通过多数投票法，综合n个二分类器的输出，从而获得每个输入样本来自不同城市的的概率。

步骤2、未知样本在训练集城市上的预测

步骤2.1、对未知样本进行数据预处理

对未知样本的的宏基因测序数据进行预处理，包括质控、丰度定量、丰度分级，获得未知样本在选取的特征菌株上面的丰度分级三元值。具体处理过程与步骤1.1相同，此处不再赘述。

步骤2.2、特征选择

从未知样本的所有菌株的集合中，挑选出一部分有区分鉴定能力的菌株，作为特征菌株。具体处理过程与步骤1.2相同，此处亦不再赘述。

步骤2.3、未知样本在训练集城市上的概率预测

将未知样本的所有特征菌株的丰度分级三元值做为输入，输入至预测模型中，得到该样本来自n个城市的概率z_i，(i＝1，2，...，n)。未知样本在城市i上的预测结果为：

其中，z_i为模型预测的未知样本在城市i上的概率，x代表模型的输入特征，w和θ为模型通过训练所学习得到的参数。

步骤2.4、地理定位

如果通过其他技术手段已经可以确定该未知样本来自这n个城市，那么，我们就可以简单的将预测概率最高的城市当作是该未知样本的来源城市。

由于预测器仅能输出未知样本在训练集上面的n个城市的概率，而不能预测出没有训练样本的城市。因此，为了预测未知样本在训练集外城市的概率，采用数值插值的办法，利用未知样本在训练集内城市的概率，在地图上进行插值，从而获得未知样本在训练集外的城市上面的概率。

若不能确定未知样本是否是来自训练集的n个城市，则地图上所有的城市均有可能成为目标城市。对于训练集外的目标城市，未知样本在其上的概率计算如下。假设有n个训练集内城市在指定坐标系上的坐标为(x_i,y_i),(i＝1,2,...,n),然后，未知样本在这些城市的概率分别为z_i,(i＝1,2,...,n)。

假设目标城市在同一指定坐标系上的坐标为(x,y)。

采用Kriging插值的方法(参考资料：Jean-Paul Chiles and Nicolas Desassis,‘Fifty Years of Kriging’,Handbook of Mathematical Geosciences,pp 589-612.)，得到该未知样本来自目标城市的概率：

这里的

是未知样本来自目标城市(x,y)的概率估计值，λ_i是城市i的权重系数，z_i是该样本来自城市i的概率。在Kriging插值法里，权重系数是能够满足点(x,y)处的估计值

与真实值z的差最小的最优系数，即

约束条件为：

同时满足无偏估计的条件，即

具体做法简述如下：Kriging方法的优化目标为：

即：

为了简化公式我们定义

从而将优化目标简化为：

定义半方差函数r_ij＝σ²-C_ij，其中σ²表示方差。在Kriging假设下，z是空间均一的，即在空间任意一点(x,y)上，z都有同样的期望e和方差σ²。我们通过半方差函数可以将优化目标的最优解转换以下形式：

这里的φ是拉格朗日乘数。将上述方程组转换成矩阵形式即：

对于r_ij，Kriging插值假设空间上相近的两点的属性相近，即r_ij和(i,j)两点的距离d_ij存在着函数关系，我们通过已知两点i，j的距离d_ij和半方差r_ij，拟合一条最优拟合曲线，来刻画d和r的关系，得到函数：

r＝r(d)

通过该拟合函数，对于任意两点，我们可以通过计算其距离，得到它们的半方差r_ij，然后对上述矩阵求逆即可以得到Kriging插值系数的最优解。从计算得到的Kriging插值系数，我们可以计算出样本在训练集外城市的概率。

上述指定坐标系为地理坐标系。虽然空间距离相近的城市菌群特征会更相似，但是空间距离并无法完全反映城市间的菌群相似性。比如，东海岸线上的城市Ce和西海岸线上的城市Cw在地理距离上横跨整个大陆。但是Ce和Cw同样在海岸线上，可能具有相似的地理环境，这两个城市间的相似性会高于与其相邻的内陆城市。因此，我们将地理坐标转换成生物坐标，用生物距离代替地理距离，采用生物坐标系做为以上通过Kriging插值来预测概率值的办法的指定坐标系。

将训练样本的所有特征菌株的丰度分级三元值做为输入，通过流形学习(manifold learning)的方法TSNE进行降维，从而得到训练集中每个样本的二维坐标。对于训练集内城市，通过训练集中来自该城市的所有样本的二维坐标计算中心点坐标，将其作为该城市的生物坐标。城市所在的经纬度坐标作为其地理坐标，通过仿射变换(affinetransform)将训练集内城市的地理坐标转换为其对应的生物坐标。并将训练集外城市的地理坐标也通过此仿射变换，转化为生物坐标，如图4所示，从而可以得到所有城市的生物坐标。

最后，在生物坐标系上，基于该样本在训练样本的来源城市上的概率值，通过上述Kriging插值的方法对所有城市进行概率估计，从而得到该样本概率最大点的城市，做为该未知样本的预测城市。

在训练样本的来源城市上的概率值是通过步骤2预测得到的，而在待考核的未知城市上的概率值是通过Kriging插值得到的。从而得到该样本概率最大点的生物坐标,并得到该未知样本的城市。在本实例中，我们从在线数据库(https://simplemaps.com/data/ world-cities)下载了所有城市的经纬度地理位置信息，把地图上的所有城市都列为待考核的未知城市。在实际应用时，可以考虑一个可能的来源城市的列表，做为待考核的未知城市。

本发明基于微生物宏基因组对未知样本的来源城市进行预测定位，且在对训练样本和未知样本进行数据处理的过程中，对样本中各菌株的丰度进行分级，通过多个门限值将双精度的菌株丰度转换成离散的多元值。该分级方法是一种将连续值转换为离散值的量化方法，提取不同菌株的丰度值之间的显著差异，而忽略微小差异。通过该分级方法去噪，从而增加算法的稳定性和鲁棒性。对相较于现有的定位方法，本发明预测准确性高。

此外，本发明通过设定一指定坐标系，并将训练集内的城市和训练集外的城市均用该指定坐标系下的坐标来表示，然后采用插值法对指定坐标系上的所有城市进行概率计算，概率最高的城市即为未知样本的来源城市，而该城市可能不存在于训练样本的来源城市集合中。也就是说，本发明不仅可以预测出属于训练样本来源城市的未知样本，也可以预测出属于训练样本来源城市外的其他城市的未知样本，进一步提高了未知样本地理未知预测的准确度。

以上所述，仅是本发明实施例而已，并非对本发明的技术范围作任何限制，故凡是依据本发明的技术实质对以上实施例所作的任何细微修改、等同变化与修饰，均仍属于本发明技术方案的范围内。