一株臨床耐藥菌高通量測(cè)序數(shù)據(jù)的深度挖掘

傅行禮 鄭學(xué)明 解冰 劉永平 高玉華 吳曉蓉 孔德華

(1江蘇大學(xué)醫(yī)學(xué)部，江蘇 鎮(zhèn)江 212001；2鎮(zhèn)江市丹徒區(qū)中醫(yī)院；3南京市高淳人民醫(yī)院)

臨床上很多疾病的發(fā)生是由于致病菌的感染，而且很多致病菌已經(jīng)進(jìn)化為對(duì)多種抗生素有抗性的菌株〔1〕。傳統(tǒng)的實(shí)驗(yàn)室培養(yǎng)方法可以得到細(xì)菌的形態(tài)學(xué)特征，從而知道該細(xì)菌的類(lèi)型，抗生素敏感實(shí)驗(yàn)也能確定該菌株對(duì)哪些抗生素有抗性，這些結(jié)果對(duì)疾病的治療提供了重要的信息。但要研究該菌株的基因組特征、來(lái)源及抗性機(jī)制等就必須進(jìn)行全基因組測(cè)序。隨著高通量測(cè)序技術(shù)的發(fā)展，細(xì)菌基因組測(cè)序已經(jīng)同聚合酶鏈?zhǔn)椒磻?yīng)(PCR) 一樣方便和普及，并成為分子生物學(xué)研究的常規(guī)技術(shù)〔2〕。

高通量二代測(cè)序大大促進(jìn)了生物信息學(xué)的發(fā)展，生物信息學(xué)家已經(jīng)開(kāi)發(fā)出了各種各樣的軟件工具對(duì)測(cè)序數(shù)據(jù)進(jìn)行處理，但每一種軟件都只解決數(shù)據(jù)處理某一個(gè)環(huán)節(jié)的一個(gè)具體問(wèn)題〔3〕。怎么分析研究大量的測(cè)序原始數(shù)據(jù)，進(jìn)而挖掘出有用的生物學(xué)信息，對(duì)于生物信息學(xué)家來(lái)說(shuō)依然是一個(gè)巨大的挑戰(zhàn)。而且如果能利用好現(xiàn)有的生物信息學(xué)工具對(duì)測(cè)序原始數(shù)據(jù)進(jìn)行充分挖掘?qū)?huì)為實(shí)驗(yàn)研究提供重要的線(xiàn)索和方向，必將會(huì)大大加快研究的步伐。

本研究通過(guò)對(duì)一株臨床耐藥菌高通量測(cè)序數(shù)據(jù)的分析研究探討信息提取、數(shù)據(jù)挖掘的流程。

1 材料和方法

1.1高通量測(cè)序數(shù)據(jù)的下載和質(zhì)量控制 從美國(guó)國(guó)立生物技術(shù)信息中心(NCBI)的序列片段歸檔數(shù)據(jù)庫(kù)(SRA)下載高通量全基因組測(cè)序(WGS)的原始數(shù)據(jù)(SRR2134675)，測(cè)序樣品是從成年女性血液中分離得到的Escherichia coli 菌株，測(cè)序方法是Illumina MiSeq雙端(PE)測(cè)序?？股乜剐詫?shí)驗(yàn)表明該菌株對(duì)氨芐青霉素(ampicillin)、環(huán)丙沙星(ciprofloxacin)、磷霉素(fosfomycin)等多種抗生素有耐藥性。下載的SRA文件用sratools 的fastq-dump 命令解壓生成兩個(gè)fastq文件。用Sickle 軟件去掉平均質(zhì)量值(Phred qualitiy)低于20(1% 測(cè)序錯(cuò)誤)的序列，再去掉含有不確定堿基(N)數(shù)大于10%及長(zhǎng)度不到30 bp的序列，從而得到高質(zhì)量的數(shù)據(jù)進(jìn)行下一步的分析。

1.2基因組denovo 組裝和注釋 用華大基因開(kāi)發(fā)的SOAPdenovo軟件對(duì)處理后的fastq 數(shù)據(jù)進(jìn)行從頭組裝〔4〕。 選擇不同的kmer 優(yōu)化組裝的結(jié)果，最后得到該菌株基因組組裝的重疊群(Scaffolds)，選擇大于500 bp的重疊群序列進(jìn)行下一步的分析。用細(xì)菌基因組注釋軟件Prokka對(duì)該菌株基因組草圖進(jìn)行注釋〔5〕。

1.3COG分類(lèi)分析 從NCBI下載微生物蛋白質(zhì)相鄰類(lèi)的聚簇(COG)數(shù)據(jù)，建立本地blast數(shù)據(jù)庫(kù)，用該菌株預(yù)測(cè)的所有蛋白質(zhì)序列(fasta格式)逐一和本地COG 數(shù)據(jù)庫(kù)中的蛋白質(zhì)序列進(jìn)行blastp 比對(duì)(E 值≤10-10,同一性≥40%,覆蓋長(zhǎng)度≥80%)，選擇匹配分?jǐn)?shù)最高蛋白質(zhì)的COG 分類(lèi)注釋查詢(xún)蛋白，統(tǒng)計(jì)得到該基因組編碼蛋白質(zhì)的COG功能分類(lèi)〔6〕。

1.4抗生素抗性基因分析 將該菌株基因組注釋的蛋白質(zhì)序列上傳到抗生素抗性數(shù)據(jù)庫(kù)(ARDB)，使用blastp 程序和數(shù)據(jù)庫(kù)中的抗性蛋白進(jìn)行序列比對(duì)(E 值≤10-7,同一性≥40%,覆蓋長(zhǎng)度≥80%)，根據(jù)同源性分析該菌株基因組中的抗性基因及抗性類(lèi)型〔7〕。

1.5測(cè)序短序列(Reads)回帖(Read Mapping)和突變檢出(Variation calling) 用多位點(diǎn)序列分型(MLST)軟件MLST2.0分析該菌株基因組組裝的Scaffolds序列，從而得到該菌株的分子分型〔8〕。通過(guò)分子分型找到和該菌株親緣關(guān)系非常近的菌株作為參考基因組，用BWA( Burrows-Wheeler-Alignment Tool)軟件把質(zhì)控后的Reads Read Mapping到參考基因組上得到SAM(Sequence Alignment/Map) 文件，然后用SAMtools軟件進(jìn)行SAM-BAM格式轉(zhuǎn)換及排序(Sorting)，再用GATK軟件對(duì)INDEL附近的Reads進(jìn)行局部的重新比對(duì)，從而將比對(duì)的錯(cuò)誤率降到最低，最后用BCFtools 軟件檢出突變(Variation Calling)：?jiǎn)魏塑账岫鄳B(tài)性(SNPs)和插入缺失突變(INDELs)〔9～12〕。過(guò)濾掉等位基因頻率低于75%及測(cè)序深度(Read-Depth)小于5的突變位點(diǎn)，生成包含高可信度突變位點(diǎn)信息的VCF(Variant Call Format)文件。

1.6突變位點(diǎn)的分析 從NCBI 下載大腸桿菌T131菌株EC958基因組的序列文件(fasta)及基因組注釋文件(GFF3)，用snpEff 軟件建立該基因組的注釋數(shù)據(jù)庫(kù)〔13〕。 用snpEff軟件對(duì)上一步BCFtools軟件生成的VCF 文件進(jìn)行注釋?zhuān)和蛔兾稽c(diǎn)所在基因，突變對(duì)基因表達(dá)的影響，突變的危害程度等。

2 結(jié) 果

2.1基因組的一般特征分析 下載的全基因組測(cè)序數(shù)據(jù)共包含2 359 983對(duì)測(cè)序片段Reads，共1.1 G個(gè)堿基的序列數(shù)據(jù)。FastQC軟件分析發(fā)現(xiàn)測(cè)序數(shù)據(jù)不含接頭(Adaptor)和條形碼(Barcode)序列，但其中一個(gè)文件的測(cè)序數(shù)據(jù)的3′末端序列平均質(zhì)量不高，質(zhì)量差異大。經(jīng)過(guò)Sickle 軟件質(zhì)控后得到2 290 048 對(duì)高質(zhì)量的測(cè)序片段。 SOAPdenovo 軟件從頭組裝該菌株基因組共產(chǎn)生584 個(gè)重疊群，N50為 230 351 bp，GC含量為50.80%，預(yù)測(cè)基因組大小為5 278 663 bp。 組裝時(shí)kmer 的選擇會(huì)大大影響組裝的質(zhì)量，最好的辦法就是測(cè)試不同的kmer,然后看N50的大小,選擇N50 最高的組裝結(jié)果。測(cè)序的數(shù)據(jù)量和測(cè)序質(zhì)量是決定基因組組裝好壞的最根本的因素，在測(cè)序數(shù)據(jù)量一定的條件下，選擇高質(zhì)量的序列進(jìn)行組裝會(huì)大大提高N50值。

Prokka 軟件對(duì)組裝的重疊群序列(Fasta文件)進(jìn)行分析，結(jié)果表明該菌株基因組共有4 861個(gè)編碼蛋白質(zhì)的開(kāi)發(fā)閱讀框(ORFs)。進(jìn)一步對(duì)這些蛋白質(zhì)進(jìn)行COG 分類(lèi)分析，結(jié)果見(jiàn)圖1。該菌株基因組編碼的蛋白質(zhì)中最多的是糖轉(zhuǎn)運(yùn)與代謝蛋白質(zhì)(Carbohydrate Transport and Metabolism)，467個(gè)ORFs；其次是氨基酸轉(zhuǎn)運(yùn)與代謝蛋白質(zhì)(Amino Acid Transport and Metabolism)，423個(gè) ORFs及一般功能預(yù)測(cè)蛋白質(zhì)(General Function Prediction Only)，351個(gè) ORFs。通過(guò)對(duì)基因組的組裝和分析，可以得到該菌株基因組的整體特征，但由于基因組的復(fù)雜性和二代測(cè)序數(shù)據(jù)讀長(zhǎng)的限制，無(wú)法得到該菌株基因組的精細(xì)圖，那些沒(méi)有覆蓋的序列(gap)大部分是基因組中的重復(fù)序列，不屬于蛋白質(zhì)編碼區(qū)和基因表達(dá)調(diào)控區(qū)。

圖1 預(yù)測(cè)蛋白質(zhì)的COG 分類(lèi)圖

2.2抗生素抗性基因分析結(jié)果 抗生素抗性基因數(shù)據(jù)庫(kù)ARDB預(yù)測(cè)出該菌株基因組中含有26 種抗性基因，見(jiàn)表1。該菌株基因組編碼的蛋白質(zhì)GHCCNJDC_03951 屬于mdtg 抗性基因類(lèi)型，GHCCNJDC_03939 屬于mdth 抗性基因類(lèi)型，這和該菌株對(duì)磷霉素有抗性的實(shí)驗(yàn)結(jié)果一致；GHCCNJDC_04736屬于bl2be_ctxm 類(lèi)型，這一預(yù)測(cè)結(jié)果也和該菌株對(duì)氨芐青霉素有抗性的實(shí)驗(yàn)結(jié)果一致。此外還有屬于acra 抗性基因類(lèi)型的GHCCNJDC_00668，屬于acrb抗性基因類(lèi)型的GHCCNJDC_00669等。在這些預(yù)測(cè)的抗生素抗性基因中，除了ACD 類(lèi)β-內(nèi)酰胺酶外，絕大多數(shù)屬于多藥耐藥外排泵 (Multidrug Resistance Efflux Pump)。

表1 抗生素抗性基因預(yù)測(cè)結(jié)果

續(xù)表1 抗生素抗性基因預(yù)測(cè)結(jié)果

2.3菌株的鑒定和比較基因組學(xué)研究 通過(guò)從頭組裝，雖然得到了該菌株基因組的草圖，但仍然無(wú)法回答該菌株是從什么菌株進(jìn)化而來(lái)，其基因組發(fā)生了哪些突變從而產(chǎn)生耐藥性。首先，用分子分型軟件MLST對(duì)該菌株的重疊群序列進(jìn)行分析，結(jié)果顯示該菌株為ST131。MLST軟件是根據(jù)七個(gè)保守的基因序列(adk，fumC，gyrB，icd，mdh，purA及recA)把大腸桿菌分成不同的序列類(lèi)型(ST)，相同ST的菌株之間基因組差異小，親緣關(guān)系近。 由于大腸桿菌不同菌株基因組的差異，如果不進(jìn)行分子分型找到參考基因組而直接用大腸桿菌MG1655菌株作為參考基因組，將會(huì)找到幾十萬(wàn)乃至上百萬(wàn)的突變位點(diǎn)，而這些突變位點(diǎn)絕大多數(shù)都是由于錯(cuò)誤的比對(duì)導(dǎo)致的。

通過(guò)分子分型，筆者知道了該菌株的進(jìn)化地位，并找到了進(jìn)行比較基因組研究的參考基因組。接下來(lái)，把該菌株Reads Mapping到大腸桿菌T131菌株EC958 基因組上，從而找到該菌株和參考基因組不一致的突變序列。BWA 回帖結(jié)果顯示該菌株和大腸桿菌EC958 菌株基因組有很高的同源性，平均測(cè)序深度(Depth)為97倍,覆蓋度幾乎百分之百(僅僅100個(gè)位點(diǎn)的測(cè)序深度為零)。突變檢出(Variation Calling)共找到1 562個(gè)突變位點(diǎn)，包括1 536個(gè)高可信度SNPs和 26個(gè)插入缺失突變(INDELs)。

2.4突變位點(diǎn)snpEff 注釋結(jié)果 為了進(jìn)一步分析這些突變，采用snpEff 軟件結(jié)合參考基因組的注釋數(shù)據(jù)對(duì)突變數(shù)據(jù)(VCF) 進(jìn)行注釋。注釋結(jié)果如表2所示。 絕大多數(shù)突變屬于SNP(1 536個(gè))，少數(shù)屬于INDEL(26個(gè))。SNP突變中位于基因編碼區(qū)的有1 187個(gè)，位于基因表達(dá)調(diào)控區(qū)349個(gè)。位于編碼區(qū)的SNP突變中，859個(gè)是不改變蛋白質(zhì)序列的同義突變(synonymous mutation)，錯(cuò)義突變(missense mutation)有326個(gè)， 此外還有2個(gè)終止突變(stop_gained mutation)。兩個(gè)終止突變將使相應(yīng)的蛋白質(zhì)翻譯提前終止：OmpF 蛋白(WP_000977908.1)基因編碼區(qū)第250位的鳥(niǎo)嘌呤轉(zhuǎn)換成胸腺嘧啶 (c.250G>T) 從而導(dǎo)致該蛋白84位的谷氨酸變成終止位點(diǎn)(p.Glu84*)；Ves蛋白(WP_000455604.1)基因編碼區(qū)的突變(c.460C>T)，導(dǎo)致終止密碼子的產(chǎn)生(p.Gln154*)， 從而使該蛋白的翻譯提前終止。大多數(shù)插入刪除(INDELs)突變都會(huì)導(dǎo)致讀碼框的改變，以致不能產(chǎn)生正常的蛋白質(zhì)。深入分析這些突變位點(diǎn)將有可能發(fā)現(xiàn)導(dǎo)致抗生素抗性的基因和突變。

表2 snpEff注釋結(jié)果

3 討 論

隨著第二代高通量測(cè)序技術(shù)的發(fā)展，數(shù)據(jù)庫(kù)中各種生物序列數(shù)據(jù)呈爆炸式增長(zhǎng)，這些生物大數(shù)據(jù)中蘊(yùn)含豐富的信息，怎么去分析挖掘這些信息是目前面臨的急需解決的問(wèn)題。本研究通過(guò)對(duì)一株臨床多藥抗性菌株的高通量二代測(cè)序原始數(shù)據(jù)進(jìn)行充分挖掘，最終找到可能的抗藥性突變位點(diǎn)，為后續(xù)的研究提供了重要的線(xiàn)索。

通過(guò)對(duì)該菌株基因組測(cè)序數(shù)據(jù)的處理和從頭組裝，我們得到了該菌株的基因組草圖。由于基因組中很難組裝的部分都位于高度重復(fù)區(qū)域，這些區(qū)域一般不是蛋白編碼基因或表達(dá)調(diào)控序列，所以根據(jù)基因組草圖，我們能預(yù)測(cè)出ORFs并進(jìn)行多藥抗性基因分析、COG分析等。再根據(jù)該菌株的基因組草圖進(jìn)行分子分型，從而確定了該菌株的進(jìn)化地位及參考基因組。通過(guò)把測(cè)序片段回帖到參考基因組、突變檢出、突變注釋等進(jìn)行一系列的生物信息學(xué)分析，最終得到了該菌株和參考基因組不一致的所有突變位點(diǎn)，這些突變將為研究該菌株的耐藥機(jī)制提供重要的線(xiàn)索。

隨著越來(lái)越多的細(xì)菌全基因組精細(xì)圖譜的產(chǎn)生，分子分型將越來(lái)越精準(zhǔn)，我們將更準(zhǔn)確地確定菌株的進(jìn)化地位，檢出突變，為實(shí)驗(yàn)科學(xué)提供更有價(jià)值的線(xiàn)索。隨著測(cè)序數(shù)據(jù)的爆炸式增長(zhǎng)，生物信息學(xué)家需要開(kāi)發(fā)出能大規(guī)模分析處理二代測(cè)序原始數(shù)據(jù)的優(yōu)秀挖掘算法和分析流程，這將使我們能從臨床致病菌的測(cè)序數(shù)據(jù)中挖掘出更多更有價(jià)值的信息。