TransFlow基于全基因组测序结核分枝杆菌近期传播自动化分析流程

　　TransFlow: a Snakemake workflow for transmission analysis of Mycobacterium tuberculosis whole-genome sequencing data
　　Pan J，Li X，Zhang M，Lu Y，Zhu Y，Wu K，Wu Y，Wang W，Chen B，Liu Z，Wang X，Gao J.
　　Bioinformatics，2023，39(1): btac785.
　　doi: 10.1093/bioinformatics/btac785.
　　PMID: 36469333.
　　近日，《生物信息学》（Bioinformatics）杂志发表了浙江省疾病预防控制中心结核病预防控制所、浙江省医学精准检验与监测研究重点实验室和浙江省中医药大学等单位合作开发的基于全基因组测序数据的结核分枝杆菌传播自动化分析流程TransFlow。TransFlow是一个基于Snakemake流程管理工具所开发的自动化分析流程，具有免费开源、安装简单、使用方便、运行高效等优点，为结核病的近期传播研究和科学防控提供了重要的工具。
　　背 景
　　结核病是由结核分枝杆菌（Mycobacterium tuberculosis，MTB）引起的感染性疾病，仍然是我国乃至全球面临的重大公共卫生问题。作为一种传染性疾病，掌握其传播规律是防控的关键。近年来，随着高通量测序技术的发展和测序成本的骤减，全基因组测序（whole-genome sequencing，WGS）越来越多地被运用于研究结核病的传播。基于WGS的分子流行病学分析比IS6110-RFLP和 MIRU-VNTR等具有更高的分辨率，同时可以通过描述核苷酸的替换顺序来追踪菌株的传播方向和传播链，为流行病学调查提供更可靠的线索。
　　但是，对WGS产生的大量数据进行分析和解释是目前该方法应用的主要瓶颈之一，极大地限制了WGS在临床上的广泛推广应用。虽然已经有一些公开发表的生物信息学分析流程和网站，可以在一定程度上实现MTB WGS数据的操作和分析，例如TB-Profiler、Mykrobe、MTBSeq 和 SAM-TB，可以检测耐药突变并鉴定MTB的谱系，有些还可以进行系统发育关系和传播成簇分析。然而，目前对于大规模MTB样本测序数据的传播探测研究，包括但不限于传播成簇和近期传播率分析、传播方向和传播网络构建，传播风险因素推测等，仍然缺乏一套标准化的数据分析流程。此外，近年来结核病和生物信息学研究领域涌现出一些新颖的分析算法和工具。因此，我们基于Snakemake流程管理工具，将MTB研究领域最新发展的先进工具组合成一套免费的、快速的和易用的自动化分析流程——TransFlow。
　　设计流程
　　TransFlow 是一个免费开源的MTB WGS分析流程，可在Windows、Linux和Mac OS等多个操作系统上安装，用户仅需要输入测序所得到的样本双端测序FASTQ数据和包含样本对应流行病学特征数据（如采样时间、患者年龄、性别和家庭住址地理坐标等）的元数据文件即可运行流程进行分析（图1）。整个 TransFlow的框架由以下五个独立且连贯的分析模块组成：
　　（1）质量控制：对于测序所得的原始数据进行质量控制，去除在测序和建库过程中人为添加的引物、接头，以及测序产生的低质量序列等；
　　（2）MTB过滤：采用比对人类和其他微生物基因组的方式去除可能的宿主和非MTB的序列；
　　（3）序列比对和突变体识别：将获得的纯净序列与参考基因组进行比对，检测全基因组的单核苷酸多态性（SNPs）；
　　（4）近期传播探测：分析菌株之间SNP数量的差异，评估菌株间的亲缘关系（遗传距离），探测具有近期传播关系的菌株簇，对每个传播簇重新构建其传播关系网络；
　　（5）传播风险因素推测：根据元数据文件中输入的病例流行病学特征数据与菌株成簇结果进行单因素回归分析，推测传播风险因素。
　　TransFlow是高度自动化的，所有模块会按照先后顺序自动运行完成。此外，每个模块都可以独立运行，用户可以通过调整相关参数获得满意的结果，如手动过滤低质量样本，尝试不同的传播检测方法或基因组成簇阈值等。TransFlow 的一个重要且独特的优势是其底层框架可以实现断点续投。如果出现错误，或者需要调整数据和参数，流程会直接重新执行相应的模块，而无需从头开始。TransFlow 是完全开源的，主要以 Python 和 R 两种编程语言实现，使用 Conda软件和环境管理系统来自动按照所需的各种依赖。TransFlow提供了完整的软件安装和用户使用的说明，通用参数设置提供了预设的默认参数和详细的说明，为用户提供参考。此外，随软件一起附带了一个示例数据集，包括 FASTQ 和元数据文件，用于快速体验软件的全部功能。
　　数据质控和过滤
　　TransFlow首先采用FastQC软件检查测序数据的质量，并为每个FASTQ文件生成质控报告，然后使用MultiQC软件将所有结果整合成一个HTLM网页交互式报告（图2）。TransFlow使用Trimmomatic和fastp软件去除在测序和建库过程中人为添加的引物接头以及测序过程中产生的低质量序列等。此外，TransFlow还提供了去除可能的宿主或非MTB的序列污染的功能。
　　基于泛基因组的SNP差异探测
　　为了克服单一参考基因组对不同谱系菌株变异检测可能产生的偏差，TransFlow采用PANPASCO软件进行配对SNP距离计算，基于一个由146个覆盖四种主要谱系（第1到第4谱系）的MTB基因组所组成的泛基因组。TransFlow会过滤掉在参考泛基因组中的PE/PPE基因家族、其他重复基因和可移动遗传元件等变异检测错误率较高的区域的SNP。最终，TransFlow会输出所有菌株之间的SNP差异数量（遗传距离）的对称矩阵，并以一个聚类热图来可视化菌株之间亲缘关系的远近（图 3A）。同时，TransFlow绘制出一个直方图来显示所有菌株之间遗传距离的频率分布，其中突出显示了0到12个SNPs距离的菌株对数量（图 3B）。这些结果可以给用户提供判断测序菌株之间是否存在一定近期传播关系的初步证据。
　　传播探测分析
　　TransFlow可通过分析菌株间的遗传距离进行传播探测分析，重建可能有传播关系成簇菌株间的传播网络，为传播溯源提供线索。TransFlow提供了两种不同的传播聚类方法，分别是基于SNP的方法和基于传播的方法，以推断可能具有近期传播关系的样本。基于SNP的方法的原理是如果两个样本的SNP距离小于或等于一个固定的阈值，则判断它们属于同一个传播簇。相较于基于SNP的方法仅考虑SNP距离，基于传播的方法则进一步考虑了采样日期、分子钟速率和传播过程等先验信息。基于传播的方法是，如果样本对之间以给定概率所估计的传输事件数量低于阈值，则将它们判断为具有近期传播关系。这样获得的传播簇不仅包含了研究群体间直接的传播事件，还包括了它们之间未抽样的隐匿传播事件。TransFlow输出一个 TSV 文件，包含所有样本的成簇结果和传播簇的编号，各传播簇按其成员数量排序。此外，TransFlow 输出两个饼图来分别显示成簇样本占比（图 4A）和所有传播簇成员数量分布（图 4B）的统计结果。接下来，TransFlow进一步使用SeqTrack 算法对包含至少三个样本的传播簇推测样本之间的传播先后顺序，重建传播关系网络（图 4C）。除了样本间的SNP距离和样本采样日期，用户还可以输入样本的地理坐标来表示它们的空间连通性，以改进局部传播关系的推断。
　　传播危险因素推断
　　TransFlow 进一步提供了推断与传播相关的流行病学风险因素的功能，以期为结核病社区精准防控提供参考。用户需要在元数据文件中提供所有需要检测的流行病学特征数据，例如年龄、性别、居住地、既往结核病治疗史以及糖尿病或艾滋病病毒感染状况等。TransFlow使用R包gtsummary对配置文件中指定的流行病学特征进行传播聚类单变量回归分析，自动检测数据集中连续的、分类的和二分的变量，执行适当的描述性统计，还包括每个变量的缺失量。最后，TransFlow生成一个可供发表的统计分析汇总表（图5）。
　　创建总结报告
　　TransFlow 的结果最终显示在一份用户友好的交互式 HTML 报告中。该报告包括统计数据的汇总与全部可视化图表，并附有数据质控、变异检测、传播成簇探测、传播网络重建和传播风险因素推断的详细方法及参数的说明。
　　真实数据集测试结果
　　为了验证TransFlow的分析效果，我们将其应用于来自一项上海地区基于人群的回顾性结核病研究的真实数据集（SRA数据库编号：SRP058221）。该研究共收集了 324 株耐多药结核病患者的 MTB 分离株。作者首先通过 VNTR 基因分型筛选了 125 个样本，然后在其中的122个样本中成功地进行了WGS。此外，我们还从原文作者那里取得了相关的流行病学数据。对于该数据集，我们测试了基于传播的方法，共鉴定出了103个测序菌株分布在36个传播簇中，成簇比例为84%（图4），这与原始研究的结果几乎相同（103株菌和38个传播簇，成簇比例为84%)。例如，图4C显示了传播簇2的重建传播网络，它与原始研究中的集群9相同，并补充了假定的传播轨迹。我们可以进一步地整合原论文中推断的流行病学联系，以复现可能的传播场景。如图4D所示，假定的指示病例是一位丈夫 (12_1614)，他随后将 MTB 传染给了他的妻子 (12_0659)。其后，小区游戏室发生传播事件，感染了其他3例患者（10_0183、10_2010及12_1050）。此外，可以确定另两名与游戏室没有任何流行病学联系的患者（10_1007 和 11_0426）与患者 10_0183 和 10_2010 有关联，也分别处于传播链中。为了确定与传播相关的危险因素，还纳入了177例在原始研究人群中被VNTR基因分型确定为独立的病例。在总共 299 例具有可用流行病学调查结果的病例中，评估了聚集病例和独特病例之间六项流行病学特征（年龄、性别、治疗史、痰涂片结果、治疗结果和北京血统）的差异。与原论文一致，结果表明年龄是耐多药结核病传播的一个假定危险因素（图5），这意味着45岁或以上的患者比其他患者更有可能处于耐多药结核病的传播集群中。
　　意义及展望
　　我们提出了一种新的基于WGS的结核病传播分析流程TransFlow，该流程快速、高效、可定制和易于使用，是研究人员有效和现代化的工具。完整的工作流程从原始读取的质量控制和MTB序列过滤开始。该流程包括数据质控和过滤、序列比对、变异检测、遗传距离计算、传播成簇探测、传播网络重建和传播风险因素推断等步骤，最后生成一个详细的可交互式的网页版总结报告，可为结核病的传播监测及防控研究提供有用的信息。
　　基于全基因组测序的结核病分子流行病学生物信息学技术仍在快速发展中。然而，关于参考基因组的选择和近期传播的SNP差异阈值仍然存在争议。为了克服这些挑战，TransFlow 采用了两种新开发的软件，分别是PANPASCO 和 TransCluster。PANPASCO 利用包含四个主要谱系（第1~4谱系）全部序列的泛基因组和成对距离算法来减少遗传距离计算的偏差。TransCluster是一种新型的传播聚类识别工具，它将采样时间、SNP距离、传播速率和分子钟速率等纳入其传播概率模型，以提高传播成簇分析的识别率和样本的适应性。我们欢迎用户对本流程的各种反馈和建议，并不断改进和更新模块，提高传播分析的可信度，以期推动WGS技术在结核病防控中的应用。