人体体内和体表均定植着大量微生物，其数量与人体细胞数量相近，而其功能基因的数量更是人体功能基因的百倍有余，因此人体微生物组也被视为人类基因组信息的重要补充，与人体健康息息相关[1-2]。此外，众多研究表明，人体微生物群落(尤其是肠道微生物群落)的紊乱，会影响宿主的免疫、代谢、神经、内分泌等多系统的功能，也是自身免疫性疾病、肥胖、糖尿病、自闭症、抑郁症、肠炎和癌症等多种慢性疾病的重要因素之一[3-5]。由此可见，微生物的构成(Composition)多样性和功能(Function)多样性，对宿主(人体)的表型有着重要的影响。如何能深入研究微生物群落中的构成单元(菌株，Strain)和功能单元(基因，Gene)与表型(Phenotype)之间的相互关系，是当前微生物组领域的一个核心问题。基于宏基因组数据(Metagenomics data)的有参考基因组(以下简称“有参”)的菌株分析算法有望成为回答该问题的最佳高通量手段。

菌株是微生物分类和研究中的最基础生命实体[6]。同一菌种(Species)的不同菌株之间在代谢能力、定植能力、致病能力等表型上可能存在很大的差异。例如，大肠杆菌中只有特定的菌株才会产生致癌物[7]，而其他的则是非致病性的[8]。反过来，同一菌种中也只有特定的菌株被证明对抑制肿瘤有效[9]。这些性状/功能的多样性，往往是基因多样性(Gene diversity)和单核苷酸多态性(Single nucleotide polymorphism，SNP)所导致的[10-11] (图 1)。其中，SNP可能会影响同源基因的转录和翻译等过程，从而影响蛋白的功能并导致性状表型的变化。而基因多样性，则反映了同源基因的有无，直接导致了不同菌株间特定功能的有无。而为了更好地展示这一概念，2005年Tettelin等提出了微生物泛基因组(Pan-genomoe)的概念[12]，用以指代某一物种全部同源基因的组成及序列信息，从而清晰反映某个同源基因在不同菌株间的出现情况。根据在菌株中出现的频次可分为：1)核心基因组(Core genome)，指所有菌株中都存在的基因；2)非必需基因组(Dispensable genome)，指只在部分菌株中存在的基因。

而一个微生物群落的某项表型，往往是群落里各菌株的相关功能的集合的体现。因此，群落里的菌株构成与其功能多样性，对整个群落的表型甚至宿主的表型都有着重大的影响(图 1)。

菌株代谢能力多样性会导致不同微生物群落对相同代谢底物产生不同的代谢产物。而随着肠道微生物与药物代谢研究的深入，菌株的代谢性状对药效的影响日益明确[13]。这使得微生物组个性化有望成为继基因组个性化后的第二类个性化医疗指标。因此，从样品中高通量地判别菌株，在未来的个性化治疗领域中将有巨大的价值和需求[14]。

菌株定植能力的多样性会导致群落移植到新的生境时，可能出现只有部分菌株能定植的情况。此外，研究发现，菌群在受到扰动的情况下(如抗生素处理)可能出现不同菌株的演替[15]。越来越多的研究关注微生物制剂在临床使用中的定植情况[16]以及使用后原有微生物群落的变化[17]。而对于复杂群落的应用(如粪菌移植)，更是需要对复杂菌株构成进行监控，从而确保其定植、功效、风险等[18]。

菌株致病能力的多样性会直接影响宿主的疾病表型。把某些致病菌株移植到小鼠模型上，可以观测到疾病的表型；而把同菌种的非致病菌株进行移植，则小鼠不会发病[19]。而在传染病中找到真正的致病菌株更是控制疫情的重要基础[20]。此外，研究表明，即使疾病组和对照组的菌种比例没有差异，但菌株仍可能存在差异[21]。并且菌株的选择性在疾病干预过程中有着重要作用[22]，因此菌株是一类重要的疾病判别特征[23]。而随着微生物与疾病的关系研究逐步从研究单菌株传染病向研究多菌株复杂群落的失调与非传染性疾病发生关系转移[24]，微生物构成和功能的多样性对宿主表型的影响显得更为常见和复杂，使得获取整个复杂群落的构成的需求变得日益重要。

综上所述，高通量、高分辨率(菌株水平)地获取微生物群落的构成和功能信息，在科学研究、工业生产和临床应用中都有着重要的价值和前景。

宏基因组测序是通过对样品中全部微生物的DNA进行提取后建立文库，并对文库内的随机片段进行序列测序和比对分析的一种高通量研究方法。它能很大程度上弥补培养组学的不足，也为更高通量地研究菌株的组成和功能信息提供了技术可能。宏基因组数据分析，主要可分为无参考基因组(无参)的组装分析和有参考基因组(有参)的比对分析两大类。

无参组装分析通过对测序片段进行组装(Assembly)并通过比较基因组学方法完成分析、解读。高质量的组装有助于获取完整微生物单倍型(Haplotype)，不但可以获取SNP信息用于菌株的追踪[25]，也可用于基因水平转移、进化选择的推演[25-26]，更能挖掘出未被发现的新菌株。然而由于基因组高度相似，直接从宏基因组测序数据完成同物种的不同菌株的组装十分困难[27]。尽管研究人员近年来开发出了一些运用多样本的复杂分箱方法[28]或辅以长读长测序数据的组装算法[29]，但这些方法往往需要大量计算资源，而且对目标物种的丰度及测序深度有一定要求。还会引入嵌合体和不完整基因组，并不适用于大量样品或者低丰度样品的高通量研究[8-30]。

有参比对分析受限于参考菌株基因组的数量[8]。随着微生物培养组学和全基因组测序的不断发展，利用参考基因组进行有参比对分析已经成为现实[8, 30-31]。基于宏基因组数据使用有参比对的方法，对特定的微生物菌种进行菌株构成分析和功能分析是完全可行的。本文将重点介绍有参比对的分析方法。

对应菌株的两类多样性(图 1)，宏基因组数据的菌株分析方法按使用及反馈的差异信息可以分为两大类(表 1)：基于SNP进行识别和基于基因组成(Genetic constitution)差异进行识别。SNP反映的是菌株在不同位点上的单碱基多样性；基因组成差异反映的是菌株的基因多样性。

在基于SNP识别的方法中，根据使用的参考基因组范围不同，又可分为两大类：使用全部基因组信息和仅使用标记基因(Marker gene)上的信息。使用全部的基因序列信息时，信息最全面，理论上菌株的辨识程度比用标记基因的方法高；但是因为比对的区域大，对计算资源要求很高，导致一次运行只能针对有限的菌株对象，常见情况是一次只能对一个菌种下的菌株进行区分，因此目前多用于“精细”追踪目标菌株在不同样品中的出现情况。使用标记基因时，因为只使用区分能力最好的部分基因作为标志基因，所以不容易受重组(Recombination)或前噬菌体(Prophage)等非保守区域的干扰，并且计算时可考虑更多菌株，甚至可以一次分析大部分需要考虑的多个菌种的菌株；但对于在标志基因区域高度相似的菌株，这类方法无法分辨。此外，因为不同方法使用的标志基因有所不同，所以其数据库与样品中菌株的匹配程度会对分析结果产生一定影响。

此外，虽然表 1中的方法都可以一定程度上发现样品中存在多个菌株的情况，但并不是所有方法都能够反馈多个菌株的比例，其中PanPhlAn和StrainPhlAn都假设每个菌种中只有一个优势菌株并只反馈该菌株的信息，而StrainSifter只会反馈样品间的进化距离，并无具体菌株比例。且PanPhlAn反馈的不是组成信息，而是各样品中目标菌种的泛基因组分布情况。

基于SNP的分析方法是菌株分析算法的主流。利用的是宏基因组数据中单碱基水平的信息，通过使用宏基因组片段序列和参考序列(全基因组或标记基因)的比对结果进行SNP位点的等位基因相对比例分析，最后根据这些比例信息，推算菌株的组成。因此，这类方法只能反馈菌株构成丰度，无法直接提供功能基因相关信息。

StrainSifter[20]和StrainEst[32]属于以全基因组为对象的方法。其中，StrainSifter完全是用菌株的全基因组信息，一次只能追踪一个菌株及其旁系变化，并以此信息构建样品间的进化关系图，但不反馈具体的菌株丰度。它主要用于快速并准确地检测样品中是否存在特定菌株及其旁系。而StrainEst则是以菌种数据库为单位，对该菌种数据库内的所有已知菌株的全基因组序列进行分析，找出差异代表区域(相当于标记基因)及SNP位点，然后再把宏基因组序列片段和这些区域进行比对，一定程度上结合了标记基因方法的思想，其与标记基因方法的主要差别是并非固定使用一套标记基因。此外，StrainEst会反馈具体的菌株丰度。因此，StrainEst适用于对特定菌种的已知菌株进行相对丰度的分析，研究菌株相对丰度和表型之间的关系。

ConStrains[33]、Strain Finder[31]、StrainPhlAn[8]都属于以标记基因为对象的方法。ConStrains和StrainPhlAn使用的都是目标菌种中的特异性标记基因库(MetaPhlAn[34])，主要差别在推算具体菌株丰度的算法上：ConStrains使用蒙特卡罗马尔可夫模型来推算菌株的丰度并结合修正后的赤池信息量准则(Corrected Akaike information criterion)来选择最优模型与结果；StrainPhlAn则侧重利用SNP信息计算样品间的进化关系，并不反馈各样品内的具体菌株丰度。此外StrainPhlAn有一个固定假设：一个样品中一个菌种只有一个代表菌株，因此该代表菌株的丰度等于该菌种丰度。Strain Finder使用通用标记基因为参考[35]，其默认数据库只针对人肠道微生物，但用户可以重新自定义构建。它使用的是最大期望(Expectation maximization，EM)算法，通过拟合不同样品中的等位基因频率变化来获取不同菌株的代表SNP序列及各样品中对应的丰度组成。Strain Finder被设计用于同时追踪粪菌移植前后的复杂群落样品中所有人肠道微生物各菌种内的菌株的丰度变化情况，也可用于对相近样品进行人肠道菌种的菌株普查。因为其方法允许识别数据库内不存在的菌株，所以可以一定程度上识别新的菌株。但由于该方法提取的SNP位置依赖于输入的样本，在不同的样本上获取的菌株信息并不能直接进行比较，因此它是整体样品分析的方法，需要相互比较的样品必须统一进行分析。此外，因为EM算法的特性，其结果都是局部最优，因此需要多次拟合来尽可能得到全局最优解[31]。

基于基因组成差异的菌株分析算法，利用的是宏基因组数据中基因水平的信息(如各基因的丰度与各参考菌株中对应基因的存在或缺失情况的拟合)。PanPhlAn[30]是目前笔者所知唯一的该类方法。PanPhlAn对目标菌种的各菌株全基因组序列中注释的基因进行聚类以合并为基因家族(Gene family)，选取每个基因家族的代表序列构建该菌种的泛基因组图谱作为参考数据库。通过宏基因组样本数据中序列与参考数据库的比对，根据计算出的每个基因家族的覆盖度判断其在该样本中的存在或缺失情况。该方法可直接提供每个基因家族在样品各菌株中的存在与缺失情况，从而反映该菌株的功能特征信息，因此可用于菌株功能与表型关系的研究(如菌株的毒力、抗药性等)。但该方法只提供优势菌株的基因家族组成信息，并无法提供菌株的构成丰度。而且由于复杂群落样品中物种间相似基因的干扰，样品中目标菌株的某些本应缺失的基因家族可能被误判为存在。

除了上述两大类菌株分析方法外，还有一些基于特定假设的菌株分析算法。例如：metaMLST[36]是基于细菌鉴定中的多位点分型(Multi-locus sequence typing，MLST)原理，通过从宏基因组数据中获取等位子组成信息，然后与PubMLST数据库[37]进行比对，从而实现对菌株的鉴定。Sigma[38]是一款对病原菌菌株进行鉴定和定量的工具。在假设该菌株已确定存在的前提下，通过将宏基因组数据匹配到已知的参考基因组上，实现对最优基因组的判定、定量和变异检测。此外，MetaPhlAn2[34]和mOTU2[39]都是较为常用的有参菌种分析方法，并提供菌株分析相关信息。其中，StrainPhlAn是在MetaPhlAn的基础上进行菌种分析；mOTU2则支持在菌种分析的基础上进行SNP分析，从而推断菌株组成。针对三代测序长读长数据，Dilthey等[40]开发了MetaMaps进行菌种菌株的鉴定，并可以进行基因的存在/缺失判别。但是该类数据和宏基因组测序数据差异较大，不在此进行讨论。

随着宏基因组测序应用的快速增长，宏基因组数据分析方法的评估备受重视，CAMI (Critical Assessment of Metagenome Interpretation)组织先后发布了其模拟数据体系[41]和评估平台[42]，并对7个分析方法进行了评估。随后，Ye等对20个宏基因组有参比对分析方法进行了详细评估[43]。但这些评估都只精细到菌种层次而非菌株层次，因此目前仍缺乏一个针对菌株分析方法的评估体系。而现有的评估体系在用于评估菌株时，会存在下面3个潜在问题。

首先，已有评估的菌株分析算法都只使用了模拟标准数据作为评估数据，且只考虑了以下两个因素：多菌株混合结合浓度梯度混合[31-33]。其中，ConStrains展示了其结果与模拟标准相比，在JS离散量(Jensen-shannon divergence, JSD)上优于随机结果。Strian Finder和StrainEst则分别显示其结果优于Constrain，且StrainEst还优于PanPhlAn和Sigma。但这些评估在选择菌株时，均未考虑菌株之间的序列相似度，无法深入了解和评判序列相似度对各菌株分析算法的影响，因此并不能真切反映各菌株分析算法的特点。

其次，通过实验制备的真实标准样品(如HMP mock[44]、ILPP mock[45]和MBARC[46])，都只是不同菌种的菌株混合，并无同菌种多菌株混合，只能反映菌种间的差异，无法用于评估菌株分析算法在同种多菌株存在场景下的精准度。

最后，目前的6个主流方法，数据反馈的侧重点差异较大，只有3个方法能反馈样品内的各菌株丰度信息，其余3个方法均只反馈各样品间的菌株基因组的差异信息，因此能通用的评估指标极其有限。而在反馈菌株丰度的方法中，只有StrainEst反映的菌株可直接与已知全基因组的参考菌株对应，而Strain Finder和ConStrains均只提供样品中代表性的SNP序列，导致分析结果和真实标准难以统一比较，大大提高了评估的难度。

由于菌株在微观上具有宿主特异性，在宏观上反映样本来源的地域、饮食等环境因素的影响，故菌株信息可用于追踪个体间传播途径，揭示环境因素的关联关系[25]。此外，由于不同菌株的基因及功能存在差异，菌株与宿主之间存在着不同的作用机制。菌株的这一特点可用于菌株作用机制研究以及关键菌株的发现筛选[47]。宏基因组菌株分析方法的出现，使得对菌株的研究脱离了实验培养的限制，为菌株研究和应用提供了新的分析思路和方法。

近两年来众多大型微生物组研究项目[48]相继启动。这些项目产生的大量宏基因组数据, 有助于科学家研究地域、饮食、年龄、疾病等因素与人体微生物群落的关系。通过对公共宏基因组数据集进行基于SNP的菌株分析，研究人员能够构建出不同菌株的进化树，并将进化树的结构与环境因素对应起来。通过对这些菌株进行泛基因组分析，研究人员能够进一步了解外界因素或者表型对应的菌株的生理功能特性。Truong等使用StrainPhlAn对多个国家的人体肠道宏基因组样本进行分析，发现微生物群体结构与群体的地理结构相关联[8]。De Filippis等使用PanPhlAn对意大利地区人的粪便样本进行宏基因组分析，发现素食或者杂食个体的普氏菌Prevotella copri中的菌株在基因层面存在明显差别[49]。Tett等发现P. copri存在4个主要的进化支，并且它们在西方人群中出现较少[50]。Vatanen等对DIABIMMUNE研究中来自芬兰、爱沙尼亚和俄罗斯的早期儿童数据集进行SNP以及宏基因组组装的菌株分析, 发现国家之间在拟杆菌属Bacteroides与双歧杆菌属Bifidobacterium中的菌种菌株上表现出明显的地域性特点，菌株的差别会对早期婴儿肠道微生物结构带来影响[51]。Zolfo等使用MetaMLST、StrainPhlAn以及PanPhlAn三种工具对1 614例城市环境样本的宏基因组数据进行了菌株分析，并从进化和功能的角度对结果进行了讨论[52]。

病原菌的监控与溯源是环境公共卫生领域面临的一个重要问题。许多研究开始使用宏基因组数据进行高效快捷的病原菌菌株筛查，以更好地支持政府和医疗机构的决策。Fresia等使用MetaMLST对城市的海滨及污水系统的样本进行宏基因组菌株分析，发现了与临床相关的一些病原菌的存在[53]。Hamner等使用Nanopore MinION测序仪以及CosmosID软件对河流样本进行宏基因组分析，检测出肠出血性和肠致病性大肠杆菌Escherichia coli菌株，以及其他病原性E. coli菌株[54]。一些研究尝试将菌株筛查用于食品安全领域。Walsh等使用MetaMLST、StrainPhlAn以及PanPhlAn对刚果地区的发酵乳制品进行检测，在一些样本中发现了可能的病原性E. coli以及肺炎克雷伯氏菌Klebsiella pneumoniae菌株[55]。

此外，由于菌株在宿主个体内有高度特异性，一些研究者尝试通过宏基因组数据研究菌株的传播途径，尤其是母婴间的传递[56-59]。Tamburini等使用菌株溯源工具StrainSifter，并用其比较血液中分离的病原菌菌株序列与肠道宏基因组数据，从而判断这些菌株是否来源于病人肠道[20]。Olm等使用菌株分析工具InStrain研究剖腹产与顺产的新生儿菌株，发现剖腹产的婴儿菌株可能来源于医院[60]。这些溯源工作对指导医院针对性地开展病原预防工作有着重要意义。

随着越来越多的细菌在健康和适应症改善中的促进作用被发现，对益生菌菌株以及粪菌移植(FMT)的效果研究也越来越受到重视。由于益生菌菌株的基因组序列往往已经明确，菌株分析方法可用于追踪益生菌菌株在肠道的留存状况。Zmora等通过检查宏基因组样本中是否包含益生菌菌株特异性基因来判断相关益生菌菌株在粪便样本中的有无，并由此分析益生菌菌株在不同个体肠道中的定植情况，然后从机制上阐明特定益生菌菌株是如何通过与肠道微生物互作来发挥效果的[17]。另外，尽管许多益生菌研究已经开始使用宏基因组研究益生菌对于肠道微生物群落的影响[61-62]或者进行益生菌的成分鉴定[63]，但是宏基因组菌株构成分析方法还没有被普遍使用。

粪菌治疗(FMT)是近年来兴起的肠道微生物治疗方法，该方法将健康人肠道微生物移植入患者体内，从而达到系统性改变患者肠道微生态的目的。在FMT治疗过程中, 医疗人员可以通过追踪供体粪便样本中的菌株在患者体内的留存和变化情况[20]，进行粪菌移植治疗效果的评估以及治疗方案的改进和优化[31]。此外，如果已知有益菌株的关键SNP或者基因，可以对供体样本进行评分和排序[64]。

随着微生物在疾病诊断和治疗方面的研究进展，不断有微生物被发现与人体疾病存有密切关联。一些微生物的组合已经用于体外诊断。但是，即便对于同一适应症(如肿瘤靶向药物治疗[65]和自闭症[66])，不同研究得出的疾病相关细菌可能并不一致，有的甚至相互矛盾。对于这种现象，一个可能的情况是当前主流的宏基因组分析工具主要在菌种水平进行分析，并没有考虑到同一菌种下菌株水平的差异。研究报道过许多菌种的菌株的确存在较大差异，例如P. copri作为人体常见的一类肠道微生物，与人体葡萄糖以及胰岛素耐受既表现出正相关[67]，也表现出负相关[68]关系。因此，对宏基因组数据从菌株水平进行进一步的挖掘，有助于发现与表型高度相关的特定菌株或者菌株组合。比如，Fang等基于MIDAS对一名克罗恩病患者的样本进行了菌株水平的时间序列分析，发现炎症状态与特定的E. coli菌株相关联[69]。

菌株分析方法可以从新的角度(特定SNP、基因)对病人表型进行分层。从机器学习和模型训练的角度看，这种处理方法相当于增加了宏基因组数据的特征数目, 而且新增加的特征更加反映真实世界的情况。对于疾病/健康样本、治疗响应/不响应样本的宏基因组菌株分析，有助于进行以菌株为靶点的治疗方案的开发或者将菌株本身作为潜在药物进行开发[70-71]。

研究表明，不同肠道微生物菌株会参与药物的代谢，从而对药物的药效及毒性产生影响[72-74]。宏基因组菌株分析技术的进步，可以帮助制药公司更好地对适用人群进行评估，并采取一定的微生物调控的措施，以提高治疗的效果和减少副作用。

随着测序技术的发展，菌株的分析方法和结果的质量也会迎来革新式的发展。首先，长读长第三代测序技术(PacBio、Oxford nanopore)及高通量染色质构象捕获技术(Hi-C)在宏基因组数据的应用将大幅提高组装的完整性和准确性[29, 75-76]，这些从菌群样本直接获得的高质量的基因组不仅可以用于无参的菌株分析，也可丰富补充现有基因组数据库，为有参的菌株分析提供更完整的参考数据。其次，这些新技术使还原质粒等可移动基因元件的菌株归属和研究其动态变化成为可能，为研究菌株的生态功能和进化发展提供了新的维度。届时，选择性地采用高质量的三代或Hi-C测序组装，并结合更廉价的二代测序提升研究规模，或将成为主流的微生物组研究策略，随之也将涌现出一批整合多种测序技术数据的生物信息算法及工具。另外，培养组学的快速发展，正极大地丰富着菌株基因组信息数据库。又因为通过对纯培养物进行的全基因组测序所获得的菌株基因组信息比通过宏基因组数据进行拼装所得的信息更为准确和可靠[77]，所以，培养组学将成为大量高质量菌株基因组的重要来源并将助力宏基因组数据分析迈入菌株水平时代。而随着计算能力的进一步发展，有参菌株分析方法可能会往基于全基因组数据信息方向发展，并且能同时进行大量菌种的菌株分析，进一步提高算法的通用性和结果的信息含量。

随着应用案例的增长，基于宏基因组测序的有参菌株分析算法在科研、临床等方面都证明了其重要价值。但菌株分析算法的开发，目前仍处于起步阶段，各方法都有主要的针对性使用场景和对应的基本假设。根据应用场景，大体可分为两类：第一类分析方法针对目的菌株进行分析。研究者需要提供菌株的有关信息，如全基因组序列，然后从样本中提取该菌株的定性和定量信息。如果已知相关菌株的背景信息, 第一类方法可以很好地为临床或者政策安排提供决策依据。最常见的应用是从环境或者宿主样本中对特定菌株(如病原菌菌株)的溯源和筛查(如StrainSifter)。这类方法对准确性要求高，但对结果的全面性要求较低。此外，由于肠道微生物菌株积极参与许多药物的代谢，这些预设的菌株信息可用于对相关适应症患者的分层，以优化临床治疗方案，实现精准治疗。如果菌株在适应症中的作用明确，第一类分析方法也可以在粪菌移植治疗中作为供体的指标对供体进行选择。第二类分析方法则能够从样本宏基因组数据中反映菌株的微观多样性(Microdiversity)，强调样本中有什么样的菌株构成，各自的比例是多少。这类方法可以监测菌株构成在多样本间的动态变化，应用于分析环境或者时间对菌株构成带来的影响(如母体到新生儿环境改变，抗生素干预或者FMT治疗对菌株的影响等)。通过与样本相关的表型信息相关联，第二类方法可用于对表型相关的菌株的筛选，在益生菌组合工艺优化、关键菌(如生物标记物、病原菌株、菌株药物)发现上具有潜在的应用价值。总而言之，在选用方法时，一定要注意算法的基本假设是否符合自身课题的特点，然后根据课题的需求进行针对性选择。

然而，目前仍缺少对宏基因组数据的菌株分析方法的详细评估，导致对算法的优缺点、精准度和可靠性了解甚少，对选择合适的方法造成了阻碍。建立专门针对菌株水平的评估体系和平台，对现有的菌株分析方法进行系统评估，将有助于研究者针对数据选择合适的方法。此外，由于现有方法对SNP和基因组成差异的割裂使用，单一算法无法同时获取菌株水平的丰度构成信息和功能基因信息。开发新的算法实现菌株构成和基因谱的关联分析，将进一步提高研究的可靠性和应用价值。