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

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

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

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

高通量二代測序大大促進了生物信息學的發(fā)展，生物信息學家已經(jīng)開發(fā)出了各種各樣的軟件工具對測序數(shù)據(jù)進行處理，但每一種軟件都只解決數(shù)據(jù)處理某一個環(huán)節(jié)的一個具體問題〔3〕。怎么分析研究大量的測序原始數(shù)據(jù)，進而挖掘出有用的生物學信息，對于生物信息學家來說依然是一個巨大的挑戰(zhàn)。而且如果能利用好現(xiàn)有的生物信息學工具對測序原始數(shù)據(jù)進行充分挖掘?qū)閷嶒炑芯刻峁┲匾木€索和方向，必將會大大加快研究的步伐。

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

1 材料和方法

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

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

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

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

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

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

2 結 果

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

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

圖1 預測蛋白質(zhì)的COG 分類圖

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

表1 抗生素抗性基因預測結果

續(xù)表1 抗生素抗性基因預測結果

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

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

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

表2 snpEff注釋結果

3 討 論

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

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

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