点击蓝字 关注我们
StrainPanDA：使用泛基因组分解技术对宏基因组数据进行菌株组成和关联的基因成分谱的重建
● StrainPanDA使用多宏基因组样本的泛基因组覆盖分布进行分析，同时获取微生物群落中共存菌株的组成和基因成分谱信息
● StrainPanDA能准确且可靠地推断模拟数据集中各菌株的组成比例和基因成分谱
● StrainPanDA所提供的菌株组成和基因成分谱的关联信息，进一步帮助我们深入理解微生物群落的适应性变化和菌株特异性功能（例如营养利用和致病性）之间关系
视频解读
Bilibili：https://www.bilibili.com/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]。
图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
为了评估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：婴儿肠道微生物组中长双歧杆菌亚种的演替与母乳喂养模式和营养利用选择有关
本案例中我们专注于长双歧杆菌（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]。
图4. 肠道宏基因组分析揭示了FMT后个体化的亚种组成以及亚种特异功能和表型之间的关联
肠道宏基因组分析揭示了FMT后个体化的亚种组成以及亚种特异功能和表型之间的关联。（A）常见肠道菌种（表S8）的亚种组成预测结果的聚类显示出明显的个体特征。从原始论文[37]中收集的受试者ID按不同的颜色标记。（B）样本之间亚种组成和物种组成的差异（Bray–Curtis差异）。配对类型被分为三组进行比较：个体间（来自不同个体的样本）、个体内（来自相同个体的样本）和供体-受体（FMT供体与受体的FMT后样本）。（C）卵形拟杆菌及其亚种的相对丰度之间的关系（通过中心对数比转换进行归一化处理）。直线表示拟合线性回归（阴影区域：95%置信区间）。侧面的密度图显示了相应变量的分布。在菌种水平上，卵形拟杆菌在复发组中富集。（D）卵形拟杆菌各亚种的泛基因组预测结果统计。（E）预测的卵形拟杆菌亚种中碳水化合物活性酶（CAZy）的基因家族谱。所有亚种都共享的CAZy基因未被显示。FMT，粪菌移植。
讨 论
为了研究种内基因成分谱差异，当前基于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(https://github.com/xbiome/StrainPanDA)和Zenodo(https://doi.org/10.5281/zenodo.6668661)发布. 其中，和真实实验数据分析相关的源代码在https://github.com/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）。
置信分数被用于对基因出现预测进行排序，并且在评估实验中用于精度召回曲线的生成。
为了选择合适的菌株数目(即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]（版本2016.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(https://github.com/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_201904），以确保大肠杆菌的丰度下限符合要求。选择菌株时，挑选具有95%-99%的全基因组配对ANI的大肠杆菌菌株来代表不同的亚种（表S3）。在每个模拟数据集中，我们通过Dirichlet分布随机生成20个比例组合（即20个样品）（表S2）。所有模拟数据集均由SimStr流程生成(https://github.com/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（2021 5月23日从GitHub下载）分别使用其预建数据库和默认参数运行（StrainEst：ftp://ftp.fmach.it/metagenomics/strainest/ref/; PStrain：https://github.com/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）：
● 未知君首席生信科学家，深圳市海外高层次人才“孔雀计划”获得者
● 研究方向涉及菌株分析算法，微生物组代谢组等多组学的生物信息分析，疾病知识图谱，以及微生物在自闭症改善中的作用机制
● 麻省综合医院及博德研究所博士后
● 深圳未知君生物科技有限公司CEO、创始人，并受聘担任中国科学院微生物研究所生物工程专业硕士研究生企业导师等 会职务
● 深圳市海外高层次人才“孔雀计划”获得者，《麻省理工科技评论》“35岁以下科技创新35人”。
更多推荐
（▼ 点击跳转）
iMeta文章中文翻译+视频解读
iMeta封面 | 宏蛋白质组学分析一站式工具集iMetaLab Suite(加拿大渥太华大学Figeys组)
????
iMeta | 华南农大陈程杰/夏瑞等发布TBtools构造Circos图的简单方法
????
iMeta | 南农沈其荣团队发布微生物 络分析和可视化R包ggClusterNet
????
iMeta | 南科大夏雨组纳米孔测序揭示微生物可减轻高海拔冻土温室气体排放
????
iMeta | 华南农大曾振灵/熊文广等-家庭中宠物犬与主人耐药基因的共存研究
????
iMeta | 南医大陈连民等综述从基因组功能角度揭示肠菌对复杂疾病的潜在影响
联系我们
iMeta主页：http://www.imeta.science
出版 ：https://onlinelibrary.wiley.com/journal/2770596x
投稿：https://mc.manuscriptcentral.com/imeta
邮箱：office@imeta.science
iMeta
微微
文章知识点与官方知识档案匹配，可进一步学习相关知识OpenCV技能树首页概览11206 人正在系统学习中
声明：本站部分文章及图片源自用户投稿，如本站任何资料有侵权请您尽早请联系jinwei@zod.com.cn进行处理,非常感谢！