点击蓝字 关注我们
StrainPanDA：使用泛基因组分解技术对宏基因组数据进行菌株组成和关联的基因成分谱的重建
/10.1002/imt2.41
RESEARCH ARTICLE
●8月1日，中国科学院深圳先进技术研究院戴磊团队在iMeta在线发表了题为“StrainPanDA: linked reconstruction of strain composition and gene content profiles via pangenome-based decomposition of metagenomic data”的文章。
●该研究开发可同时提取共存菌株的组成和基因成分谱的菌株分析工具StrainPanDA，其所提供的菌株组成和基因成分谱的关联信息，进一步帮助我们深入理解微生物群落的适应性变化和菌株特异性功能之间的关系。
● 第一作者：胡函、谭宇翔、李陈浩
● 通讯作者：戴磊(lei.dai@);
谭验(yant@)
● 其他作者：陈俊宇、寇岩、徐振江、刘洋彧
●主要单位：中国科学院深圳先进技术研究院，深圳未知君生物科技有限公司，美国麻省综合医院及哈佛医学院
摘 要
具有功能多样性的种内不同菌株可以在微生物群落中共存。然而，目前的菌株分析方法无法通过分析宏基因组数据提供菌株组成与其基因成分谱之间的直接联系。在此，我们提出了StrainPanDA（菌株水平的泛基因组分解分析）。它使用多个宏基因组样品的泛基因组覆盖分布信息，同时获取微生物群落中共存菌株的组成和基因成分谱多样性。我们使用模拟数据集系统地验证了StrainPanDA的准确性和鲁棒性。为了证明以基因为核心的菌株分析方法的效果，我们应用StrainPanDA分析了婴儿以及接受粪菌移植治疗患者的肠道宏基因组样本。我们展示了菌株组成和基因成分谱的关联信息对于理解微生物群落的适应性变化与菌株特异性功能（例如营养利用和致病性）之间关系至关重要。最后，StrainPanDA对计算资源的要求极低，可以批量并行处理群落中的多个物种。简而言之，StrainPanDA可以应用于宏基因组数据以检测分子功能与微生物/宿主表型之间的关联，从而产生可供验证的假设，并在菌株或亚种水平上获得新的生物学见解。
亮 点
●StrainPanDA使用多宏基因组样本的泛基因组覆盖分布进行分析，同时获取微生物群落中共存菌株的组成和基因成分谱信息
●StrainPanDA能准确且可靠地推断模拟数据集中各菌株的组成比例和基因成分谱
●StrainPanDA所提供的菌株组成和基因成分谱的关联信息，进一步帮助我们深入理解微生物群落的适应性变化和菌株特异性功能（例如营养利用和致病性）之间关系
视频解读
Bilibili：/video/BV1fd4y1278T/
Youtube：https://youtu.be/4AiHIjZDrl4
中文翻译、PPT、中/英文视频解读等扩展资料下载
请访问期刊官网：http://www.imeta.science/
全文解读
引 言
大量证据表明同一菌种内可以有多个菌株在微生物群落中共存[1，2]。同一菌种共存的菌株存在基因成分的差异（即非核心基因组），这大多是由水平基因转移（HGT）产生的[3-5]。而非核心基因组的种内多样性可能导致显著的表型差异（例如在营养物质利用、致病性和抗生素耐药性等方面表现出不同），并在微生物跨环境适应过程中发挥重要作用[6-9]。此外，许多与宿主微生物组相关的健康状况已被发现是由个别菌株功能引起的[8，10-14]。
宏基因组测序提供了一种不需培养却仍能研究复杂微生物群落的组成和功能的方法，给微生物组研究带来了革命性的变化。当前常用的宏基因组数据定量分析工具主要提供物种层级的分类组成信息[15-18]。在已知的微生物分离株基因组迅速增加的同时[1，9，19，20]，宏基因组数据的高分辨率分析揭示了菌种内部显著的变化[21，22]。宏基因组菌株分析方法已被用于追踪菌株传播或扩散[23，24]，研究微生物菌株的群体遗传学[25]，以及对特定感兴趣的菌株进行分型[26-32]。
微生物菌株的基因成分谱决定了其生物学功能。迄今为止，大多数的菌株分析方法使用单核苷酸多态性（SNV）来鉴定菌株组成[25，28，29，33，34]。通过假设SNV单倍型与基因成分谱[4]之间存在关联，基于SNV的方法可以间接对菌种内基因成分的变化进行估算。然而，对于许多菌种而言，已有研究显示SNV单倍型并不能反映HGT引起的微生物遗传多样化[4]。相应的，目前基于泛基因组的方法可以推断出宏基因组样本中最优势菌株的基因成分谱[35]，但无法提供样品中共存菌株的组成丰度和基因成分谱信息。建立菌株组成和菌株基因成分谱的联系可以为研究微生物群落的适应性变化以及与宿主互作等方面提供关键线索，但目前基于参考基因组的生物信息学工具并没有解决这个问题。
为了满足日益增长的从宏基因组数据中推断菌株功能的需求，我们开发了一种名为StrainPanDA（菌株水平泛基因组分解分析）的新方法，通过使用宏基因组数据的泛基因组覆盖分布作为输入来同时获取共存菌株的组成和基因成分谱信息（图1）。我们通过利用多维度模拟数据集验证了StrainPanDA的性能。结果表明StrainPanDA能够从宏基因组数据中准确推断菌株的组成和基因成分谱。为了证明StrainPanDA在人类微生物组研究中的实际应用价值，我们分析了配对母婴的[36]以及接受粪菌移植（FMT）治疗的患者的[29，37]纵向肠道微生物组样本。我们发现StrainPanDA能够鉴别菌株特异性功能与微生物群落适应性变化（或宿主表型）之间的关联，从而在种内的水平上产生新的生物学见解和可供验证的分子机制假设。
结 果
通过泛基因组覆盖分布的分解推导菌株组成和基因成分
宏基因组数据中特定菌种的泛基因组覆盖分布由所有共存菌株的基因成分组成。如果多个宏基因组样品中菌株组成存在不同，原则上可以从泛基因组覆盖分布信息推导出样品内菌株的组成以及每个菌株的基因成分谱[38]。基于这种思路，StrainPanDA的主要算法旨在将基因家族丰度矩阵D分解为两个矩阵的乘积，即基因成分谱矩阵P和菌株组成矩阵S（图1A，详见方法）。此处，来自宏基因组数据的泛基因组覆盖分布由矩阵D表示，其中Dij是宏基因组样本j中基因家族i的计数的归一化结果。共存菌株的基因成分谱由二元矩阵P表示，其中元素Pij表示菌株j中存在/不存在基因家族i。样品中共存菌株的组成丰度由S矩阵表示，其中，Sij是菌株i在样本j中的相对丰度（和）。在StrainPanDA的具体实现中，基因家族丰度矩阵D通过非负矩阵分解（NMF）[39-41]求解矩阵P和S。这种处理使StrainPanDA能够同时描述共存菌株的组成和基因成分谱。
StrainPanDA软件提供了一套完整的菌株分析的工作流（图1B），包含从多个宏基因组样本中导入原始测序数据，进行短读序列映射，菌株分解以及下游注释（详见方法）。为了支持对菌株的基因成分变化的解读，StrainPanDA整合了几个数据库的功能注释，包括但不限于京都基因和基因组百科全书KEGG[42], 碳水化合物活性酶CAZy[43]和毒力因子数据库VFDB[44]。
图1. StrainPanDA工作流程图
（A）基因家族丰度矩阵D通过非负矩阵分解为P和S两个矩阵的乘积（见方法部分）。基因成分谱矩阵P是一个二元矩阵，表示每个菌株中基因家族的存在/不存在状态。菌株组成矩阵S表示每个样品中共存菌株的相对丰度。在示例中，宏基因组样品的大小S = 15，菌株数K = 6。（B）StrainPanDA分析的工作流程以逐个菌种分析的方式进行。每个菌种的分析包括将宏基因组短读序列映射到泛基因组数据库，菌株分解和基因家族谱的功能注释。CAZy，碳水化合物活性酶; GF，基因家族; KEGG，京都基因和基因组百科全书; StrainPanDA，菌株水平泛基因组分解分析; VFDB，毒力因子数据库
StrainPanDA在模拟数据中能准确预测菌株组成和基因成分谱
我们使用模拟的宏基因组数据（见方法）验证StrainPanDA的性能。对基于大肠杆菌的菌株模拟混合数据（从2到8个菌株，见方法），StrainPanDA预测的菌株组成总体上接近实际组成（Ground Truth）（图2A），并且在共存菌株数量较少时（2和4个菌株）性能更好。为了进行定量比较，我们计算了预测结果和实际菌株组成之间的Jensen-Shannon散度（JSD）和Matthews相关性系数（MCC），其中JSD和MCC已被广泛用于菌株分析工具的评估[28，30，33]。在共存菌株数量较少时（2和4个菌株），StrainPanDA预测的菌株组成优于主流的基于SNV的方法，包括StrainEst[30]和PStrain[33]（后者根据ConStrains进行了优化[28]）。虽然StrainEst的结果在共存菌株数量较少时更容易有假阳性，其性能在6和8个菌株时比其他方法更好。此外，我们生成了具有不同测序误差水平（图S2），不同测序深度（图S3）和不同背景噪声（与不同宏基因组数据集混合，表S1，图S4）的大肠杆菌多菌株混合的模拟数据集。与基于SNV的方法相比，StrainPanDA在预测菌株组成方面具有较好的鲁棒性。
图2.使用模拟宏基因组数据验证StrainPanDA的性能
（A） 在大肠杆菌的多菌株混合模拟数据（pWGS数据集，1×测序深度，见方法部分）中，把实际菌株组成（Ground Truth）与由StrainPanDA和现有工具（StrainEst和PStrain）预测的菌株组成进行比较。实际大肠杆菌菌株的数量（n = 2、4、6和8）按行显示。每个堆叠条对应一个模拟样本。菌株按相对丰度排序显示。如果预测菌株的数量超过实际菌株的数量，则将额外的菌株都归为“Extras”。（B）实际菌株和预测菌株组成之间的Jensen–Shannon散度（JSD）。No BG，没有背景；BG1/BG2/BG3，混合了三种不同宏基因组数据集作为背景的大肠杆菌的模拟数据（WGSBG数据集，100倍背景，见方法部分）。每个点代表一个模拟样本（n = 20）。（C） 评估了不同菌种的实际菌株组成和预测菌株组成之间的JSD。每个点代表一个模拟样本（n = 24）。缺失的结果标记为“NA”。（D）大肠杆菌菌株的实际和预测基因家族谱（使用的模拟数据与图A相同）。每行是一个基因家族，每列是一个菌株。聚类基于欧几里得距离。（E）评估了不同菌种共存菌株基因家族谱的精确召回曲线下面积（AUPRC）。每个点代表一个菌株（n = 4个菌株）的AUPRC。褐色点代表基因家族谱的随机猜测结果（见方法部分）。配对t检验得出的p值：*p < 0.05，**p < 0.01和***p < 0.0 1。长双歧杆菌、艰难梭菌、粪肠球菌、普拉梭菌和普雷沃氏菌。pWGS，纯全基因组的模拟数据；StrainPanDA，菌株水平泛基因组分解分析；WGSBG，含背景的全基因组的模拟数据。
为了评估StrainPanDA在不同菌种中的表现，我们生成了常见的人类肠道菌种（长双歧杆菌，艰难梭菌，粪肠球菌，普拉梭菌和普雷沃氏菌）的混合模拟数据（表S2，见方法）。与其他方法相比，在4菌株条件下，StrainPanDA对所有物种的菌株组成的预测都最准确（JSD为0.021±0.006）（图2C）。
虽然目前基于SNV的方法可以从宏基因组样品中获取共存菌株的组成，但它们不能直接预测菌株所对应的基因成分。在此，我们展示StrainPanDA能够同时获取每个宏基因组样品中的菌株组成以及菌株之间的基因成分差异。在大肠杆菌的模拟数据中，StrainPanDA预测的基因家族组成接近实际组成（图2D，图S5;在4个菌株的模拟数据中，精确度在0.91-0.96之间，召回率在0.87-0.96之间）。值得关注的是，StrainPanDA能够推断出未包含在参考基因组数据库中的菌株的基因家族谱。所有大肠杆菌菌株的精确召回曲线下面积（AUPRC）超过0.95，表明StrainPanDA能够以高灵敏度和高精度重建菌株的基因成分谱（图S6）
我们进一步评估了六个人类肠道菌种的模拟数据的预测基因成分谱结果。平均的AUPRC高于0.9，明显优于随机猜测（图2E）。预测的基因家族谱对不同的测序误差、测序深度和真实宏基因组数据背景都具有鲁棒性（表S4）。此外，为了证明StrainPanDA鉴定菌株特异性基因的能力，我们在模拟数据中将致病性大肠杆菌爆发菌株O104[45]与其他大肠杆菌菌株进行混合（表S5）。所有与疫情爆发相关的基因家族都被StrainPanDA成功预测重现（图S7）。最后，StrainPanDA和PanPhlAn/PanPhlAn3在推断基因成分谱方面的表现相当（图S8）。值得注意的是，PanPhlAn和PanPhlAn3只能报告每个宏基因组样品中最优势菌株（即相对丰度最高的菌株）的基因成分谱;相比之下，StrainPanDA可以鉴定所有共存菌株的基因成分谱。
综上所述，我们的评测结果表明，StrainPanDA可以准确预测宏基因组样品中共存菌株的组成谱和基因成分谱。在接下来的部分，我们将展示StrainPanDA在两套纵向宏基因组数据集中的应用，以阐明肠道微生物组在亚种水平上的多样性。
应用案例1：婴儿肠道微生物组中长双歧杆菌亚种的演替与母乳喂养模式和营养利用选择有关
在菌株水平上直接推断种群结构和菌株间基因成分谱差异对于理解微生物群落的生态学至关重要。在此，我们应用StrainPanDA来研究婴儿肠道微生物群落中共存亚种的适应性变化。我们分析了前人发表的数据集，其中包括来自约100对母婴的肠道宏基因组样本（婴儿在3个时间点采样：新生儿，4个月和12个月）[36]（表S6）。在物种水平上，作者发现婴儿肠道微生物群落的组成在每个采样时间点都具有独特的特征，并且母乳喂养的停止显然与从婴儿肠道微生物群落变化为类似成人的群落的成熟过程有关[36]。
本案例中我们专注于长双歧杆菌（Bifidobacterium longum）的亚种分析，因其在婴儿肠道微生物群落的发育中发挥重要作用[36，46-50]，并且在本研究中发现它在4个月的婴儿样本中富集（图S9）。有趣的是，我们发现长双歧杆菌的亚种组成随着时间的推移存在明显的演替模式（即优势亚种的变化）（图3A）。在StrainPanDA预测的3个长双歧杆菌亚种中，长双歧杆菌亚种2在母亲的肠道微生物群落中占主导地位。在婴儿的肠道微生物群落中，长双歧杆菌亚种3在新生儿中最为普遍，而亚种1在中间时间点（4个月）短暂增加，然后被亚种2（12个月）超越。根据原始研究中提供的饮食记录，我们进一步将婴儿样本分为两个不同的类别，即连续时间点间的“停止母乳喂养”和“继续母乳喂养”（见方法）。可以明显看出，长双歧杆菌亚种1在继续母乳喂养的婴儿中相对富集（图3B）。相比之下，一旦停止母乳喂养，长双歧杆菌亚种1的优势地位就被亚种2接替。
长双歧杆菌亚种组成与母乳喂养模式之间的关联提示菌种内可能存在营养利用功能上的竞争（图3C）。根据KEGG[42]和CAZy[43]的功能注释，我们发现各个长双歧杆菌亚种的营养利用基因成分谱存在明显差异。长双歧杆菌亚种1具有独特的基因家族（图3C中标记为红色）。这些基因家族与人乳寡糖（HMO）代谢的关键酶相关，包括半乳糖苷酶，α-L-岩藻糖苷酶，唾液酸酶及其相应的碳水化合物活性酶（GH33，GH29和GH95）（图S12-S13和表S7）。此外，尿素代谢相关的基因家族仅在亚种1中发现。因此，长双歧杆菌亚种1在利用母乳中的HMO和尿素[51]方面的独特功能潜力可以在母乳喂养条件下体现竞争优势，这与我们观察到亚种1在4个月婴儿肠道宏基因组中的短暂优势一致。
我们的发现与过往报道中关于断奶后婴儿的长双歧杆菌菌株中HMO利用基因特征[52]以及亚种频率的变化趋势[4，46，47，53，54]一致。与长双歧杆菌参考基因组相比，亚种1的功能特征表明它可能与婴儿长双歧杆菌亚种相对应（图S14），后者是从婴儿中分离出来的已知与母乳喂养有关的菌株[46，55，56]。
图3. 婴儿肠道微生物群中长双歧杆菌亚种的演替可归因于营养利用的选择
（A）来自多个时间点的母亲和婴儿（新生儿、4个月和12个月）的三个长双歧杆菌亚种组成结果的三元图。每个点代表一个样本。（B）连续时间点之间长双歧杆菌亚种相对丰度的变化。根据随后时间点的母乳喂养状态，将婴儿分为两组（紫色，停止母乳喂养，N = 86；红色，持续母乳喂养，N = 71）。Student t检验：****p<0.00005；ns，不显著；（C）预测长双歧杆菌亚种的基因家族谱。与宿主聚糖、尿素酶和CRISPR代谢相关的基因家族（KEGG注释）被挑选出来进行展示。左边的颜色条表示基因功能簇的类别。每行是基因家族的一个子类，亚种1中特有的基因家族用红色标记。热图中的色阶表示特定子类中的基因家族覆盖情况的归一化结果（即属于该子类的基因家族的测出率）。HMO，人乳寡糖；KEGG，京都基因和基因组百科全书。
总之，我们的结果表明StrainPanDA能够识别菌株特异性功能（通过基因成分谱的预测结果）和适应性变化（通过菌株组成的预测结果）之间的关联，从而在亚种水平上对微生物生态学产生新的生物学见解和可供验证的假设。
应用案例2：对FMT后的肠道宏基因组分析揭示个体化的亚种特征和亚种特有的功能
粪菌移植（FMT）将健康供体的细菌菌株引入受体，并对受体肠道微生物群落的结构和功能产生深远影响[57，58]。在此，我们应用StrainPanDA分析FMT治疗克罗恩病的临床试验中的宏基因组样本[37]，其中包括患者17例（FMT组8例，假手术组9例）和每个患者的多个样本（4-8个时间点）。我们分析了人类肠道宏基因组中常见细菌物种的亚种组成（表S8）。预测的亚种组成的分层聚类表现出明显的个体特征，并且这样的状况持续稳定了24周内（图4A）。同一个体（在多个时间点采样，即“个体内”组）的样品之间的亚种组成的两两成对距离显著低于不同个体（即“个体间”组）样品之间的成对距离（图4B，p值<10-15），在种水平上也可以看到类似的模式。此外，我们发现配对的FMT供体和受体（即“供体-受体”组）之间亚种组成的差异显著低于“个体间”组（p值 < 0.001），表明供体菌株的定植以及供体和受体菌株有共存现象[37]。我们还发现供体肠道细菌的定植在亚种水平（效应大小=1.4）比在菌种水平（效应大小= 0.78）更明显（图4B）。同样，我们应用StrainPanDA分析了一套独立的来自艰难梭菌感染患者FMT的宏基因组样本数据集[29]。我们在FMT后肠道宏基因组中观察到清晰的亚种定植的模式，与原始研究中基于SNV的菌株分析结果一致[29]（图S15）。总之，我们展示了StrainPanDA能够描述个体之间亚种组成的差异并跟踪菌株的传播。
为了阐明肠道微生物组在维持克罗恩病患者缓解中的潜在作用，我们进一步研究了与FMT后临床结果相关联的菌株水平遗传特征。原始研究指出，在FMT后复发的患者中拟杆菌属的菌种存在富集情况[37]。我们将分析重点放在卵形拟杆菌（Bacteroides ovatus）上，因为我们发现卵形杆菌在复发个体中富集（FDR校正后的p值= 0.2）（图4C-E）。在预测所得的卵形拟杆菌亚种中，亚种2的相对丰度与菌种整体的丰度变化呈正相关（斯皮尔曼相关系数 = 0.64，FDR校正的p值 < 10-6，图S16）。我们发现不同卵形拟杆菌亚种之间的基因成分存在很大差异（图4D）。有趣的是，我们发现卵形拟杆菌亚种2比其他亚种具有更多的CAZy家族基因，表明其利用多种碳源的功能潜力和潜在的竞争优势（图4E）。因此，卵形拟杆菌亚种2的菌株特异性代谢功能也许可以解释为什么它在菌种内占主导地位（~20%）及其与物种丰度的正相关（图4C和图S16）。此外，卵形拟杆菌亚种2携带几种菌株特异性毒力因子基因（例如，IV型分泌系统，胆固醇依赖性细胞溶解素），它们可能有助于解释卵形拟杆菌与FMT后复发之间的正相关关系（图S16）。例如，胆固醇依赖性细胞溶解素是一种能形成孔隙的毒素，可破坏宿主细胞质膜[59]。而宿主细胞质膜的完整性与炎症性肠病有关[60]。此外，我们注意到卵形拟杆菌亚种4在缓解组中更丰富（图4C，FDR校正的p值<10-5）;该卵形拟杆菌亚种是否有助于FMT后的缓解仍有待在未来的研究中进行验证。同样，我们对普通拟杆菌（Bacteroides vulgatus）进行了StrainPanDA分析。该菌种在复发个体中更为丰富，并且其亚种之间存在明显的功能差异（图S17）。
综上所述，StrainPanDA提供的菌株组成与基因成分谱的联系，可以极大地促进我们对种以下水平的微生物生态学的理解。对于与宿主健康密切相关的微生物，这种联系有助于形成关于分子功能（例如致病性）与临床结果之间的可供验证的假设，而这些验证工作可以在分离株的相关实验中进行。
图4. 肠道宏基因组分析揭示了FMT后个体化的亚种组成以及亚种特异功能和表型之间的关联
肠道宏基因组分析揭示了FMT后个体化的亚种组成以及亚种特异功能和表型之间的关联。（A）常见肠道菌种（表S8）的亚种组成预测结果的聚类显示出明显的个体特征。从原始论文[37]中收集的受试者ID按不同的颜色标记。（B）样本之间亚种组成和物种组成的差异（Bray–Curtis差异）。配对类型被分为三组进行比较：个体间（来自不同个体的样本）、个体内（来自相同个体的样本）和供体-受体（FMT供体与受体的FMT后样本）。（C）卵形拟杆菌及其亚种的相对丰度之间的关系（通过中心对数比转换进行归一化处理）。直线表示拟合线性回归（阴影区域：95%置信区间）。侧面的密度图显示了相应变量的分布。在菌种水平上，卵形拟杆菌在复发组中富集。（D）卵形拟杆菌各亚种的泛基因组预测结果统计。（E）预测的卵形拟杆菌亚种中碳水化合物活性酶（CAZy）的基因家族谱。所有亚种都共享的CAZy基因未被显示。FMT，粪菌移植。
讨 论
本文报告了一种名叫StrainPanDA（菌株水平泛基因组分解分析）的新方法，能从宏基因组数据中同时获取共存菌株的组成谱及菌株相应的基因成分谱。我们的评估结果表明，StrainPanDA对模拟数据进行了准确而稳健的预测。StrainPanDA预测的菌株组成与主流方法相比结果相近或更优；同时，即使对于不在参考基因组数据库中的菌株，其预测的基因成分谱也很接近真实基因成分谱。此外，我们将StrainPanDA应用于真实的宏基因组数据集，以了解这些真实案例中细菌种内多样性与生物学意义的关联。例如，我们发现，婴儿肠道微生物群落中长双歧杆菌亚种的组成变化与饮食变化有关，且特定的长双歧菌亚种在利用母乳营养方面的独特功能潜力可能会为其带来竞争优势。我们还展示了通过对菌株组成和菌株基因成分谱的关联分析能够形成直接的功能解释以及可供验证的假设。
为了研究种内基因成分谱差异，当前基于SNV的方法假设SNV谱与基因成分谱之间存在直接关联。然而，许多在核心基因组中高度相似的微生物基因组的共有基因不到70%[3]，这表明基于SNV推断基因成分谱可能并不可靠。相反，StrainPanDA采用基于泛基因组的方法能直接推断物种内多个共存菌株的基因成分谱。我们目前的方法依赖于从指定微生物物种的基因组集合中构建泛基因组，因此它不考虑物种之间的基因转移。为了解释种间基因转移，一个可能的解决方案是包括一个“假定的移动元素”库，以扩展每个物种的泛基因组。StrainPanDA的预测依赖于泛基因组数据库，但不限于库内可用的参考基因组；因此，只要泛基因组中包含相关基因家族，StrainPanDA也能识别新菌株中相应的基因家族分布。尽管我们重点比较了StrainPanDA与其他基于参考数据库的方法，但需要指出的是，基于宏基因组进行基因组组装（MAG）的方法也可以从宏基因组数据中进行菌株分析。例如，DESMAN[61]可以预测每个菌株的基因组草图；其他最近开发的基于MAG的菌株分析方法包括mixtureS [62]和STRONG [63]等。与基于参考数据库的方法相比，基于MAG的方法可以识别新的物种和基因，但MAG的质量将极大地影响结果。此外，基于MAG的方法比基于参考数据库的方法需要更高的测序深度，这阻碍了其应用于低丰度的物种。相比之下，测序深度不是StrainPanDA的限制因素（图S3-S4）。
StrainPanDA最适合分析具有共享菌株的多个宏基因组样本，例如纵向数据研究。虽然本研究中的分析侧重于人类肠道微生物，但只要目标物种的泛基因组可用，StrainPanDA也能拓展到不同环境中的微生物群落。StrainPanDA的性能，包括预测菌株组成和基因成分谱的准确性，随着样本量（图S18）和测序深度（图S3）的增加而提高。要在典型的宏基因组数据集上应用StrainPanDA，需要至少有10个样本，并且有关菌种的相对丰度高于1%。由于StrainPanDA算法的特性，可能很难区分开具有遗传嵌合性的或基因成分谱高度相似（即缺乏菌株特异性特征）的菌株，因此StrainPanDA最适合在亚种水平上进行分析[3]。最后，与基于MAG的方法相比，StrainPanDA对计算资源的要求极低（图S19），并且可以扩展至并行处理多个物种。
结 论
StrainPanDA能够从宏基因组数据中准确分析菌株组成和基因成分谱。我们预期，将StrainPanDA应用于快速增长的宏基因组数据集，特别是在微生物组的时空特征[64-67]的背景下，将有助于阐明分子功能与微生物/宿主表型之间的新关联，以及种内水平的微生物生态学。
代码和数据可用性
StrainPanDA的源代码已在GitHub(/xbiome/StrainPanDA)和Zenodo(/10.5281/zenodo.6668661)发布. 其中，和真实实验数据分析相关的源代码在/xbiome/StrainPanDA-data/tree/main/example#readme. 所有的补充材料（文本、图、表、中文翻译版本或视频）也可从线上获取。
方法细节
泛基因组数据库生成以及宏基因组数据的映射
本研究中使用的菌种的泛基因组采用PanPhlAn（版本1.2.8）推荐的步骤来生成[35]。来自每个菌种的基因组从美国国立生物技术信息中心（NCBI）下载。基因组之间的平均核苷酸相似度（ANI）使用mash工具进行计算 （版本1.1）[68]。代表菌株使用基因组间两两成对ANI ≤ 99%作为标准进行挑选，并用作参考基因组来进行泛基因组的构建。已经注释的基因从参考基因组中提取出来，并使用usearch（v7）按照95%相似性聚类成基因家族，由此构建泛基因组数据库。鸟枪法宏基因组数据使用PanPhlAn（版本1.2.8）[35]映射到泛基因组数据库上，其中PanPhlAn分别采用Bowtie2（版本2.4.1）[70]和SAMtools（版本0.1.19）[71]工具对短读序列进行映射和计数。将属于同一基因家族的基因的读数（采用每百万短读序列中来自于某基因每千碱基长度的短读序列数RPKM进行归一化）进行加和而得到基因家族谱，然后将所有宏基因组样本的基因家族谱进一步合并以得到统一的基因家族谱矩阵。为了解决序列匹配中的潜在噪音问题，当RPKM值低于阈值的时候（默认值是10），基因家族丰度截取为0。截取之后，所有样本中都不存在的基因家族被移除，不再用于未来分析。此外，如果检测到的基因家族数目低于0.9 × gmin (gmin为所有参考基因组的基因家族数目中最小的值)，该样本被过滤掉。
StrainPanDA算法
StrainPanDA核心算法工作是将有关菌种的基因家族丰度矩阵D分解为两个矩阵的乘积（图1）。
其中基因家族丰度矩阵D是一个N×S非负矩阵，而Dij是宏基因组样本j中基因家族i的计数的归一化值。基因成分谱矩阵P是一个N×K的二元矩阵，当基因家族i存在于菌株j时，Pij为1，否则为0。菌株组成矩阵S是K×S矩阵，而Sij是菌株i在样本j中的相对丰度（并且）. N是有关菌种的泛基因组的基因家族数目，而S是宏基因组样本的数目，K是菌株数目（即因式分解的秩）。
为了估算P和S，我们使用非负矩阵分解（NMF）得到近似解P'和S'，这样就考虑了两个矩阵的非负限制条件（优化过程采用了R包NMF[40,41]版本0.21.0中的snmf/r算法）。稀疏性限制(即目标函数中的正则化项)的加入保障了因子分解结果的唯一性[41,72]。矩阵S'随后被转化为相对丰度。我们对矩阵P'进行二元化处理，基于的假设是“存在”的基因家族的矩阵元素的值应该高于“不存在”的基因家族的元素的值，而且由于P是二元矩阵，矩阵元素应该有一个比较紧凑的分布（图S1B）。简言之，对每个菌株j, 我们寻找矩阵元素形成的概率密度曲线的峰(pmax），使得峰右侧(Pij > pmax)的矩阵元素数目等同于有关菌种的基因家族数目的期望值 (即泛基因组数据库中种下面各参考基因组的基因家族数目的平均)。我们随后在0与峰值之间选取一个θ (θ 默认值为0.5 × pmax)将密度曲线一分为二，其中基因家族权值若大于θ则被视为存在。样本j中的基因家族i采用如下置信分数Cij进行计算：
置信分数被用于对基因出现预测进行排序，并且在评估实验中用于精度召回曲线的生成。
为了选择合适的菌株数目(即NMF的秩)，我们在默认1-12的范围内拣选满足以下条件的最小的菌株数目：（1）任一菌株在所有样本中的平均相对丰度应该大于t2（t2默认值为0.1），（2）所有菌株的基因家族的数目应该大于t3 × gmin （t3默认值为0.5），gmin 是所有参考基因组的基因家族数目的最小值，（3）两两菌株间的基因家族谱的Jaccard距离应当大于t1（t1默认值是0.1）。程序同时提供了一个选项允许用户自行指定菌株数目。在本研究中，我们并没有在StrainPanDA的评估实验与案例分析中预设菌株的数目。
通过模拟数据评估StrainPanDA
基于大肠杆菌菌株的模拟数据。我们生成了四种类型的模拟数据：1）无错误（ErrFree）：通过ART模拟器[73]（版本.06.05；参数：-ef -ss HS25 -l 150 -m 270 -s 27）从大肠杆菌的参考基因组中随机选取片段。2）ART模拟的测序错误（ARTErr）：使用ART在无错误数据的基础上添加序列错误（参数：-ss HS25 -l 150 -m 270 -s 27）。3）WGS的真实测序错误（pWGS）：通过seq-tk(/lh3/seqtk，版本1.3；默认参数）从选定菌株的全基因组测序数据中随机抽取片段；4） pWGS数据与真实宏基因组数据集混合的模拟数据（WGSBG）：使用三种不同的宏基因组数据集（见表S4。BG1:IBD、BG2:FMT、BG3:MI，如图2所示）以不同比例（1、5、10、25和100倍）与大肠杆菌的pWGS数据混合。通过Kraken2分析宏基因组样本（版本2.1.1；数据库：minikraken2_ v2_8GB_04），以确保大肠杆菌的丰度下限符合要求。选择菌株时，挑选具有95%-99%的全基因组配对ANI的大肠杆菌菌株来代表不同的亚种（表S3）。在每个模拟数据集中，我们通过Dirichlet分布随机生成20个比例组合（即20个样品）（表S2）。所有模拟数据集均由SimStr流程生成(/xbiome/StrainPanDA/tree/main/SimStr). 对于每个菌株，其基因组大小被视为1×测序深度，并用于计算产生的短读序列数量。将菌株的最小相对丰度（即出现频率）设置为5%做为基本单位。例如，对于我们在本研究中提到的1×测序深度的大肠杆菌模拟数据，频率为5%的菌株的数据大小约为4.5兆碱基（MB），而该数据集中每个样本的总深度始终为20×, 大小约90MB（1×测序深度为一个单位，共20个单位）。为了评估样本量的影响，我们还通过Dirichlet分布生成了一个包含4个菌株共400个组合样品的模拟数据集。将该模拟数据分10次运行（每次40个样本），并在每次运行中分别进一步降低采样数量到20、15、10和5个样本，从而判断不同样品量对结果的影响。
不同肠道菌种的模拟数据。我们分别生成了6个菌种的模拟数据，包括长双歧杆菌、艰难梭菌、大肠杆菌、粪肠球菌、普拉梭菌和普雷沃氏菌（Sync6，5×测序深度的pWGS）（表S2）。对于每个物种，将4个菌株按5%、10%、25%和60%的相对丰度进行排列组合，总共生成24个样本（表S2）。所有6个物种都在预先建立的StrainPanDA和PStrain数据库中。StrainEst的预建数据库中只有长双歧杆菌、大肠杆菌和粪肠球菌，因此与StrainEst相比时，其他3个物种未被考虑（图2C）。
菌株组成结果的评估。对StrainEst（从docker获取的v1.2.4版本）和PStrain（ 5月23日从GitHub下载）分别使用其预建数据库和默认参数运行（StrainEst：ftp://ftp.fmach.it/metagenomics/strainest/ref/; PStrain：/wshuai294/PStrain). 通过SimStr比较了不同方法的菌株组成分布结果。对于图2A（堆叠条形图）所示的菌株组成，相对丰度低于0.01的菌株结果被认为是噪音进行过滤，剩余菌株按最终的相对丰度（重新归一化为1后）按降序排序。排序后，超过模拟菌株数量的被预测菌株会被归类为“Extras”。
我们使用两种常用指标来评估不同方法的菌株组成的结果：
（1）Jensen-Shannon散度（JSD）[74]：通过phyloseq[75]（R包）计算菌株组成结果和实际菌株组成之间的JSD。其中，菌株的匹配关系按相对丰度的降序进行排列。如果预测的菌株数量与实际菌株数量不同，则向维数较低的向量添加零。JSD是对称的，并且在[0,1]的区间内。它反映了菌株组成分布之间的非相似性，这意味着JSD=0时表示预测的数值完全准确。
（2）Matthews相关系数[76]（MCC）：
其中TP是真阳性数，TN是真阴性数，FP是假阳性数，FN是假阴性数。MCC的范围为-1到1，其中1表示预测结果完全一致，0表示预测结果是随机的，而-1表示完全不一致。
由于PStrain的结果没有菌株注释，我们计算MCC时仅计算由StrainPanDA和StrainEst预测的菌株组成。对于StrainEst，预测菌株的注释由参考基因组直接产生。对于StrainPanDA，我们计算预测菌株的基因谱和所有参考基因组基因谱之间的Jaccard距离(JD(A,B)=1-|A∩B|÷|A∪B|)，并使用与预测菌株的Jaccard距离最小的参考基因组做为注释。如果模拟数据中的菌株在参考基因组数据库中存在，则直接将参考基因组的ID与实际菌株进行比较，以确定预测菌株是否为真阳性。如果模拟数据中的菌株不在参考基因组数据库中，我们使用系统发育树来确定预测菌株是否为真阳性（表S4）。简单来说，我们使用parsnp[77]（版本1.5.1，默认参数）生成一个系统发育树，包括了模拟数据中使用的菌株的基因组和所有参考基因组。如果预测菌株的注释与实际菌株的系统发育树距离在阈值内（对大肠杆菌而言，阈值 = 0.05，对应于99% ANI），则认为预测结果为真阳性。
评估基因家族谱的预测结果。我们对于模拟数据集中的每个菌株，根据从NCBI下载的参考基因组，通过ART模拟器[73]（Sync-Single数据集）生成5×测序深度的无错误短读序列。Sync-Single数据集中每个菌株都有三个重复，然后通过PanPhlAn和PanPhlAn3-v3.1（敏感模式的默认参数）生成每个菌株的实际基因家族谱。在两个或多个重复中发现的基因家族被视为“存在”。将每个菌株的实际基因家族谱（参考谱）与预测的基因家族谱进行比较（图2D），计算菌株的预测基因家族谱与其参考谱（或参考基因组的基因家族谱的随机抽样）之间的Jaccard距离（图S5）。每个菌株的基因家族谱的精确召回曲线由R包PRROC[78]生成，其中使用置信分数对预测的基因家族进行排序（图S6）。对于随机猜测组，通过从泛基因组中抽取N个基因家族作为“存在”来生成1000个随机基因家族谱，其中N是参考基因组中每个菌株的基因家族数的平均数。为了证明StrainPanDA识别菌株特异性基因的能力，将致病性大肠杆菌菌株O104（GCF_002983645）引入4个菌株的模拟数据集中（Sync O104数据集，pWGS，5×测序深度）。我们使用Scholz等人[35]整理的疫情爆发相关基因来评估StrainPanDA的基因成分预测（表S5）。我们还将StrainPanDA预测的基因家族图谱与PanPhlAn和PanPhlAn3的预测进行了比较（图S8）
运行时间的评估。StrainPanDA的运行时间使用Linux中的time命令来测量。所有这些测试都在Intel（R）Xeon（R）Gold 6238 CPU@2.10GHz和16 GB内存的工作站上运行。在4株大肠杆菌菌株的模拟数据（见方法部分中的大肠杆菌菌株模拟数据）的抽样数据上运行StrainPanDA，以计算样本量与运行时间（秒）的关系。在pWGS和25倍WGSBG数据集上运行StrainPanDA，以计算菌株数与运行时间（秒）的关系。
StrainPanDA在真实案例上的应用
母婴肠道宏基因组数据。ERP005989的所有可用样本均从EBI下载（其中8个样本失败，表S6）。根据饮食信息排除无饮食信息的婴儿。将剩下的84名婴儿分为3个不同组（B_F_F：4个月时停止母乳喂养；B_B_F:12个月时终止母乳喂养；B_B_B：持续母乳喂养；不包括F_B_M样本）（表S6）。长双歧杆菌基因家族覆盖度不足的样本在预处理步骤被StrainPanDA过滤，并从下游分析中排除。“持续母乳喂养”组包括连续时间点（即新生儿和4个月之间，或4个月和12个月之间）保持母乳喂养的婴儿。“停止母乳喂养”组包括在4个月或12个月前停止母乳喂养的婴儿。为了对亚种进行功能解释，我们将注释到同一KEGG家族（9月1日下载）或CAZy家族（7月31日下载）的基因家族进行归组。为了进一步分析与母乳喂养相关的关键功能，我们从相关参考文献[51,79]中整理了一组KEGG直系同源家族。根据文献[51,79]（表S7），所有KEGG直系同源家族可以进一步分为子类和类。基因家族覆盖值按每个子类中的被检出基因占总基因数的比例进行计算（图3C）。
FMT供体-受体宏基因组数据。原始数据从ENA下载（注册号：PRJNA625520源于克罗恩病研究[37]，PRJEB23524源于艰难梭菌感染研究[29]）。对于克罗恩病数据集，使用Kraken2[18]（版本2.0.8-beta）和miniKraken数据库（v2_8GB_04_beta）估计物种相对丰度。为了确定与缓解或复发相关联的物种，我们计算每个受试者不同时间点的平均相对丰度，并使用Wilcoxon秩和检验来选择差异菌种（FMT之前或复发后收集的样本被丢弃）。随后使用StrainPanDA预测亚种的相对丰度，并通过中心对数比变换进行归一化，以计算与物种丰度的Spearman相关性。对于亚种的功能注释，毒力因子和CAZy注释取自前述构建的菌种特异的数据库。
基因家族的功能注释
为了获得KEGG注释的基因家族，使用usearch[69]（v11.0.667；'-ublast'）把基因家族代表序列和KEGG直系同源信息（-07-20版）进行比对。保留相似度>50%和查询覆盖度>50%的比对结果。为了获得CAZy注释的基因家族，使用seqkit[80]（0.15.0）将基因家族序列翻译成六个开放阅读框。翻译后的氨基酸被run_dbcan[81]（2.0.11）用作输入。run_dbcan使用DIAMOND[82]（2.08）、HMMer[83]（3.3.2）和Hotpep[84]（2.08）以及默认参数来预测CAZy注释。只有当至少两个程序预测出CAZy注释时，该注释才会被选择。如果一个基因家族被多个不同的CAZy注释匹配，则只使用第一个注释。为了注释毒力因子，通过DIAMOND[82]（2.0.8）（blastp，查询覆盖率>50%，相似度>50%）将基因家族中心序列和毒力因子数据库（VFDB， 4月9日）进行匹配。
引文格式：
Hu, Han, Yuxiang Tan,Chenhao Li, Junyu Chen, Yan Kou, ZhenjiangZech Xu, Yang‐Yu Liu, Yan Tan, and Lei Dai. .“StrainPanDA: Linked reconstruction of straincomposition and gene content profiles viapangenome‐based decomposition of metagenomicdata.” iMeta. e41
作者简介
胡函（第一作者）
●未知君首席生信科学家，深圳市海外高层次人才“孔雀计划”获得者
●研究方向涉及菌株分析算法，微生物组代谢组等多组学的生物信息分析，疾病知识图谱，以及微生物在自闭症改善中的作用机制
谭宇翔（第一作者）
●中国科学院深圳先进技术研究院助理研究员。广东省“珠江人才计划”海外青年引进计划和深圳市海外高层次人才“孔雀计划”获得者。
●研究方向为针对微生物组的生物信息学与大数据科学。当前研究主要结合宏基因组测序和培养组技术，进行深入至菌株层级的生物信息学分析。相关学术成果已发表于iScience、iMeta、Journal of Integrative Plant Biology、Journal of Functional Foods等期刊。
李陈浩（第一作者）
●麻省综合医院及博德研究所博士后
●主要研究方向：计算生物学、宏基因组、单细胞及空间转录组。以第一作者或通讯作者在Nature Medicine, Microbiome, PLOS Computational Biology等期刊发表论文
戴磊（通讯作者）
●中国科学院深圳先进技术研究院研究员、博士生导师，合成微生物组研究中心主任
●国家重点研发计划青年项目负责人，入选《麻省理工科技评论》中国区“35 岁以下科技创新 35 人”。戴磊实验室在定量生物学、合成生物学与微生物组学的交叉领域进行原创性研究，致力于对宿主共生微生物组的结构和功能进行精准表征和调控，解决人体健康、农业生产等重大问题
谭验（通讯作者）
●深圳未知君生物科技有限公司CEO、创始人，并受聘担任中国科学院微生物研究所生物工程专业硕士研究生企业导师等社会职务
●深圳市海外高层次人才“孔雀计划”获得者，《麻省理工科技评论》“35岁以下科技创新35人”。
更多推荐
（▼ 点击跳转）
iMeta文章中文翻译+视频解读
iMeta封面 | 宏蛋白质组学分析一站式工具集iMetaLab Suite(加拿大渥太华大学Figeys组)
▸▸▸▸
iMeta | 东农吴凤芝/南农韦中等揭示生物炭抑制作物土传病害机理
▸▸▸▸
iMeta | 华南农大陈程杰/夏瑞等发布TBtools构造Circos图的简单方法
▸▸▸▸
iMeta | 叶茂/时玉等综述环境微生物组中胞内与胞外基因的动态穿梭与生态功能
▸▸▸▸
iMeta | 南农沈其荣团队发布微生物网络分析和可视化R包ggClusterNet
▸▸▸▸
iMeta | 华南师大王璋组综述人体肺部微生物组与人类健康和疾病之间的隐秘关联
▸▸▸▸
iMeta | 南科大夏雨组纳米孔测序揭示微生物可减轻高海拔冻土温室气体排放
▸▸▸▸
iMeta | 北大陈峰/陈智滨等发表口腔微生物组研究中各部位取样的实验方法(Protocol)
▸▸▸▸
iMeta | 华南农大曾振灵/熊文广等-家庭中宠物犬与主人耐药基因的共存研究
▸▸▸▸
iMeta | 深圳先进院马迎飞组开发基于神经网络分析肠道菌群的方法
▸▸▸▸
iMeta | 南医大陈连民等综述从基因组功能角度揭示肠菌对复杂疾病的潜在影响
期刊简介
“iMeta” 是由威立、肠菌分会和本领域数百位华人科学家合作出版的开放获取期刊，主编由中科院微生物所刘双江研究员和荷兰格罗宁根大学傅静远教授担任。目的是发表原创研究、方法和综述以促进宏基因组学、微生物组和生物信息学发展。目标是发表前10%(IF > 15)的高影响力论文。期刊特色包括视频投稿、可重复分析、图片打磨、青年编委、前3年免出版费、50万用户的社交媒体宣传等。2月正式创刊发行！
联系我们
iMeta主页：http://www.imeta.science
出版社：/journal/2770596x
投稿：/imeta
邮箱：office@imeta.science
微信公众号
iMeta
责任编辑
微微
如果觉得《iMeta | 深圳先进院戴磊组开发可同时提取共存菌株的组成和基因成分谱的菌株分析工具...》对你有帮助，请点赞、收藏，并留下你的观点哦！
