牙齦卟啉單胞菌編碼基因重注釋研究

徐曉捷，計得偉，張欣悅，張無忌，張會雄*

(1.電子科技大學生命科學與技術學院，神經信息教育部重點實驗室，成都610054；

2.電子科技大學信息醫學中心，成都 610054；

3. 成都中醫藥大學針灸推拿學院，成都610000)

牙齦卟啉單胞菌編碼基因重注釋研究

徐曉捷1,2，計得偉1,2，張欣悅3，張無忌1,2，張會雄1,2*

(1.電子科技大學生命科學與技術學院，神經信息教育部重點實驗室，成都610054；

2.電子科技大學信息醫學中心，成都 610054；

3. 成都中醫藥大學針灸推拿學院，成都610000)

摘要：為了確保牙齦卟啉單胞菌生物大分子信息的準確性，對NCBI數據庫中的3株牙齦卟啉單胞菌的注釋信息進行研究。首先，準備好蛋白質編碼與非編碼序列正負樣本，用基于Z曲線理論的Fisher判別法對正負樣本集進行訓練，確定一個判斷ORF編碼或非編碼的閾值t0，由閾值作為判別條件來識別所有的ORFs，判斷基因片段是否具有編碼蛋白質的功能，由此閾值為判別標準排除掉3株牙齦卟啉單胞菌基因組中錯誤的基因注釋信息。然后，用Prodigal基因預測軟件對牙齦卟啉單胞菌進行基因預測，基因預測結果與原始功能已知基因進行比對，挑選出具有不同5’終端的ORFs，將這些具有不同5’終端的ORFs與功能已知的基因片段進行比對，找到重疊率小于20%的候選基因。最后，對這些候選基因用Blast進行序列比對找到滿足條件的新基因，并為這些新基因添加功能注釋信息?；谝陨戏椒ü才懦?17個非編碼的開放式閱讀框，并找到了30個NCBI數據庫中缺失的編碼蛋白質的新基因。

關鍵詞：牙周??；牙齦卟啉單胞菌；基因重注釋；新基因

牙周疾病是常見的危害人類牙齒的主要口腔疾病。而牙齦卟啉單胞菌被認為是牙周疾病最重要的致病菌之一，與多種牙周疾病有密切關系。牙周炎是一種慢性口腔疾病，破壞牙齒支持組織，包括膠原蛋白、纖維和骨骼。牙周疾病是由細菌引起的一類感染性疾病，而牙齦卟啉單胞菌(Porphyromonasgingivalis,P.gingivalis)被認為是牙周疾病最重要的致病菌之一。且與成年人、青少年的牙周炎、牙周膿腫、牙槽骨膿腫、牙髓感染以及難治性牙周炎有關。牙齦卟啉單胞菌是牙周病細菌病因學研究的熱點[1]。牙齦卟啉單胞菌不僅可以引起發炎，它還與動脈粥樣硬化以及肥胖病的發生有關[2-5]，且牙齦卟啉單胞菌引起的口腔感染能夠通過侵犯主動脈的組織循環加速內皮細胞凋亡[5]，造成內皮功能紊亂，許多研究描述了牙周炎導致內皮功能障礙，可通過牙周治療來改善內皮功能[6]。Curtis等發現，在牙齦卟啉單胞菌W50菌株的55-kDa大外膜上存在著一個由重組活化基因(Recombination activation gene, rag)B編碼的相對分子質量為免疫顯性表面抗原，與牙周病患者的免疫球蛋白G抗體能否發揮作用有密切關系[7]。通過揭示牙齦卟啉單胞菌生物大分子(如核酸、蛋白質等)的結構，并探索其在遺傳信息和細胞信息的傳遞方式，有助于研究牙齦卟啉單胞菌的致病機理，為研究牙周疾病提供依據。

在基因組公共數據庫中已有牙齦卟啉單胞菌基因組的功能注釋信息，但是由于很多原因，都有可能造成基因組注釋出現有蛋白質功能編碼基因被丟棄，或非編碼蛋白質功能編碼基因被錯誤標記為功能編碼部分的情況出現?？赡墚敃r基因組數據庫數據量的局限性，或相似基因注釋存在錯誤等，導致基因預測軟件會產生一部分錯誤注釋的基因，即非編碼的開放式閱讀框被預測為編碼基因。這就需要研究人員定期對基因組注釋信息進行更新。如Bocs等就在26個原核生物全基因組中就發現34%的基因是被錯誤注釋的[8]。還有一種情況是一些真正編碼蛋白質的基因，由于種種原因卻被丟棄掉了，可以通過一些從頭預測的基因查找工具結合基因相似性比對來探測這些基因并為它們添加正確的生物功能信息。近幾年，隨著基因測序技術的快速發展，尤其是第二代基因測序技術的出現，越來越多的微生物基因組完成了測序，并被上傳至公共核苷酸數據庫。大量的基因序列數據為人們挖掘更多的生物信息提供了絕佳的機會。與此同時，這也對基因注釋信息的準確性提出了更高的要求[9]。如果一個物種的基因組注釋出現了錯誤，那么不僅會影響基于此基因組的后續研究工作，還可能導致與此基因組具有親緣關系的其他基因組的相關研究工作出現問題，因此為了保證基因注釋信息的準確性，需要對數據庫中已測序基因組的注釋信息進行定期的檢查[10]。

針對以上問題，下載了NCBI數據庫中最新的牙齦卟啉單胞菌全基因組的注釋信息，用基于Z曲線理論的Fisher判別法識別假設基因，排除3株牙齦卟啉單胞菌數據庫中被錯誤注釋的假陽性的開放式閱讀框(Open reading frames, ORFs)，共排除了117個非編碼ORFs。增加新基因，即一些真正的能編碼蛋白質的基因，由于種種原因被丟棄掉了，需要用基因預測工具并結合基因相似性比對，或通過實驗手段探測這些數據庫中丟失的基因并為它們添加正確的生物功能注釋信息。如Zhou等就通過轉錄分析和相似性搜索相結合的方法為野油菜黃單胞菌(Xanthomonascampestris)添加了306個新蛋白編碼基因[11]。用Prodigal基因預測軟件對3株牙齦卟啉單胞菌進行基因預測，把預測基因與原始基因注釋信息進行比對，保留重疊率低于20%的預測基因為候選基因，并通過Blast對候選基因進行比對，滿足條件的則被認為是要找的新基因，共找到了30個NCBI數據庫中缺失的新基因。

1材料和方法

1.1　數據來源

本研究所用的數據主要由兩部分組成，一部分是牙齦卟啉單胞菌的全基因組各染色體DNA序列文件(文件擴展名為fna)，另一部分是該物種對應的基因在染色體上的位置分布及編碼蛋白質功能信息等基因注釋數據(文件擴展名為ptt)。這兩部分數據都可以從美國國家生物技術信息中心(NCBI)所提供的核酸序列公開數據庫(GenBank)的Ftp下載中心(ftp://ftp.ncbi.nlm.nih.gov/genomes/Bacteria/)獲得。牙齦卟啉單胞菌全基因組總共包括3個，均是完全測序且在2014年7月之前下載的，它們的全名依次是：PORPHYROMONAS_GINGIVALIS_ATCC_33277_UID58879，PORPTYROMONAS_GINGIVALIS_TDC60_UID67407，PORPHYROMONAS_GINGIVALIS_W83_UID57641，對應的參考序列號為：NC_010729，NC_015571，NC_002950。

基因組注釋文件中包含基因片段編碼蛋白質功能的描述信息，根據這些描述信息把基因分為三類。第一類是具有明確功能描述的基因，此類基因一般會有確定的基因名稱，如gyrB表示DNA旋轉酶B亞單位的編碼蛋白質。第三類是功能描述為Hypothetical Protein的基因，即在基因注釋中不能確定功能信息的假設基因。余下的基因歸為第二類基因，一般是在注釋文件中具有Family、Putative、Domain等描述詞的基因。而第三類基因中還不確定哪些基因真正具有蛋白質編碼功能，哪些不具有蛋白質編碼功能。因此本文將重點關注第三類基因。

1.2　ORFs判定

要排除基因注釋中的非編碼ORFs，關鍵在于建立一個模型和識別方法對所有需要驗證的ORFs進行判定。Z-fisher是基于Z曲線理論對假設基因進行檢驗并排除非編碼ORFs[12, 13]。在任意一個基因序列片段或ORF中，把基因序列分為3個相位，第1相位對應第1、4、7、…個堿基所在的位置；第2相位對應第2、5、8、…個堿基所在的位置；第3相位對應第3、6、9、…個堿基所在的位置。根據基因序列的Z變換原理，任意一個基因片段或ORF可由33位空間中的一個點來標識，這33個分量將用作基因編碼區的識別變量。具體理論基礎和實現過程可參考文獻[12-13]。

1.3　去除過注釋基因的過程

在重注釋過程中首先要排除錯誤注釋的基因信息。基于從頭預測的基因預測軟件(Gene finder)會產生一部分錯誤注釋的基因，即非編碼ORFs被預測為編碼基因，這部分基因需要從注釋文件中刪除。對于本步驟過程的討論可以參考文獻[9]。Zfisher是專業為檢查和排除細菌或古細菌非編碼ORFs而設計的開源服務系統，可在http://147.8.74.24/Zfisher/獲得[9]，步驟見圖1。

圖1　判斷第三類基因中的基因序列是否編碼蛋白質的流程圖Fig.1　The flowchart of judging the gene sequence whetherencoding the protein or not

1.4　查找新基因的過程

在對已測序的基因組進行注釋的過程中，為了保證較低的假陽性，一些真正編碼蛋白質的基因可能會被遺漏。本研究中使用Blast在線服務中的Blastx程序對所有候選基因的核苷酸序列進行查詢。如果一個候選基因的Blast結果同時滿足以下4個條件：(1) Evalue<1×10-20，(2) Query Cover>60%，(3) Ident>50%，(4) 候選基因與同源相似基因的長度差<20%，則此候選基因是要找的新基因[9]，并為這些新基因添加正確的基因功能信息，具體實現步驟見圖2。

圖2　用Prodigal基因預測軟件對牙齦卟啉單胞菌的基因預測及發現新基因的過程Fig.2　The process of predicting the candidate genesfrom P.gingivalis uesed Prodigal gene predictionsoftware and discovery new genes

2結果與討論

2.1　基因組大小與基因數量的線性關系

在對牙齦卟啉單胞菌基因組進行重注釋之前，先對基因組大小與基因數目之間的關系進行統計分析，本文中用到了2 638個細菌或古細菌的全基因序列及對應的基因注釋信息(包括3個牙齦卟啉單胞菌)作為統計分析對象，根據物種的基因組注釋信息可以統計出每個染色體的大小及注釋的基因數目，并繪制二者的散點分布圖(見圖3)。圖中x軸表示基因組的大小(單位為kb)，y軸表示基因數目，從圖中可以發現這2 638個細菌或古細菌的基因組大小與基因數目之間具有很強的正相關性(相關系數R=0.994)，這說明隨著物種基因組的增大，其包含的基因數目也應該隨之增多。Mira等也提出，與真核生物相比，大部分原核生物(包括細菌和古細菌等)的編碼蛋白質基因緊密的分布在染色體上[14]。此外，由于原核生物中缺少內含子，所以其基因結構比真核生物要簡單?？赡苷沁@種緊密的染色體結構以及簡單的基因結構，使得細菌或古細菌的基因組大小與基因數目間具有強征相關性。

圖3　基因組大小與基因數目關系分布圖Fig.3　Linear correlation between genome size and gene number

通過繪制基因組大小與基因數目的線性擬合線(圖中黑色虛線)，我們發現大部分細菌或古細菌分布在擬合線附近，有部分物種的注釋基因數目遠多于(或少于)擬合值。針對本文的研究對象，3個牙齦卟啉單胞菌(圖中實心圓點)，也有類似的規律。由于3個牙齦卟啉單胞菌的基因組大小比較相近(約2 300 K)，所以它們在圖中幾乎分布在同一垂直線上。我們可以發現3個牙齦卟啉單胞菌的注釋基因數目分布在擬合性兩側，在基因組大小與基因數量關系方面，這3個牙齦卟啉單胞菌未顯示出任何異常。

2.2　去除非編碼的ORFs

以P.gingivalisATCC33277為例，基于Fisher判別模型，對正負樣本集進行訓練，得到判別的閾值，然后比對所有第三類基因，根據閾值判別每一個基因片段是否真正編碼蛋白質。在P.gingivalisATCC33277中，有36個假設基因判定為非編碼ORFs(見表1)。P.gingivalisW83沒有排除的非編碼ORFs。P.gingivalisTDC60排除81個非編碼ORFs(見表2)。

表1 　P.gingivalisATCC33277中排除的36個非編碼ORFs基因片段同義號

在一個指定的細菌基因組中，所有的蛋白質編碼基因都應該有相似的核苷酸組成結構[15]，也就是說P.gingivalisATCC33277中的假設基因需要與其功能已知基因具有相似的核苷酸結構，否則將被判定為非編碼ORFs。相似性核苷酸結構的判定，正是通過判別模型來確定的，在判別模型中會根據33個識別變量確定此核苷酸序列的閾值，通過此閾值判定是否編碼蛋白質，排除這36個假設基因正是基于此判別方法[12]。下圖是P.gingivalisATCC33277菌株1 125個功能已知基因(藍色*圓點標記)和36個非編碼ORFs(黑色*圓點標記)的核苷酸散點分布圖(見圖4)。注:*圖中顏色標注見電子版(http://swxxx.alljournals.cn/index.aspx)(2015年第4期)。

表2 　P.gingivalisTDC60中排除的81個非編碼ORFs基因片段同義號

圖4　P.gingivalis ATCC33277全基因組G+C含量散點分布圖Fig.4　P.gingivalis ATCC33277genome G+C content scatter distribution

從圖中可以觀察到絕大部分的功能已知基因與非編碼ORFs相分離。而且幾乎所有的功能已知基因都位于45度對角線上方，這說明其第二相位G+C含量要低于第三相位G+C含量。而36個非編碼ORFs中絕大部分分布在45度對角線附近，這表明其第二、三相位的G+C含量基本相同。由此可見編碼功能蛋白質將會影響基因的核苷酸結構分布[13, 16, 17]。因此，由于這36個假設基因與功能已知基因具有不同的核苷酸結構，在判別模型中得到的判別值不滿足編碼蛋白質的Z曲線閾值，導致其被排除為非編碼ORFs。

2.3 找到新基因，添加功能信息

使用Blast在線服務對所有候選基因的核苷酸序列進行查詢。如果一個候選基因的Blast結果同時滿足4個條件：(1) Evalue<1×10-20，(2) Query Cover>60%，(3) Ident>50%，(4)候選基因與同源相似基因的長度差<20%，我們就認為此候選基因是要找的新基因。通過以上方法，從3株牙齦卟啉單胞菌中分別找到了不同數量的新基因。在P.gingivalisTDC60中找到了6個新基因(見表3)。這6個新基因的基因位置與原注釋中的基因位置重疊率很低，全部小于0.05%，其中還包括5個重疊率幾乎為0的新基因，即原注釋信息中幾乎沒有覆蓋到的基因。根據同源基因的功能描述確定新基因的功能信息，同時這6個新基因也被賦予各自同源基因的功能注釋信息，如新基因348 817-348 960(+)則被注釋為轉座酶(Transposase)。

表4和表5分別是P.gingivalisATCC33277和P.gingivalisW83中發現的新基因以及其相應的功能注釋信息。

表3　P.gingivalisTDC60中發現的6個新基因信息

表4　P.gingivalisATCC33277中發現的5個新基因信息

表5　P.gingivalisW83中發現的19個新基因信息

3結論與展望

基因組重注釋方法是根據Fisher判別法識別3株牙齦卟啉單胞菌所有第三類基因(假設基因)，判定基因片段是否具有編碼蛋白質功能?；诖朔椒◤?株牙齦卟啉單胞菌中共排除了117個非編碼ORFs。對牙齦卟啉單胞菌使用基于從頭預測方法的基因識別工具Prodigal查找候選新基因，并以最新的基因數據庫為基礎進行Blast在線相似性比對查找同源基因，最后根據設定的參數閾值對結果進行過濾篩選，確定滿足條件的新基因并添加對應的基因功能信息，在本文中為牙齦卟啉單胞菌共添加了30個新基因。經過本文的重注釋，可能仍然還存在未排除的非編碼ORFs和未找到的新基因。為保證結果的可靠性，使用特異性較低的方法排除非編碼ORFs(低至54%)，同時在查找新基因的過程中只保留高相似度的結果(高達99%)。隨著這兩個參數的變化，發現新基因的數量和排除的非編碼基因的ORF的數量都有可能會變化。本研究中，用Prodigal基因預測軟件識別基因位置，后續可以擴展使用更多其他的基因預測軟件對假設基因進行驗證，以確保結果的可靠性。

參考文獻

[1]黃定明, 吳亞菲. 牙齦卟啉單胞菌的分型及其致病作用[J].國外醫學: 口腔醫學分冊,2002, 29(4): 213-215.

HUANG Dingming, WU Yafei. Typing and pathogenic role of porphyromonas gingivalis aeromonas[J]. Foreign Medical: Stomatology Volume, 2002,29(4): 213-215.

[2]SHAH P K. Plaque disruption and thrombosis: potential role of inflammation and infection[J].Cardiology in Review, 2000,8(1): 31-39.

[3]KUVIN J T,KIMMELSTIEL C D.Infectious causes of atherosclerosis[J].American Heart Journal,1999,137(2):216-226.

[4]CAI Y, KOBAYASHI R, HASHIZUME-TAKIZAWA T, et al. Porphyromonas gingivalis infection enhances Th17 responses for development of atherosclerosis[J].Archives of Oral Biology, 2014, 59(11): 1183-1191.

[5]AO M , MIYAUCHI M , INUBUSHI T, et al. Infection with porphyromonas gingivalis exacerbates endothelial Injury in obese mice[J].PloS One,2014,9(10): e110519-e110519.

[6]GURAV A N. The implication of periodontitis in vascular endothelial dysfunction[J].European Journal Of Clinical Investigation, 2014,44(10): 1000-1009.

[7]HANLEY S A , ADUSE-OPOKU J , CURTIS M A . A 55-Kilodalton immunodominant antigen of porphyromonas gingivalis W50 Has arisen via horizontal gene transfer[J].Infection and Immunity, 1999, 67(3): 1157-1171.

[8]BOCS S,DANCHIN A,MéDIGUE C.Re-annotation of genome microbial coding-sequences:finding new genes and inaccurately annotated genes[J].BMC Bioinformatics,2002,3(1):1-10.

[9]GUO F B , XIONG L , TENG L , et al. Re-annotation of protein-coding genes in 10 complete genomes of Neisseriaceae family by combining similarity-based and composition-based methods[J].DNA Research,2013,20(3):273-286.

[10]CAMUS J C,PRYOR M J ,MéDIGUE C,et al.Re-annotation of the genome sequence of mycobacterium tuberculosis H37Rv[J].Microbiology,2002,148(10):2967-2973.

[11]ZHOU L,VORH?LTER F J,HE Y Q ,et al.Gene discovery by genome-wide CDS re-prediction and microarray-based transcriptional analysisinphytopathogenXanthomonascampestris[J].BMC Genomics,2011, 12(1):359.

[12]ZHANG C T , ZHANG R . Analysis of distribution of bases in the coding sequences by a diagrammatic technique[J].Nucleic Acids Research, 1991, 19(22): 6313-6317.

[13]ZHANG C T , CHOU K C . A graphic approach to analyzing codon usage in 1562 Escherichia coli protein coding sequences[J].Journal of Molecular Biology,1994, 238(1): 1-8.

[14]MIRA A , OCHMAN H , MORAN N A . Deletional bias and the evolution of bacterial genomes[J].Trends Genet,2001, 17(10): 589-596.

[15]ZHANG C T , WANG J . Recognition of protein coding genes in the yeast genome at better than 95% accuracy based on the Z curve[J].Nucleic Acids Research,2000, 28(14): 2804-2814.

[16]GUO F B. The distribution patterns of bases of protein-coding genes, non-coding ORFs, and intergenic sequences in pseudomonas aeruginosa PA01 genome and its implications[J].Journal of Biomolecular Structure and Dynamics,2007,25(2):127-133.

[17]CHEN L L , ZHANG C T . Seven GC-rich microbial genomes adopt similar codon usage patterns regardless of their phylogenetic lineages[J].Biochemical And Biophysical Research Communications,2003,306(1): 310-317.

Re-annotation ofPorphyromonasgingivaliscoding-sequences

XU Xiaojie1,2, JI Dewei1,2, ZHANG Xinyue3, ZHANG Wuji1,2, ZHANG Huixiong1,2*

(1.SchoolofLifeScienceandTechnology,KeyLabofNeuroinformationofMinistryofEducation,

UniversityofElectronicScienceandTechnology(UESTC),Chengdu610054,China;

2.MedicalInformaticsCenter,UESTC,Chengdu610054,China;

3.SchoolofAcupunctureandMassage,ChengduUniversityofTCM,Chengdu610054,China)

Abstract：To ensure accuracy ofP.gingivalisbiological macromolecules information,we investigated the annotations of the 3P.gingivalisbased on NCBI database. Firstly, we prepared protein-coding and non-coding sequences as positive and negative samples,respectively,and used Fisher Discriminant which was designed based on Z curve theory to determine the threshold t0,which was used as the criterion to determine whether the gene encoding the protein or not. We firstly excluded the wrong annotation information from three stains ofP.gingivalisbased on the threshold. Secondly, theP.gingivaliswere predicted with the prodigal gene prediction software. We used the predicted genes compared to the original known-function genes and selected the ORFs with different 5’terminals, identified the candidate genes with overlapping rate of less than 20% from the ORFs with different 5’terminals.Finally, we used the sequence alignment software Blast to find the candidate genes that meet the conditions. We excluded 117 non-coding open reading frames, and found 30 new protein-coding genes that were not annotated in the NCBI database.

Keywords：Periodontal disease;Porphyromonasgingivalis; Re-annotation; New genes

中圖分類號：Q343.1+2

文獻標志碼：A

文章編號：1672-5565(2015)04-205-07

doi:10.3969/j.issn.1672-5565.2015.04.01

作者簡介：徐曉捷,女,碩士研究生,研究方向：生物醫學工程;E-mail:517170490@qq.com.*通信作者：張會雄，副教授，研究方向：移動互聯與公眾健康;E-mail:940351908@qq.com.

基金項目：中央高?；究蒲袠I務費(ZYGX2013J100)；2014年非全日制專業學位研究生教研教改項目(ZY2014009)。

收稿日期：2015-07-19;修回日期：2015-09-10.