沈應(yīng)博，史曉敏，沈建忠，汪洋，王少林

?

全基因組測序與生物信息學(xué)分析在細(xì)菌耐藥性研究中的應(yīng)用

沈應(yīng)博*，史曉敏*，沈建忠，汪洋，王少林

中國農(nóng)業(yè)大學(xué) 動物醫(yī)學(xué)院，北京 100193

沈應(yīng)博, 史曉敏, 沈建忠, 等. 全基因組測序與生物信息學(xué)分析在細(xì)菌耐藥性研究中的應(yīng)用.生物工程學(xué)報(bào), 2019, 35(4): 541–557.Shen YB, Shi XM, Shen JZ, et al. Application of whole genome sequencing technology and bioinformatics analysis in antimicrobial resistance researches. Chin J Biotech, 2019, 35(4): 541–557.

隨著耐藥細(xì)菌的大量出現(xiàn)及廣泛傳播，細(xì)菌耐藥性成為全球備受關(guān)注的問題。耐藥細(xì)菌的特征如耐藥基因、毒力因子、質(zhì)粒分型等以及不同菌株間親緣關(guān)系對于細(xì)菌耐藥性流行病學(xué)及分子生物學(xué)的研究有著十分重要的意義。但是傳統(tǒng)的技術(shù)手段如聚合酶鏈?zhǔn)椒磻?yīng) (Polymerase chain reaction，PCR) 和脈沖場凝膠電泳 (Pulsed field gel electrophoresis，PFGE) 等得到的結(jié)果不夠全面且精確度低，對于現(xiàn)有的研究存在很大的局限性。全基因組測序技術(shù) (Whole genome sequencing，WGS) 和生物信息學(xué)分析 (Bioinformatics analysis) 由于能夠快速詳盡地得到耐藥細(xì)菌的特征，也能更加精細(xì)地判斷不同菌株間的進(jìn)化關(guān)系，逐漸成為更加有效的技術(shù)手段，為耐藥性研究提供了有效的幫助。因此，文中系統(tǒng)地介紹全基因組測序分析流程中的各個步驟，主要包括文庫構(gòu)建、平臺測序以及后期數(shù)據(jù)分析三大方面的不同方法和其相應(yīng)的特點(diǎn)，期望相關(guān)研究人員對此能夠有更全面的了解，并得到一定的幫助。

耐藥性，全基因組測序，生物信息學(xué)

二十世紀(jì)末期人們對細(xì)菌耐藥性的忽略及抗生素的大量使用導(dǎo)致耐藥細(xì)菌的廣泛流行，對臨床抗感染治療造成了極大的威脅。而且耐藥細(xì)菌能夠在動物、環(huán)境、食品等環(huán)節(jié)間相互傳播，使細(xì)菌耐藥性逐漸成為全世界醫(yī)學(xué)、政界、媒體等各界廣泛關(guān)注的重大公共衛(wèi)生安全問題[1-3]。細(xì)菌耐藥性問題日趨嚴(yán)重，目前每年約有70萬人死于耐藥菌造成的感染性疾病，英國經(jīng)濟(jì)學(xué)家Jim O’Neill預(yù)測2050年該數(shù)據(jù)可能會上升到每年1 000萬人，并累計(jì)給全世界GDP (Gross domestic product) 造成高達(dá)100萬億美元的損失[4]。隨著多重耐藥菌甚至泛耐藥菌的廣泛流行，尤其是對目前臨床治療細(xì)菌感染的兩類“最后一道防線”藥物——碳青霉烯和多黏菌素耐藥的菌株，嚴(yán)重威脅了動物及人類的健康。其中最重要的兩個基因碳青霉烯耐藥基因NDM-1[5]和多黏菌素耐藥基因[6]均是由質(zhì)粒介導(dǎo)的可轉(zhuǎn)移耐藥基因，可在不同種屬細(xì)菌或不同媒介間相互傳播擴(kuò)散，同時還可與其他耐藥基因共存而成為“超級細(xì)菌”[7-12]，使得治療這些細(xì)菌造成的感染性疾更加困難。面對越來越嚴(yán)峻的細(xì)菌耐藥性形勢，加強(qiáng)耐藥菌株的監(jiān)控以及耐藥菌株特征的分析顯得尤為重要。精確、快速、便捷地獲得耐藥菌株的詳盡信息 (包括耐藥基因、毒力因子、質(zhì)粒分型等特征) 以及菌株間親緣關(guān)系等數(shù)據(jù)對于對抗細(xì)菌耐藥性將會有很大的幫助。隨著科學(xué)技術(shù)的發(fā)展，全基因組測序技術(shù) (Whole genome sequencing，WGS) 或稱為二代測序技術(shù) (Next-generation sequencing，NGS)，同時結(jié)合生物信息學(xué)分析技術(shù) (Bioinformatics analysis) 已逐漸成為科學(xué)界研究細(xì)菌耐藥性的重要技術(shù)手段[13]。

全基因組測序技術(shù)不僅可獲得單一菌落的基因組信息，還可獲得混合基因組的信息 (Metagenomics，宏基因組)，即包括非常規(guī)培養(yǎng)和不可培養(yǎng)的細(xì)菌DNA信息。相較于Sanger測序 (一代測序)，WGS無需針對不同DNA片段或細(xì)菌種屬設(shè)計(jì)特定引物，而是測序獲得隨機(jī)的序列后裝配成相對完整的基因組。然而不同測序平臺對每段序列的讀取長度的差異，導(dǎo)致測序結(jié)果的不同[14-15]。通過WGS不僅可獲得近乎完整的細(xì)菌DNA信息，包括種屬 (Species)、耐藥基因 (Antimicrobial resistance genes)、毒力因子 (Virulence-associated genes)、轉(zhuǎn)移元件 (Mobile elements) 等信息，還可對多個細(xì)菌間的基因組信息進(jìn)行比較，對于耐藥菌株的分子流行病學(xué)和傳播機(jī)制研究至關(guān)重要。

WGS自問世以來因其價格昂貴，一直未被廣泛使用，近年來，各大測序平臺的發(fā)展致使價格持續(xù)下降，該技術(shù)被廣泛運(yùn)用于基礎(chǔ)和臨床細(xì)菌耐藥性的研究中[16-17]。雖然WGS技術(shù)在通量方面有優(yōu)勢，但因其測序讀長有一定的局限性 (小于600 bp)，而無法獲得細(xì)菌完整的DNA信息。然而三代測序技術(shù) (Third-generation sequencing) 則突破了WGS在讀長方面的壁壘，實(shí)現(xiàn)了無需PCR擴(kuò)增對每條DNA分子進(jìn)行單獨(dú)測序，平均讀長可達(dá)10 kb甚至更長，可越過一些二代測序技術(shù)難以測通的重復(fù)序列[18]。對于二代及三代測序技術(shù)而言，從樣品DNA提取到數(shù)據(jù)分析(圖1)，每一個步驟使用的方法不同獲得的數(shù)據(jù)質(zhì)量也大相徑庭，最終都將對結(jié)果造成不同程度的影響。

近年來，隨著WGS技術(shù)在細(xì)菌耐藥性領(lǐng)域的應(yīng)用，科學(xué)家們借助該技術(shù)取得了重大的研究進(jìn)展。例如在質(zhì)粒介導(dǎo)的多黏菌素耐藥基因–[6,19-23]以及碳青霉烯耐藥基因NDM-17[24]等新型耐藥基因的發(fā)現(xiàn)過程中，WGS技術(shù)發(fā)揮著不可替代的作用。Hadziabdic等[25]、Kr?ger等[26]、Porse等[27]利用WGS技術(shù)與生物信息學(xué)分析方法闡明了耐藥基因、耐藥質(zhì)粒的進(jìn)化過程；Holt等[28-29]、Wong等[30]、Duchêne等[31]通過該技術(shù)對耐藥菌或病原菌 (志賀桿菌、多重耐藥沙門菌) 在某個地區(qū)或全球范圍內(nèi)的分子進(jìn)化過程進(jìn)行了詳細(xì)的研究，分析了其起源及進(jìn)化因素，為細(xì)菌耐藥性的控制提出了新的科學(xué)依據(jù)。另一方面，WSG數(shù)據(jù)不僅可以在微觀世界給予我們很多信息，而且可以結(jié)合宏觀數(shù)據(jù)在耐藥性的進(jìn)化以及其相關(guān)風(fēng)險因素的分析中發(fā)揮著重要作用。例如，Shen等[32]、Wedley等[33]將WGS數(shù)據(jù)與宏觀數(shù)據(jù)進(jìn)行聯(lián)合分析，獲得了耐藥菌/基因的產(chǎn)生或傳播的相關(guān)風(fēng)險因素。綜上所述，WGS技術(shù)助力了細(xì)菌耐藥性領(lǐng)域的重大發(fā)展，已逐漸成為推動該研究領(lǐng)域發(fā)展的重要技術(shù)手段，但是由于該技術(shù)存在較高的技術(shù)壁壘，導(dǎo)致其受眾面相對較窄。因此，我們將針對二代全基因組測序技術(shù) (由于三代測序技術(shù)價格頗高暫未廣泛應(yīng)用，所以本文僅作簡單介紹) 從文庫構(gòu)建、平臺測序、數(shù)據(jù)分析 (組裝拼接、基本特征分析、核苷酸多樣性分析等) 三大方面系統(tǒng)地評估不同環(huán)節(jié)中不同方法的特點(diǎn)，為細(xì)菌耐藥性中單一細(xì)菌的研究提供幫助。此外，宏基因組技術(shù)目前在耐藥性研究方面的應(yīng)用相對較少，且與單基因組測序在文庫構(gòu)建、平臺測序兩方面基本相同，僅在數(shù)據(jù)分析方面有所不同，因此本文在此不再單獨(dú)詳細(xì)介紹。

1 文庫構(gòu)建

對于WGS技術(shù)而言，雖然對最終的結(jié)果影響最大的是樣品DNA的提取質(zhì)量[34]，但DNA文庫構(gòu)建的影響也不可忽視。隨著技術(shù)的發(fā)展，二代測序DNA文庫的構(gòu)建方法越來越多，雖各有千秋，但其原理和流程基本一致 (圖2)。

DNA文庫構(gòu)建的核心步驟是DNA片段化，常用的方法有物理破碎 (Physical)、酶切 (Enzymatic) 和化學(xué)打斷 (Chemical)。其中應(yīng)用最廣的是物理破碎 (如超聲破碎法) 或酶切法 (如非特異性核酸內(nèi)切酶和轉(zhuǎn)座酶Tagmentation)[35]。超聲破碎通常使用Covaris儀器 (Covaris，Woburn，MA) 和Bioruptor儀器 (Diagenode，Belgium)。酶切法主要包括使用DNaseⅠ或者Fragmentase或者兩者混合物 (New England Biolabs，Ipswich，MA) 將DNA片段化。雖然兩種方法均可以將DNA片段化并獲得有效的測序結(jié)果，但相較于物理破碎方法，酶切法會產(chǎn)生更多人為的DNA片段插入或缺失。新的酶切法——Illumina’s Nextera tagmentation (Illumina，San Diego，CA) 的出現(xiàn)，不僅減少了人為引入的錯誤，而且大大縮短了DNA文庫構(gòu)建的時間[35]。然而文庫質(zhì)量與DNA的GC含量和片段化方法有著很大關(guān)系，對于GC含量高或低的樣品，使用物理片段化方法效果更好。

圖1 全基因組測序分析簡易流程圖[36]

圖2 DNA文庫構(gòu)建流程圖

目前市面上存在很多商業(yè)化的DNA文庫構(gòu)建試劑盒，隨著試劑盒的快速更新?lián)Q代，其價格下降的同時文庫質(zhì)量在提高，并且初始DNA量的要求也在降低。盡管如此，一個重要的原則是DNA初始量越大，后續(xù)需要擴(kuò)增的循環(huán)越少，得到基因組信息越接近真實(shí)。KAPA Hyper Prep Kit Illumina?platforms (Kapa Biosystems，Wilmington，MA) 試劑盒在DNA初始量足夠的情況下可實(shí)現(xiàn)無PCR擴(kuò)增完成DNA文庫的構(gòu)建。眾多的商業(yè)化試劑盒中均有詳細(xì)的流程，在具體過程中有些許差別，在此我們僅對幾款運(yùn)用廣泛的試劑盒在關(guān)鍵的幾個步驟上進(jìn)行整體的比較(表1)。末端修復(fù)、加A尾、接頭連接以及PCR擴(kuò)增是DNA文庫構(gòu)建中關(guān)鍵的步驟，其中主要的優(yōu)化程序是將末端修復(fù)和加A尾這兩步進(jìn)行合并以縮短總體的時間。然而只有少部分如Truseq?DNA PCR-free和KAPA的兩款試劑盒可以做到PCR-free，這不僅能減少文庫構(gòu)建時間，也能最大程度地保證DNA的真實(shí)性，減少PCR擴(kuò)增帶來的偏差。此外，磁珠純化是各類試劑盒中不可缺少的步驟，該步驟費(fèi)時費(fèi)力，且磁珠的價格昂貴，若能在此步驟上進(jìn)行優(yōu)化將會對DNA文庫的構(gòu)建效率有很大的提高。

表1 DNA文庫構(gòu)建試劑盒的特征比較

a: end repair & A-tailing are performed in the same reaction system; b: size selection after end repair; c: the PCR amplification step is selective.

目前，三代全基因組測序技術(shù)還未得到廣泛的應(yīng)用，其相應(yīng)的文庫試劑盒較少，主要包括PacBio公司和Oxford Nanopore Technology公司推出的一系列試劑盒 (SMRTbell Barcoded Adapter Prep Kit、SMRTbell Damage Repair Kit-SPv3、Rapid Barcoding Kit、Rapid Sequencing Kit、Ligation Sequencing Kit、1D^2 Sequencing Kit等)。PacBio公司推出試劑盒的建庫流程主要包括DNA片段化、末端修復(fù)、接頭連接、文庫片段純化、雜交引物和聚合酶綁定，而Oxford Nanopore Technology公司的試劑盒建庫流程更加簡化。對于兩款快速試劑盒 (Rapid Barcoding Kit/Rapid Sequencing Kit) 而言其工作流程主要包括酶切和接頭連接，對于Ligation Sequencing Kit和1D^2 Sequencing Kit這兩款試劑盒而言，其工作流程主要包括可選片段化、末端制備、接頭連接、磁珠吸附。相較于二代全基因組測序文庫試劑盒，三代測序文庫試劑盒最大的特點(diǎn)是無需PCR擴(kuò)增，文庫制備耗時較短。相較于PacBio公司，Oxford Nanopore Technology公司推出的試劑盒實(shí)現(xiàn)了片段化可選，且無需對文庫進(jìn)行純化，耗時更少，例如使用Rapid Barcoding Kit/Rapid Sequencing Kit 制備測序文庫僅需要10 min。

2 測序平臺

目前市場上的全基因組測序平臺主要由Illumina (San Diego)、ThermoFisher (Waltham)、PacBio (Pacific Biosciences) 和Oxford Nanopore Technology開發(fā)的一系列平臺如Miseq、Hiseq、X Ten、Ion、PacBio RSⅡ、PacBio Sequel以及MinION等，不同平臺的測序技術(shù)大相徑庭，Illumina測序平臺主要通過讀取不同色熒光標(biāo)記的可逆終止核苷酸的圖像得到最終的序列結(jié)果[37]；ThermoFisher測序平臺主要是使用半導(dǎo)體芯片計(jì)算每添加一個核苷酸pH值的改變而推算得到最后的序列結(jié)果[38]；PacBio測序平臺主要通過納米技術(shù)和現(xiàn)代光學(xué)系統(tǒng)對單分子合成中的堿基磷酸基團(tuán)上的熒光信號進(jìn)行識別，并將熒光信號轉(zhuǎn)化為序列結(jié)果[39]；Oxford Nanopore Technology測序平臺主要通過將單分子堿基穿過納米孔蛋白的電流信號轉(zhuǎn)換為序列結(jié)果[40]。在數(shù)據(jù)產(chǎn)出、運(yùn)行時間、序列讀長等指標(biāo)上不同公司的測序平臺甚至同一公司不同系列的平臺均存在較大的差別，可根據(jù)具體的情況選擇最適合的平臺 (表2)。

3 組裝拼接

每一條讀長 (Reads) 的DNA信息將存儲于FASTQ格式的文件中，并通過組裝拼接算法將原始序列 (Raw reads) 拼接成更長的片段 (Contig或者Node)。盡管我們希望拼接得到的序列可以代表完整的基因組信息，但由于二代測序技術(shù)在測序讀長上的劣勢，而細(xì)菌基因組上又存在許多比單一讀長還長的重復(fù)序列無法被測序拼接，因此導(dǎo)致基因組序列拼接后出現(xiàn)多個裂縫 (Gap) 而不完整。三代測序技術(shù)具有讀長上的優(yōu)勢，可以不被重復(fù)序列所限制，所以能夠得到完整的基因組信息。對于太小的質(zhì)粒 (小于10 kb)，Pacbio平臺在文庫構(gòu)建時反而容易將其忽略[39,41]，所以對于含有小質(zhì)粒的菌株需要注意所構(gòu)建的文庫的片段大小。由于序列拼接需要處理龐大的數(shù)據(jù)占用較多的計(jì)算機(jī)資源，因此對計(jì)算機(jī)的硬件要求相對較高，且多數(shù)軟件屬于命令類型 (Command line) 需要通過終端（Terminal）代碼運(yùn)行，不易操作。此外，商業(yè)化的軟件可提供可視化的多功能操作平臺，但并非所有實(shí)驗(yàn)室能夠承受其昂貴的價格。

表2 不同測序平臺的相關(guān)參數(shù)

aManufacture’s data.bPreviously reported data.cThe number of chip/run.dThe number of SMRT cells.eIon 314TMChip v2 or Ion 314TMChip v2 BC.fIon 510Chip.gFor 16 SMRT cells.

對于二代測序短讀長 (Short reads) 的結(jié)果，最常用的拼接方法是基于德布魯因圖 (de Bruijn Graph，DBG) 算法開發(fā)的軟件，而該算法最常見的問題則是當(dāng)不同短序列間有重復(fù)時很難分辨其中錯誤的堿基，這將導(dǎo)致重復(fù)部分的序列在拼接過程中被排除[42]。為解決這個問題，在該算法的基礎(chǔ)上進(jìn)一步將原始短序列分割成更小的序列，即為-mers，隨后被降低到-1 mers (圖3)，并通過歐拉算法 (Eulerian walk)獲得最短的-1 mers的可能路徑，從而將序列拼接起來，減少重復(fù)區(qū)域的錯誤拼接。

Velvet是一款基于DBG算法的軟件，用于重測序 () 數(shù)據(jù)的組裝拼接[43]。該軟件包括velveth和velvetg兩個組件，前者用于-mer的構(gòu)建，后者用于-mer陣列的圖形搭建。VelvetOptimizer是一個Perl腳本由Simon Gladman和Torsten Seemann開發(fā)，用于參數(shù)的自動優(yōu)化 (http://www.vicbioin formatics.com/software.velvetoptimiser.shtml)。通過-mers和DBG算法的合用，Velvet可增加拼接富含重復(fù)序列菌株的可能性。

SPAdes也是一款基于DBG算法開發(fā)的適用于多種測序平臺重測序數(shù)據(jù)的拼接軟件[44]。該算法通過以下4個步驟：1) 組裝圖通過錯誤校準(zhǔn)算法后被構(gòu)建為多片段大小的圖；2) 通過-mer和DBG兩種方法進(jìn)行估計(jì)；3) 構(gòu)建配對的組裝圖被構(gòu)建；4) 得到拼接序列的合集后通過原始序列與該結(jié)果比對校正后得到最終結(jié)果。此外，該軟件還推出了一種新的針對質(zhì)粒拼接的算法[49]。

圖3 DBG算法與k-mers的簡易示意圖

IDBA-UD同樣是基于DBG算法針對短讀長序列拼接而開發(fā)的軟件[45]。該軟件使用相鄰序列的測序深度 (Depth) 來改進(jìn)相關(guān)序列測序深度的閾值，然后通過本地配對序列的組裝來減少重復(fù)序列造成的間隔 (Gap)，通過這樣的方法可以減少在短讀長序列中遺失的-mers中的序列信息。但由于該程序只提供命令類型的版本，用戶界面不夠友好，操作較困難而不推薦使用。

RAY依舊是一款基于DBG算法的程序，但其特殊之處在于依賴歐拉距離算法，這種算法定義某一特殊序列子集為種子 (Seeds) 并通過添加這些種子對拼接的序列進(jìn)行延伸[42]。這樣的延伸過程被啟發(fā)法 (Heuristics) 或命令法 (Commands) 控制，將在種子與序列無交叉重疊時停止。這樣拼接出來的序列長度相對較短，但是錯誤較少。表3列出了幾種常用的序列拼接軟件及其各自的特點(diǎn)。

對于三代測序結(jié)果而言，拼接方法多是基于Overlap Layout Consensus (OLC) 算法或DBG算法[50]。三代測序技術(shù)的超長reads，導(dǎo)致其單堿基隨機(jī)錯誤率較高，所以三代測序的原始reads不推薦直接用于組裝。為解決該問題，三代測序結(jié)果在組裝前需對reads進(jìn)行校正。部分三代測序結(jié)果的組裝程序需要用二代測序結(jié)果對三代測序結(jié)果進(jìn)行清洗，例如PacBioToCA[51]、LSC[52]等，而另一部分程序則可通過三代測序結(jié)果中的短序列對長序列進(jìn)行校正，如HGAP (Hierarchical genome-assembly process) 在小型基因組的校正和組裝中表現(xiàn)良好[46]。HGAP首先對三代測序序列按照一定規(guī)則進(jìn)行有序排序，選出其中較長的序列作為“種子”序列，然后用三代數(shù)據(jù)中的較短序列校正較長序列的錯誤，將校正后的“種子”序列進(jìn)行組裝[46]。

由于MinION測序的發(fā)展晚于Pacbio，所以關(guān)于MinION測序數(shù)據(jù)的組裝軟件多數(shù)是在PacBio組裝軟件的基礎(chǔ)上發(fā)展而來或是同時適用于兩者，例如Canu[47]和unicycler[48]等。

Canu是專為單分子序列拼接而設(shè)計(jì)，由Celera Assembler優(yōu)化而來，通過三代測序數(shù)據(jù)進(jìn)行自身校正拼接[47]。Canu的優(yōu)化源于新的重疊和拼接算法，包括一種基于tf-idf加權(quán)MinHash的適宜重疊策略以及一種避免折疊分叉重復(fù)和單倍型的裝配算法[47]。Canu運(yùn)行包括3個階段：reads校正、reads修剪和單序列構(gòu)建 (Unitig construction)。Unitig指許多短片段交叉重疊后裝配的正確有效的長片段。在所有階段中，第一步均是建立一個短序列重疊數(shù)據(jù)庫 (Read and overlap database)。1) 校正階段主要為估計(jì)正確的reads，生成正確的reads和選出用于校正的最佳短序列重疊群；2) 修剪階段主要是標(biāo)識出輸入序列重疊群中重疊區(qū)域，對其進(jìn)行修剪，獲取基于重疊群延伸開的最長序列；3) 裝配階段首先進(jìn)行錯誤序列的識別，然后構(gòu)建重疊序列的最佳overlap直方圖 (Best overlap graph，BOG)，最終輸出組裝好的序列[47]。

表3 不同序列拼接軟件的特點(diǎn)比較

Unicycler是將二代測序技術(shù)與三代測序技術(shù)的數(shù)據(jù)結(jié)合起來，生成更加精確、完整的基因序列[48]。首先，將短reads進(jìn)行拼接獲得精確、連續(xù)的拼接序列。然后，通過長reads確定短reads的最佳拼接方式[48]。主要分為以下7個步驟：1) SPAdes (v3.6.2 or later) 構(gòu)建de Bruijn圖形程序集，然后通過Unicycler的相關(guān)算法，平衡contigs數(shù)與死端 (Dead ends) 數(shù)，從而獲得最優(yōu)的裝配圖；2) 將測序深度和contigs依據(jù)貪心 (Greedy) 算法確定contigs的多重性；3) 通過SPAdes獲得短reads間的鏈接序列，從而鏈接序列將單拷貝contigs連起來；4) 將長reads與多個單拷貝contigs進(jìn)行比對，可以獲得單拷貝contigs之間的鏈接序列，并用單拷貝contigs對長reads進(jìn)行矯正；5) 確定單拷貝contig之間的鏈接序列的可信度，根據(jù)其可信度高低獲得短contigs的最佳排列方式；6) 基于高質(zhì)量的鏈接序列，將contigs進(jìn)行拼接獲得長contigs；7) 通過短contigs對長序列進(jìn)行比對，以減少堿基錯配、小序列插入或缺失[48]。

4 數(shù)據(jù)分析軟件

通過拼接好的序列可獲得細(xì)菌的許多信息，對于細(xì)菌耐藥性研究而言，以下幾方面的信息是我們所需要的：1) 細(xì)菌種屬 (Species)；2) 細(xì)菌攜帶耐藥基因 (Antimicrobial resistance genes，ARGs)、毒力因子 (Virulence-associated genes，VAGs) 以及插入序列 (Insert sequence，IS) 的情況等；3) 細(xì)菌攜帶質(zhì)粒 (Incomplete types，Inc types) 類型；4) 細(xì)菌的多位點(diǎn)序列分型 (Multi-locus sequence type，MLST) 等。為得到菌株的這些信息，通常是將拼接好的序列或原始數(shù)據(jù)用不同的軟件與特定的數(shù)據(jù)庫進(jìn)行比對，尋找對應(yīng)的基因信息，相應(yīng)的軟件可以按照操作方式劃分為網(wǎng)頁類型 (Web-based tools) 和命令類型 (Command-line)。

網(wǎng)頁類型工具可提供更為直觀的用戶界面和簡便的操作環(huán)境。對于細(xì)菌種屬鑒定來說，KmerFinder是一個很好的工具，如果使用其命令版本處理拼接好的序列 (Contigs) 只需約9 s，而對于原始數(shù)據(jù) (Raw reads) 則需要大約190 s53]。且該工具已在Center for Genomic Epidemiology (CGE，https://cge.cbs.dtu. dk/services/KmerFinder/) 網(wǎng)站上開放使用。另一個使用較多的網(wǎng)頁工具是NCBI (National Center for Biotechnology Information) 提供的比對工具BLAST (Basic Local Alignment Search Tool，https://blast.ncbi. nlm.nih.gov/Blast.cgi)。該工具由NCBI提供，可使用NCBI上所有的數(shù)據(jù)庫，具有較全面的信息，但解釋結(jié)果時需要通過對各種參數(shù)和結(jié)果嚴(yán)格地篩選并結(jié)合相應(yīng)的背景知識去判斷，相對繁瑣。另一款程序Rapid Annotation using Subsystem Technology (RAST，http://rast.nmpdr.org)[54]能夠快速注釋基因組信息，但由于參數(shù)的設(shè)置不同且算法的相對固定導(dǎo)致結(jié)果的準(zhǔn)確性有所下降。此外，在耐藥基因、毒力因子、質(zhì)粒分型以及菌株分子分型等特征的鑒定方面可分別使用CGE網(wǎng)站中的ResFinder[55]、VirulenceFinder[56]、PlasmidFinder[57]和MLST[58]等工具完成。該網(wǎng)站擁有方便快捷的獨(dú)立數(shù)據(jù)分析工具而且還有一個完整的批量數(shù)據(jù)分析流程 (Bacterial Analysis Pipeline)[59]，但是有些工具的數(shù)據(jù)庫并不全，例如VirulenceFinder中的數(shù)據(jù)庫僅有李斯特菌、金黃色葡萄球菌、大腸桿菌和腸球菌四個數(shù)據(jù)庫。

在命令類型的工具方面，PathoScope可以直接通過原始數(shù)據(jù)鑒別細(xì)菌的種屬而不需要對序列進(jìn)行拼接組裝[60]，該分析軟件多用于宏基因組數(shù)據(jù)的菌群結(jié)構(gòu)分析，如腸道微生物[61-62]、皮膚菌群結(jié)構(gòu)[63]的研究。Clinical PathoScope作為該工具為適應(yīng)臨床樣本分析而延伸出的附加程序，可在25 min內(nèi)從多種屬的樣品中檢測到致病菌，且準(zhǔn)確率可達(dá)94.7%[64]。PROKKA是一款快速注釋原核生物基因組信息的軟件，該軟件可以預(yù)測基因的位置和其相應(yīng)的功能[65]。該軟件在4核的電腦上完成大腸桿菌K-12的基因組注釋僅需6 min，且準(zhǔn)確率可達(dá)99.63%，此外一項(xiàng)試驗(yàn)表明PROKKA在預(yù)測基因的數(shù)量上優(yōu)于RAST[65]。在Stoesser等[66]、Yang等[67]的研究中，利用該軟件對NDM基因定位，并對其所在的質(zhì)粒進(jìn)行了注釋。此外，在耐藥基因、毒力因子、質(zhì)粒類型及分子分型方面，一款快速高效的軟件SRST2 (Short read sequencing typing) 可通過原始短序列對任何序列的數(shù)據(jù)庫進(jìn)行比對計(jì)算并得到準(zhǔn)確的結(jié)果[68]。SRST2推薦使用耐藥基因數(shù)據(jù)庫ARDB[69]、毒力因子數(shù)據(jù)庫VFDB[70]以及分子分型數(shù)據(jù)庫PubMLST (https://pubmlst.org/databases. shtml)。由于該軟件可利用不同數(shù)據(jù)庫完成耐藥菌的多種分子特征分析，從而使其成為細(xì)菌耐藥性研究中應(yīng)用最廣泛的軟件之一。例如Shen等[32]、Wang等[71]在其各自的研究中 (細(xì)菌耐藥性的分子流行病學(xué)、耐藥基因的傳播途徑等) 中均使用了該軟件進(jìn)行分子特征的定義。各類軟件信息見表4。

5 SNP分析及進(jìn)化樹構(gòu)建

獲得菌株內(nèi)在信息 (基本特征) 后，進(jìn)一步研究不同菌株之間的外部聯(lián)系 (親緣關(guān)系)，對于是否存在克隆傳播或是判斷流行株的暴發(fā)來源有著重要的意義。上文提及的MLST可作為菌株間關(guān)聯(lián)性的一種指標(biāo)，但因?yàn)檫@種方法只針對菌株的7個基因進(jìn)行種類的劃分，在精確性上有著一定的局限性，所以我們將針對其他精確性更高的基于細(xì)菌間單核苷酸多態(tài)性 (Single nucleotide polymorphism，SNP) 的方法進(jìn)行討論。不同的方法在精確性、操作性等方面有著較大的差異，但基本都需要通過代碼進(jìn)行操作且對計(jì)算機(jī)硬件要求高。

CSI Phylogeny是一款在CGE網(wǎng)站上可用的基于參考序列檢測SNP的工具，具有高保守性和準(zhǔn)確性高的特點(diǎn)[72]。這款軟件使用Burrows-Wheeler Aligner (BWA) 將目標(biāo)序列與參考序列進(jìn)行比對后通過用戶設(shè)置的參數(shù)提取出SNPs，并可檢查SNP是否在所有序列中均存在，最終通過FastTree軟件[73]構(gòu)建最大釋然法樹。

NDtree是另一款在CGE網(wǎng)站上可用的工具，它可將原始序列生成-mers后與參考序列進(jìn)行比對，然后通過公式計(jì)算菌株間SNP的數(shù)量。這些數(shù)據(jù)將生成一個矩陣 (Matrix) 文件并使用Phylip (http://evolution.genetics.washington.edu/phylip.html)計(jì)算各菌株間的進(jìn)化關(guān)系。值得一提的是該方法具有保守程度的參數(shù)設(shè)定，可能會導(dǎo)致不正確的結(jié)果。

kSNP3是一款可不使用參考序列并不需提供多序列校準(zhǔn)文件的工具[74]。該軟件使用-mer分析去推斷SNPs并適用于任何種屬細(xì)菌。kSNP3可選擇對核心SNPs進(jìn)行分析，生成多個文件包括核心和非核心SNPs信息的文件、Newick格式的簡約法 (Parsimony)、鄰接法 (Neighbor-joining) 和最大釋然法樹的文件。Wilson等利用該軟件繪制了來源于澳大利亞食品生產(chǎn)鏈條中耐藥李斯特菌的系統(tǒng)發(fā)育樹，確定了其耐藥性的關(guān)鍵遺傳標(biāo)記[75]。

Roary是一款命令類型工具，可快速檢測各菌株泛基因組 (Pan-genome) 中SNPs的情況。由于該軟件每個樣本需使用注釋的拼接序列，所以所有的分析菌株必須來自同一個種[76]。序列中的編碼區(qū)將會被CD-HIT工具轉(zhuǎn)換為蛋白序列，然后使用BLASTP工具對所有菌株中的蛋白序列進(jìn)行搜索比對，最終通過Markov cluster algorithm (MCL) 算法[77]分為不同組別后與CD-HIT先前轉(zhuǎn)換的序列合并為最后的結(jié)果。來自1 000株鼠傷寒沙門菌序列的分析在單核CPU的計(jì)算機(jī)上只運(yùn)行了4.3 h，并且泛基因組的正確率達(dá)到了100%[76]。因其快速高效，該軟件多用于大量數(shù)據(jù)的分析，例如Moradigaravand等通過該軟件對10年內(nèi)分離的205株粘質(zhì)沙雷氏菌進(jìn)化分析后發(fā)現(xiàn)，該菌在進(jìn)化過程中多次獲得不同耐藥基因，揭示了多重耐藥粘質(zhì)沙雷氏菌的陽性率逐漸上升的原因[78]。

4 各類分析軟件的特征

Table 4 Features, advantages and disadvantages of various analysis software

Pan-Seq是另一款比較序列間泛基因組差異的軟件，包括3個組件：Novel region finder (NRF)、Core and accessory genome finder (CAGF) 和Locus selector (LS)[79]。NRF工具使用MUMmer[80]鑒別序列間的差異位點(diǎn)后CAGF通過MUMmer alignment將這些差異位點(diǎn)序列添加到初始的泛基因組中，再把泛基因組分成碎片與原始序列進(jìn)行校對，通過BLASTn算法將不同的碎片分為核心和附加基因組，最終通過LS工具可識別輸入序列間不同基因的SNPs。López-Camacho等人利用該軟件分析了肺炎克雷伯菌在燒傷重癥監(jiān)護(hù)病房內(nèi)暴發(fā)期間的耐藥性演變，確定了抗生素選擇壓力在此次暴發(fā)的出現(xiàn)和進(jìn)化過程中發(fā)揮了重要作用[81]。

Lyve-SET是最近新報(bào)道的基于參考序列SNP的高質(zhì)量SNP分析工具。該工具通過序列最小最大覆蓋度 (Coverage) 分辨堿基一致性，并舍棄那些只在單向序列里出現(xiàn)的SNPs，同時可選擇性排除特殊的噬菌體 (Phage-specific regions) 和重復(fù)區(qū)域 (Repeat regions) 以提高SNPs的準(zhǔn)確性和可靠性[82]。

Harvest是一款基于核心基因組SNP序列分析的軟件，可快速提取菌株間核心基因組的SNPs信息并構(gòu)建進(jìn)化樹，同時包含插件Gingr可動態(tài)的可視化數(shù)據(jù)。該軟件既可指定參考序列，也可隨機(jī)選擇參考序列，但是兩個菌株間差異過大的話可能會被排除[83]。Shen等基于該軟件的系統(tǒng)發(fā)育分析，結(jié)合種群結(jié)構(gòu)貝葉斯分析，發(fā)現(xiàn)來源于中國30個省市287株-positive測序分離株具有4個明顯的譜系，且4個譜系與省份無明顯相關(guān)性[32]。

獲取SNPs信息后可通過不同的算法對菌株間的進(jìn)化關(guān)系構(gòu)建進(jìn)化樹，其中貝葉斯 (Bayesian) 和最大釋然法是最常用的方法，它們相較于鄰接法和簡約法具有更高的精確性，但是對于大量的樣本分析需要耗費(fèi)過多的計(jì)算機(jī)資源，存在一定的缺陷[84]。Holt等在2013年利用該軟件對越南地區(qū)索氏志賀菌的系統(tǒng)發(fā)育進(jìn)行了分析，為病原體在新宿主群體中的微進(jìn)化提供了一種獨(dú)特的，高分辨率的研究思路[29]。

Randomized Axelerated Maximum Likelihood (RAxML) 是一款以最大釋然法為基礎(chǔ)的進(jìn)化樹構(gòu)建軟件[85]。該算法首先生成一個假定的進(jìn)化樹，然后經(jīng)過幾步的優(yōu)化調(diào)整后得到新的進(jìn)化樹，當(dāng)經(jīng)過優(yōu)化后的進(jìn)化樹不再增加進(jìn)化關(guān)系的合理性時便停止重復(fù)該過程。當(dāng)用戶使用該軟件處理核苷酸和SNP數(shù)據(jù)時，必須使用GTRCAT (GTR，General time-reversible) 模型通過添加命令“-m GTRCAT”來校準(zhǔn)進(jìn)化樹，但當(dāng)樣本數(shù)小于50時，不推薦使用該模型。該軟件主要用于大樣本數(shù)據(jù)的系統(tǒng)發(fā)育分析，例如Casali等利用該軟件完成了1 035株耐藥結(jié)核桿菌的系統(tǒng)發(fā)育分析，闡明了耐藥結(jié)核桿菌在俄羅斯人群中的演變過程[86]。

FastTree是另一款基于最大釋然法開發(fā)的軟件[73]，通過4個階段完成進(jìn)化樹的構(gòu)建。1) 創(chuàng)建一個起始的進(jìn)化樹和存儲內(nèi)部節(jié)點(diǎn)的文件；2) 初始進(jìn)化樹的長度通過調(diào)換相鄰節(jié)點(diǎn)和重排分枝后逐漸減少；3) 通過數(shù)學(xué)模型得到最大釋然樹(CAT)；4) 進(jìn)化樹的可信度需通過Shimodaira-Hasegawa (SH) 測試[87]。同樣類似于RAxML軟件，F(xiàn)astTree在處理核苷酸和SNP數(shù)據(jù)時也需要使用GTR+CAT模型通過添加命令“-gtr”來完成進(jìn)化樹的校準(zhǔn)。該軟件常用于宏基因組數(shù)據(jù)的“系統(tǒng)發(fā)育重建”，對揭示微生態(tài)菌群結(jié)構(gòu)的進(jìn)化過程有一定的作用。

相較于最大釋然法，MrBayes則是基于貝葉斯算法設(shè)計(jì)的一款軟件，該軟件使用Markov chain Monte Carlo (MCMC) 算法，自定義選項(xiàng)和參數(shù)設(shè)定較多，因此難以操作[88]。

進(jìn)化樹的文件最常用的是NEXUS和Newick格式，許多圖形用戶界面 (GUI，Graphical User Interface) 的軟件比如FigTree (http://tree.bio.ed. ac.uk/software/figtree/)、MEGA[89]、Archaeopteryx[90]以及網(wǎng)頁類型的工具iTOL[91]等均可對進(jìn)化樹進(jìn)行可視化的分析。相關(guān)基因組比對軟件以及進(jìn)化樹構(gòu)建軟件見表5。

6 討論

隨著細(xì)菌耐藥性問題的日益嚴(yán)重，針對致病性耐藥菌暴發(fā)甚至共生菌的耐藥性監(jiān)控等相關(guān)研究顯得越來越重要，尤其是這些研究可作為提出暴發(fā)控制、耐藥性逆轉(zhuǎn)等理論的基礎(chǔ)。值得關(guān)注的是，隨著全基因組測序技術(shù)的迅速發(fā)展，測序質(zhì)量不斷提高的同時，測序費(fèi)用呈顯著下降趨勢。該技術(shù)相較于傳統(tǒng)的技術(shù)手段，可獲得更全面、更準(zhǔn)確的結(jié)果，將會在未來研究中得到廣泛的應(yīng)用[92]。

除了DNA提取質(zhì)量上的差別外，文庫質(zhì)量對于測序結(jié)果的影響也十分顯著，不同DNA文庫構(gòu)建試劑盒有著不同的特點(diǎn)[93]。很顯然物理破碎DNA的方法相較于化學(xué)法能帶來更少的誤差對樣品GC含量的嗜好也相對較少。令人振奮的是，目前已有許多公司支持機(jī)器自動化的操作來構(gòu)建DNA文庫，大量減少了人力和物力，同時減少了人工操作處理樣本時帶來的實(shí)驗(yàn)誤差。此外，將文庫構(gòu)建步驟簡化、趨向無PCR步驟以及減少磁珠純化等均可減少文庫構(gòu)建過程中所耗費(fèi)的時間，也可獲得更加趨近于真實(shí)的基因組信息。

目前主流的測序平臺還是由Illumina公司推出的幾款平臺，很顯然，在數(shù)據(jù)的產(chǎn)出量上這幾款平臺有著明顯的優(yōu)勢，并且測序結(jié)果中的堿基錯誤率也相對較低，對于大量單菌或是宏基因組的測序推薦選擇這類平臺。相反的ThermoFisher公司開發(fā)的幾款測序平臺在測序讀長上有著很大的優(yōu)勢，但是相對的產(chǎn)出量上則低了很多，這對于緊急情況單菌的測序有著潛在的應(yīng)用前景。目前，二代測序技術(shù)的測序讀長僅在300–600 bp之間，對于富含重復(fù)序列的細(xì)菌來說，難以得到一個完整的結(jié)果。針對這樣的情況，三代測序技術(shù)——PacBio平臺[39]以及MinION平臺[94]在其測序讀長 (平均>10 kb) 上有著很大的優(yōu)勢，能夠解決上述二代測序技術(shù)所遇到的問題，可以得到細(xì)菌的完整基因組信息即完成圖。然而測序讀長長的同時單堿基的錯誤率便會升高，且測序深度和產(chǎn)出量均不高，所以若能結(jié)合二代測序的數(shù)據(jù)進(jìn)行校對，那么得出的結(jié)果將會十分可靠。目前Pilon軟件即可用二代數(shù)據(jù)對三代數(shù)據(jù)進(jìn)行拼接結(jié)果的校準(zhǔn)[94]。與此同時，在PacBio平臺測序價格高居不下的情況下，測序成本也將增加，雖然MinION平臺的價格相對低一些，但Nanopore測序技術(shù)還處于初級階段，該技術(shù)及后續(xù)的生物信息學(xué)分析還需進(jìn)一步的開發(fā)。目前，已有不少研究者選用MinION平臺開展耐藥性的研究，Li等利用該測序技術(shù)結(jié)合二代測序技術(shù)不僅獲得了質(zhì)粒完整序列[96]，而且對插入序列IS在沙門菌中的質(zhì)粒異質(zhì)性進(jìn)行了研究[97]。Ludden等利用同樣的技術(shù)證明了污水中攜帶碳青霉烯耐藥基因質(zhì)粒可能在不同細(xì)菌之間交換，為環(huán)境中耐藥基因傳播機(jī)制的研究提供了一種新的研究方法[87]。

表5 基因組比較分析工具以及進(jìn)化樹構(gòu)建軟件特點(diǎn)

雖然目前有許多商業(yè)化的軟件可對WGS數(shù)據(jù)進(jìn)行分析，例如BioNumerics (Applied Maths，Biomérieux)、CLC Genomic Workbench (Qiagen) 和SeqSphere (Ridom) 等，這些軟件雖可完成多種生物信息學(xué)分析，如序列組裝拼接、SNPs提取以及進(jìn)化樹構(gòu)建等，但每個模塊需要單獨(dú)購買使用權(quán)。雖然商業(yè)化軟件價格高昂，但它們具有友好的用戶界面，加上有非常詳盡的用戶手冊可以幫助用戶去了解和使用它們的功能，降低了使用者的門檻。此外，本文也介紹了少數(shù)幾款基于網(wǎng)頁類型的工具，如包含了眾多分析模塊的CGE網(wǎng)站，其開源免費(fèi)的特點(diǎn)造成了巨大的訪問量，導(dǎo)致其效率偏低，且由于其屬于公共網(wǎng)絡(luò)，數(shù)據(jù)的保密措施還有待商榷。與上述工具不同的是，我們主要介紹了多款命令類型的軟件，這些軟件均為開源免費(fèi)的，大多數(shù)都需要在Linux、Ubuntu、Mac OS X等系統(tǒng)環(huán)境下運(yùn)行，少數(shù)可在Window系統(tǒng)下運(yùn)行，因此限制了這些軟件在非生物信息學(xué)背景的研究人員中使用。另外在批量處理數(shù)據(jù)時對系統(tǒng)硬件要求相對較高，所有操作均需使用代碼完成，且無可視化的圖形操作界面，給無生物信息學(xué)背景的分析人員又增添了許多困難。

7 展望

全基因組測序技術(shù)尤其是長片段測序技術(shù)的發(fā)展，極大地推動了細(xì)菌耐藥性的研究進(jìn)程，尤其是在新耐藥基因的發(fā)現(xiàn)、耐藥基因的傳播機(jī)制、耐藥菌/耐藥基因的進(jìn)化分析以及細(xì)菌耐藥性的風(fēng)險因素分析等方面都起到了至關(guān)重要的作用?；凇癘ne Health”的理念，耐藥性在人、動物與環(huán)境之間相互傳播的研究還需要更進(jìn)一步更深入的證據(jù)，而WGS技術(shù)已逐漸展露出其在細(xì)菌耐藥性研究中的潛力，可提供全面精確的數(shù)據(jù)，并將會成為未來多年內(nèi)最為重要的耐藥性研究手段之一。希望通過本文對文庫構(gòu)建、測序平臺、數(shù)據(jù)分析完整流程的詳細(xì)介紹，能夠讓更多的研究者對該技術(shù)有更全面的了解，并推進(jìn)此項(xiàng)技術(shù)在耐藥性領(lǐng)域的應(yīng)用。

[1] World Health Organization. Worldwide country situation analysis: response to antimicrobial resistance. Switzerland: World Health Organization, 2015.

[2] Laxminarayan R, Sridhar D, Blaser M, et al. Achieving global targets for antimicrobial resistance. Science, 2016, 353(6302): 874–875.

[3] Woolhouse M, Ward M, van Bunnik B, et al. Antimicrobial resistance in humans, livestock and the wider environment. Philos Trans Roy Soc B Biol Sci, 2015, 370(1670): 20140083.

[4] O’Neil J. Antimicrobial resistance: tackling a crisis for the health and wealth of nations. London: The Review on Antimicrobial Resistance, 2014.

[5] Kumarasamy KK, Toleman MA, Walsh TR, et al. Emergence of a new antibiotic resistance mechanism in India, Pakistan, and the UK: a molecular, biological, and epidemiological study. Lancet Infect Dis, 2010, 10(9): 597–602.

[6] Liu YY, Wang Y, Walsh TR, et al. Emergence of plasmid-mediated colistin resistance mechanism MCR-1 in animals and human beings in China: a microbiological and molecular biological study. Lancet Infect Dis, 2016, 16(2): 161–168.

[7] Dalmolin TV, Castro L, Mayer FQ, et al. Co-occurrence ofandKPC-2in a clinical isolate ofin Brazil. J Antimicrob Chemother, 2017, 72(8): 2404–2406.

[8] Sun P, Bi ZW, Nilsson M, et al. Occurrence ofKPC-2,CTX-M, andinfrom well water in rural China. Antimicrob Agents Chemother, 2017, 61(4): e02569-16.

[9] Tac?o M, dos Santos Tavares R, Teixeira P, et al.andKPC-3insequence type 744 after meropenem and colistin therapy, Portugal. Emerg Infect Dis, 2017, 23(8): 1419–1421.

[10] Wang Y, Tian GB, Zhang R, et al. Prevalence, risk factors, outcomes, and molecular epidemiology ofpositivein patients and healthy adults from China: an epidemiological and clinical study. Lancet Infect Dis, 2017, 17(4): 390–399.

[11] Wang Y, Zhang RM, Li JY, et al. Comprehensive resistome analysis reveals the prevalence of NDM and MCR-1 in Chinese poultry production. Nat Microbiol, 2017, 2: 16260.

[12] Zhong LL, Zhang YF, Doi Y, et al. Coproduction of MCR-1 and NDM-1 by colistin-resistantisolated from a healthy individual. Antimicrob Agents Chemother, 2017, 61(1): e01962-16.

[13] Bertelli C, Greub G. Rapid bacterial genome sequencing: methods and applications in clinical microbiology.Clin Microbiol Infect, 2013, 19(9): 803–813.

[14] Jünemann S, Sedlazeck FJ, Prior K, et al. Updating benchtop sequencing performance comparison. Nat Biotechnol, 2013, 31(4): 294–296.

[15] Loman NJ, Misra RV, Dallman TJ, et al. Performance comparison of benchtop high-throughput sequencing platforms. Nat Biotechnol, 2012, 30(5): 434–439.

[16] Hasnain SE, O’Toole RF, Grover S, et al. Whole genome sequencing: a new paradigm in the surveillance and control of human tuberculosis. Tuberculosis (Edinb), 2015, 95(2): 91–94.

[17] Lecuit M, Eloit M. The diagnosis of infectious diseases by whole genome next generation sequencing: a new era is opening. Front Cell Infect Microbiol, 2014, 4: 25.

[18] Choi SC. On the study of microbial transcriptomes using second- and third-generation sequencing technologies. J Microbiol, 2016, 54(8): 527–536.

[19] Xavier BB, Lammens C, Ruhal R, et al. Identification of a novel plasmid-mediated colistin-resistance gene,, in, Belgium, June 2016. EuroSurveill, 2016, 21(27): pii=30280.

[20] Yin WJ, Li H, Shen YB, et al. Novel plasmid-mediated colistin resistance genein. mBio, 2017, 8(3): e00543-17.

[21] Carattoli A, Villa L, Feudi C, et al. Novel plasmid-mediated colistin resistancegene inand, Italy 2013, Spain and Belgium, 2015 to 2016. Eurosurveillance, 2017, 22(31): 30589.

[22] Borowiak M, Fischer J, Hammerl JA, et al. Identification of a novel transposon-associated phosphoethanolamine transferase gene,, conferring colistin resistance in-tartrate fermentingenterica subsp. enterica serovar Paratyphi B. J Antimicrob Chemother, 2017, 72(12): 3317–3324.

[23] Wang XM, Wang Y, Zhou Y, et al. Emergence of a novel mobile colistin resistance gene,, in NDM-producing. Emerg Microbes Infect, 2018, 7: 122.

[24] Liu ZH, Wang Y, Walsh TR, et al. Plasmid-mediated novelNDM-17gene encoding a carbapenemase with enhanced activity in a sequence type 48strain. Antimicrob Agents Chemother, 2017, 61(5): e02233–16.

[25] Hadziabdic S, Fischer J, Malorny B, et al.transfer and microevolution of avian native IncA/C2NDM-1-Carrying Plasmid pRH-1238 during a broiler chicken infection study. Antimicrob Agents Chemother, 2018, 62(4): e02128–17.

[26] Kr?ger C, Kary SC, Schauer K, et al. Genetic regulation of virulence and antibiotic resistance in. Genes, 2017, 8(1): 12.

[27] Porse A, Sch?nning K, Munck C, et al. Survival and evolution of a large multidrug resistance plasmid in new clinical bacterial hosts. Mol Biol Evol, 2016, 33(11): 2860–2873.

[28] Holt KE, Baker S, Weill FX, et al.genome sequencing and phylogenetic analysis indicate recent global dissemination from Europe. Nat Genet, 2012, 44(9): 1056–1059.

[29] Holt KE, Thieu Nga TV, Thanh DP, et al. Tracking the establishment of local endemic populations of an emergent enteric pathogen. Proc Natl Acad Sci USA, 2013, 110(43): 17522–17527.

[30] Wong VK, Baker S, Pickard DJ, et al. Phylogeographical analysis of the dominant multidrug-resistant H58 clade ofidentifies inter- and intracontinental transmission events. Nat Genet, 2015, 47(6): 632–639.

[31] Duchêne S, Holt KE, Weill FX, et al. Genome-scale rates of evolutionary change in bacteria. Microb Genom, 2016, 2(11): e000094.

[32] Shen YB, Zhou HW, Xu J, et al. Anthropogenic and environmental factors associated with high incidence ofcarriage in humans across China. Nat Microbiol, 2018, 3(9): 1054–1062.

[33] Wedley AL, Dawson S, Maddox TW, et al. Carriage of antimicrobial resistantin dogs: prevalence, associated risk factors and molecular characteristics. Vet Microbiol, 2017, 199: 23–30.

[34] Costea PI, Zeller G, Sunagawa S, et al. Towards standards for human fecal sample processing in metagenomic studies. Nat Biotechnol, 2017, 35(11): 1069–1076.

[35] Marine R, Polson SW, Ravel J, et al. Evaluation of a transposase protocol for rapid generation of shotgun high-throughput sequencing libraries from nanogram quantities of DNA. Appl Environ Microb, 2011, 77(22): 8071–8079.

[36] https://image.baidu.com

[37] Bentley DR, Balasubramanian S, Swerdlow HP, et al. Accurate whole human genome sequencing using reversible terminator chemistry. Nature, 2008, 456(7218): 53–59.

[38] Yergeau E, Lawrence JR, Sanschagrin S, et al. Next-generation sequencing of microbial communities in the Athabasca River and its tributaries in relation to oil sands mining activities. Appl Environ Microb, 2012, 78(21): 7626–7637.

[39] Rhoads A, Au KF. PacBio sequencing and its applications. Genomics, Proteomics & Bioinformatics, 2015, 13(5): 278–289.

[40] Jain M, Olsen HE, Paten B, et al. The oxford nanopore minION: delivery of nanopore sequencing to the genomics community. Genome Biol, 2016, 17: 239.

[41] Ku CS, Roukos DH. From next-generation sequencing to nanopore sequencing technology: paving the way to personalized genomic medicine. Expert Rev Med Devices, 2013, 10(1): 1–6.

[42] Boisvert S, Laviolette F, Corbeil J. Ray: simultaneous assembly of reads from a mix of high-throughput sequencing technologies. J Comput Biol, 2010, 17(11): 1519–1533.

[43] Zerbino DR, Birney E. Velvet: algorithms forshort read assembly using de Bruijn graphs. Genome Res, 2008, 18(5): 821–829.

[44] Bankevich A, Nurk S, Antipov D, et al. SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J Comput Biol, 2012, 19(5): 455–477.

[45] Peng Y, Leung HCM, Yiu SM, et al. IDBA-UD: aassembler for single-cell and metagenomic sequencing data with highly uneven depth. Bioinformatics, 2012, 28(11): 1420–1428.

[46] Chin CS, Alexander DH, Marks P, et al. Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data. Nat Methods, 2013, 10(6): 563–569.

[47] Koren S, Walenz BP, Berlin K, et al. Canu: scalable and accurate long-read assembly via adaptive-mer weighting and repeat separation. Genome Res, 2017, 27(5): 722–736.

[48] Wick RR, Judd LM, Gorrie CL, et al. Unicycler: resolving bacterial genome assemblies from short and long sequencing reads. PLoS Comput Biol, 2017, 13(6): e1005595.

[49] Antipov D, Hartwick N, Shen M, et al. plasmidSPAdes: assembling plasmids from whole genome sequencing data. Bioinformatics, 2016, 32(22): 3380–3387.

[50] De Lannoy C, De Ridder D, Risse J. The long reads ahead:genome assembly using the MinION. F1000Res, 2017, 6: 1083.

[51] Karalius J. Resources for advanced bioinformaticians working in plant and animal genomes with SMRT sequencing. Menlo Park, CA: Pacific Biosciences, 2015.

[52] Au KF, Underwood JG, Lee L, et al. Improving PacBio long read accuracy by short read alignment. PLoS ONE, 2012, 7(10): e46679.

[53] Larsen MV, Cosentino S, Lukjancenko O, et al. Benchmarking of methods for genomic taxonomy. J Clin Microbiol, 2014, 52(5): 1529–1539.

[54] Aziz RK, Bartels D, Best AA, et al. The RAST server: rapid annotations using subsystems technology. BMC Genomics, 2008, 9: 75.

[55] Zankari E, Hasman H, Cosentino S, et al. Identification of acquired antimicrobial resistance genes. J Antimicrob Chemother, 2012, 67(11): 2640–2644.

[56] Joensen KG, Scheutz F, Lund O, et al. Real-time whole-genome sequencing for routine typing, surveillance, and outbreak detection of verotoxigenic. J Clin Microbiol, 2014, 52(5): 1501–1510.

[57] Carattoli A, Zankari E, García-Fernández A, et al.detection and typing of plasmids using PlasmidFinder and plasmid multilocus sequence typing. Antimicrob Agents Chemother, 2014, 58(7): 3895–3903.

[58] Larsen MV, Cosentino S, Rasmussen S, et al. Multilocus sequence typing of total-genome-sequenced bacteria. J Clin Microbiol, 2012, 50(4): 1355–1361.

[59] Thomsen MCF, Ahrenfeldt J, Cisneros JLB, et al. A bacterial analysis platform: an integrated system for analysing bacterial whole genome sequencing data for clinical diagnostics and surveillance. PLoS ONE, 2016, 11(6): e0157718.

[60] Hong CJ, Manimaran S, Shen Y, et al. PathoScope 2.0: a complete computational framework for strain identification in environmental or clinical sequencing samples. Microbiome, 2014, 2: 33.

[61] Faith JJ, Guruge JL, Charbonneau M, et al. The long-term stability of the human gut microbiota. Science, 2013, 341(6141): 1237439.

[62] Cox LM, Yamanishi S, Sohn J, et al. Altering the intestinal microbiota during a critical developmental window has lasting metabolic consequences. Cell, 2014, 158(4): 705–721.

[63] Oh J, Byrd AL, Deming C, et al. Biogeography and individuality shape function in the human skin metagenome. Nature, 2014, 514(7520): 59–64.

[64] Byrd AL, Perez-Rogers JF, Manimaran S, et al. Clinical PathoScope: rapid alignment and filtration for accurate pathogen identification in clinical samples using unassembled sequencing data. BMC Bioinformatics, 2014, 15: 262.

[65] Seemann T. Prokka: rapid prokaryotic genome annotation. Bioinformatics, 2014, 30(14): 2068–2069.

[66] Stoesser N, Giess A, Batty EM, et al. Genome sequencing of an extended series of NDM-producingisolates from neonatal infections in a Nepali hospital characterizes the extent of community- versus hospital-associated transmission in an endemic setting. Antimicrob Agents Chemother, 2014, 58(12): 7347–7357.

[67] Yang P, Xie Y, Feng P, et al.NDM-5carried by an IncX3 plasmid insequence type 167. Antimicrob Agents Chemother, 2014, 58(12): 7548–7552.

[68] Inouye M, Dashnow H, Raven LA, et al. SRST2: rapid genomic surveillance for public health and hospital microbiology labs. Genome Med, 2014, 6: 90.

[69] Liu B, Pop M. ARDB—antibiotic resistance genes database. Nucleic Acids Res, 2009, 37(S1): D443–D447.

[70] Chen LH, Zheng DD, Liu B, et al. VFDB 2016: hierarchical and refined dataset for big data analysis—10 years on. Nucleic Acids Res, 2016, 44(D1): D694–D697.

[71] Wang Y, Zhang RM, Li JY, et al. Comprehensive resistome analysis reveals the prevalence of NDM and MCR-1 in Chinese poultry production. Nat Microbiol, 2017, 2: 16260.

[72] Kaas RS, Leekitcharoenphon P, Aarestrup FM, et al. Solving the problem of comparing whole bacterial genomes across different sequencing platforms. PLoS ONE, 2014, 9(8): e104984.

[73] Price MN, Dehal PS, Arkin AP. FastTree: computing large minimum evolution trees with profiles instead of a distance matrix. Mol Biol Evol, 2009, 26(7): 1641–1650.

[74] Gardner SN, Slezak T, Hall BG. kSNP3.0: SNP detection and phylogenetic analysis of genomes without genome alignment or reference genome. Bioinformatics, 2015, 31(17): 2877–2878.

[75] Wilson A, Gray J, Chandry PS, et al. Phenotypic and genotypic analysis of antimicrobial resistance amongisolated from Australian food production chains. Genes, 2018, 9(2): 80.

[76] Page AJ, Cummins CA, Hunt M, et al. Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics, 2015, 31(22): 3691–3693.

[77] Enright AJ, van Dongen S, Ouzounis CA. An efficient algorithm for large-scale detection of protein families. Nucleic Acids Res, 2002, 30(7): 1575–1584.

[78] Moradigaravand D, Boinett CJ, Martin V, et al. Recent independent emergence of multiple multidrug-resistantclones within the United Kingdom and Ireland. Genome Res, 2016, 26(8): 1101–1109.

[79] Laing C, Buchanan C, Taboada EN, et al. Pan-genome sequence analysis using Panseq: an online tool for the rapid analysis of core and accessory genomic regions. BMC Bioinformatics, 2010, 11: 461.

[80] Kurtz S, Phillippy A, Delcher AL, et al. Versatile and open software for comparing large genomes. Genome Biol, 2004, 5(2): R12.

[81] López-Camacho E, Gómez-Gil R, Tobes R, et al. Genomic analysis of the emergence and evolution of multidrug resistance during aoutbreak including carbapenem and colistin resistance. J Antimicrob Chemother, 2014, 69(3): 632–636.

[82] Katz LS, Griswold T, Williams-Newkirk AJ, et al. A comparative analysis of the Lyve-SET phylogenomics pipeline for genomic epidemiology of foodborne pathogens. Front Microbiol, 2017, 8: 375.

[83] Treangen TJ, Ondov BD, Koren S, et al. The Harvest suite for rapid core-genome alignment and visualization of thousands of intraspecific microbial genomes. Genome Biol, 2014, 15(11): 524.

[84] Strimmer K, von Haeseler A. Likelihood-mapping: a simple method to visualize phylogenetic content of a sequence alignment. Proc Natl Acad Sci USA, 1997, 94(13): 6815–6819.

[85] Stamatakis A, Ludwig T, Meier H. RAxML-III: a fast program for maximum likelihood-based inference of large phylogenetic trees.Bioinformatics, 2005, 21(4): 456–463.

[86] Casali N, Nikolayevskyy V, Balabanova Y, et al. Evolution and transmission of drug-resistant tuberculosis in a Russian population. Nat Genet, 2014, 46(3): 279–286.

[87] Price MN, Dehal PS, Arkin AP. FastTree 2-approximately maximum-likelihood trees for large alignments. PLoS ONE, 2010, 5(3): e9490.

[88] Ronquist F, Teslenko M, van der Mark P, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol, 2012, 61(3): 539–542.

[89] Kumar S, Stecher G, Tamura K. MEGA7: molecular evolutionary genetics analysis version 7.0 for bigger datasets. Mol Biol Evol, 2016, 33(7): 1870–1874.

[90] Han MV, Zmasek CM. phyloXML: XML for evolutionary biology and comparative genomics. BMC Bioinformatics, 2009, 10: 356.

[91] Letunic I, Bork P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic Acids Res, 2016, 44(W1): W242–W245.

[92] Land M, Hauser L, Jun SR, et al. Insights from 20 years of bacterial genome sequencing. Funct Integr Genomics, 2015, 15(2): 141–161.

[93] Aigrain L, Gu Y, Quail MA. Quantitation of next generation sequencing library preparation protocol efficiencies using droplet digital PCR assays—a systematic comparison of DNA library preparation kits for Illumina sequencing. BMC Genomics, 2016, 17: 458.

[94] van der Helm E, Imamovic L, Hashim Ellabaan MM, et al. Rapid resistome mapping using nanopore sequencing. Nucleic Acids Res, 2017, 45(8): e61.

[95] Walker BJ, Abeel T, Shea T, et al. Pilon: an integrated tool for comprehensive microbial variant detection and genome assembly improvement. PLoS ONE, 2014, 9(11): e112963.

[96] Li R, Xie M, Dong N, et al. Efficient generation of complete sequences of MDR-encoding plasmids by rapid assembly of MinION barcoding sequencing data. GigaScience, 2018, 7(3): 1–9.

[97] Li R, Chen K, Chan E W-C, et al. Resolution of dynamic MDR structures among the plasmidome of Salmonella using MinION single-molecule, long-read sequencing. J Antimicrob Chemother, 2018, 73(10): 2691–2695.

Application of whole genome sequencing technology and bioinformatics analysis in antimicrobial resistance researches

Yingbo Shen*, Xiaomin Shi*, Jianzhong Shen, Yang Wang, and Shaolin Wang

College of Veterinary Medicine, China Agricultural University, Beijing 100193, China

The emergence and spread of antimicrobial resistance has become a serious global issue. Bacterial characteristics, such as antimicrobial resistance genes, virulence-associated genes, plasmid types, and phylogenetic relationship among different strains, are the keys to unravel the occurrence and dissemination of antimicrobial resistance. However, the accuracy and efficiency of the traditional techniques, such as polymerase chain reaction and pulsed field gel electrophoresis is insufficient to underlying the mystery of antimicrobial resistance. Recently, the whole genome sequencing and high-throughput bioinformatics analysis have been successfully used in antimicrobial resistance studies, helping scientists to obtain the nature of antimicrobial resistance bacteria quickly, and more precisely to paint the evolutionary relationship among different strains. Therefore, in this study, we aim to systematically introduce the recent development of whole genome sequencing analysis, including different methods and corresponding characteristics of library preparation, platform sequencing, data analysis, and the latest application of the technology in the antimicrobial resistance research. We hope that this review can provide more comprehensive knowledge about whole genome sequencing and bioinformatic analysis for antimicrobial resistance research.

antimicrobial resistance, whole genome sequencing, bioinformatics

10.13345/j.cjb.180350

August 30, 2018;

December 17, 2018

National Key Research and Development Program of China (No. 2016YFD0501301), National Natural Science Foundation of China (No. 31572568).

Shaolin Wang. Tel: +86-10-62734255; E-mail: shaolinwang@cau.edu.cn

*These authors contributed equally to this study.

科技部重點(diǎn)研發(fā)計(jì)劃(No. 2016YFD0501301)，國家自然科學(xué)基金(No. 31572568) 資助。

(本文責(zé)編 郝麗芳)