結核分枝桿菌全基因組測序數據分析方法與應用進展

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

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

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

一、MTB基因組的基本特征

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

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

MTB的基因組研究主要分為以下5個基本步驟，主要涉及MTB樣品的制備、WGS數據產出、數據質控與預處理、變異檢測，以及數據分析和可視化等內容。

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

2.WGS數據產出：抽提樣品的DNA后，通過構建測序文庫和進行高通量測序來獲得WGS數據。

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

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

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

三、WGS數據分析常用的環境

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

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

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

四、結核病WGS研究常用的分析軟件

近年來，隨著高通量測序技術的成熟和應用，結核病WGS研究領域的相關分析方法和工具也取得了快速發展，大量優秀的軟件、流程、在線分析平臺相繼發布，對推動本領域的研究做出了貢獻。

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

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

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

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

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

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

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

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

(三)菌株傳播檢測

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

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

(四)菌株進化分析

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

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

(五)其他分析工具

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

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

五、測序數據和分析代碼的分享

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

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

六、結語與展望

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

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

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