结核分枝杆菌全基因组测序数据分析方法与应用进展

一、MTB基因组的基本特征

H37Rv标准株的基因组全长约400万个碱基，包含3906个蛋白质编码基因，可编码参与脂质代谢的各种酶类，以及2个具有重复结构的富含甘氨酸的蛋白质家族PE和PPE，后两者是MTB与其他细菌的区别之处。人感染的MTB具有高度的克隆性，根据单核苷酸多态性（single nucleotide polymorphism, SNP）的差异和缺失可以将感染人类的结核分枝杆菌复合群分为7个主要的系统发育谱系，即第1至第7谱系。MTB在人体内的基因组突变速率大约为0.04~2.2突变·基因组-1·年-1，不同谱系间有明显的差异。由于MTB基因组的单克隆性，并且不同菌株间很难发生重组或者基因水平转移，因此，其主要通过核心基因或启动子的自发变异获得耐药性。目前已证实的耐药靶基因有rpoB（利福平）、katG和inhA（异烟肼）、rpsl和rrs（链霉素）、embB（乙胺丁醇）、gyrA和gyrB（氟喹诺酮类）和pncA（吡嗪酰胺）等。

二、WGS在MTB基因组研究中的基本步骤

MTB的基因组研究主要分为以下5个基本步骤，主要涉及MTB样品的制备、WGS数据产出、数据质控与预处理、变异检测，以及数据分析和可视化等内容。

1.MTB样品的制备：从痰液样本培养物中提取MTB的DNA。

2.WGS数据产出：抽提样品的DNA后，通过构建测序文库和进行高通量测序来获得WGS数据。

3.数据质控与预处理：对测序所得原始数据（raw data）进行质量控制，重点关注测序数据的总测序数据量、高质量测序数据比例（Q30）和GC含量等指标，以满足下游分析的要求。数据预处理包括去除在测序和建库过程中人为添加的引物、接头，以及测序过程中产生的低质量序列等。建议采用比对人类和其他微生物基因组的方式去除可能的宿主和非MTB序列，再将获得的纯净序列（clean data）与参考基因组进行比对，并主要使用比对率（测序数据中成功比对到参考基因组的比例）、覆盖率（参考基因组被成功比对的比例）以及平均测序深度这三个指标对比对结果进行质控。

4.变异检测：基于比对结果进行SNP、插入/缺失（insertion-deletion, indel）和结构变异（structure variation, SV）等基因组变异的检测，并基于参考基因组对变异进行注释。PE/PPE基因家族、其他重复基因和可移动遗传元件等区域的变异检测错误率较高，通常在后续分析中被排除。

5.数据分析和可视化：基于基因组变异信息，可以进一步鉴定MTB的谱系或亚种、预测菌株的耐药性、监测MTB的传播等。并可以选择合适的图形将数据可视化，提高结果的可读性，有利于生物学规律的观察和总结。

三、WGS数据分析常用的环境

WGS数据分析需要在专门的软件环境下开展，熟悉常用的编程语言能够帮助研究者更好地利用现有工具分析数据。目前，本领域的分析工具主要集中在Shell和Python这两种语言环境下运行。这两类语言环境下有很多可利用的生物信息学软件，研究者只需要通过极少的代码串联现有的工具就可以实现数据分析的自动化。对于高通量测序数据的处理则需要使用高性能的服务器，Linux是其常用的服务器操作系统。

Shell语言是Linux操作系统的命令和程序设计语言，几乎所有的生物信息学分析工具都可以在Linux服务器的Shell环境下运行，而在其他系统环境中搭建分析流程则非常困难。如果研究者的电脑运行的是Windows操作系统，则需要安装远程访问Linux服务器的软件，如Xshell或PuTTY等。如果是Mac OS系统，研究者就需要使用系统自带的Terminal程序实现远程访问Linux服务器。

Snakemake是基于Python的一款流程搭建工具，继承了Python语言简单易读、逻辑清晰、便于维护的特点，同时还支持Python语法，非常适合新的使用者。Snakemake的基本组成单位叫“规则”，即rule；每个rule里面又有多个元素（input、output、run等）。它的执行逻辑就是将各个rule利用input/output连接起来，形成一个完整的工作流，即当检测到input，就执行相应rule；检测到output，就跳过相应rule，根据这一规则，Snakemake还可以实现断点续投。结合Conda软件包管理工具，Snakemake可以轻松解决各种软件安装的依赖问题。Visual Studio Code是一款免费跨平台的代码编辑软件，支持使用SSH连接Linux服务器进行远程开发，保持开发与分析工作环境的一致性。

四、结核病WGS研究常用的分析软件

近年来，随着高通量测序技术的成熟和应用，结核病WGS研究领域的相关分析方法和工具也取得了快速发展，大量的软件、流程、在线分析平台相继发布，对推动本领域的研究做出了贡献。

（一）数据处理和变异检测

原始数据需要进行数据质量过滤，包括过滤测序接头、低质量序列、低复杂度序列、重复序列等，常用的质控和过滤软件有fastp和Trimmomatic等。原始测序数据经低质量序列过滤后，可用Kraken软件去除来源于人和非MTB物种的序列。测序数据经过清理，下一步是将序列定位或比对到参考基因组上，序列比对常用BWA和Bowtie等工具，输出的标准定位文件格式为SAM/BAM。可使用SAMtools和Picard软件来处理和分析SAM/BAM文件。常用的基因组变异检测工具有SAMtools/BCFtools、GATK和freebayes等软件。检测到的所有变异结果存储在VCF格式文件中，需要进一步结合质量值、测序深度、重复性等参数进行过滤，得到可信度高的变异数据集。此外，还可以整合多种工具进行变异检测，保留具有高度一致性的变异结果以进一步提高可信度。为了从检测到的变异中获得生物学功能等方面的信息，可使用SnpEff软件进行变异注释。可以基于参考基因组通过SAMtools构建多个菌株全部变异的一致性序列，用于后续的遗传距离计算和系统发育树构建。

（二）MTB谱系鉴定和耐药性检测

相比于标准基因分型技术，WGS具有更高的鉴别能力，可以根据SNP的差异和缺失来识别MTB菌株的谱系/亚型。同时，WGS可以在全基因组水平上检测所有已知耐药基因的变化信息，其效能已获得世界卫生组织的肯定。国内外研究者开发了几款自动化分析工具，只需导入原始测序数据即可获得菌株的基因组变异检测、谱系鉴定和耐药性预测结果。本文将重点介绍以下3款近期发表并被广泛引用的软件平台。

1.TB-Profiler分析软件：该软件由伦敦卫生与热带医学院的Taane G. Clark教授团队在2015年发布，同时提供了网页版在线工具以及开源的可本地化的命令行版本，可通过Conda软件包管理器快速安装。此外，研究者可根据需求个性化地修改TB-Profiler使用的突变数据库，使之纳入新发现的耐药突变来提高耐药检测的准确性。版本的TB-Profiler还进一步支持了Oxford Nanopore MinION三代测序平台产生的长片段序列的分析。

2.Mykrobe分析软件：同样在2015年，欧洲生物信息中心Zamin Iqbal教授团队发布了基于Kmer算法的MTB分析软件Mykrobe，提供了Windows和Mac OS系统的安装版本，可轻松部署在PC和笔记本电脑上。该软件同样免费开源并且自带图形化操作界面，软件分析速度快且易用性强，但下游分析功能略少。

3.SAM-TB分析平台：该平台是由复旦大学基础医学院高谦教授团队与深圳市慢性病防治中心合作建立的一个MTB综合数据分析平台，具有易于访问、界面友好、操作简单、功能丰富等优点。该平台在MTB谱系鉴定和耐药性预测的基础上，还提供了系统发育树重建、菌株间SNP距离计算和非结核分枝杆菌混合感染鉴定等功能。SAM-TB测序数据分析平台的建立为我国WGS技术在结核病耐药和传播监测网络上的建设提供了重要基础。

上述工具对耐药性的检测采用的是直接关联法，即通过判断是否存在数据库中的已知耐药相关变异来判断是否耐药。虽然其对一线抗结核药物有很好的预测效果，但对预测二线抗结核药物则不太理想。近年来，一些基于WGS数据的机器学习类耐药预测方法被证明能够快速且准确地预测MTB的耐药性，同时能够发现新的耐药位点并有助于解释耐药机制。如GenTB是哈佛医学院Maha R. Farhat教授团队2021年发布的一种基于神经网络的结核病耐药在线预测工具，相较于TB-Profiler和Mykrobe软件在一线和二线抗结核药物耐药性预测效果，其基准测试的结果均有所提升。

（三）菌株传播检测

高通量测序技术的发展使得快速监测MTB传播成为可能。WGS技术可通过检测菌株间SNP差异并结合分子进化算法鉴定其传播方向和传播链，识别传染源和传播缺失环节。鉴于MTB的遗传多样性非常低，通常使用5或12个SNP的差异阈值来表明流行病学联系。除此之外，研究人员近期还陆续开发了一些方法来改进WGS技术对MTB传播的探测效果。PANPASCO软件是一种基于线性泛基因组图谱的遗传距离计算方法，能够有效减少不同谱系菌株测序数据比对的损失率，提高SNP检测的分辨率，在多个数据集测试中表现出比传统方法更好的传播探测效果，具有较好的普适性。PANPASCO也是基于Snakemake软件的开发，适用于大规模样本的传播检测。除了基于单一SNP差异阈值的菌株分型之外，Transcluster软件是一种通过推测新的传播事件来识别近期传播簇的方法，其综合考虑了菌株的传播速率、可能发生传播的病例采样时间的间隔和基因组间SNP的差异数，用以估计菌株间传播事件发生的概率和次数，以此判断是否具有流行病学联系。

这些基于WGS的方法已被证明比接触追踪表现更好，并且较经典分型方法（例如可变数目串联重复序列分型）具有更高的分辨率。在准确识别近期传播簇的基础上，可以进一步结合分子进化方法推测其内部的传播网络（传播链），常用的软件有SeqTrack和TransPhylo等。SeqTrack是早的使用整体传播树的构建对研究的样本群体进行传播网络推断的工具之一，TransPhylo则是在此基础上加入了对流调信息的分析，综合考虑菌株在宿主体内的进化情况，从而对传播网络进行推断。因此，TransPhylo对样本数据中的流调信息具有更高的要求。传播网络可以通过Cytoscape和igraph等软件进行可视化，并结合病例之间时空分布和接触情况进一步分析传播顺序和传播源。

（四）菌株进化分析

国内外已有较多研究运用系统发育理论并结合复杂的进化模型与方法，从MTB的遗传序列中提取流行病学信息，进而重建结核病流行过程中病原体时间、空间甚至表型范围上的进化过程。系统发育树是进化研究的核心，主流建树软件众多，其中MEGA属于图形化软件，因界面友好而被广泛使用，方法包括距离法、大简约法、大似然法和贝叶斯法，其中距离法又包括少进化法和邻接法。由于邻接法建树极快，通常用于建树尝试阶段，而正式建树常选用可靠性高的大似然法。其他常用的进化树构建软件还有RAxm[x]l、IQ-TREE和FastTree。这三款软件都是基于大似然法进行系统发育树的构建，RAxm[x]l和IQ-TREE可以构建出更优似然值的系统发育树，但是需要消耗更多的计算资源和时间，而FastTree则可以更加快速地完成系统发育树构建，但性能与稳定性不如前者。

近年来，随着新发和再发传染病事件的上升趋势，一种新型的带有时间戳的贝叶斯进化树正在兴起，其节点和分支带有病原体可能被引入当地传播的时间，有助于在结核病暴发和流行期间实时进行疫情管理。BEAST是目前常用的贝叶斯物种分化时间估计软件之一。通过软件的图形界面导入序列、设置分类群、序列收集日期、核苷酸替代模型、分子钟类型、树先验模型并调整参数的权重，结合马尔科夫链蒙特卡罗算法采样，收敛后得到高可靠性的带分歧时间的群体进化树以及分子钟速率的估计。可利用Evolview或iTOL网站在线进行进化树的可视化和美化。

（五）其他分析工具

一些其他领域的分析方法在MTB基因组学中的研究也得到了推广和应用，如全基因组关联分析（genome-wide association study, GWAS）在人类疾病相关基因的鉴定中发挥了巨大作用。由于已知的耐药突变位点不能解释所有耐药表型，近期GWAS分析也应用于MTB研究中，用于大规模探索SNP与表型之间的关系。事实上，关联研究可以使用各种遗传数据类型，包括SNP、indel和SV等，以及不同的表型，如菌株毒性和传播性等。近年来陆续有众多的细菌GWAS分析工具公布，如基于回归算法的pyseer和基于收敛算法的hogwash等。基于收敛算法的GWAS分析对于小样本数据可以发现更优的结果，但是对于克隆群体的效果不佳，同时对于大样本数据需要较多的计算资源。

此外，目前已发布的用于混合感染的检测软件如MixInfect、QuantTB和SplitStrains等均可用于分辨由多菌株混合感染引发结核病的情况。其中，MixInfect是在SNP检测结果的基础上使用贝叶斯模型进行混合感染的分析；QuantTB是将待检测样本的测序数据与已经搭建好的混合感染代表菌株数据库进行比对来评估样本的混合感染情况，评估结果的准确性高度依赖数据库中的数据；而SplitStrains则是使用更复杂的期望大化算法来分析样本的混合感染情况，对低深度的测序数据以及遗传距离较近的菌株混合感染具有更优的检测性能。

五、测序数据和分析代码的分享

通常WGS测序的原始数据应在文章发表时上传至NCBI、EBI和DDBJ等国际数据中心。中国科学院建立的组学原始数据归档库（Genome Sequence Archive, GSA）是国内被国际期刊认可的组学数据发布平台，填补了我国相关数据库的空白，极大地便利了国内研究者测序数据的递交、管理和分享。

越来越多的结核病研究文章在发表的同时，会在Github之类的平台公开分析代码和测试数据，并在研究者的反馈下不断优化和升级，使后续相关分析结果更加合理。笔者基于Snakemake工具开发了一套MTB全基因组测序数据自动化分析流程——TBSeqPipe（https://github.com/KevinLYW366/TBSeqPipe）。该流程从原始测序数据出发，可对MTB样本进行谱系鉴定、耐药性预测、遗传距离计算、群体进化分析以及混合感染检测，并生成可视化的分析总结报告，方便了国内外研究者的使用。

六、结语与展望

高通量测序技术通量的提高和价格的下降，极大地推动了WGS技术在结核病分子流行病学研究中的应用。WGS技术可检测完整的MTB基因组，既可以迅速获得全面、详细、准确的临床菌株及其耐药性信息，及时为临床用药及个体化治疗提供指导，还可以进一步应用于MTB微进化和分子流行病学研究，为结核病防控策略提供重要依据。尽管如此，现在仍缺乏一致的、国际的WGS数据分析金标准，难以将不同流程间产生的异质性很高的检测结果进行相互比较。目前，虽然已经有一些专门为MTB基因组学分析开发的算法和软件，但仍处于发展的初级阶段，还有很多有待改进的方向，如常用耐药性检测工具仅限于少数已知的编码药物靶向蛋白质基因上的关键突变，对在耐药机制中研究较少的二线抗结核药物的预测结果较差。提高MTB基因组数据分析的效率，熟练掌握主流分析的基本思路和常用工具是基础，通过编程来实现分析自动化可以极大地提高工作效率和可重复性。此外，对已发表的数据进行归纳整理、提高可用性以及进一步挖掘也十分必要。