李相辰 柳正衛(wèi) 陸燁瑋 朱業(yè)蕾 張明五 蔣錦琴 彭小軍 王煒欣 高俊順 王曉萌

據(jù)世界衛(wèi)生組織估計，2020年全世界約有990萬例新發(fā)結(jié)核病患者，因結(jié)核病死亡的數(shù)量更是高達150萬例[1]。為了應(yīng)對嚴峻的結(jié)核病疫情，世界衛(wèi)生組織發(fā)布了“終止結(jié)核病戰(zhàn)略(2016—2035年)”的目標，聯(lián)合國召開了結(jié)核病高級別會議，形成了旨在加強終止結(jié)核病的行動和投資的結(jié)核病政治宣言。但目前距離該目標仍有很大差距，亟需研發(fā)新的技術(shù)和方法[2]。

自1998年首個結(jié)核分枝桿菌(Mycobacteriumtuberculosis, MTB)標準菌株(H37Rv)的基因組完成圖公布至今[3]，全基因組測序(whole genome sequencing，WGS)技術(shù)經(jīng)歷了從第一代的Sanger測序到現(xiàn)在的第二代高通量短讀長測序和第三代的單分子長片段測序的快速發(fā)展[4]。隨著高通量測序技術(shù)的成熟和測序成本的降低，WGS技術(shù)已被廣泛應(yīng)用于MTB的研究中，產(chǎn)生了海量的基因組數(shù)據(jù)，但對數(shù)據(jù)的深入分析面臨技術(shù)挑戰(zhàn)[5]。生物信息學作為一門將生物學、計算機科學及統(tǒng)計學結(jié)合起來的交叉學科，在生物數(shù)據(jù)的獲取、管理、分析和解釋方面都具有巨大優(yōu)勢[6]。為此，筆者對MTB基因組的生物信息學分析方法和應(yīng)用進行綜述，為同領(lǐng)域研究者更方便、更靈活的開展數(shù)據(jù)分析和快速選擇研究分析工具提供參考。

一、MTB基因組的基本特征

H37Rv標準株的基因組全長約400萬個堿基，包含3906個蛋白質(zhì)編碼基因，可編碼參與脂質(zhì)代謝的各種酶類，以及2個具有重復結(jié)構(gòu)的富含甘氨酸的蛋白質(zhì)家族PE和PPE，后兩者是MTB與其他細菌的區(qū)別之處[3]。人感染的MTB具有高度的克隆性，根據(jù)單核苷酸多態(tài)性(single nucleotide polymorphism，SNP)的差異和缺失可以將感染人類的結(jié)核分枝桿菌復合群分為7個主要的系統(tǒng)發(fā)育譜系，即第1至第7譜系[7]。MTB在人體內(nèi)的基因組突變速率大約為0.04～2.2突變·基因組-1·年-1，不同譜系間有明顯的差異[8]。由于MTB基因組的單克隆性，并且不同菌株間很難發(fā)生重組或者基因水平轉(zhuǎn)移，因此，其主要通過核心基因或啟動子的自發(fā)變異獲得耐藥性[9]。目前已證實的耐藥靶基因有rpoB(利福平)、katG和inhA(異煙肼)、rpsl和rrs(鏈霉素)、embB(乙胺丁醇)、gyrA和gyrB(氟喹諾酮類)和pncA(吡嗪酰胺)等[10]。

二、WGS在MTB基因組研究中的基本步驟

MTB的基因組研究主要分為以下5個基本步驟，主要涉及MTB樣品的制備、WGS數(shù)據(jù)產(chǎn)出、數(shù)據(jù)質(zhì)控與預(yù)處理、變異檢測，以及數(shù)據(jù)分析和可視化等內(nèi)容。

1.MTB樣品的制備：從痰液樣本培養(yǎng)物中提取MTB的DNA。

2.WGS數(shù)據(jù)產(chǎn)出：抽提樣品的DNA后，通過構(gòu)建測序文庫和進行高通量測序來獲得WGS數(shù)據(jù)。

3.數(shù)據(jù)質(zhì)控與預(yù)處理：對測序所得原始數(shù)據(jù)(raw data)進行質(zhì)量控制，重點關(guān)注測序數(shù)據(jù)的總測序數(shù)據(jù)量、高質(zhì)量測序數(shù)據(jù)比例(Q30)和GC含量等指標，以滿足下游分析的要求。數(shù)據(jù)預(yù)處理包括去除在測序和建庫過程中人為添加的引物、接頭，以及測序過程中產(chǎn)生的低質(zhì)量序列等。建議采用比對人類和其他微生物基因組的方式去除可能的宿主和非MTB序列，再將獲得的純凈序列(clean data)與參考基因組進行比對[5]，并主要使用比對率(測序數(shù)據(jù)中成功比對到參考基因組的比例)、覆蓋率(參考基因組被成功比對的比例)以及平均測序深度這三個指標對比對結(jié)果進行質(zhì)控。

4.變異檢測：基于比對結(jié)果進行SNP、插入/缺失(insertion-deletion，indel)和結(jié)構(gòu)變異(structure variation，SV)等基因組變異的檢測，并基于參考基因組對變異進行注釋。PE/PPE基因家族、其他重復基因和可移動遺傳元件等區(qū)域的變異檢測錯誤率較高，通常在后續(xù)分析中被排除[5]。

5.數(shù)據(jù)分析和可視化：基于基因組變異信息，可以進一步鑒定MTB的譜系或亞種、預(yù)測菌株的耐藥性、監(jiān)測MTB的傳播等。并可以選擇合適的圖形將數(shù)據(jù)可視化，提高結(jié)果的可讀性，有利于生物學規(guī)律的觀察和總結(jié)。

三、WGS數(shù)據(jù)分析常用的環(huán)境

WGS數(shù)據(jù)分析需要在專門的軟件環(huán)境下開展，熟悉常用的編程語言能夠幫助研究者更好地利用現(xiàn)有工具分析數(shù)據(jù)。目前，本領(lǐng)域的分析工具主要集中在Shell和Python這兩種語言環(huán)境下運行。這兩類語言環(huán)境下有很多可利用的生物信息學軟件，研究者只需要通過極少的代碼串聯(lián)現(xiàn)有的工具就可以實現(xiàn)數(shù)據(jù)分析的自動化。對于高通量測序數(shù)據(jù)的處理則需要使用高性能的服務(wù)器，Linux是其最常用的服務(wù)器操作系統(tǒng)。

Shell語言是Linux操作系統(tǒng)的命令和程序設(shè)計語言，幾乎所有的生物信息學分析工具都可以在Linux服務(wù)器的Shell環(huán)境下運行，而在其他系統(tǒng)環(huán)境中搭建分析流程則非常困難。如果研究者的電腦運行的是Windows操作系統(tǒng)，則需要安裝遠程訪問Linux服務(wù)器的軟件，如Xshell或PuTTY等。如果是Mac OS系統(tǒng)，研究者就需要使用系統(tǒng)自帶的Terminal程序?qū)崿F(xiàn)遠程訪問Linux服務(wù)器。

Snakemake是基于Python的一款流程搭建工具，繼承了Python語言簡單易讀、邏輯清晰、便于維護的特點，同時還支持Python語法，非常適合新的使用者[11]。Snakemake的基本組成單位叫“規(guī)則”，即rule；每個rule里面又有多個元素(input、output、run等)。它的執(zhí)行邏輯就是將各個rule利用input/output 連接起來，形成一個完整的工作流，即當檢測到input，就執(zhí)行相應(yīng)rule；檢測到output，就跳過相應(yīng)rule，根據(jù)這一規(guī)則，Snakemake還可以實現(xiàn)斷點續(xù)投。結(jié)合Conda軟件包管理工具，Snakemake可以輕松解決各種軟件安裝的依賴問題。Visual Studio Code是一款免費跨平臺的代碼編輯軟件，支持使用SSH連接Linux服務(wù)器進行遠程開發(fā)，保持開發(fā)與分析工作環(huán)境的一致性。

四、結(jié)核病WGS研究常用的分析軟件

近年來，隨著高通量測序技術(shù)的成熟和應(yīng)用，結(jié)核病WGS研究領(lǐng)域的相關(guān)分析方法和工具也取得了快速發(fā)展，大量優(yōu)秀的軟件、流程、在線分析平臺相繼發(fā)布，對推動本領(lǐng)域的研究做出了貢獻。

(一)數(shù)據(jù)處理和變異檢測

原始數(shù)據(jù)需要進行數(shù)據(jù)質(zhì)量過濾，包括過濾測序接頭、低質(zhì)量序列、低復雜度序列、重復序列等，常用的質(zhì)控和過濾軟件有fastp[12]和Trimmomatic[13]等。原始測序數(shù)據(jù)經(jīng)低質(zhì)量序列過濾后，可用Kraken軟件去除來源于人和非MTB物種的序列[14]。測序數(shù)據(jù)經(jīng)過清理，下一步是將序列定位或比對到參考基因組上，序列比對常用BWA[15]和Bowtie[16]等工具，輸出的標準定位文件格式為SAM/BAM?？墒褂肧AMtools[17]和Picard[18]軟件來處理和分析SAM/BAM文件。常用的基因組變異檢測工具有SAMtools/BCFtools[17]、GATK[18]和freebayes[19]等軟件。檢測到的所有變異結(jié)果存儲在VCF格式文件中，需要進一步結(jié)合質(zhì)量值、測序深度、重復性等參數(shù)進行過濾，最終得到可信度高的變異數(shù)據(jù)集。此外，還可以整合多種工具進行變異檢測，保留具有高度一致性的變異結(jié)果以進一步提高可信度。為了從檢測到的變異中獲得生物學功能等方面的信息，可使用SnpEff軟件進行變異注釋[20]。最后可以基于參考基因組通過SAMtools構(gòu)建多個菌株全部變異的一致性序列，用于后續(xù)的遺傳距離計算和系統(tǒng)發(fā)育樹構(gòu)建。

(二)MTB譜系鑒定和耐藥性檢測

相比于標準基因分型技術(shù)，WGS具有更高的鑒別能力，可以根據(jù)SNP的差異和缺失來精準識別MTB菌株的譜系/亞型[21]。同時，WGS可以在全基因組水平上檢測所有已知耐藥基因的變化信息，其效能已獲得世界衛(wèi)生組織的肯定[22]。國內(nèi)外研究者開發(fā)了幾款自動化分析工具，只需導入原始測序數(shù)據(jù)即可獲得菌株的基因組變異檢測、譜系鑒定和耐藥性預(yù)測結(jié)果。本文將重點介紹以下3款近期發(fā)表并被廣泛引用的軟件平臺。

1.TB-Profiler分析軟件：該軟件由倫敦衛(wèi)生與熱帶醫(yī)學院的Taane G. Clark教授團隊在2015年發(fā)布[23]，同時提供了網(wǎng)頁版在線工具以及開源的可本地化的命令行版本，可通過Conda軟件包管理器快速安裝。此外，研究者可根據(jù)需求個性化地修改TB-Profiler使用的突變數(shù)據(jù)庫，使之納入新發(fā)現(xiàn)的耐藥突變來提高耐藥檢測的準確性。最新版本的TB-Profiler還進一步支持了Oxford Nanopore MinION三代測序平臺產(chǎn)生的長片段序列的分析。

2.Mykrobe分析軟件：同樣在2015年，歐洲生物信息中心Zamin Iqbal教授團隊發(fā)布了基于Kmer算法的MTB分析軟件Mykrobe[24]，提供了Windows和Mac OS系統(tǒng)的安裝版本，可輕松部署在PC和筆記本電腦上。該軟件同樣免費開源并且自帶圖形化操作界面，軟件分析速度快且易用性強，但下游分析功能略少。

3.SAM-TB分析平臺：該平臺是由復旦大學基礎(chǔ)醫(yī)學院高謙教授團隊與深圳市慢性病防治中心合作建立的一個MTB綜合數(shù)據(jù)分析平臺[25]，具有易于訪問、界面友好、操作簡單、功能豐富等優(yōu)點。該平臺在MTB譜系鑒定和耐藥性預(yù)測的基礎(chǔ)上，還提供了系統(tǒng)發(fā)育樹重建、菌株間SNP距離計算和非結(jié)核分枝桿菌混合感染鑒定等功能。SAM-TB測序數(shù)據(jù)分析平臺的建立為我國WGS技術(shù)在結(jié)核病耐藥和傳播監(jiān)測網(wǎng)絡(luò)上的建設(shè)提供了重要基礎(chǔ)[26]。

上述工具對耐藥性的檢測采用的是直接關(guān)聯(lián)法，即通過判斷是否存在數(shù)據(jù)庫中的已知耐藥相關(guān)變異來判斷是否耐藥。雖然其對一線抗結(jié)核藥物有很好的預(yù)測效果，但對預(yù)測二線抗結(jié)核藥物則不太理想。近年來，一些基于WGS數(shù)據(jù)的機器學習類耐藥預(yù)測方法被證明能夠快速且準確地預(yù)測MTB的耐藥性，同時能夠發(fā)現(xiàn)新的耐藥位點并有助于解釋耐藥機制[27-29]。如GenTB是哈佛醫(yī)學院Maha R. Farhat教授團隊2021年發(fā)布的一種基于神經(jīng)網(wǎng)絡(luò)的結(jié)核病耐藥在線預(yù)測工具，相較于TB-Profiler和Mykrobe軟件在一線和二線抗結(jié)核藥物耐藥性預(yù)測效果，其基準測試的結(jié)果均有所提升[30]。

(三)菌株傳播檢測

高通量測序技術(shù)的發(fā)展使得快速監(jiān)測MTB傳播成為可能。WGS技術(shù)可通過檢測菌株間SNP差異并結(jié)合分子進化算法鑒定其傳播方向和傳播鏈，識別傳染源和傳播缺失環(huán)節(jié)[31]。鑒于MTB的遺傳多樣性非常低，通常使用5或12個 SNP的差異閾值來表明流行病學聯(lián)系[31]。除此之外，研究人員近期還陸續(xù)開發(fā)了一些方法來改進WGS技術(shù)對MTB傳播的探測效果。PANPASCO軟件是一種基于線性泛基因組圖譜的遺傳距離計算方法，能夠有效減少不同譜系菌株測序數(shù)據(jù)比對的損失率，提高SNP檢測的分辨率，在多個數(shù)據(jù)集測試中表現(xiàn)出比傳統(tǒng)方法更好的傳播探測效果，具有較好的普適性[32]。PANPASCO也是基于Snakemake軟件的開發(fā)，適用于大規(guī)模樣本的傳播檢測。除了基于單一SNP差異閾值的菌株分型之外，Transcluster軟件是一種通過推測新的傳播事件來識別近期傳播簇的方法，其綜合考慮了菌株的傳播速率、可能發(fā)生傳播的病例采樣時間的間隔和基因組間SNP的差異數(shù)，用以估計菌株間傳播事件發(fā)生的概率和次數(shù)，以此判斷是否具有流行病學聯(lián)系[33]。

這些基于WGS的方法已被證明比接觸追蹤表現(xiàn)更好，并且較經(jīng)典分型方法(例如可變數(shù)目串聯(lián)重復序列分型)具有更高的分辨率[34]。在準確識別近期傳播簇的基礎(chǔ)上，可以進一步結(jié)合分子進化方法推測其內(nèi)部的傳播網(wǎng)絡(luò)(傳播鏈)，常用的軟件有SeqTrack[35]和TransPhylo[36]等。SeqTrack是最早的使用整體傳播樹的構(gòu)建對研究的樣本群體進行傳播網(wǎng)絡(luò)推斷的工具之一，TransPhylo則是在此基礎(chǔ)上加入了對流調(diào)信息的分析，綜合考慮菌株在宿主體內(nèi)的進化情況，從而對傳播網(wǎng)絡(luò)進行推斷。因此，TransPhylo對樣本數(shù)據(jù)中的流調(diào)信息具有更高的要求[37]。傳播網(wǎng)絡(luò)可以通過Cytoscape[38]和igraph[39]等軟件進行可視化，并結(jié)合病例之間時空分布和接觸情況進一步分析傳播順序和傳播源。

(四)菌株進化分析

國內(nèi)外已有較多研究運用系統(tǒng)發(fā)育理論并結(jié)合復雜的進化模型與方法，從MTB的遺傳序列中提取流行病學信息，進而重建結(jié)核病流行過程中病原體時間、空間甚至表型范圍上的進化過程[40-41]。系統(tǒng)發(fā)育樹是進化研究的核心，主流建樹軟件眾多，其中MEGA屬于圖形化軟件，因界面友好而被廣泛使用[42]，方法包括距離法、最大簡約法、最大似然法和貝葉斯法，其中距離法又包括最少進化法和鄰接法。由于鄰接法建樹極快，通常用于建樹嘗試階段，而正式建樹常選用可靠性高的最大似然法。其他常用的進化樹構(gòu)建軟件還有RAxML[43]、IQ-TREE[44]和FastTree[45]。這三款軟件都是基于最大似然法進行系統(tǒng)發(fā)育樹的構(gòu)建，RAxML和IQ-TREE可以構(gòu)建出更優(yōu)似然值的系統(tǒng)發(fā)育樹，但是需要消耗更多的計算資源和時間，而FastTree則可以更加快速地完成系統(tǒng)發(fā)育樹構(gòu)建，但性能與穩(wěn)定性不如前者[46]。

近年來，隨著新發(fā)和再發(fā)傳染病事件的上升趨勢，一種新型的帶有時間戳的貝葉斯進化樹正在興起，其節(jié)點和分支帶有病原體可能被引入當?shù)貍鞑サ臅r間，有助于在結(jié)核病暴發(fā)和流行期間實時進行疫情管理[47]。BEAST是目前最常用的貝葉斯物種分化時間估計軟件之一[48]。通過軟件的圖形界面導入序列、設(shè)置分類群、序列收集日期、核苷酸替代模型、分子鐘類型、樹先驗?zāi)Ｐ筒⒄{(diào)整參數(shù)的權(quán)重，結(jié)合馬爾科夫鏈蒙特卡羅算法采樣，收斂后得到高可靠性的帶分歧時間的群體進化樹以及分子鐘速率的估計。最后可利用Evolview[49]或iTOL[50]網(wǎng)站在線進行進化樹的可視化和美化。

(五)其他分析工具

一些其他領(lǐng)域的分析方法在MTB基因組學中的研究也得到了推廣和應(yīng)用，如全基因組關(guān)聯(lián)分析(genome-wide association study，GWAS)在人類疾病相關(guān)基因的鑒定中發(fā)揮了巨大作用[51]。由于已知的耐藥突變位點不能解釋所有耐藥表型，近期GWAS分析也應(yīng)用于MTB研究中，用于大規(guī)模探索SNP與表型之間的關(guān)系[52]。事實上，關(guān)聯(lián)研究可以使用各種遺傳數(shù)據(jù)類型，包括SNP、indel和SV等，以及不同的表型，如菌株毒性[53]和傳播性[54]等。近年來陸續(xù)有眾多的細菌GWAS分析工具公布，如基于回歸算法的pyseer[55]和基于收斂算法的hogwash[56]等?；谑諗克惴ǖ腉WAS分析對于小樣本數(shù)據(jù)可以發(fā)現(xiàn)更優(yōu)的顯著性結(jié)果，但是對于克隆群體的效果不佳，同時對于大樣本數(shù)據(jù)需要較多的計算資源[53]。

此外，目前已發(fā)布的用于混合感染的檢測軟件如MixInfect[57]、QuantTB[58]和SplitStrains[59]等均可用于分辨由多菌株混合感染引發(fā)結(jié)核病的情況。其中，MixInfect是在SNP檢測結(jié)果的基礎(chǔ)上使用貝葉斯模型進行混合感染的分析[57]；QuantTB是將待檢測樣本的測序數(shù)據(jù)與已經(jīng)搭建好的混合感染代表菌株數(shù)據(jù)庫進行比對來評估樣本的混合感染情況，評估結(jié)果的準確性高度依賴數(shù)據(jù)庫中的數(shù)據(jù)[58]；而SplitStrains則是使用更復雜的期望最大化算法來分析樣本的混合感染情況，對低深度的測序數(shù)據(jù)以及遺傳距離較近的菌株混合感染具有更優(yōu)的檢測性能[59]。

五、測序數(shù)據(jù)和分析代碼的分享

通常WGS測序的原始數(shù)據(jù)應(yīng)在文章發(fā)表時上傳至NCBI、EBI和DDBJ等國際數(shù)據(jù)中心。中國科學院建立的組學原始數(shù)據(jù)歸檔庫(Genome Sequence Archive, GSA)是國內(nèi)首個被國際期刊認可的組學數(shù)據(jù)發(fā)布平臺，填補了我國相關(guān)數(shù)據(jù)庫的空白，極大地便利了國內(nèi)研究者測序數(shù)據(jù)的遞交、管理和分享[60]。

越來越多的結(jié)核病研究文章在發(fā)表的同時，會在Github之類的平臺公開分析代碼和測試數(shù)據(jù)，并在研究者的反饋下不斷優(yōu)化和升級，使后續(xù)相關(guān)分析結(jié)果更加合理。筆者基于Snakemake工具開發(fā)了一套MTB全基因組測序數(shù)據(jù)自動化分析流程——TBSeqPipe(https://github.com/KevinLYW366/TBSeqPipe)。該流程從原始測序數(shù)據(jù)出發(fā)，可對MTB樣本進行譜系鑒定、耐藥性預(yù)測、遺傳距離計算、群體進化分析以及混合感染檢測，并最終生成可視化的分析總結(jié)報告，方便了國內(nèi)外研究者的使用。

六、結(jié)語與展望

高通量測序技術(shù)通量的提高和價格的下降，極大地推動了WGS技術(shù)在結(jié)核病分子流行病學研究中的應(yīng)用[5]。WGS技術(shù)可檢測完整的MTB基因組，既可以迅速獲得全面、詳細、準確的臨床菌株及其耐藥性信息，及時為臨床用藥及個體化治療提供指導，還可以進一步應(yīng)用于MTB微進化和分子流行病學研究，為結(jié)核病精準防控策略提供重要依據(jù)。盡管如此，現(xiàn)在仍缺乏一致的、國際公認的WGS數(shù)據(jù)分析金標準，難以將不同流程間產(chǎn)生的異質(zhì)性很高的檢測結(jié)果進行相互比較[61]。目前，雖然已經(jīng)有一些專門為MTB基因組學分析開發(fā)的算法和軟件，但仍處于發(fā)展的初級階段，還有很多有待改進的方向，如常用耐藥性檢測工具僅限于少數(shù)已知的編碼藥物靶向蛋白質(zhì)基因上的關(guān)鍵突變，對在耐藥機制中研究較少的二線抗結(jié)核藥物的預(yù)測結(jié)果較差[62]。提高MTB基因組數(shù)據(jù)分析的效率，熟練掌握主流分析的基本思路和常用工具是基礎(chǔ)，通過編程來實現(xiàn)分析自動化可以極大地提高工作效率和可重復性。此外，對已發(fā)表的數(shù)據(jù)進行歸納整理、提高可用性以及進一步挖掘也十分必要。

利益沖突所有作者均聲明不存在利益沖突

作者貢獻李相辰、柳正衛(wèi)、陸燁瑋和朱業(yè)蕾：文獻檢索、資料收集整理和文稿撰寫；張明五、蔣錦琴和彭小軍：資料收集整理、文稿修訂和編輯；王煒欣和高俊順：數(shù)據(jù)分析、方法應(yīng)用、資源支持、文稿修訂和編輯；王曉萌：文章概念提出、項目分析、資金支持、初稿撰寫、文稿修訂、審核和編輯