iMeta | 中科院微生物所王军组-MetaSVs：结合长、短reads用于宏基因组结构变异分析和可视化流程

MetaSVs:一个结合长、短reads用于宏基因组结构变异分析和可视化流程

iMeta主页：http://www.imeta.science

结构变异(structural variants, SVs，包括大片段的插入、缺失、倒置和易位)显著影响微生物基因组中基因的功能，微生物组中的SVs与多种生物过程和人类疾病相关。随着测序和生物信息学技术的进步，越来越多的测序数据和分析工具已经被广泛用于微生物组SVs分析，导致对更专用的SVs分析工作流程的更高需求。此外，由于各种测序技术独特的检测偏差，包括短reads测序(如Illumina平台)和长reads测序(如Oxford Nanopore和PacBio)，基于多个平台的SV发现对于全面识别各种SVs是必要的。在这里，我们建立了一个整合了Nanopore长reads和Illumina短reads的MetaSVs流程来分析微生物基因组中的SVs，并进一步确定可以反映代谢差异的差异SVs。我们的流程为研究人员提供了在不需要特定技术专长的情况下轻松获得微生物基因组中的SVs和相关代谢通路，这对于对宏基因组中的SVs感兴趣但缺乏复杂生物信息学知识的研究人员特别有用。

引  言

微生物组是指在人体器官和许多其他类型的环境生态位中的复杂微生物群落，在改善人类和植物的健康状况、营养获取和免疫反应、影响动植物健康和改变环境生化特性等方面具有重要作用。其中，SVs是微生物基因组中高度可变的片段，影响基因组分子和细胞过程、调控功能、3D结构以及基因表达。SVs会使来自同一物种的细胞表现出复杂的表型和不同的生长速率，进而导致不同的适应和占据不同的生态位。此外，多项研究表明，在宏基因组基因组中SVs表现出显著的个体间变异性和高度的个体内稳定性，因此可以作为个性化的“微生物组指纹”来区分属于相同或不同个体的宏基因组样本。

尽管以上的这些证据都可以看出SVs的重要性，但目前对SVs的鉴定主要是通过基于短reads的宏基因组分析来完成的，这限制了它们检测的准确性。然而，长reads测序技术的快速发展使得产生数千个碱基对的reads成为可能，牛津纳米孔技术(Oxford Nanopore Technology, ONT)甚至可以达到2 Mbp的reads长度。特别是对于不能被短reads覆盖的SVs，可以通过覆盖相应基因组区域的ONT reads更直接地进行验证。而且ONT reads是足够长的，可以跨越细菌基因组上的大部分重复区域，从而显著提高基因组组装的连续性。但是由于ONT reads的高错误率导致无法区分<10%差异的重复拷贝。因此，它很可能会影响基因组组装质量，导致一些潜在的错误。近年来，越来越多的微生物组研究人员将ONT长reads和Illumina短reads结合起来，得到了更完整的微生物基因组，解决了串联重复序列和SVs等复杂区域，为SVs分析提供了前所未有的潜力。

宏基因组SVs的分析通常需要一系列复杂的分析，这可能会让很多想要进行此类分析的研究人员感到困惑。目前，有几种合适的生物信息学工具可以执行这些步骤，但是对于如何将多个分析软件的结果整合成一个功能集成的工作流程，还没有系统的研究。设计工作流程需要大量的生物信息学知识、专业知识、资源和工具。特别是对于没有编程或没有生物信息学专业知识的科学家存在许多障碍，例如，前期工具安装，确定各种参数值，后期如何在有效地将多个工具整合在一起等等。再加上数以百万计的测序数据被用作源输入，许多分析步骤，如序列比对，基因组组装和bins (high‐quality draft genomes)分箱，既耗时又参数繁重，急需一个高效，可追溯而且灵活的工作流程。在这里，我们已经建立了一个集成的流程MetaSVs，用于微生物组SVs分析，使用ONT长reads和Illumina短reads的混合数据来综合分析宏基因组SVs (图1)。通过在Linux命令行运行一个Python脚本，就可以获得高质量的宏基因组组装基因组( high‐quality metagenome‐assembled genomes，MAGs)，检测基因组SVs的数量及分布，并进一步确定与潜在代谢功能相关的差异SVs。

图 1. MetaSVs的分析流程图

结  果

MetaSVs的配置

MetaSVs的源代码和示例结果已存放在GitHub。该流程的配置教程可以在这里找到: https://github.com/Wlab518/SV_procedure。

MetaSVs的输入文件

图2. MetaSVs的主程序

(A)主程序的唯一输入参数(config.ini文件)。红色字体表示需要针对不同项目修改的参数，绿色字体表示行首用#标记的注释信息；(B)主程序(call_SVs_procedure.py)；(C)主程序生成的shell子脚本，运行过程以及子脚本运行结果。以step12 (KEGG富集分析)为例，用红色框标记。

MetaSVs的主程序

主程序是用Python 3编写的一个Python脚本(名为call_SVs_ procedure.py) (图2B)，它使用基于配置的方式传递参数，并在Linux下作为命令行脚本代码运行。这个Python脚本只需要一个名为config.ini的文件作为输入，上面详细描述了这个参数(参见图2A)。该程序的核心任务包括13个步骤(表1)，解决具体任务:创建软链接(step0)，质量控制和序列统计(step1)，去除宿主reads (step2)，宏基因组组装和评估(step3)，提取高质量的基因组bins和去复制(step4和5)，物种分类和bins的基因模型(step6和7)，SVs的检测和可视化(step8-10)，KEGG (Kyoto Encyclopedia of Genes and Genomes)对SVs相关基因的富集分析(step11和12)。

在执行核心任务之前，程序首先读取“config.ini”文件中提供的参数，并对这些参数进行初步检查，特别是需要的软件和Conda环境。然后，自动生成一组shell脚本(表1中以“.sh”结尾)，其中包含执行每个步骤所需的所有命令，并依次执行。在每个shell脚本的执行过程中，有关每个样本的命令被写入不同的子脚本，从而可以同时并行处理多个样本，其中并行处理的样本数量由配置文件的“multiprocessing”参数决定。

表1. MetaSVs主程序中每个步骤及其功能描述

同时，程序能够有效处理数据还具备了两个关键特征:“循环执行”和“断点运行”。循环执行意味着，当由于各种损坏或错误(如文件、网络或磁盘问题)而导致程序失败时，程序执行不是简单地退出，而是使用相同的脚本代码进入下一次执行。在循环执行三次之后，如果失败仍然存在，程序将退出。断点运行意味着程序能够从最近的检查点恢复，而不是在中断时覆盖原本可用的中间文件，从而避免了从进程开始重新启动的需要。

MetaSVs的依赖脚本

MetaSVs的结果输出

为了方便用户，我们提供了一个来自人类肠道的示例数据集，它可以在“test”文件夹中找到，包括其源代码、运行过程和脚本的运行结果 (图2C)。为了支持过程可追溯性，首先生成一个“shell/”目录来保存脚本文件，然后将脚本执行的具体过程存储在当前文件夹的“run_*.sh”子文件夹中的“step*.sh.*e*”(存储错误信息)和“step*.sh.*o*”(存储正常运行过程)文件中。以第12步 KEGG富集分析为例，“step12.KEGG_enrichment.sh”脚本的运行过程存放在“run_step12.KEGG_enrichment.sh”文件夹中。如果脚本运行正常，则生成“step12.KEGG_enrichment.1.sh.Check”文件。如果运行失败，将在循环执行三次后退出。此外，如果用户希望从特定步骤重新启动程序，他们可以删除“shell/”路径下该步骤运行的整个“run_step*”文件夹，然后重新运行主程序，主程序将从该步骤开始，而不是从流程的开始。从step0到step12的结果输出文件分别存储在以下目录中:“rawdata_10ge/”、“qc/”、“cleandata/”、“assembly/”、“binning/”、“drep_bins/”、“gene_model/”、“taxonomy/”、“SV/”、“KEGG_enrichment/”和“result_stat/”。

图3. 人类肠道微生物组的结构变异(SVs)

(A) Anaerobutyricum hallii 的T1和T2组中SVs数量(包括插入、缺失、重复和易位)。Wilcoxon检验，*** p < 0.05；(B)SVs在A. hallii的参考MAG基因组上的分布。灰色圆圈为参考MAG, T1组在圈内，T2组在圈外；发生SVs的基因进行KEGG功能富集分析的结果，包括映射到每个功能通路的基因数量(C)和对应的ko ID (D)。p值来自Fisher检验。KEGG, Kyoto Encyclopedia of Genes and Genomes; MAG, metagenome‐assembled genome。

案例研究及结果

这我们提供了一个人类肠道微生物组的简化示例，展示本方法的参数配置、程序运行以及结果输出过程，进而证明MetaSVs的可行性。在这个例子中，我们对5名健康个体的10份粪便样本(根据采样时间分为两组:第一次采样的T1组和10天后第二次采样的T2组)进行测序，得到的Illumina和ONT测序数据用于分析微生物基因组中发生在基因上的SVs，并进一步比较了SVs的组间差异代谢物。所有步骤都在默认的参数设置下运行，并记录每个步骤的总运行时间(表S1)。整个示例项目运行了在6天18小时46分钟，其中大部分时间主要花费在了基因组组装和bins分箱这两步。首先我们评估了原始测序数据的质量，其中Q20和Q30 (phred值为20或30的碱基的百分比，代表高质量读取的百分比)作为Illumina reads的质量指标(表S2)，序列长度分布作为ONT reads的质量指标(图S2)。同时，质控和去除宿主数据后的序列统计也在表S2中给出。然后，去除宿主后的reads被用于混合组装策略(默认软件MetaSpades)进行宏基因组组装。每个样本使用平均1 × 10⁶个ONT reads(超过90%的reads长度为3000 bp)和6 × 10⁷个Illumina 150 bp的对端reads，我们的流程组装结果平均包含9 × 10⁴个contigs，总计3.3 × 10⁸ bp，其中最大contigs达到2192261 bp, N50值为35,646 bp(表S3)。随后，我们将混合组装得到的contigs分成代表单个细菌物种的MAGs，总共获得447个MAGs(每个样品38-56个MAGs)，其中MAGs的序列完整性>70%，污染度<10%。去除重复的MAGs (即属于同一菌种)后，共获得158种细菌，其中20种细菌存在于5个以上的样本中(表S4)。对于这20个物种，我们选择在每个物种中得分最高的MAG作为参考(详细信息列于表2)，然后检测SVs。以A. hallii为例，我们的流程共鉴定出243个插入、121个缺失、16个重复和197个易位(图3A)。在A. hallii中，所有SVs类型在T1组和T2组之间差异显著(Wilcoxon检验和T检验，p < 0.05)。此外，我们检测并比较了T1组和T2组中在参考MAG基因上SVs的分布，发现SVs之间存在显著差异，特别是插入和易位(图3B和表S5)。我们进一步基于KEGG数据库对发生SVs的基因进行了功能富集分析。结果显示，共有64条通路富集(表S6)，其中一条(“Starch and sucrose metabolism”通路)在T2组显著富集(图3C,D)。这表明SVs通过影响基因的完整性，进一步导致菌株水平代谢功能的差异。SVs的功能后果导致基因功能和关键代谢物与细菌之间的独特关系的潜在破坏，这最终转化为细菌-宿主表型关联的差异。

表2. 示例项目中基于细菌宏基因组组装得到的高质量基因组(MAGs)的详细信息

方  法

我们的流程是从Illumina原始测序reads和ONT reads开始(在用Guppy进行based-calling期间，ONT reads已经进行了质控)，然后是Illumina测序数据的质控，去除宿主reads，宏基因组组装，bins分箱，bins去复制，基因预测，物种注释，SVs检测以及可视化，最后是进行KEGG富集分析(图1)。

首先，使用MetaWRAP软件的“read_qc”模块对Illumina原始序列进行质量控制。然后是分别用bowtie2和Minimap2对宿主基因组进行定位，从而识别Illumina和ONT数据中的宿主reads，并去除识别得到的宿主reads。使用SeqKit统计质控以及去除宿主序列后的reads数。对于基因组组装，我们提供了五种技术供用户选择，包括只用ONT数据进行组装的Canu和Flye，只用Illumina数据进行组装的Spades，以及使用OPERA - MS和MetaSpades的混合组装，其中基因组组装的质量通过Quast进行评估。同时，我们比较了这5种组装的组装效率、基于这些组装得到的MAGs数量和MAGs质量，以及从MAGs中识别的SVs数目。结果表明，基于MetaSpades的混合组装实现了高的SVs发现(图S3)， 而且获得了更高完整性、更低污染率的MAGs (表S7)。因此，MetaSpades软件是我们流程基因组组装策略的默认选项。

然后使用MetaWRAP软件将组装得到的contigs打包成MAGs。在bins分箱期间，细化和重组的默认阈值是完整性>70%和污染率<10%。然后使用dRep对所有的MAGs进行合并以及去重复。对于去重后的MAGs，物种注释信息通过gtdbtk获得，基因信息由Prokka预测。每个物种根据以下公式计算得到Score值，选出Score值最高的MAG用作后续SVs检测的参考：

Score = 1*Completeness − 5*Contamination + 0.5 log(N50)

然后使用MUM&- Co软件检测在50%样本中存在的物种基因组上的SVs事件(插入、删除、易位和倒位)。值得注意的是，这50%的阈值是为了后续统计分析:如果一个物种是罕见的，并且只出现在非常有限的样本中，那么很可能不会进行任何有意义的统计比较。此外，考虑到所有的SVs事件总是相对的，我们选择在所有样本中相同物种的所有MAGs中完整性最高的MAG作为参考，这不仅使所有工作更容易，而且允许在MAGs比对期间尽可能多的恢复来自同一物种基因组上的所有SVs。为了进一步降低SVs假阳性发现的可能性，不考虑MAGs 的contigs起止点10 bp内的SVs事件。此外，我们使用Minimap2软件将ONT reads重新映射到MAG参考序列，然后使用Integrative Genomics Viewer将映射结果可视化，以手动验证SVs的可靠性。我们随机选择的一组SVs，结果证实了当前方法检测到的SVs中有97%是正确的(图S1)。最终的SVs数量和在基因上的分布分别通过R软件包ggplot2和Circos (www.circos.ca)软件可视化。最后使用diamond将所有MAGs的所有基因序列映射到KEGG Orthology (KO)数据库中，将其id转换为KEGG id。然后用Kobas对制图结果进行注释。利用R包clusterProfiler中的enricher函数，其中将预测到发生SVs的基因作为前景基因，将上述基因预测结果作为背景进行KEGG富集分析。

讨  论

高灵活性和可扩展性

在这项研究中，我们开发了一个整合短(Illumina)和长(ONT)reads的流程MetaSVs，用于宏基因组SVs分析。通过该流程，图1中描述的微生物SVs分析可以很容易地实现。最重要的是，我们通过使用一个对用户友好型的配置文件来增加SVs分析的灵活性和可扩展性。该配置文件提供了流程中用到的关键参数，并为每个参数提供了多个选项，从而可以轻松控制整个流程。根据分析需求和用户偏好，可以灵活修改上述配置文件提供的默认分析参数。此外，鉴于高质量的宏基因组组装是宏基因组SVs分析的关键组成部分，为了适应不同的样本集，我们提供三种软件程序供选择，以获得高质量的基因组。因此，我们的流程具有高度的灵活性和可扩展性，允许用户选择所需的输入数据，分析工具和参数，并可应用于人体肠道和其他器官，水，土壤和其他环境中的各种微生物组。

运行高效性

随着Illumina和ONT/PacBio技术的蓬勃发展，目前在宏基因组研究中产生数百万个测序数据是常见的，生物信息学分析越来越依赖于高效的框架和工作流程来处理大量的序列数据, 这也推动了对宏基因组SVs研究的需求。因此，为了实现高效率，断点续跑和多样本并行等功能是有效处理数据所必需的。特别是，流程可能会由于各种原因而失败，例如网络或磁盘问题、文件损坏或错误；因此，能够从最近的检查点恢复而不是覆盖现有文件，为此类任务提供了所需的附加功能。此外，实现并行对于多样本处理具有重要价值，大大缩短了生物信息学分析的执行时间，提高了执行效率。在这里，我们加入了上面提到的两个关键特性，以支持MetaSVs流程的有效执行。如果中断，我们的流程将从中断的地方继续，从而消除了从脚本开始重新启动的需要。同时，我们提供了在一次运行中并行处理多个样本的能力，但用户需要注意可用资源和运行任务。最后，一个值得注意的优点是，由于资源抢占、任务阻塞或网络问题导致的故障可以通过我们的流程在一定时间间隔后使用执行循环有效地解决。

可追溯性和透明性

鉴于宏基因组SVs的发现涉及一系列复杂的过程，除了上述高效和灵活的特点外，我们还致力于提高可追溯性和透明度。为了确保数据处理步骤的透明性，可以在执行路径下找到每个核心步骤的所有代码、参数配置和函数(如表1所示)。这增加了分析的透明度，并帮助用户重用流程/分析，从而使数据分析过程更容易遵循。出于可追溯性的目的，准确的输入、关键的中间产物、最终输出和相关的外部数据都保存在指定的文件中。特别是，我们保留了脚本执行的具体过程(包括正常和错误信息)，以提供透明的过程描述，从而支持研究结果的可追溯性。然而，即使可以执行流程，其操作的正确性也只能通过额外的统计测试来验证。我们的示例数据集可用于测试流程的正确性和可执行性，用户还可以使用自己的数据进行基准测试。

局限性和前景

虽然MetaSVs为研究宏基因组SVs在功能上的影响提供了一个急需的解决方案，但仍然存在一些局限性。其中一个挑战是推动标准化流程描述语言的努力。由于我们的流程依赖于多个复杂的软件，因此在将所有脚本和依赖项打包到Conda包中时，存在一些环境依赖冲突方面的问题。例如，MetaWRAP软件需要Python 2，而Prokka需要Python 3。目前，我们使用Python 2和Python 3将所有依赖分别打包到两个Conda包中，以解决依赖冲突和版本迭代问题。但是，考虑到用户的安装体验，我们进一步使用Docker软件对Conda环境进行容器化，最终提供一个可以方便直接安装在Linux系统上的metasv容器。流程的另一个挑战是允许调度器更好地分配内存和计算资源。众所周知，组装和bins分箱是高内存任务，目前只能通过并行减少样本数量来缓解。然而，这也意味着时间消耗的增加，因为在某些情况下，减少任务时间消耗和平衡内存负载是不兼容的。除了内存和时间使用之外，未来还需要解决进一步的问题，例如使用专门为宏基因组数据设计的软件实现更好的组装，降低长读测序的成本，以及找出应用程序中长读和短读的最佳比例。我们的流程代表了当前最先进的混合组装方法和可用工具的现状，已经实现了足够高质量的MAGs，并随后生成了SVs的配置文件，可以分析其功能含义。与使用Illumina‐only组装或直接将reads映射到参考数据库的常见方法相比，混合流程提高了宏基因组组装的质量，覆盖了常见的分类群，并大大扩展了SV的检测范围(表S7和图S4)。此外，值得注意的是，ONT reads更容易验证SV结果的可信度。总之，我们的流程通过将ONT reads纳入宏基因组分析，扩大了遗传变异的检测范围，并进一步研究了它们与代谢组的复杂关系，这为宏基因组研究和深入的功能分析开辟了新的途径和机会。

（▼ 点击跳转）

“iMeta” 是由威立、肠菌分会和本领域数百千华人科学家合作出版的开放获取期刊，主编由中科院微生物所刘双江研究员和荷兰格罗宁根大学傅静远教授担任。目的是发表所有领域高影响力的研究、方法和综述，重点关注微生物组、生物信息、大数据和多组学等。目标是发表前10%(IF > 20)的高影响力论文。期刊特色包括视频投稿、可重复分析、图片打磨、青年编委、前3年免出版费、50万用户的社交媒体宣传等。2022年2月正式创刊发行！发行后相继被Google Scholar、ESCI、PubMed、DOAJ、Scopus等数据库收录！2024年6月获得首个影响因子23.7，位列全球SCI期刊前千分之五(107/21848)，微生物学科2/161，仅低于Nature Reviews，同学科研究类期刊全球第一，中国大陆11/514！

“iMetaOmics” 是“iMeta” 子刊，主编由中国科学院北京生命科学研究院赵方庆研究员和香港中文大学于君教授担任，是定位IF>10的高水平综合期刊，欢迎投稿！