劉碩，曾志，曾凡才，杜萌澤

研究報告

基于序列相似性和Z曲線方法重注釋原核生物蛋白編碼基因

劉碩1，曾志1，曾凡才2，杜萌澤2

1. 電子科技大學(xué)生命科學(xué)與技術(shù)學(xué)院，成都 611731 2. 西南醫(yī)科大學(xué)基礎(chǔ)醫(yī)學(xué)院，分子生物與生物化學(xué)教研室，瀘州 646000

隨著測序技術(shù)的不斷發(fā)展，產(chǎn)生了海量的基因組測序數(shù)據(jù)，極大地豐富了公共遺傳數(shù)據(jù)資源。同時為了應(yīng)對大量基因組數(shù)據(jù)的產(chǎn)生，基因組比較和注釋算法、工具不斷更新，使得聯(lián)合多種注釋工具得到更準(zhǔn)確的蛋白編碼基因的注釋信息成為可能。目前公共數(shù)據(jù)庫的原核生物基因組測序和裝配有些是10多年前的，存在大量預(yù)測的功能未知的編碼基因。為了提升美國國家生物信息中心(National Center for Biotechnology Information, NCBI)數(shù)據(jù)庫中基因組的注釋質(zhì)量，本研究聯(lián)合使用多種原核基因識別算法/軟件和基因表達(dá)數(shù)據(jù)重注釋1587個細(xì)菌和古細(xì)菌基因組。首先，利用Z曲線的33個變量從177個基因組原注釋中識別獲得3092個被過度注釋為蛋白編碼基因的序列；其次，通過同源比對為939個基因組中的4447個功能未知的蛋白編碼基因注釋上具體功能；最后，通過聯(lián)合采用ZCURVE 3.0和Glimmer 3.02以及Prodigal這3種高精度的、廣泛使用且基于算法不同而互補(bǔ)的基因識別軟件來尋找漏注釋基因。最終，從9個基因組中找到了2003個被漏注釋的蛋白編碼基因，這些基因?qū)儆诙鄠€蛋白質(zhì)直系同源簇(clusters of orthologous groups of proteins, COG)。本研究使用新的工具并結(jié)合多組學(xué)數(shù)據(jù)重新注釋早期測序的細(xì)菌和古細(xì)菌基因組，不僅為新測序菌株提供注釋方法參考，而且這些重注釋后得到的細(xì)菌基因序列也會對后續(xù)基礎(chǔ)研究有所幫助。

細(xì)菌；重注釋；Z曲線；假定ORFs；非蛋白編碼ORFs

基因組分析是生命科學(xué)研究中非常關(guān)鍵的基礎(chǔ)工作，只有獲得準(zhǔn)確的基因組注釋才能為后續(xù)的分析和研究提供有力支撐并得到有意義的結(jié)果。目前，美國國家生物信息中心(National Center for Biotech-nology Information, NCBI)積累了大量10多年前完成測序的基因組。受限于早期有限的基因注釋工具[1]，這些基因組注釋的精確度有待于進(jìn)一步提升。例如，在原核生物基因組中尋找小的基因時容易出現(xiàn)此類基因被漏注釋而丟失的問題[2]。隨著測序效率的提升，幾何級數(shù)增長的測序量意味著這些注釋錯誤將會通過同源搜索等方式被放大[3]。Breitwieser等[4]檢查了所有的細(xì)菌和古細(xì)菌基因組后發(fā)現(xiàn)2250個原核生物基因組居然包含了人類基因組的信息。為了確?；蚪M信息的準(zhǔn)確性[3,5]，相關(guān)數(shù)據(jù)庫需要時常更新原始的基因組注釋信息。當(dāng)前隨著各種組學(xué)數(shù)據(jù)的日益完善，通過聯(lián)合多工具多組學(xué)的方式對基因組精確重注釋也具備可行性。

自1995年第1個原核生物流感嗜血桿菌()基因組被注釋以來，大量的原核生物基因組陸續(xù)完成測序和注釋[5]。根據(jù)GenBank[6]的記錄，自1999年4月至2019年10月，被測序序列的數(shù)量從3,525,418增長至216,763,706，增長倍數(shù)約為60倍。近年來基因識別軟件，尤其是識別原核生物中基因的新方法和軟件在被陸續(xù)開發(fā)和不斷升級[7~14]，其中較廣泛使用的有IPred[9]、Gene-MarkS[11]、Glimmer[12]、EasyGene[13]和ZCURVE[14]等。這些新的基因探測工具不僅可以用于注釋基因組從而獲得更精確的結(jié)果，而且還可賦予假定蛋白詳盡的功能[15~19]。這些工具中基于Z-curve理論的軟件ZCURVE 3.0[8]具有93.7%的準(zhǔn)確率，與具有93.0%準(zhǔn)確率的Glimmer 3.02[20]相當(dāng)，而且Z-curve通過把單核苷酸、雙核苷酸和三核苷酸的堿基組成轉(zhuǎn)換成33個空間變量來表示DNA序列的特征，并通過對正樣本和負(fù)樣本的學(xué)習(xí)來對未知序列進(jìn)行編碼蛋白能力的預(yù)測[15,16]。該方法在原核生物的編碼蛋白的預(yù)測上準(zhǔn)確度高、應(yīng)用廣[21]。這兩個軟件聯(lián)合使用可以獲得互相補(bǔ)充的結(jié)果。更準(zhǔn)確的細(xì)菌和古細(xì)菌基因組信息有助于更精確地推導(dǎo)生命的起源和演化，如Weiss等[22]對原核生物基因組的610萬個蛋白編碼基因的全部簇和進(jìn)化樹進(jìn)行研究以推斷最后的共同祖先(the last universal common ancestor, LUCA)，而不準(zhǔn)確的注釋會影響該研究結(jié)論的可靠性。

本研究對1587個測序和注釋較早的細(xì)菌和古細(xì)菌菌株的基因組進(jìn)行了重新注釋，這些基因組包含功能已知的開放閱讀框(open reading frame, ORF)和假定ORFs (功能未知的ORFs)，而假定ORFs中有些實(shí)際上是非蛋白編碼基因的序列，此外基因組還包含被漏注釋掉但實(shí)際上是可編碼蛋白基因的那些序列。參考基因表達(dá)數(shù)據(jù)并聯(lián)合3個注釋準(zhǔn)確率較高的原核生物識別軟件(ZCURVE 3.0[8]、Glimmer 3.02[20]以及Prodigal[23])，移除掉一些之前被過度注釋為基因的序列，同時給功能未知的ORFs注釋了新的功能，并獲得了之前漏注釋的基因。

1 材料與方法

1.1 數(shù)據(jù)的選取

選用來自30個門(phylum)1587個20年內(nèi)被陸續(xù)測序的細(xì)菌和古細(xì)菌的基因組來進(jìn)行重注釋(附表1)，基因組來自992個物種(附表2)。基因組的裝配序列和注釋信息均于2019年1月獲取自NCBI數(shù)據(jù)庫，它們的測序時間為1999年~2019年。

在使用33個Z-curve變量去除1587個基因組的過度注釋基因時，選取它們的一些具有已知功能的ORFs作為正樣本，并選取每個基因組的此類ORFs隨機(jī)打亂1000次后產(chǎn)生的完全隨機(jī)的序列作為負(fù)樣本。

為判斷計算識別的序列為真正的非蛋白編碼序列，從GEO[24]和paxdb[25]下載了蛋白豐度或者RNA豐度數(shù)據(jù)。原基因組中第二類功能不確定的ORFs即假定ORFs如果被Z-curve 33變量方法識別為可能的非蛋白編碼序列，同時其沒有對應(yīng)的蛋白或者mRNA被檢測到，即認(rèn)定其為過度注釋為基因的序列，繼而從基因組基因列表中移除。

本研究對1587個基因組中的9個模式物種基因組進(jìn)行新基因的注釋。它們分別是嗜酸氧化亞鐵硫桿菌(ATCC 23270)、炭疽芽孢桿菌(str. Ames)、枯草芽孢桿菌(subsp. subtilis str. 168)、大腸桿菌(str. K-12 substr. MG1655)、流感嗜血桿菌(Rd KW20)、腦膜炎奈瑟球菌(MC58)、腸道沙門氏菌(subsp. enterica serovar Typhi str. CT18)、金黃色葡萄球菌(subsp. aureus NCTC 8325)和釀膿鏈球菌(SF370)。

1.2 重注釋流程

整個重注釋的流程分為3個部分：(1)識別過度注釋為基因的序列；(2)未知功能ORFs的功能注釋；(3)識別欠注釋的基因(圖1)，其具體的方法學(xué)細(xì)節(jié)部分見1.3、1.4、1.5和1.6。

1.3 尋找過度注釋為基因的序列

Z-curve把DNA序列中出現(xiàn)在第一(1、4、7...)、第二(2、5、8...)和第三(3、6、9...)密碼子位的A、G、C、T的堿基的頻率分別用a1、g1、c1、t1；a2、g2、c2、t2；a3、g3、c3、t3來表示，根據(jù)Z曲線理論，一條DNA序列可以用三維平面中的點(diǎn)P(=1、2、3)來表示，而1、2、3所對應(yīng)的x、y、z可以通過下邊的公式來計算得出，據(jù)此得出Z-curve的單核苷酸的9個變量(9=3×3)。

圖1 重注釋流程圖

當(dāng)用12()或23()來表示DNA序列的第一第二密碼子位或第二第三密碼子位時，上邊的公式可以衍生出以下公式，其中代表A、G、C或者T中的一種堿基。表示密碼子相位組合，=12表示第一第二相位，同理=23表示第二第三相位，可以得到24個變量(24=3×4×2)。

得到以上Z-curve 33變量后，文章使用每個基因組對應(yīng)的正樣本和負(fù)樣本進(jìn)行訓(xùn)練。從正樣本和負(fù)樣本中隨機(jī)抽取1/10的樣本作為測試集，并進(jìn)行十次交叉驗(yàn)證。最終平均預(yù)測準(zhǔn)確率均大于99%，證明該預(yù)測方法穩(wěn)定可靠。

將沒有提供明確功能的ORFs作為測試集可尋找過度注釋為基因的序列。預(yù)測為非蛋白編碼基因的序列且同時沒有在公共數(shù)據(jù)庫中找到表達(dá)數(shù)據(jù)的假定ORFs將被歸為過度注釋為基因的序列?；虻谋磉_(dá)數(shù)據(jù)來自不同的RNA測序集和GEO數(shù)據(jù)集(附表3)，這些數(shù)據(jù)可以幫助排除掉具有表達(dá)數(shù)據(jù)的基因。

1.4 給功能未知ORFs注釋功能

假如一些假定/功能未知的ORFs沒有被識別為非蛋白編碼序列而是被預(yù)測為蛋白編碼基因，研究就會通過同源搜索賦予其功能，使用的工具為BLAST[26]，賦予基因功能的同源搜索的參數(shù)設(shè)置為值(-value) < 1e-20，覆蓋率(Coverage) > 80%以及一致性(Identity) > 70%。而假如被預(yù)測為非蛋白編碼基因的序列其RNA-seq的數(shù)據(jù)顯示為Count/TPM/ RPKM/CPM > 0，即存在，可以認(rèn)為該序列具有表達(dá)數(shù)據(jù)，則需通過BLAST來注釋其功能，根據(jù)經(jīng)驗(yàn)參數(shù)，把賦予基因功能的同源搜索的參數(shù)同樣設(shè)置為值 < 1e-20，覆蓋率> 80%以及一致性 > 70%。

1.5 補(bǔ)充漏注釋基因

基因組注釋經(jīng)常會丟掉一些真正的基因[2,27]，使用多種基因搜索工具可以盡可能減少該錯誤的發(fā)生。此處聯(lián)合ZCURVE 3.0[8]、Glimmer 3.02[20]和Prodigal[23]來注釋9種重要的模式生物的代表株。這3個軟件識別出的不在原基因組基因列表中的ORFs將被全部加入候選基因序列集。對此集合中的序列用BLAST[26](https://blast.ncbi.nlm.nih.gov/Blast. cgi)進(jìn)行同源比對后只保留那些與公共數(shù)據(jù)庫中功能已知的基因有顯著相似性的候選基因序列。由于此處目的為注釋被遺漏的基因，研究根據(jù)經(jīng)驗(yàn)參數(shù)把同源比對的參數(shù)設(shè)置為比1.4部分的閾值稍寬松的值，即< 1e-20，覆蓋率 > 60%，一致性 > 60%，并用2個閾值對模式生物大腸桿菌的代表株str. K-12 substr. MG1655被注釋后得到的246個新基因進(jìn)行了2次同源比對，通過對比，發(fā)現(xiàn)寬松的閾值更能確保注釋得到的新基因的完整性，詳細(xì)結(jié)果參見2.4。

1.6 同源簇注釋

結(jié)構(gòu)和功能相似的蛋白編碼序列屬于同一個蛋白簇。在此，為注釋得到的新基因進(jìn)行同源簇COG注釋可以幫助更好去理解基因組如何正常行使功能。研究同時聯(lián)合使用了EggNOG[28]、WebMGA[29]和NCBI中的COG數(shù)據(jù)庫[30]來注釋這些基因序列的同源簇，其對應(yīng)的COG編號可在附表7中找到。EggNOG包含2031個物種和19萬個同源簇和對應(yīng)功能注釋數(shù)據(jù)。序列分析器WebMGA[29]和NCBI上的COG數(shù)據(jù)庫[30]用于進(jìn)一步完善同源簇功能注釋。

2 結(jié)果與分析

2.1 1587個物種的注釋現(xiàn)狀

根據(jù)這1587個物種的具體測序時間，繪制了同一年份被測序的基因組數(shù)目占全部基因組比例的柱形圖(圖2，附表1)，其中2012年之前被測序的基因組占整體基因組的97%以上。根據(jù)CVTree[31]列出的門(phylum)的數(shù)量，與細(xì)菌和古細(xì)菌相關(guān)的門(phylum)為41個，選取的細(xì)菌和古細(xì)菌的基因組在門(phylum)水平的覆蓋度為73%，研究選取的基因組盡可能包含了模式微生物的代表菌株如大腸桿菌和枯草芽孢桿菌等屬于同一物種的注釋時間不同的多種血清型，整體來說，選取的基因組在Refseq數(shù)據(jù)庫中裝配和注釋的時間較早。

圖2 1587個基因組的測序時間分布

2.2 原注釋中3092個假定ORFs為非蛋白編碼序列

在這1587個基因組的原始注釋中，平均有1.26%的假定基因重注釋后信息發(fā)生改變(圖3A)。使用功能已知的ORFs作為正樣本，選用對正樣本序列隨機(jī)打亂1000次后產(chǎn)生的序列作為負(fù)樣本，進(jìn)行訓(xùn)練后通過十重交叉驗(yàn)證后得到預(yù)測序列編碼能力的平均準(zhǔn)確率為0.9985 (圖3B)。接著，研究依次對1587個細(xì)菌基因組中的功能未知ORFs進(jìn)行預(yù)測，其中56,462個ORFs被預(yù)測為可能的非蛋白編碼序列(附表4)。

由于計算方法本身的偏差，本研究使用實(shí)驗(yàn)數(shù)據(jù)進(jìn)一步確定這些序列是否為非蛋白編碼序列。當(dāng)且僅當(dāng)這些序列在公共數(shù)據(jù)庫中查詢不到表達(dá)數(shù)據(jù)，才認(rèn)為它們確實(shí)不具有編碼功能。研究查詢了GEO等數(shù)據(jù)庫并最終確定177個基因組中共3092個ORFs為過度注釋為基因的序列(附表5)。需要重點(diǎn)提及的是其中43個基因組有20個以上的假定ORFs被識別為這類序列(表1)。

圖3 假定ORFs的比率及不同準(zhǔn)確率對應(yīng)的基因組數(shù)量及大豆根瘤菌(B. japonicum USDA 110)3類序列兩個密碼子位GC含量

A：不同比率的重注釋假定ORFs對應(yīng)基因組的數(shù)目；B：不同識別非編碼ORFs的準(zhǔn)確率對應(yīng)基因組的數(shù)目；C：大豆根瘤菌3類序列對應(yīng)的第2和第3位密碼子GC含量。

之前的研究已經(jīng)報道過編碼的ORFs在第二和第三密碼子位的位置上有著不同的GC含量分布[16]，即大多數(shù)功能已知的ORFs的第三密碼子位的GC含量比其第二密碼子位的GC含量都高，而對于非編碼ORFs，第二密碼子位的GC含量則是接近第三密碼子位的GC含量。在43個基因組中，大豆根瘤菌(USDA 110)基因組的假定ORFs被識別到了最多的非蛋白編碼序列(147條)，其第二密碼子位和第三密碼子位的GC含量的分布和以上的論斷是類似的(圖3C)，這驗(yàn)證了本研究鑒定出的非編碼ORFs具有可信性。

表1 識別出過度注釋的ORFs多于20的菌株基因組的信息

2.3 4447個功能未知ORFs被注釋上確定的功能

本研究用同源比對的方法為那些功能未知的ORFs注釋上功能。最終，939個基因組中共計4447個ORFs被注釋了確切功能(附表6)。其中在33個基因組中，有超過20個假定ORFs被注釋上具體功能(表2)。

這些新注釋的功能對生命活動非常重要。例如本研究發(fā)現(xiàn)1587個基因組中的601個功能未知的ORFs與一些質(zhì)膜蛋白基因序列相似。質(zhì)膜蛋白控制細(xì)胞內(nèi)外物質(zhì)的交換和信息的傳遞，參與信號通路的調(diào)控[32]。其中沙眼衣原體(D/UW-3/CX)中共有5個功能未知的蛋白被注釋為質(zhì)膜蛋白，沙眼衣原體可以引發(fā)炎癥病變[33]，如它與子宮頸癌直接或間接相關(guān)。這些被注釋上新功能的假定ORFs的詳細(xì)信息請見附表6。

2.4 對新基因選取寬松閾值和嚴(yán)格閾值進(jìn)行同源比對的結(jié)果

鑒于1.4和1.5部分進(jìn)行同源比對的閾值不同，研究對模式生物大腸桿菌的代表株str. K-12 substr. MG1655注釋后得到的246個新基因進(jìn)行了第二次同源比對，選擇的閾值為值 < 1e-20，覆蓋率 > 80%以及一致性 > 70%，對兩次比對做比較后發(fā)現(xiàn)有9個基因滿足寬松的閾值(值 < 1e-20，覆蓋率 > 60%，一致性 > 60%)，而當(dāng)設(shè)置嚴(yán)格的閾值(值 < 1e-20，覆蓋率 > 80%以及一致性 > 70%)時會被丟失掉。這9個基因的功能和對應(yīng)值在表3中被列出，這些功能沒有出現(xiàn)在str. K-12 substr. MG1655的現(xiàn)有注釋中，因此認(rèn)為它們可能對菌株發(fā)揮功能起著不可或缺的作用，而這恰恰是對該菌株注釋新基因的意義所在。

表2 含有20個以上的假定ORFs被注釋上準(zhǔn)確功能的基因組的信息

2.5 9個基因組中新識別出2003個新基因

ZCURVE 3.0具有93.7%的準(zhǔn)確率[8]，Glimmer 3.02具有93.0%的準(zhǔn)確率[20]，這兩個軟件和假陽性率較低的Prodigal[23]被一起用來識別漏注釋的基因，預(yù)測出的新的編碼ORFs形成的并集被用來同源比對。最終，在9個基因組中識別到了2003個新基因(表4，附表7)。在最初的注釋中基因數(shù)量相對較少的釀膿鏈球菌和流感嗜血桿菌分別得到了104和123個新基因，新基因的數(shù)量大約是原始基因數(shù)目的6%。腦膜炎奈瑟球菌和腸道沙門氏菌獲得了很多新注釋到的基因，它們占據(jù)了其原始基因數(shù)目的14%。這些被新識別的基因的名字和其起始位置的信息存放在附表7中。

注釋同源蛋白簇對研究具有相似功能的基因幫助頗多。通過綜合EggNOG[28]、WebMGA[29]和COG數(shù)據(jù)庫[30]，共有1073個新識別的基因被注釋到對應(yīng)的COGs(圖4)，其豐度就是某蛋白質(zhì)直系同源簇所包含的新識別的基因占全部新識別的基因的比例。這1073條序列的詳細(xì)COG注釋可以在附表7中查詢到，而所有1587個基因組重注釋的具體信息則可以在附表8中被查到，該表列出的具體信息包含其基本信息即NC_ID和基因組大小(bp)，還包含其原始注釋信息即蛋白質(zhì)數(shù)目、功能已知的基因和假定基因，并列出了通過重注釋得到的注釋信息即假定ORFs中被注釋上新功能的ORFs數(shù)目和非編碼ORFs的數(shù)目，還提供了9個基因組新基因的數(shù)量，并提供了通過Z-curve 33變量方法識別非編碼ORFs的準(zhǔn)確率。

3 討論

本研究結(jié)合3種原核基因組注釋軟件和表達(dá)數(shù)據(jù)重注釋了1587個細(xì)菌和古細(xì)菌的基因組。首先，識別了過度注釋為基因的序列，隸屬于177個基因組的3092個假定ORFs被識別為非蛋白編碼序列；接下來，還賦予了假定ORFs以功能，939個基因組的4447個基因被通過相似性搜索而注釋獲得準(zhǔn)確功能，并且在9個基因組中識別了2003個屬于多個蛋白質(zhì)直系同源簇的新基因。

重注釋可以保證進(jìn)一步分析的數(shù)據(jù)的準(zhǔn)確性。在真核生物中有研究表明多拷貝的基因會出現(xiàn)功能的分化[34]，原核生物也存在多拷貝的基因，它們是否也存在這種分化，可以通過開展重注釋后獲得準(zhǔn)確的信息來進(jìn)一步研究。而正如前文所提及的那樣，伴隨著測序數(shù)據(jù)的增加和技術(shù)的不斷進(jìn)步，對已測序菌株的基因組進(jìn)行重新注釋的必要性也在增大，同時根據(jù)Breitwieser的研究[4]，微生物基因組可能會被人類基因組所污染而導(dǎo)致基因的丟失，這就更加要求研究者們用多種方法來更新注釋。本次進(jìn)行重新注釋的1587個菌株可歸為992個物種，每個物種涵蓋多種血清型，例如沙門氏菌根據(jù)抗原的不同可以分為不同的血清型[35]，其基因組層面可能會存在差異。Liu等[33]通過全基因組序列對6個沙門氏菌(Salmonella)菌株建立進(jìn)化樹顯示亞種的相同血清型的C和A進(jìn)化距離較遠(yuǎn)，卻與不同血清型的進(jìn)化距離較近。研究通過CVTree[31]對15個亞種構(gòu)建進(jìn)化樹發(fā)現(xiàn)血清型之間會呈現(xiàn)不同的聚集，與Liu等[33]繪制的進(jìn)化樹一致(附圖1)。因此，通過重注釋不同血清型的菌株來使得基因組信息更完善，從而幫助深入探索血清型之間更準(zhǔn)確的演化關(guān)系。而研究參考的數(shù)據(jù)集為1587個基因組中已經(jīng)被注釋為具有明確功能的ORFs，為了驗(yàn)證這些基因具有表達(dá)數(shù)據(jù)，針對大腸桿菌的功能已知的基因進(jìn)行驗(yàn)證，研究下載了對應(yīng)該物種的表達(dá)數(shù)據(jù)并尋找這些功能已知的基因中有多少具有表達(dá)數(shù)據(jù)，通過與GEO數(shù)據(jù)庫中GSE56133[36]和GSE118058[37]的表達(dá)數(shù)據(jù)比較，發(fā)現(xiàn)2835個功能已知的基因中有2806個在這兩個數(shù)據(jù)集中有表達(dá)數(shù)據(jù)，這進(jìn)一步證實(shí)了現(xiàn)有明確功能注釋的基因可以作為參考集來進(jìn)行預(yù)測。

表3 大腸桿菌(E.coli str. K-12 substr. MG1655)滿足寬松閾值的新基因

表4 9個菌株的名稱、NC序列號、基因組大小、基因總數(shù)和新注釋基因的數(shù)目

圖4 特定同源簇對應(yīng)的新基因占全部新基因的比例

A：RNA加工和修飾；B：染色質(zhì)結(jié)構(gòu)和動力學(xué)；Y：核結(jié)構(gòu)；Z：細(xì)胞骨架；W：細(xì)胞外結(jié)構(gòu)；D：細(xì)胞周期控制，細(xì)胞分裂，染色體分裂；O：翻譯后修飾，蛋白反轉(zhuǎn)，伴侶；Q：次生代謝產(chǎn)物的生物合成、運(yùn)輸和分解代謝；I：脂質(zhì)轉(zhuǎn)運(yùn)與代謝；J：翻譯，核糖體結(jié)構(gòu)和生物發(fā)生；H：輔酶運(yùn)輸與代謝；F：核苷酸轉(zhuǎn)運(yùn)與代謝；V：防衛(wèi)機(jī)制；N：細(xì)胞遷移；C：能量產(chǎn)生和轉(zhuǎn)化；U：細(xì)胞內(nèi)傳輸，分泌和囊泡轉(zhuǎn)運(yùn)；P：無機(jī)離子轉(zhuǎn)運(yùn)與代謝；T：信號轉(zhuǎn)導(dǎo)機(jī)制；K：轉(zhuǎn)錄；E：氨基酸轉(zhuǎn)運(yùn)和代謝；M：細(xì)胞壁和細(xì)胞膜的生物發(fā)生；G：碳水化合物運(yùn)輸和代謝；S：功能未知；R：只能預(yù)測大致功能；L：復(fù)制重組和修復(fù)。

得益于眾多基因識別軟件的升級如ZCURVE 3.0[8]和不斷更新的公共數(shù)據(jù)庫如NCBI的GEO[24]，使得尋找新的編碼蛋白ORFs和用表達(dá)數(shù)據(jù)驗(yàn)證成為可能。在識別新的基因方面，本研究聯(lián)合ZCURVE 3.0[8]、Glimmer 3.02[20]和Prodigal[23]對原核基因組重新注釋，其中Glimmer 3.02是基于隱馬爾可夫模型，識別準(zhǔn)確率高達(dá)93.0%，但是它對序列的局部特征有著強(qiáng)烈的依賴性，而Prodigal是通過對現(xiàn)有的細(xì)菌和古細(xì)菌基因組注釋信息進(jìn)行訓(xùn)練，因此對已知基因/保守基因的識別效果優(yōu)于其余的軟件，但預(yù)測未報道的基因時會有些許偏差。而ZCURVE 3.0是基于DNA序列的全局統(tǒng)計特征，它將Fisher線性判別替換為支持向量機(jī)(SVM)[38]以提高靈敏度，并且鑒于在預(yù)測過程中由于負(fù)樣本的隨機(jī)性容易產(chǎn)生假陽性的基因，ZCURVE 3.0依據(jù)了ORFs核酸分布的花瓣模型[39]在訓(xùn)練集中產(chǎn)生5類負(fù)樣本并逐次分類，然后保留在多次分類中均被預(yù)測為基因的那些ORFs，并且算法還把33個包含零階和一階的Z曲線變量增加為額外包含二階和三階Z曲線變量的765個變量，可以最大化地優(yōu)化程序的預(yù)測效果。選取三個程序混合預(yù)測更能保證預(yù)測的準(zhǔn)確率。

總之，本研究基于序列的全局特征和局部特征，結(jié)合同源比對對多個物種的多種菌株的基因組進(jìn)行了全面的重新注釋。未來，伴隨著基因組多組學(xué)數(shù)據(jù)的增加以及相應(yīng)的注釋工具和數(shù)據(jù)庫的完善[40~43]，公共數(shù)據(jù)庫里基因組的注釋將會更加準(zhǔn)確。

致謝

感謝電子科技大學(xué)生命科學(xué)與技術(shù)學(xué)院的郭鋒彪教授在研究開展中給予的幫助。

附錄：

附圖和附表詳見文章電子版www.chinagene.cn。

[1] M?rk S, Holmes I. Evaluating bacterial gene-finding HMM structures as probabilistic logic programs., 2012, 28(5): 636–642.

[2] Warren AS, Archuleta J, Feng WC, Setubal JC. Missing genes in the annotation of prokaryotic genomes., 2010, 11(1): 131.

[3] Salzberg SL. Next-generation genome annotation: we still struggle to get it right., 2019, 20(1): 92.

[4] Breitwieser FP, Pertea M, Zimin AV, Salzberg SL. Human contamination in bacterial genomes has created thousands of spurious proteins., 2019, 29(6): 954–960.

[5] Fleischmann RD, Adams MD, White O, Clayton RA, Kirkness EF, Kerlavage AR, Bult CJ, Tomb JF, Dougherty BA, Merrick JM, Mckenney K, Sutton G, FitzHugh W, Fields C, Gocayne JD, Scott J, Shirley R, Liu LL, Glodek A, Kelley JM, Weidman JF, Phillips CA, Spriggs T, Hedblom E, Cotton MD, Utterback TR, Hanna MC, Nguyen DT, Saudek DM, Brandon RC, Fine LD, Fritchman JL, Fuhrmann JL, Geoghagen NSM, Gnehm CL, McDonald LA, Small KV, Fraser CM, Smith HO, Venter JC. Whole-genome random sequencing and assembly ofRd., 1995, 269(5223): 496–512.

[6] Benson DA, Karsch-Mizrachi I, Lipman DJ, Ostell J, Sayers EW. Genbank., 2009, 37(Data-base issue): D26–D31.

[7] Yu JF, Xiao K, Jiang DK, Guo J, Wang JH, Sun X. An Integrative method for identifying the over-annotated protein-coding genes in microbial genomes., 2011, 18(6): 435–449.

[8] Hua ZG, Lin Y, Yuan YZ, Yang DC, Wei W, Guo FB. ZCURVE 3.0: identify prokaryotic genes with higher accuracy as well as automatically and accurately select essential genes., 2015, 43(W1): W85– W90.

[9] Zickmann F, Renard BY. IPred-integrating ab initio and evidence based gene predictions to improve prediction accuracy., 2015, 16(1): 134.

[10] Keilwagen J, Wenk M, Erickson JL, Schattat MH, Grau J, Hartung F. Using intron position conservation for homology-based gene prediction., 2016, 44(9): e89.

[11] Besemer J, Lomsadze A, Borodovsky M. GeneMarkS:a self-training method for prediction of gene starts in microbial genomes. Implications for finding sequence motifs in regulatory regions., 2001, 29(12): 2607–2618.

[12] Kelley DR, Liu B, Delcher AL, Pop M, Salzberg SL. Gene prediction with Glimmer for metagenomic sequences augmented by classification and clustering., 2012, 40(1): e9.

[13] Larsen TS, Krogh A. EasyGene-a prokaryotic gene finder that ranks ORFs by statistical significance., 2003, 4(1): 21.

[14] Guo FB, Ou HY, Zhang CT. ZCURVE: a new system for recognizing protein-coding genes in bacterial and archaeal genomes., 2003, 31(6): 1780–1789.

[15] Du MZ, Guo FB, Chen YY. Gene re-annotation in genome of the extremophileby using bioinformatics methods., 2011, 29(2): 391–401.

[16] Guo FB, Xiong LF, Teng JL, Yuen KY, Lau SK, Woo PC. Re-annotation of protein-coding genes in 10 complete genomes of Neisseriaceae family by combining similarity- based and composition-based methods., 2013, 20(3): 273–286.

[17] Lei Y, Kang SK, Gao JX, Jia XS, Chen LL. Improved annotation of a plant pathogen genomepv. oryzae PXO99A., 2013, 31(3): 342–350.

[18] Mao Y, Yang X, Liu Y, Yan Y, Du Z, Han Y, Song Y, Zhou L, Cui Y, Yang R. Reannotation ofstrain 91001 Based on Omics Data., 2016, 95(3): 562–570.

[19] Pfeiffer F, Bagyan I, Alfaro‐Espinoza G, Zamora‐Lagos MA, Habermann B, Marin‐Sanguino A, Oesterhelt D, Kunte HJ. Revision and reannotation of theDSM 2581T genome., 2017, 6(4): e00465.

[20] Delcher AL, Bratke KA, Powers EC, Salzberg SL. Identifying bacterial genes and endosymbiont DNA with Glimmer., 2007, 23(6): 673–679.

[21] Zhang R, Zhang CT. A Brief Review:The Z-curve theory and its application in genome analysis., 2014, 15(2): 78–94.

[22] Weiss MC, Sousa FL, Mrnjavac N, Neukirchen S, Roettger M, Nelson-Sathi S, Martin WF. The physiology and habitat of the last universal common ancestor., 2016, 1(9): 16116.

[23] Hyatt D, Chen GL, Locascio PF, Land ML, Larimer FW, Hauser LJ. Prodigal: prokaryotic gene recognition and translation initiation site identification., 2010, 11(1): 119.

[24] Barrett TT, Troup DB, Wilhite SE, Ledoux P, Rudnev D, Evangelista C, Kim IF, Soboleva A, Tomashevsky M, Marshall KA, Phillippy KH, Sherman PM, Muertter RN, Edgar R. NCBI GEO: archive for high-throughput functional genomic data., 2009, 37(Database issue): D885–D890.

[25] Wang M, Weiss M, Simonovic M, Haertinger G, Schrimpf SP, Hengartner MO, von Mering C. PaxDb, a database of protein abundance averages across all three domains of life., 2012, 11(8): 492–500.

[26] McGinnis S, Madden TL. BLAST: at the core of a powerful and diverse set of sequence analysis tools., 2004, 32(Suppl.2): W20–W25.

[27] Wood DE, Lin H, Levy-Moonshine A, Swaminathan R, Chang YC, Anton BP, Osmani L, Steffen M, Kasif S, Salzberg SL. Thousands of missed genes found in bacterial genomes and their analysis with COMBREX., 2012, 7(1): 37.

[28] Huerta-Cepas J, Szklarczyk D, Forslund K, Cook H, Heller D, Walker MC, Rattei T, Mende DR, Sunagawa S, Kuhn M, Jensen LJ, Mering CV, Bork P. eggNOG 4.5: a hierarchical orthology framework with improved functional annotations for eukaryotic, prokaryotic and viral sequences., 2016, 44(Database issue): D286–D293.

[29] Wu ST, Zhu ZW, Fu LM, Niu BF, Li WZ. WebMGA: a customizable web server for fast metagenomic sequence analysis., 2011, 12(1): 444.

[30] Tatusov RL, Galperin MY, Natale DA, Koonin EV. The COG database: a tool for genome-scale analysis of protein functions and evolution., 2000, 28(1): 33–36.

[31] Qi J, Luo H, Hao BL. CVTree: a phylogenetic tree reconstruction tool based on whole genomes., 2004, 32: W45–W47.

[32] Hockenbery D, Nu?ez G, Milliman C, Schreiber RD, Korsmeyer SJ. Bcl-2 is an inner mitochondrial membrane protein that blocks programmed cell death., 1990, 348(6299): 334–336.

[33] Liu WQ, Feng Y, Wang Y, Zou QH, Chen F, Guo JT, Peng YH, Jin Y, Li YG, Hu SN, Johnson RN, Liu GR, Liu SL.C: Genetic divergence fromand pathogenic convergence with., 2009, 4(2): e4510.

[34] Vankuren NW, Long M. Gene duplicates resolving sexual conflict rapidly evolved essential gametogenesis functions., 2018, 2(4): 705–712.

[35] Minor LL, Bockemühl J. 1987 supplement (no.31) to the schema of Kauffmann-White., 1988, 139(3): 331–335.

[36] Dwyer DJ, Belenky PA, Yang JH, Macdonald IC, Martell JD, Takahashi N, Chan CT, lobritz MA, Braff D, Schwarz EG, Ye JD, Pati M, Vercruysse M, Ralifo PS, Allison KR, Khalil AS, Ting AY, Walker GC, Collins JJ. Antibiotics induce redox-related physiological alterations as part of their lethality., 2014, 111(20): E2100–E2109.

[37] Hadjeras L, Poljak L, Bouvier M, Morin-Ogier Q, Canal l, Cocaign-Bousquet M, Girbal L, Carpousis AJ. Detachment of the RNA degradosome from the inner membrane ofresults in a global slowdown of mRNA degradation, proteolysis of RNase E and increased turnover of ribosome-free transcripts., 2019, 111(6): 1715–1731.

[38] Kim S, Yu Z, Kil RM, Lee M. Deep learning of support vector machines with class probability output networks., 2015, 64: 19–28.

[39] Guo FB. The distribution patterns of bases of protein- coding genes, non-coding ORFs, and intergenic sequences inPA01 genome and its implication., 2007, 25(2): 127–133.

[40] Uyar B, Yusuf D, Wurmus R, Rajewsky N, Ohler U, Akalin A. RCAS: an RNA centric annotation system for transcriptome-wide regions of interest., 2017, 45(10): e91.

[41] Huang Y, Liu Q, Chi LJ, Shi CM, Wu Z, Hu M, Shi H, Chen H. Application of BIG-Annotator in the genome sequencing data functional annotation and genetic diagnosis., 2018, 40(11): 1015–1023.黃瑩, 劉琪, 池連江, 石承民, 吳禎, 胡敏, 石宏, 陳華. BIG-Annotator: 基因組測序數(shù)據(jù)高效功能注釋及其在遺傳診斷中的應(yīng)用. 遺傳, 2018, 40(11): 1015–1023.

[42] Bick JT, Zeng SQ, Robinson MD, Ulbrich SE, Bauersachs S. Mammalian Annotation Database for improved annotation and functional classification of Omics datasets from less well-annotated organisms., 2019, 2019: 1–16.

[43] Ravindran SP, Lüneburg J, Gottschlich L, Tams V, Cordellier M. Daphnia stressor database: Taking advantage of a decade of Daphnia ‘-omics’ data for gene annotation., 2019, 9(1): 11135.

附圖1亞種不同血清型的系統(tǒng)發(fā)育樹

Supplementary Fig. 1 The phylogenetic tree of different serological types ofsub-species

Comprehensive re-annotation of protein-coding genes for prokaryotic genomes by Z-curve and similarity-based methods

Shuo Liu1, Zhi Zeng1, Fancai Zeng2, Mengze Du2

The development of sequencing technology has generated huge genomicsequencing information and largely enriched public genetic resources. To analyze such big data, the algorithms and tools for comparison and annotation of genomes are updated continually, enabling genome annotation with higher accuracyvarious annotation tools. Many prokaryotic genomes in public database were sequenced and assembled more than a decade ago, and they contained multiple genes with unknown functions. To improve the current annotation for those genomes in NCBI, we re-annotate 1587 bacterial and archaeal genomes using multiple prokaryotic gene recognition algorithms/softwares and gene expression data.The 33 Z-curve variableswere appliedto recognize sequences that were over-annotated to genes of 1587 bacterial and archaeal genomes deposited in public databases, anda total of 3092sequences belonging to 177 genomes were recognized as sequences over-annotated as protein-coding genes.Next, 4447 protein-coding genes with unknown functions from 939 genomes were annotated with definite functions by similarity search. Finally, we recognized 2003 missed protein-coding genesthat belong to known COG (clusters of orthologous groups of proteins) of nine genomes using three methods (ZCURVE 3.0, Glimmer3.02 and Prodigal), which are accurate and frequently used for gene finding. Their algorithms are different and complementary. This is a comprehensive study for re-annotation of bacterial and archaeal genomes with new tools combining multi-omics data, which should provide a reference for annotation of newly sequenced strains, and also benefit further fundamental researches with the bacterial gene sequences obtained after re-annotation.

bacteria; re-annotation; Z-curve; hypothetical ORFs; non-coding ORFs

2020-02-20;

2020-05-11

電子科技大學(xué)理科實(shí)力提升計劃項(xiàng)目(編號：Y0301902610100202)資助[Supported by Science Strength Improvement Plan of University of Electronic Science and Technology of China (No. Y0301902610100202)]

劉碩，在讀博士研究生，專業(yè)方向：微生物基因組學(xué)。E-mail: liushuo20022020@gmail.com

杜萌澤，博士，講師，研究方向：生物信息學(xué)。E-mail: du_mengze@foxmail.com

10.16288/j.yczz.20-022

2020/5/28 11:15:04

URI: http://kns.cnki.net/kcms/detail/11.1913.R.20200528.0949.001.html

(責(zé)任編委: 包其郁)