本发明属于宏基因组分类分析，更具体地，涉及一种基于泛基因组图对宏基因组数据进行菌株水平分类的方法。
背景技术：
1、微生物广泛存在于地球上的各种环境中，如海洋、土壤和胃肠道，并在维持生态平衡和宿主健康方面发挥着不可或缺的作用。微生物群落通常由多种微生物组成，具有不均匀的物种丰度和高度的多样性。由于微生物的快速突变和水平基因转移，同一物种内可能出现不同菌株，这进一步增加了微生物群落的生物多样性和复杂性。先前的研究表明，不同的菌株在环境中可能表现出不同的表型或执行不同的生物学功能。例如，大肠杆菌菌株e.coli o104具有致病性，而普通大肠杆菌菌株则无害。此外，菌株在人体不同部位的分布也存在差异。因此，对微生物群落进行菌株水平分类对于理解微生物群落动态至关重要。
2、为了揭示微生物群落的结构和功能，需要在物种水平和菌株水平估计其组成和丰度。宏基因组测序数据保留了微生物群落大多数遗传信息，使其能够进行菌株水平分类；宏基因组测序数据包括二代(ngs)短读长和三代(tgs)长读长两种类型。
3、用于ngs的物种以上水平的分类方法可以分为三类：第一类是基于标记基因(marker)的方法，如metaphlan4；第二类是dna-to-protein的方法，如kaiju；第三类是dna-to-dna的方法，如kraken2、clark等；尽管前述方法是基于ngs开发的，但如kraken2、clark这些方法也可用于tgs；专门用于ngs的单物种菌株水平分类方法是例如基于dna-to-dna的strainscan；专门用于tgs的物种水平以上的分类方法同样可以分为三类：基于标记基因的melon方法、基于dna-to-protein的megan-lr方法、以及基于dna-to-dna的metamaps方法；其中metamaps也用于tgs菌株水平分类；值得注意的是，centrifuge是一种同时适用于ngs和tgs，并且能够进行物种水平和菌株水平分类的方法。
4、然而，上述几种现有方法普遍存在一些不可忽略的技术问题：
5、第一、上述几种方法通常都选择一个完整的基因组作为每个物种的代表，由于这种单一参考基因组的策略无法捕捉一个物种内部不同菌株的基因组多样性，因此会限制这些方法进行菌株水平分类的能力。
6、第二、这些方法如metaphlan4、kaiju、kraken2、clark、strainscan、melon、megan-lr、metamaps都不能同时高效处理ngs和tgs两种类型的宏基因组数据；基于ngs数据开发的方法在处理tgs数据时，由于算法和参数的不完全适配，往往会导致性能下降；此外，基于k-mer的方法无法充分利用tgs数据的优势。
7、第三、这些方法如metaphlan4、kaiju、kraken2、clark、melon、megan-lr都只能在物种水平上进行分类，而难以实现菌株水平分类甚至多物种菌株水平分类。另外stranscan方法只能在单物种菌株水平上进行分类，同样存在局限性。
技术实现思路
1、针对现有技术的以上缺陷或改进需求，本发明提供了一种基于泛基因组图对宏基因组数据进行菌株水平分类的方法(pangenome graph based taxonomic classifier，简称pantax)。其目的在于，解决现有方法由于单一参考基因组的策略，导致菌株水平分类能力受到限制的技术问题；以及由于ngs和tgs两种类型数据特性不同，导致难以同时高效处理两种类型数据的技术问题；以及由于开发力度不够和菌株高相似性，难以进行菌株水平分类甚至是多物种菌株水平分类的技术问题。
2、为实现上述目的，按照本发明的一个方面，提供了一种基于泛基因组图对宏基因组数据进行菌株水平的方法，包括如下步骤：
3、(1)从公共数据库获取多个物种及其对应的多个完整基因组，并根据获取的所有完整基因组构建参考数据库。
4、(2)针对步骤(1)得到的参考数据库中的每个物种而言，使用该物种在该参考数据库中对应的所有参考基因组构建该物种对应的泛基因组图，并将所有物种对应的泛基因组图进行合并，以得到参考泛基因组图。
5、(3)获取用户输入的宏基因组数据，将步骤(2)得到的参考泛基因组图与该宏基因组数据进行比对，以得到该宏基因组数据中每条测序读段的多条比对结果。
6、(4)根据步骤(3)得到的宏基因组数据中每条测序读段的所有比对结果进行物种水平分类分箱处理，以得到该宏基因组数据对应的最终物种组成。
7、(5)针对步骤(4)得到的宏基因组数据对应的最终物种组成中的每个物种而言，对该物种进行物种水平分类谱分析处理，以获取该物种的相对丰度。
8、(6)针对步骤(4)得到的宏基因组数据对应的最终物种组成中的每个物种而言，根据步骤(2)得到的该物种对应的泛基因组图对该物种进行路径丰度优化pao处理，以得到宏基因组数据的最终菌株组成、以及该最终菌株组成中每个菌株的绝对丰度。
9、(7)根据步骤(4)得到的宏基因组数据对应的最终物种组成对步骤(6)得到的宏基因组数据的最终菌株组成、以及该最终菌株组成中每个菌株的绝对丰度进行丰度调整和归一化计算，以得到宏基因组数据的最终菌株组成、以及该最终菌株组成中每个菌株的相对丰度。
10、优选地，步骤(1)包括以下子步骤：
11、(1-1)从公共数据库下载获取多个物种、以及与每个物种对应的多个完整基因组；
12、(1-2)对步骤(1-1)得到的所有完整基因组进行去除质粒处理，以得到多个去除质粒后的完整基因组；
13、(1-3)针对步骤(1-1)获取的每个物种而言，使用fastani工具计算该物种对应的任意一对完整基因组之间的平均核苷酸相似性ani，作为该物种对应的该对完整基因组之间的距离得分；
14、(1-4)针对步骤(1-1)获取的每个物种而言，使用基于连通图的启发式聚类算法对步骤(1-3)得到的该物种对应的所有对完整基因组之间的距离得分进行聚类处理，以得到多个完整基因组集合，从所有完整基因组集合中选择包含完整基因组最多的完整基因组集合，并从选择的该完整基因组集合中选取预定数量的完整基因组作为该物种对应的多个参考基因组，所有物种对应的所有参考基因组构成参考数据库。
15、优选地，步骤(1-4)中，如果聚类结果小于预设的第一聚类阈值，则说明该聚类结果对应的一对完整基因组是不同物种的菌株；如果聚类结果大于等于该预设的第一聚类阈值且小于预设的第二聚类阈值，则说明该聚类结果对应的一对完整基因组是同一物种的不同菌株；如果聚类结果大于预设的第二聚类阈值，则说明该聚类结果对应的一对完整基因组是同一个菌株；
16、步骤(1-4)的聚类处理是获取大于等于该预设的第一聚类阈值且小于预设的第二聚类阈值的聚类结果作为聚类处理得到的多个完整基因组集合。
17、优选地，步骤(2)包括以下子步骤：
18、(2-1)针对步骤(1)得到的参考数据库中的每个物种而言，如果该物种在参考数据库中存在大于一个参考基因组，则使用pggb工具为该物种对应的所有参考基因组构建该物种对应的泛基因组图；如果该物种在参考数据库中仅存在一个参考基因组，则将该参考基因组分割成长度为1024bp的片段，然后将每个片段作为一个节点构建泛基因组图，从而得到该物种对应的泛基因组图。
19、(2-2)将步骤(2-1)得到的所有物种对应的泛基因组图进行合并，以得到参考泛基因组图。
20、步骤(3)中，对于二代短读长的宏基因组数据而言，是使用giraffe工具进行二者的比对；具体而言，首先使用vg工具为参考泛基因组图构建索引，以加速比对过程，然后，使用giraffe工具将短读长的宏基因组数据与参考泛基因组图进行比对，以得到短读长的宏基因组数据中每条测序读段的多条比对结果；
21、步骤(3)中，对于三代长读长的宏基因组数据而言，是使用graphaligner工具进行二者的比对，以得到长读长的宏基因组数据中每条测序读段的多条比对结果。
22、步骤(4)包括以下子步骤：
23、(4-1)针对步骤(3)得到的宏基因组数据中的每条测序读段而言，根据该测序读段的所有比对结果得到该测序读段的最优比对结果。
24、(4-2)针对步骤(3)得到的宏基因组数据中的每条测序读段而言，根据步骤(4-1)得到的该条测序读段的最优比对结果获取该条测序读段对应的物种，该宏基因组数据中所有测序读段对应的物种构成该宏基因组数据对应的第一物种组成(其具有潜在假阳性)。
25、(4-3)对步骤(4-2)得到的宏基因组数据对应的第一物种组成进行假阳性物种过滤处理，以得到第二物种组成作为该宏基因组数据对应的最终物种组成。
26、优选地，步骤(4-1)中，如果一条测序读段具有多个比对结果，只保留最优比对结果，以确保每条测序读段对应一个物种；具体而言，在短读长的比对结果中，每条测序读段只有一个比对结果，因此该比对结果都是最优比对结果；而在长读长的比对结果中，通过比对质量mapq值、匹配的碱基数、同一性和测序读段长度从所有比对结果中确定最优比对结果；
27、步骤(4-3)具体为，针对第一物种组成中的每一条物种而言，从该物种对应的多个测序读段的最优比对结果获取多个测序读段对应的mapq值，在该物种对应的所有测序读段中，如果其中至少有一条测序读段的mapq值达到最高分60；并且至少十分之一的测序读段的mapq值大于2，则该物种被鉴定为真，否则该物种被鉴定为假阳性而后过滤。
28、优选地，步骤(5)包括以下子步骤：
29、(5-1)针对步骤(4)得到的宏基因组数据对应的最终物种组成中的每个物种而言，根据步骤(4-1)得到的该物种对应的每个测序读段的最优比对结果获取该物种的绝对丰度。
30、本步骤是采用以下公式计算最终物种组成中第i个物种的绝对丰度ci：
31、
32、其中，i∈[1，最终物种组成中的物种总数t]，ni是第i个物种对应的测序读段总数，j∈[1，ni]，nij是第i个物种对应的第j个测序读段，|nij|是第i个物种对应的第j个测序读段的碱基数，mi是第i个物种的对应的参考基因组总数，k∈[1，mi]，gik是第i个物种对应的第k个参考基因组，|gik|是第i个物种对应的第k个参考基因组的碱基数。
33、(5-2)对步骤(5-1)得到的最终物种组成中每个物种的绝对丰度进行归一化计算，以得到该物种的相对丰度；
34、本步骤具体是采用以下公式计算最终物种组成中第i个物种的相对丰度ai：
35、
36、优选地，步骤(6)包括以下子步骤：
37、(6-1)设置计数器i＝1；
38、(6-2)判断i是否大于宏基因组数据对应的最终物种组成的物种总数t，如果是则过程结束，否则转入步骤(6-3)。
39、(6-3)针对最终物种组成中的第i个物种而言，根据步骤(4-1)得到的该物种对应的每条测序读段的最优比对结果计算该物种对应的泛基因组图中每个节点的观测丰度。
40、(6-4)针对最终物种组成中的第i个物种而言，根据步骤(2-1)得到的该物种对应的泛基因组图和步骤(4-1)得到的该物种对应的每条测序读段的最优比对结果过滤该物种中的假阳性菌株，以得到该物种对应的、第一次假阳性过滤后的菌株组成。
41、(6-5)针对最终物种组成中的第i个物种而言，根据步骤(6-3)得到的该物种对应的泛基因组图中每个节点的观测丰度，对步骤(6-4)得到该物种对应的的第一次假阳性过滤后的菌株组成进行第一次pao求解，并根据求解结果保留菌株绝对丰度大于零的菌株，以得到该物种对应的第一菌株组成、以及该第一菌株组成中每个菌株的绝对丰度。
42、(6-6)针对最终物种组成中的第i个物种而言，根据步骤(4-1)得到的该物种对应的每条测序读段的最优比对结果和步骤(6-5)得到该物种对应的第一菌株组成和第一菌株组成中每个菌株的绝对丰度过滤该物种中的假阳性菌株，以得到该物种对应的、第二次假阳性过滤后的菌株组成；
43、(6-7)针对最终物种组成中的第i个物种而言，根据步骤(6-3)得到的该物种对应的泛基因组图中每个节点的观测丰度，对步骤(6-6)得到的该物种对应的、第二次假阳性过滤后的菌株组成进行第二次pao求解(该过程与上述第一次pao求解过程完全相同，在此不再赘述)，并根据求解结果保留菌株绝对丰度大于零的菌株，以得到第二菌株组成、以及该第二菌株组成中每个菌株的绝对丰度，二者分别作为该物种的最终菌株组成、以及该最终菌株组成中每个菌株的绝对丰度。
44、(6-8)设置i＝i+1，并返回步骤(6-2)。
45、优选地，针对该物种对应的泛基因图中的每个节点而言，步骤(6-3)是采用以下公式计算该节点v的观测丰度
46、
47、其中，v∈v，v是该物种对应的泛基因组图中的所有节点，|v|是该节点v的碱基数，rv是该物种对应的覆盖该节点的所有测序读段，r∈[1，rv]，rvr是该物种对应的覆盖该节点的第r条测序读段，|rvr|是该物种对应的覆盖该节点的第r条测序读段的碱基数；
48、步骤(6-4)具体为，首先，将步骤(2-1)得到的第i个物种对应的泛基因组图中每个菌株基因组路径上的每三个连续(顺序出现)节点构建为一个三节点组vtriplet；
49、然后，获取步骤(4-1)得到的该物种对应的所有测序读段在该物种对应的泛基因组图中所覆盖的所有三节点组所构成的三节点组集合vreads；
50、随后，针对该物种对应的泛基因组图中每个菌株基因组路径而言，获取该菌株基因组路径上的所有三节点组，针对该菌株基因组路径对应每个三节点组而言，如果这个三节点组仅出现在该菌株基因组路径而没有出现在其他菌株基因组路径中，则这个三节点组为该菌株基因组路径对应的菌株特异性三节点组，该菌株基因组路径对应的所有菌株特异性三节点组构成菌株特异性三节点组集合vstrain；
51、其后，针对该物种对应的泛基因组图中的每个菌株基因组路径而言，采用以下公式计算该菌株基因组路径对应的菌株特异性三节点组集合在该物种对应的所有测序读段中出现的比例fstrain：
52、
53、其中，fstrain∈[0，1]，1表示该菌株基因组路径对应的菌株特异性三节点组集合都被该物种对应的所有测序读段所覆盖；
54、最后，如果得到的比例fstrain＜0.3，则该菌株基因组路径对应的菌株被鉴定为假阳性菌株并被过滤掉，否则该菌株被鉴定为真，并且作为第一次假阳性过滤后的菌株组成中的一个菌株。
55、优选地，步骤(6-5)中的pao目标函数为：
56、
57、其中，v是该物种对应的泛基因组图中的所有节点，p是构建该物种对应的泛基因组图的所有菌株基因组路径，是泛基因组图中节点v∈v的观测丰度，ap是菌株基因组路径p∈p所对应的菌株绝对丰度，即pao所需要求解的决策变量。
58、步骤(6-6)具体为，首先，针对步骤(6-5)得到该物种对应的第一菌株组成中的每个菌株而言，获取该菌株在该物种对应的泛基因组图中对应的菌株基因组路径，作为第一菌株组成中该菌株对应的菌株基因组路径；
59、然后，针对物种对应的第一菌株组成中每个菌株对应的菌株基因组路径而言，根据步骤(4-1)得到的该物种对应的每条测序读段的最优比对结果计算该菌株基因组路径对应的菌株特异性三节点组集合中每个菌株特异性三节点组的平均覆盖深度；本过程具体是采用以下公式计算该菌株特异性三节点组vtriplet的平均覆盖深度c(vtriplet)：
60、
61、其中，s∈s，s是该物种对应的覆盖该菌株特异性三节点组的所有测序读段，|s是该物种对应的覆盖该菌株特异性三节点组的测序读段s的碱基数，|vtriplet|是该菌株特异性三节点组vtriplet的碱基数。
62、随后，针对第一菌株组成中每个菌株对应的菌株基因组路径而言，根据该菌株基因组路径对应的菌株特异性三节点组集合中每个菌株特异性三节点组的平均覆盖深度计算该菌株基因组路径的预估菌株绝对丰度；本过程具体是采用以下公式计算该菌株基因组路径p的预估菌株绝对丰度ap,triplet：
63、
64、其中，是|vstrain|该菌株基因组路径p对应的菌株特异性三节点组集合vstrain中菌株特异性三节点组的总数。
65、其后，针对第一菌株组成中的每个菌株而言，采用以下公式计算该菌株的绝对丰度ap与该菌株对应的菌株基因组路径的预估菌株绝对丰度ap,triplet之间的差异：
66、
67、最后，如果得到的差异dstrain＞0.45，则该菌株基因组路径对应的菌株被鉴定为假阳性菌株并被过滤，否则该菌株被鉴定为真，作为第二次假阳性过滤后的菌株组成中的一个菌株。
68、优选地，步骤(7)包括以下子步骤：
69、(7-1)根据步骤(4)得到的宏基因组数据对应的最终物种组成和步骤(5-1)得到的该最终物种组成中每个物种的绝对丰度，对步骤(6)得到的宏基因组数据的最终菌株组成中每个菌株的绝对丰度进行丰度调整，以得到该菌株调整后的菌株绝对丰度；其中本步骤中，针对步骤(4)得到的宏基因组数据对应的最终物种组成中的每个物种而言，如果该物种中绝对丰度最大的菌株超过该物种对应的绝对丰度，则对该物种内所有菌株的绝对丰度进行缩放处理，以使该物种中所有菌株的绝对丰度总和等于该物种的绝对丰度；
70、(7-2)针对宏基因组数据的最终菌株组成而言，对步骤(7-1)得到的该最终菌株组成中所有菌株调整后的菌株绝对丰度进行归一化计算，以得到该最终菌株组成中每个菌株的相对丰度。
71、总体而言，通过本发明所构思的以上技术方案与现有技术相比，能够取得下列有益效果：
72、(1)本发明由于采用了步骤(2)，其通过泛基因组图的技术捕获了每个物种的多个菌株的基因组作为参考，而不仅仅是使用每个物种的代表基因组作为参考；本发明通过菌株多样性的扩展，增强了宏基因组测序读段的比对准确性和在菌株水平上组成的估计，因此能够解决现有方法由于采用单一参考基因组的策略，导致菌株水平分类能力受到限制的技术问题；
73、(2)本发明由于采用了步骤(3)，其在ngs和tgs宏基因组数据与泛基因组图的比对中，分别使用相应的短读长比对工具giraffe和长读长比对工具graphaligner，因此能够解决由于ngs和tgs两种类型数据特性不同，导致现有方法如metaphlan4、kaiju、kraken2、clark、strainscan、melon、megan-lr、metamaps都不能同时高效处理ngs和tgs两种类型的宏基因组数据的技术问题；
74、(3)本发明由于采用了步骤(4)到步骤(7)，其通过宏基因组数据与参考泛基因组比对可以获得所有物种水平分类结果，并且进一步基于泛基因组图进行路径丰度优化处理可以得到所有物种对应的菌株水平分类结果，因此能够解决现有方法如metaphlan4、kaiju、kraken2、clark、melon、megan-lr只能在物种水平上进行分类，而难以实现菌株水平分类甚至多物种菌株水平分类的技术问题。
1.一种基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，包括如下步骤：
2.根据权利要求1所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，步骤(1)包括以下子步骤：
3.根据权利要求1或2所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，
4.根据权利要求1至3中任意一项所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，步骤(2)包括以下子步骤：
5.根据权利要求4所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，
6.根据权利要求5所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，步骤(5)包括以下子步骤：
7.根据权利要求6所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，步骤(6)包括以下子步骤：
8.根据权利要求7所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，
9.根据权利要求8所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，
10.根据权利要求9所述的基于泛基因组图对宏基因组数据进行菌株水平的方法，其特征在于，步骤(7)包括以下子步骤：
