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

研究報(bào)告

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

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

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

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

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

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

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

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

1 材料與方法

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

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

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

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

本研究對(duì)1587個(gè)基因組中的9個(gè)模式物種基因組進(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 重注釋流程

整個(gè)重注釋的流程分為3個(gè)部分：(1)識(shí)別過度注釋為基因的序列；(2)未知功能ORFs的功能注釋；(3)識(shí)別欠注釋的基因(圖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所對(duì)應(yīng)的x、y、z可以通過下邊的公式來計(jì)算得出，據(jù)此得出Z-curve的單核苷酸的9個(gè)變量(9=3×3)。

圖1 重注釋流程圖

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

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

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

1.4 給功能未知ORFs注釋功能

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

1.6 同源簇注釋

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

2 結(jié)果與分析

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

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

圖2 1587個(gè)基因組的測(cè)序時(shí)間分布

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

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

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

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

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

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

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

2.3 4447個(gè)功能未知ORFs被注釋上確定的功能

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

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

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

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

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

2.5 9個(gè)基因組中新識(shí)別出2003個(gè)新基因

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

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

3 討論

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

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

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

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

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

A：RNA加工和修飾；B：染色質(zhì)結(jié)構(gòu)和動(dòng)力學(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：信號(hào)轉(zhuǎn)導(dǎo)機(jī)制；K：轉(zhuǎn)錄；E：氨基酸轉(zhuǎn)運(yùn)和代謝；M：細(xì)胞壁和細(xì)胞膜的生物發(fā)生；G：碳水化合物運(yùn)輸和代謝；S：功能未知；R：只能預(yù)測(cè)大致功能；L：復(fù)制重組和修復(fù)。

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

總之，本研究基于序列的全局特征和局部特征，結(jié)合同源比對(duì)對(duì)多個(gè)物種的多種菌株的基因組進(jìn)行了全面的重新注釋。未來，伴隨著基因組多組學(xué)數(shù)據(jù)的增加以及相應(yīng)的注釋工具和數(shù)據(jù)庫(kù)的完善[40~43]，公共數(shù)據(jù)庫(kù)里基因組的注釋將會(huì)更加準(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: 基因組測(cè)序數(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í)力提升計(jì)劃項(xiàng)目(編號(hào)：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é)任編委: 包其郁)